（计算数学专业论文）gpu通用计算在ct中的应用.pdf

摘要 伴随着p c 级微机的崛起和普及，多年来计算机图形处理器( g r a p h i c s p r o c e s s i n gu n i t ，简称g p u ) 以大大超过摩尔定律的速度高速发展。图形处理器的 发展极大地提高了计算机图形处理的速度和图形质量，并促进了与计算机图形相 关应用领域的快速发展。与此同时，图形处理器绘制流水线的高速度和并行性以 及近年来发展起来的可编程功能为图形处理以外的通用计算提供了良好的运行 平台，这使得基于g p u 的通用计算成为近两三年来人们关注的一个研究热点。 在c t 领域，g p u 通用计算技术同样有很大的应用空日j ，越来越多的人丌始 研究g p u 通用计算技术在c t 中的应用，本文利用g p u 通用计算技术实现了三 维c t 数据的仿真。用c p u 进行三维数据的仿真，其仿真速度比较慢，而且模 型越复杂所花费的时间越多。由于图像的渲染过程与c t 的数据扫描过程很相似， 受此启发，本文把c t 数据扫描的过程看作图像渲染的过程，利用g p u 通用计 算技术在g p u 上实现了c t 数据仿真。实验数据表明用g p u 进行数据仿真与传 统的c p u 进行c t 数据的仿真效率得到大大提升。同样道理把c t 重建的过程可 以看作投影的反过程，受此启发，本文利用g p u 通用计算技术实现了平行扇束 c t 图像重建，其中投影数掘是2 5 6 * 2 5 6 ，3 6 0 个角度投影，重建体大小是 2 5 6 * 2 5 6 * 2 5 6 ，测试结果表明g p u 上重建时问比c p u 下重建时i 日j 快了4 倍左右。 本文还利用g p u 通用计算技术实现了快速傅立叶变换并将结果应用到c t 滤波反投影重建算法中，将滤波反投影算法统一在g p u 进行。实验数据表明 g p u 上滤波效率不如在c p u 滤波效率，而且数据越大对比越明显，下一步将 探索这些问题产生的原因，并进一步改进g p u f f t 算法。 关键词： 图形处理器、g p u 通用计算、可编程性、c t 数据仿真 a b s t r a c t a sg r a p h i c sp r o c e s s i n gu n i t ( g p u ) h a sb e e nd e v e l o p i n gr a p i d l yr e c e n t l yw i t ha s p e e do v e rm o o r sl a w , v a r i o u sa p p l i c a t i o n sa b o u tc o m p u t e rg r a p h i c sh a v eg r o w n g r e a t l y a t t h es a m et i m e , t h eh i g h l yp r o c e s s i n gp o w e r , p a r a l l e l i s ma n d p r o g r a m m a b i l i t ya v a i l a b l en o w a d a y s o i lt h ec u r r e n tg p um a k et h eg e n e r a l 。p u r p o s e c o m p u t a t i o na v a i l a b l e t h e r ea r em a n ya p p l i c a t i o n sa b o u tg e n e r a l p u r p o s ec o m p u t a t i o no i lg p ui nc t i nt h i sp a p e r , w em a k es i m u l a t i o n so ft h et h r e e - d i m e n s i o n a lc td a t au s i n gt h e t e c h n i q u eo fg e n e r a l p u r p o s ec o m p u t a t i o no ng p u c o m p a r i n g t ot h ec o n v e n t i o n a l m e t h o d s ，e f f i c i e n c yo fd a t as i m u l a t i o nf r o mo u r m e t h o dh i g h l ya d v a n c e s t h e nw e m a d et h ep a r a l l e l 胁一b e a ms c a nu s i n gt h et e c h n i q u eo f g e n e r a l - p u r p o s ec o m p u t a t i o n o ng p u o u rn u m e r i c a le x p e r i m e n t ss h o wt h a tt h es p e e do fi m a g er e c o n s t r u c t i o ni s h i g h l ye n h a n c e dc o m p a r i n g t ot h ec o n v e n t i o n a lm e t h o d su s i n gt h ec p u i nt h i sp a p e r , w ei m p l e m e n tt h ef a s tf o u r i e rt r a n s f o r ma l g o r i t h m ( m e t h o d ) o n g p u u s i n gt h et e c h n i q u eo fg e n e r a l p u r p o s ec o m p u t a t i o no ng p u t h e n ，w ea p p l y t h i sr e s u l tt ot h ef i l t e r i n gp a r to ft h ef i l t e r i n gb a c k - p r o j e c t i o nr e c o n s t r u c t i o n a l g o r i t h m k e y w o r d s ：g r a p h i c sp r o c e s s i n gu n i t ，g e n e r a l p u r p o s ec o m p u t a t i o n ，p r o g r a m m a b i l i t y , c td a t as i m u l a t i o n h 首都师范大学学位论文原创性声明 本人郑重声明：所呈交的学位论文，是本人在导师的指导下，独立进行研 究工作所取得的成果。除文中已经注明引剧的内容外，本论文不含任何其他个 人或集体已经发表或撰写过的作品成果。对本文的研究做出重要贡献的个人和 集体，均已在文中以明确方式标明。本人完全意识到本声明的法律结果由本人 承担。 学位论文作者签名：镧f 专磊 h 期：冲卵r 首都师范大学学位论文授权使用声明 本人完全了解首都师范人学有关保留、使用学位沦文的规定，学校有权保 留学位论文并向国家主管部门或其指定机构送交论文的电子版和纸质版。有权 将学位论文用于非赢利目的的少量复制并允许论文进入学校图书馆被查阅。有 权将学位论文的内容编入有关数据库进行检索。有权将学位论文的标题和摘要 汇编出版。保密的学位论文在解密后适用本规定。 一龇荡豫肇 嗍棚闩 第1 章引言 1 1 本文的研究背景及意义 伴随着p c 缴微机的崛起和普及，多年来计算机图形处理器( ( g r a p h i c s p r o c e s s i n gu n i t ，简称g p u ) 以大大超过摩尔定律的速度高速发展。图形处理器的 发展极大地提高了计算机图形处理的速度和图形质量，并促进了与计算机图形相 关应用领域的快速发展。与此同时，图形处理器绘制流水线的高速度和并行性以 及近年束发展起来的可编程功能为图形处理以外的通用计算提供了良好的运行 平台，越来越多的人丌始用其做一些非绘制方面的计算。这些计算涉及的范围很 广，从图形输出流水线以外的非绘制处理，到几何计算m2 1 、碰撞检测3 1 、代数 运算【5 _ 7 】、优化计算【l l 】、偏微分方程p d e s ( p a r t i a ld i f f e r e n t i a le q u a t i o n s ) 数值求解等【89 _ 0 1 ，不一而足这使得基于g p u 的通用计算成为近两三年来人们 关注的一个研究热点，在c t 数据仿真以及c t 重建算法领域，g p u 通用计算技 术也逐渐引起大家的关注。 在c t 数掘仿真领域，由于图像处理中图像的渲染过程与c t 的数据扫描过 程很相似，所以与图像渲染相关的g p u 通用计算技术，在这方面会有很大的应 用空间。c t 过程可以分为三部分：数据采集、图像重建、图像显示。其中重建 算法的设计与数据采集模式直接相关，而通过重建算法重建的图像可能在不同程 度下出现各种伪影，为了更好的分析算法及改进算法，分析清楚各种伪影产生的 机理，进行数学仿真是通常必要的。通常情况下的仿真是利用数值积分方法在 c p u 上进行的，但速度慢。在c p u 上进行三维数据的仿真，2 5 6 * 2 5 6 面探测器、 3 6 0 个投影的数掘通常用需要几十分钟，而且模型越复杂所花费的时h j 越多做。 由于图像的渲染过程与c t 的数据扫描过程很相似，受此启发，本文将c t 数据 扫描的过程看作图像渲染的过程，结合g p u 通用计算技术在g p u 上实现了三维 c t 数据的仿真，与传统的c p u 进行三维c t 数据的仿真相比效率得到大大提升， 在数撕：仿真中具有较高的应用价值。 在c t 重建算法领域，由于c t 图像重建与计算机图形学的渲染过程非常相 似，c t 理论中的平行束投影和锥束投影的概念分别与计算机图形学中的正交投 影和透视投影的概念一致。c t 理论中的x 光源和探测器分别对应于计算机图形 学中的摄像机和渲染目标。c t 理论中的投影地址的概念也相应于图形学中纹理 l 坐标的概念。由于g p u 利用专用的纹理映射硬件计算纹理坐标( 投影地址) ，并 能够并行处理多个被投影的数据，利用g p u 可以将投影操作的速度提升一到两 个数量级( 与用c p u 执行投影操作相比较) 。根据这些特点本文应用g p u 通用 计算技术实现平行扇束c t 的图像重建，投影数据足2 5 6 * 2 5 6 ，3 6 0 个角度投影， 重建体大小是2 5 6 * 2 5 6 * 2 5 6 ，测试结果表明g p u 上重建时问比c p u 下重建时l 日j 快了4 倍左右。 本文也将g p u 通用计算技术应用到滤波反投影算法中。滤波反投影算法分 为两个步骤滤波、反投影，滤波占5 ，反投影占9 5 。目前反投影部分 已经可以再g p u 上实现。本文利用g p u 通用计算技术在g p u 上实现了f f t 并应用到滤波反投影算法的滤波部分，将滤波反投影算法统一在g p u 进行，实 验数据表明g p u 上滤波效率不如在c p u 滤波效率，而且数据越大对比越明显， 下一步将探索这些问题产生的原因，并迸一步改进g p u f f t 算法。 1 2 本文的主要工作 本文主要做了两部分工作。第一部分利用g p u 分别实现了多层扇柬c t 数 据的仿真和多层扇束c t 的图像重建：第二部分利用g p u 通用计算技术在g p u 上实现了快速傅立叶变换。 1 )g p u 通用计算方法分别实现了三维c t 数掘的仿真和平行扇束c t 的图像 重建。首先讲述了用g p u 进行三维c t 数掘仿真的具体过程，并与传统的方式的 仿真效率进行比较；其次讲述了平行扇束c t 的图像重建。 2 )g p u 通用计算方法实现快速傅立叶变换，首先讲述了用g p u 做快速傅立 叶变换( g p u f f t ) 具体过程，并与用c p u 做快速傅立叶变换进行比较。其次将 g p u f f t 应用到滤波反投影算法中，在g p u 上实现了滤波反投影算法中滤波部 分。 1 3 本文的结构安排 第二章介绍g p u 的通用计算的发展历史及其特点。首先介绍了图形处理器 在各个时期的发展及其应用；其次介绍g p u 通用计算的基本流程，包括g p u 绘制管道的可编程并行处理能力分布在哪些部分，它们是怎样的流处理机制； 然后介绍了g p u 通用计算在应用中的两种思路基于r e g i s t e r c o m b i n e r 和 2 基于像素程序的这两种思路。最后介绍了g p u 通用计算的优点及存在的问题。 第三章主要讲述了如何利用g p u 实现多层扇束c t 的图像重建以及如何利 用g p u 通用计算实现c t 中的数据仿真。首先讲述了c t 图像重建与计算机图 形学的渲染过程相似性，提出了在g p u 上实现c t 图像重建的思路：其次讲述 了什么是平行扇束投影，以及g p u 上c t 数据的存储方式：再次讲述了平行扇 束投影的g p u 数据仿真，并与传统的c p u 数据仿真进行比较；最后讲述了如 何利用g p u 实现c t 图像重建，并给出了实验数据与传统的在c p u 上c t 图像 重建进行效率比较。 第四章主要讲述了如何利用g p u 通用计算方法来实现快速傅立叶变换。首 先回顾了什么是傅立叶变换，它的概念以及变换公式，以及它的应用领域，其次 讲述了g p u 通用计算的方法来实现快速傅立叶变换的具体过程，最后将g p u 上做 傅立叶变换的速度，以及数据精度与在c p u 上做傅立叶变换进行比较，并将 g p u f f r 应用到滤波反投影算法中，在g p u 上实现了滤波反投影算法中滤波部 分。 第血章对垒文进行总结，并提出下一步的研究工作。 第2 章g p u 的通用计算 本章主要有三部分内容。第一部分主要是图像处理器的发展，介绍了从2 0 世纪六、七十年代到2 0 世纪九十年代图像处理器的发展过程，以及目前图像处理 器的新特征；第二部分介绍g p u 通用计算的基本流程及其应用中的两种思路，包 括g p u 绘制管道的可编程并行处理能力分布在哪些部分，它们是怎样的流处理机 制；最后介绍了g p u 通用计算的优点及存在的问题。 2 1 图形处理器的发展 我们先对图形处理器的发展历史作一个简单回顾。 2 0 世纪六、七十年代，受硬件条件的限制，图形显示器只是计算机输出的 一种工具。限于硬件发展水平，人们只是纯粹从软件实现的角度柬考虑图形用 户界面的规范问题。图形用户界面国际标准g k s ( g k s 3 d ) ，p h i g s 就是其中的 典型代表。 2 0 世纪8 0 年代初期，出现g e ( g e o m e t r ye n g i n e ) t 坦】为标志的图形处理器。 g e 芯片的出现使得计算机图形学的发展进入图形处理器引导其发展的年代。 g e 的核心是四位向量的浮点运算。它可由一个寄存器定制码定制出不同功能， 分别用于图形渲染流水线中，实现矩阵，裁剪，投影等运算。1 2 个这样的g e 单元可以完整地实现三维图形流水线的功能。芯片设计者j a m e sc l a r k 以此为核 心技术建立的s g i 公司，对图形学和图形工业发展产生巨大影响。基于s g i 图 形处理器功能的图形界面g l 及其后的o p e n g l ，成为图形用户界面事实上的 工业标准。 2 0 世纪8 0 年代和9 0 年代。g e 及其图形处理器功能不断增强和完善，使 得图形处理功能逐渐从c p u 向g p u 转移。此前计算机图形特别足j 维图形的 应用主要集中在工作站以上的处理机上，图形处理流水线的大部分功能是由 c p u 完成的。现代图形处理的流水线主要功能分为顺序处理的两个部分。第一 部分对图元实施几何变化以及对图元属性进行处理( 含部分光照计算) 。第二部 分则是扫描转换进行光栅化以后完成一系列的图形绘制处理，包含各种光照效 果和合成、纹理映射、遮挡处理、反混淆处理等。当越柬越高级的图形绘制功 4 能需要在光栅一级实现时，可以利用s i m d ( 单指令多数据) 结构进行处理。 2 0 世纪8 0 年代北卡罗莱那大学的p i x e lp l a n e 系列和p i x e lf l o w 就是使用这种 技术的典型代表。 2 0 世纪9 0 年代，n v i d i a 进入个人电脑3 d 市场。1 9 9 9 年推出具有标志意 义的图形处理器，g e f o r c e2 5 6 。第一次在图形芯片上实现了3 d 几何变换和光 照计算。此后g p u 进入高速发展时期，平均每隔6 个月就出现性能翻番的新的 g p u 。 从s g i 的g e 到n v i d i a 的g e f o r c eg p u 经历了2 0 年，芯片的线宽从3 “m 缩d , no 1 3 9 m ( 2 0 0 4 年的g c f o r c e 6 8 0 0 集成了2 2 亿晶体管) ，集成电路的逻 辑设计能力提高3 0 0 倍，但处理器数掘通道接口带宽仅提高十几倍。同时对图 形处理器计算能力的需求不断增长，出现了可编程的图形处理器。以n v i d i a 和a t l 为代表的g p u 技术i f 是适应这种趋势。 现代g p u 的发展经历一个过程。在出现g p u 概念出现以前，特殊的图 形处理硬件只出现在诸如s g i 等图形工作站上。 - 第一代g p u 出现在1 9 9 8 年，这些处理器主要处理光栅化，有些处理器 支持多纹理。主要代表有n v i d i at n t 2 ，a t ir a n g 和3 d f xv o o d 0 0 3 。 一第二代g p u 出现在1 9 9 9 年，图形处理器可以处理几何变换和光照计算， 但不具备可编程性。主要代表为n v l d l ag e f o r c e 2 5 6 ，g e f o r c e2 和a t i r a d e o n7 5 0 0 。 一第三代g p u 出现在2 0 0 1 年，图形硬件流水线可看成流处理器。可在顶 点级上支持可编程性，象素级上的可编程性则比较差。这个时候采用 g p u 进行通用计算丌始出现。主要代表有n v i d i a g e f o r c e 3 ，g e f o r c e 4 ， a t ir a d e o n8 5 0 0 。 - 第四代g p u 出现在2 0 0 3 年，以n v i d i ag e f o r c ef x ，a t ir a d e o n 9 7 0 0 9 8 0 0 为代表，象素和顶点可编程性更通用化，可以包含上千条指 令。g p u 具备浮点功能，纹理不再限制在【0 ，l 】范围，从而可以用于处 理任意数组。 _ 第五代g p u ，也就是最新的代g p u ，出现于2 0 0 4 年。功能相对以前 更为丰富，灵活。顶点着色器( v e r t e xs h a d e r ) 【玎以访问纹理，支持动态 e 分支操作；象素着色器( p i x e ls h a d e r ) 支持分支操作、子函数调用、6 4 位浮点纹理滤波和融合、多个绘制目标等。主要代表有n v i d i a g e f o r c e 6 8 0 0 。n v i d i ag e f o r c e6 8 0 0 集成了2 亿2 千2 百力晶体管，具有超标 量的1 6 条管线架构。 随着计算技术和集成电路技术的发展，图形硬件的更新速度迅猛。图2 1 g p u 、c p u 浮点运算速度对比图，g p u ( g r a p h i c s p r o c e s s i n gu n i t ) 自1 9 9 9 年 首先由n v i d i a 公司提出束后，就其发展的速度而言，是c p u 更迭速度的三 倍多。目前图形芯片的主要市场被n v i d i a 和a t i 这两家公司占领，从高端到 低端都有相应的产品来满足市场。新的图形硬件带来一些新的特征，这些特征 【1 6 1 7 】概括起来有如下几方面： 1 ) 在顶点级和像素级提供了灵活的可编程特性； 2 ) 在顶点级和像素级运算上都支持i e e e 3 2 位浮点运算； 3 ) 支持多遍绘制的操作，这样避免了多次c p u 与g p u 之日j 的数据交 换： 4 ) 支持绘制到纹理的功能( r e n d e r - t o - t e x t u r e p b u f f e r ) ，从而避免将计算 结果拷贝到纹理这一比较费时的过程； 5 ) 支持依赖纹理功能，以方便数据的索引访问，可以将纹理作为内存 来使用。 l & ， l l 舯 6 ；l 幽 _ e 墨 制傍麟年嗍 ， 占 n 事 墨 兰凹u凹 嗽 摹 幽 2 2g p u 的通用计算的基本流程及其应用中的两种思路 在图形处理器可编程技术普及以前，也出现过一些基于图形卡的数值计算。 近年来，g p u 在通用计算领域的应用越来越普遍。一般晚来，目前基于g p u 的 通用计算的基本过程如图2 2 所示。图2 2 所示最右边一列表示子素处理 ( f r a g m e n tp r o c e s s i n g ) 可以选择的方法。 | 璺i2 - 2 堆十g p u 的通用计算流程i i 5 】 这罩我们主要介绍一下基于r e g i s t e rc o m b i n e r 和基于像素程序的这两种思 路，冈为这两种方法相对于其他方法更为灵活，实现的运算也更为复杂。r e g i s t e r c o m b i n e r 最先是在g e f o r c e 2 引入的，只有2 个g e n e r a lc o m b i n e r ，目前最新的硬 件上则町以支持8 个g e n e r a lc o m b i n e r ，另外还有一个f i n a lc o m b i n e r 作为最后输 出运算的。基于r e g i s t e rc o m b i n e r 的运算需要根据其支持的运算形式和有限的 c o m b i n e r 个数来设计整个的流程，如图2 3 所示，即数据以顶点颜色的形式或者纹 理的形式传a r e g i s t e rc o m b i n e r ，利用其支持的运算可以| 一j 时完成4 个颜色通道的 计算。 幽2 3r e g i s t e rc o m b i n e l 总的流程f i 5 】 目前顶点程序不管是从支持的指令数还是寄存器的个数部较以i i 有了很大的 提高，如n v i d i ac i n e f x 目l ；i 每个顶点能支持6 4 k 条指令：常量寄存器的个数和临时 寄存器个数也都较以丽有了增加；更重要的是丌始支持动态和静态控制流的循环 和分枝操作，以及子函数操作。 i 轳”薏雾、$ e n ”“i。p 渊t 棼磷f 秘1 $ 黔e | 罐a r t 锄匏$筘0 辫馥 拣 ” 一参懒雠嘲喀 l 燕套雠薹霉蓦哟燃粥嚣l 锈 誊 j ， l ” 粼蝴撼。蝴。爹黝v 蚓e t 8 x 懈喇岫獬 黝蝴撼蝴爹毳一懈m 岫獬 。 妒 黼 “ “ ”删= = 鹣燃滋= = t e m 群a t yl v # “3 t i # l # 瓣# l 巍薛御懒德藿镶 。叛蠕燃鹱秘熊 。￥张触蛾j 辩穆譬l 辩静潞 。 。黝锵鲢$ 搭霹懿。 幽2 4 顶点程序模型【1 5 l p m ，f e b v r 0 e 键 # h r t e f 3 鬈榔；l 噻 鞠f 粼糖擂； 辩糍嚏鞴嗡l 鲢e 涔4 ， 蕊 幽2 5 像素程序模型u s 而像素程序能够访问的独立的纹理目前达到1 6 个之多；支持的指令数也多达 1 0 2 4 条；寄存器的个数增j j l 羽j 6 4 个。显然跟r e g i s t e rc o m b i n e r 相比顶点编程和像 素编程更具可编程性，有着更为灵活的操作，而且它们都提供3 2 位的浮点精度， 从而对于通用计算来说更为实用化，可以设计出更为复杂的运算。在采用顶点编 程方面，如文献 4 2 1 将矢量转换为一系列的四元组，实现了矢量、矩阵的一些运 算，但由于需要光栅化这个过程，将结果送入纹理或者帧缓存中，然后再读回到 主内存进行处理。所以如果数据量太大，可能会形成数据向显卡传输的带宽瓶颈。 图2 4 和图2 5 分别给出了顶点编程和像素编程的模型，可以看出两者很相似。虽 然目前顶点编程已支持分支操作和循环，但不能访问纹理数据，而图形硬件的像 素处理能力要比顶点处理能力高一个数量级，像素填充速度远比三角形传输速度 高，且目前的像素渲染管道也要比顶点渲染管道多，所以像素程序更适合用来做 通用计算。目前大多数文献均采用这种途径，如文献【8 ，1 0 】。 基于前人的若干工作，本文采用像素编程并结合实际中的需要在通用计算方 面做了些工作。由于基于g p u 的通用计算涉及的内容很广，而且跟具体的问题紧 密联系，本文下一章节会讨论基于像素编程快速傅立叶变换g p u 实现，以及g p u 通用计算技术在c t 数掘的仿真应用。 2 3g p u 的通用计算的优点 从图形处理机的结构看，它是利用通用机的结构实现图形并行处理的。而 g p u 却是完全专用子图形输出流水线的处理和加速。因此当g p u 的功能越来 越强时，与图形仃关的处理便自然而然地从c p u 向g p u 转移。最先发生的 9 转移是最靠近应用程序的几何变换部分，其中包括造型变换和观察变换；其次 是局部或特殊光照效果的计算和生成。当顶点级和象素级的可编程功能越来越 灵活时图形本身的处理速度和灵活性都得到了前所未有的提高。而当g p u 内部象素级的纹元达到可以参与编程的运算时，它从某种程度上模拟了类似于 p i x e lp l a n e 处理单元的部分功能，以至于向着可作通用计算的方向发展。 这 时，不可思议的事情即基于g p u 的通用计算便应运而生了。基于g p u 的通 用计算( g p g p u ：g e n e r a lp u r p o s eg p u ) 指的是利用图形卡来实现一般意义上的 计算，而不单纯是绘制。就是如何利用g p u 来实现矢量、矩阵的基本代数运 算，然后在这个基础上如何实现一些相对复杂的运算，如线性方程组的求解， 从而实现复杂应用问题的求解。相对于以l j i 采用固定渲染管道的图形硬件，上 述提到的这些新特征无疑加快了g p u 在通用计算方面的应用。而采用图形硬 件来做通用计算的主要目的是为了加速，加速的动力来自这些新硬件所具有的 以下主要优势； 1 1 一定的并行性。这一功能主要是通过多个渲染管道和r g b a 四个颜色 通道同时计算来体现的，另外在一个时钟周期内町以同时获取2 个甚至更多副 纹理。顶点程序的多个渲染管道意味着一个时钟周期可以并行处理多个顶点， 而对于像素程序同样如此。相对于并行机而言，图形卡提供的并行性虽然很弱， 但它在十分廉价的基础上为很多应用提供了一个很好的并行方案，尤其是对于 图形本身的应用来说。 2 1 高密集的运算。由于图形卡内部的内存接口位宽大于c p u 上的位宽， 如g e f o r c ef x 的内存位宽达2 5 6 位，显然高于c p u 上3 2 位的位宽，这样 整个计算的带宽大大提高。g p u 相对于c p u 来说，更适应传输大块的数据， 虽然c p u 上有c a c h e 以加速整个计算过程，但c p u 上的c a c h e 相对于图 形卡显存来说太小，一般只有6 4 k b ，而现在的显存大多都在6 4 m 以上，由 此可见一斑。 3 ) 减少了g p u 与c p u 的数据通信。尤其足当整个应用针对图形生成的 时候，不再需要在c p u 与g p u 之问进行多次数掘交换，从而可以让c p u 解 放出来做其他的事情。这些优势使得g p u 比c p u 更适用于流处理计算，因 此g p u 也被认为是一个s i m d 的并行机或者流处理器，可以用于处理大规模 数据集，使得应用得到加速。而相比之下，c p u 本质上是一个杯量计算模型， 1 0 而计算单元偏少，主要针对复杂控制和低延迟而非高带宽进行了若干优化。 第3 章利用g p u 实现多层扇束c t 的图像重建 c t 图像重建与计算机图形学的渲染过程非常相似，c t 理论中的平行束投 影和锥束投影的概念分别与计算机图形学中的正交投影和透视投影的概念一 致。c t 理论中的x 光源和探测器分别对应于计算机图形学中的摄像机和渲染 目标。c t 理论中的投影地址的概念也相应于图形学中纹理坐标的概念。由于 g p u 利用专用的纹理映射硬件计算纹理坐标( 投影地址) ，并能够并行处理多 个被投影的数据，利用g p u 可以将投影操作的速度提升- n 两个数量级( 与用 c p u 执行投影操作相比较) 。反投影操作类似幻灯机将幻灯片投影到屏幕的过 程，利用g p u 也可以方便实现。因此利用g p u 加速c t 重建受到人们的广泛 关注。在这方面，b r i a nc a b r a l ! 。引、k l a u sm u e l l e 一9 2 0 一、f a n gx u 2 0 1 、k e nc h i d l o w 等人已作了大量工作，表明利用g p u 通用计算技术一般能够将重建速度提高一 到两个数量级。 三代的扇束扫描模式是商用c t 常用的扫描模式，这种模式重建的是被测 物体断层的图像，然而在很多情况下人们希望得到物体的三维信息。因此需要 对多个断层进行扫描重建，这种扫描重建的模式称之为多层扇束的扫描重建模 式。根据这种模式中不同断层的重建具有的并行性，本章利用g p u 实现了多层 扇束下的c t 图像重建算法。 3 1 多层扇束投影 三代扇束扫描模式下探测器为线探测器或弧形探测器，在扫描过程中射线 源和探测器固定不动，转台旋转一周( 如图3 1 所示) 。为了获得重建物体的三 维信息，需要对多个断层进行扫描重建，即多层扇束扫描重建模式。在此模式 下完成一次三代扫描后，射线源和探测器同步1 二升( 或下降) 一定高度，然后 做三代扫描获得此断层的扫描数据，依次f 去获得n 层的扫描数据。这种模式 可以等价成射线源和探测器同步上升( 或下降) n 次获得在同一个角度下的不 同断层的投影数据( 见图3 2 所示) ：然后转台旋转一个角度，射线源和探测器 在同步上升( 或下降) n 次获得另一个角度下的不同断层的投影数据；依次下 去获得3 6 0 度上的不同断层的投影数据。 1 2 图3 2 是平行扇束的投影模式的示意图，计算机图形学中有正交投影和透 视投影，对平行扇束投影来说没已有的投影矩阵。因此重建过程中必须构造出 平行扇束投影矩阵。 l 鳘i3 - l 扇束c t | i 描水壶| 茔i 幽3 - 2 1 7 r 扇束投影叫：怠i 茔| 3 2 g p u 下c t 数据的存储方式 、 投影平面 空i 日j 物体 g p u 程序中的数据以纹理的格式存储，g p u 下的纹理对应c p u 程序中的 数组，二维纹理和_ 二三维纹理分别与_ ：维数组和三维数组相对应。在g p u 程序中， c t 数据需要以类似f 数身l 的纹胖形式存储。纹理的纹素相应于数组的元素。 纹素的纹理坐标相应于数组元素的指标。与数组元素指标有所不同，2 d 纹理坐 标是介于【o ，1 的浮点数，0 和1 分别相应于数组元素的最小和最大指标。矩形 纹理( m * n ) 的横坐标是是介于【o ，m 】的浮点数，纵坐标是介于【o ，n 的浮点数 ( 见图3 3 ) 由于我们重建的层数并不一定是2 的整次幂，探测器的个数也不 一定2 的整次幂，因此重建中需要用到矩形纹理。通过指定寻址模式，g p u 可 以处理 0 ，m 】范围之外的纹理坐标。 g p u 下的纹理数据格式比c p u 下的数组数据格式更加灵活，它支持单 通道、四通道等纹理格式。单通道纹理的纹素是标量；四通道纹理的纹素是4 维向量。g p u 可以并行地执行向量操作，利用这种并行性，我们可以实现并行 地对四组数据进行投影或反投影。当的的g p u 增加了对浮点纹理的支持和相应 的浮点运算单元。g p u 所支持的数据格式有：8 位整型、1 6 位整型、1 6 位浮点 和3 2 位浮点等等。 对纹理的采样相应于对数组元素的读取。不同之处在于，如果纹理坐标所 指的采样点刚好与某纹素对应，则此纹素被采样；如果不对应，则需要通过指 定某种纹理滤波方式由采样点附近的纹素值得出采样值。现代g p u 支持的滤波 方式有最近点滤波( 默认滤波方式) 、双线性滤波、三线性滤波和各向异性滤波。 值得一提的是，这种g p u 硬件加速的滤波正好可以代替c t 重建的c p u 程序 的软件插值过程。 t e x t l r e1 6 ) ( 1 。6 ) 2 ，6 )3 6 ) ( 唾6 ) ( s t t e x t u r ea d d r e s s ( 2 , 5 ，3 匈 幽3 3 矫形纹理的纹理坐标 3 3 平行扇束投影的g p u 数据仿真 由于c t 的数据模拟过程与图形学中体绘制的过程非常相似，因此考虑利 用g p u 来实现平行扇束扫描模式的数值仿真。 3 3 1 体渲染 体渲染的方法有基于二维纹理的体渲染方法和基于三维纹理的体渲染方 法。 基于二维纹理的体渲染：体数据按其空i 日j 顺序被存储于一组二维纹理中。 被重建的长方体由一组等问隔的平行的长方片代表，通过纹理映射将二维纹理 贴到这些长方片上，然后渲染这些长方片，打丌a l p h a 混合，最后将渲染结果 累加即可得到三维体数掘的投影。 这种方法的优点是插值少，它只在每个二维纹理中做插值，而二维纹理之 i h j 没有插值运算，所以缺点是绘制效果差，尤其当平行的长方片与视点方向平 行时无法应用二维纹理。所以一般情况需要三组相互垂直的二维纹理，根据与 视点方向的夹角判断用那组纹理 基于三维纹理的体渲染：体数据按其空间顺序被存储于三维纹理中。被重 建体由一组等问隔平行片组成，平行片的法线方向始终与投影方向平行。 与二维纹理的方法不同，平行片组始终与视线方向垂直，因此平行片组在 每个投影角度都需要更新，向- 且需要相应地计算出平行片顶点的三维纹理的坐 标，然后就町以按不同方向采样三维纹理了。所以这种方法需要的存储量小， 但是由于每个角度都用重新计算，夭l 此计算量大，绘制效果i ：t - - - 维纹理好， 3 3 2 数值仿真 仿真方法是基于二位纹理的仿真，与二维纹理的体绘制有所不同，体绘制 的目的是将渲染的结果最终显示在屏幕七，而我们是要得到仿真数据，因此需 要将渲染后的数掘读回内存存成数据文件。并且体渲染的结果一般是8 位的数 掘，因为人眼的分辨没有那么高，而仿真数掘要求是3 2 位的，最低也要1 6 位 的数据，因此计算量要比体绘制大的多，需要用到浮点纹理。投影矩阵也有所 不同，需要构造平行扇束投影矩阵。 如图3 - 4 所示模型足行：定义在市方体内，用g p u 通用计算方法实现数据的 is 仿真( 以下简称g p u c t ) ，首先需要将数据用纹理来表示，如图所示将立方体 划分为若干切片，每个切片绑定成不同的纹理。如图3 4 所示立方体中蓝色和 黄色切片分别与纹理t e x i 、t e x j 绑定；根据需要将立方体划分为n 个切片，并 分别与纹理t e x 0 、t e x i t e x n 1 绑定，则立方体数据就可以由这n 个有 序纹理束表示。 图3 4 中黄色的方格区域为面探测器，也可以说是投影目标区我们记做 t e x o b j 。在实现仿真的过程中纹理t e x 0 、t e x l t e x n 1 将依次渲染到t e x o b j ， 所有纹理渲染结束后t e x o b j 中的数据就是模型在某一角度下的投影，模型旋转 换一个角度后，进行同样的过程就得到另一个角度的投影。 倒3 - 4 仿真1 i 愿幽 以上方式模型旋转范围在正负4 5 度之问时得到的仿真数掘是可以的，但当 模型旋转的角度超出了萨负4 5 度范围时，t e x i 与渲染目标t e x o b j 所成角度比 较大，也就是说t e x i 整片数据只渲染到t e x o b j 很狭小的一段，如图3 5 所示 t e x i 与t e x o b j 成6 0 度央角时，t e x i 中的1 6 个数据投影到t e x o b j 中8 个探测 器上，会造成较大的误差。 图3 - 5t e x io t e x o b j 兴柏 为了减小纹理t e x i 与渲染目标t e x o b j 央角过大引起的误差，我们应当始 终控制其央角在j 下负4 5 度之白j 。当模型旋转超过这个范围时，我们可以通过沿 与原纹理垂直方向重新划分切片的方法来避免角度过大引起的误差。 g p u c t 数据仿真流程如图3 - 6 所示，结合流程图来下面探讨g p u c t 数 据仿真的具体过程。 幽36 g p u c t 数据仿真流程 1 7 具体步骤： s t e pl ：首先将三维数据模型进行划分，划分为2 5 6 片数据( 可以根据具体 的需要) ，这样每一片数据大小为2 5 6 * 2 5 6 。 s t e p2 ：将划分好的每片数掘做成2 5 6 * 2 5 6 大小的纹理，得到2 5 6 片纹理， 这样三维数据模型就由2 5 6 片纹理束表示，分别编号为t e x 0 、t e x l t e x 2 5 5 。 s t e p3 - 建立两个3 2 位f l o a t i n g - p o i n tp b u f f e r 作为纹理渲染目标记作t e x o b j ， 其中t e x o b j 中初值为0 。 s t e p4 ：将t e x 0 渲染到目标t e x o b j 做混合运算t e x o b j = t e x o b j + t e x i ；继续将 t e x i 渲染到目标作混合运算，如此循环进行直到将t e x 2 5 5 渲染到目标t e x o b j 。 此时目标中的数掘就是模型在这一角度下的平行投影数据。 仿真效果： 如图3 7 所示模型的正面透视图，一个圆柱体，圆柱体内部各个方向均匀8 个小圆柱体，如图4 - 8 所示模型的侧面透视图，8 个小圆柱体所处的深度不同。 剖3 7 模型的萨由透视图i 茔 3 - 8 模型的侧由透视幽 如图3 - 9 所示，其为模型f 面旋转3 0 度后的g p u 仿真透视图，如图3 1 0 所示为模型倾斜1 5 度仿真透视图。 幽3 - 9 模型旋转3 0 度透税i ! l 幽3 1 0 模型倾斜1 5 度透视图 每片纹理大小为2 5 6 * 2 5 6 ，投影数据大小为2 5 6 * 2 5 6 时c p u 上仿真一个投 影数掘平均所用的时自j 为5 ，6 3 5 s ，g p u 上仿真一个投影数据平均所用的时间 为0 5 0 8 s ，可以看出用g p u 做c t 数据仿真效率有了很大的提高。 3 4 利用g p u 的c t 重建方法 由于重建中9 0 的时问是在反投影部分消耗的，因此我们主要考虑投影部 j 分的加速问题，滤波的g p u 实现在下章讨论。 投影数据用四通道二维纹理存储，因此根据扫描角度的对称性，将四个对 称的角度( 屈万一屈万+ p ，2 石一p ) 数据放入四个通道中( 如图3 一1 2 所示) 。 重建过程中所有的投影数据均放入显存。 重建体的数据划分有两种方法，一种是横切片，一种是竖切片。我们在实 现中用的是竖切片的方法( 图3 1 2 所示) 。显存中只放对称的p q 个切片，四个 切片放在一张四通道的二维纹理中。 9 i 角度的投影 2 石一口角度的投影 ，一一。 蕊鬈 。、o t i i 稽；。 ，! 、，d ： j i 张 ! 。? 7 ：粼 $ 孱、：、j 7 tu i ? 。_ 式，：。、：， 、簿! 釜j 二j 万+ 口角度的 幽3 - 1 1 角度对称的俯视| 堇i 辨 投影 倒3 一1 2 数据慑切”的划分方法 首先设置多层扇束投影的投影环境，利用投影纹理技术可以计算出待重建 切片投影到投影平面的位置口，6 。然后将投影数据反投影到重建切片上。将 渲染环境设置成萨交投影的方式，运用多纹理技术与其它角度的反投影数据进 行累加。切片转动4 5 度后完成反投影部分。值得注意的是多纹理技术目前的显 卡支持8 个纹理同时渲染。 j ，- 。“ 。，- 蕊、 。 ，_ 待重建掇擘的伉置 ? j 、 ，、， 一：、 蕊； - ? 、 、 辽影 、j 7 、，j ， 、 j 二：曼重型饬层数癌、 j 蓄，7 影数据 、，7 - 一 一， 3 5 实验结果 i 垒l3 一1 3 反投影的一意幽 本实验运行平台：采用的实验平厶为i n t e l ( r ) p e n t i u m ( r ) 3 0 0 g h z ，主内存 为1 o o g b ，而显卡采用的是a t i r a d e o n x 3 0 0 x 5 5 0 s e r i e s ，显卡内存为1 2 8 m ， 显卡的核心频率为4 0 0 m h z ，驱动的版本为8 1 8 0 0 重建图像见图3 1 4 3 1 4 曲个转且叫层的电娃幽像 其中投影数据是2 5 6 * 2 5 6 ，3 6 0 个角度投影，重建体大小是2 5 6 * 2 5 6 * 2 5 6 ， 测试结果是3 2 位浮点的计算精度下需要2 4 秒的时| 日j 。比c p u 下的9 3