基于现代GPU的实时锥束重建算法研究

第３Ｏ卷第９期　２０１０年９月　核电子学与探测技术　

Ｎｕｃｌｅａｒ　Ｅｌｅｃｔｒｏｎｉｃｓ＆Ｄｅｔｅｃｔｉｏｎ　Ｔｅｃｈｎｏｌｏｇｙ　Ｖ０１．３０　Ｎｏ．９　

Ｓｅｐｔ．　２０１０　

基于现代ＧＰＵ的实时锥束重建算法研究　孙毅刚，孙修宇，张红颖　（中国民航大学航空自动化学院，天津３００３００）　

摘要：ＦＤＫ算法是一种被广泛应用的小锥角近似三维重建算法，但其运算量大，重建时间长。因　此，本文提出了一种基于现代ＧＰＵ实现的实时ＦＤＫ加速算法。该方法充分结合现代ＧＰＵ设备，通过灵　活有效的并行策略减少反投影所需要的计算量；并利用一种新的基于角度相关性的优化方案进一步提　升该算法的性能。大量实验表明该算法在不牺牲重建质量的前提下可以获得更高的加速性能。　关键词：ＦＤＫ；图形处理器；锥束重建　中图分类号：ＴＰ　３９１　文献标识码：Ａ　文章编号：０２５８－０９３４（２０１０）０９．１２６０－０６　

如何提高重建速度是当前锥束ＣＴ机研制　中亟待解决的问题之一。近年来，随着ＧＰＵ技　术的发展，基于ＧＰＵ实现的锥束重建加速算法　逐渐成为ＣＴ研究领域中的一个热门课题　ＭｕｅＵｅｒ【ｌ　等人首先给出了基于２Ｄ纹理映射的　ＳＡＲＴ加速算法，Ｘｕ　等人基于可编程管线，利　用３Ｄ纹理映射技术加速反投影操作，Ｇｕｏｒｕｉ　Ｙａｈ　等人根据ｚ轴对称性改进了文献［２］中　方法，并将其重建时间减少到５．２　Ｓ。尽管利用　传统图形ＡＰＩ加速ＣＴ重建的方法已经取得了　一定成果，但是这类方法在一定程度上增加了　算法实现的复杂性，给实际开发带来很多困难；　同时相对固定的映射模式降低了这类并行加速　算法的灵活性，无法充分发挥并行运算的优势。　Ｏ７年ＮＶＩＤＩＡ公司推出了基于统一设备构架　（Ｃｏｍｍｏｎ　Ｕｎｉｔｅｄ　Ｄｅｖｉｃｅ　Ａｒｃｈｉｔｅｃｔｕｒｅ，ＣＵＤＡ）的　ＧＰＵ。利用ＣＵＤＡ所提供的通用编程技术，重　收稿日期：２００９－０９－０３　基金项目：国家自然科学基金委员会与中国民用航　空总局联合资助项目（６０６７２１６８）；民航总局科技基　金项目（ＭＨＲＤＯ７ｚ４５）。　作者简介：孙毅刚（１９６３－），男，山东汶上人，教授，博　士，主要研究方向为ＣＴ重建，机场安检研究。　１２６０　建算法可以被当作普通的多处理器任务在　ＧＰＵ中执行。ＳｃｈｅｒｌＬ４　等人正是针对该构架提　出了ＣＵＤＡ－ＦＤＫ加速算法，并给出了实现该算　法的必要步骤。虽然文献［４］中的加速方法获　得了比传统ＧＰＵ加速算法更好的加速比，但是　其对算法本身计算复杂度的简化程度不够，并　且与ＣＵＤＡ构架的结合程度不高。因此本文　针对上述不足做了进一步改进，提出了更适用　于ＣＵＤＡ构架的并行加速策略，同时利用三角　函数关系进一步提升该加速算法的性能。与文　献［４］中ＣＵＤＡ—ＦＤＫ算法相比，本文提出的算　法在不牺牲精度的前提下获得了更显著的加速　比，可以达到实时重建的效果。　

１　ＦＤＫ算法简介　ＦＤＫ算法　是一种适用于平面圆轨道的　锥束重建算法。该算法属于滤波反投影类近似　算法，在锥角较小的情况下可以获得很好的重　建效果　。　，并且容易被改进到螺旋轨道上。图　１给出了圆弧轨迹下的全扫描示意图。　ＦＤＫ重建算法包括如下三个步骤：（１）对　投影数据进行加权处理；（２）对加权后的数据　进行横向一维滤波；（３）对滤波后的数据进行　反投影操作。公式（１）给出了ＦＤＫ算法的数学　图１　圆形扫描轨迹下的坐标系　表达式。　

ｙ，／３）：Ｊ［　

ｐ　（１３，０（　，ｙ，／３），ｂ（ｘ，，，，ｚ，／３））ｄ］　（１）　其中　

口（　，），，　）＝尺　（２）　

，ｙ，ｚ，卢）＝ｚ南（３）　Ｕ（　，Ｙ，　）＝Ｒ＋Ｙｒ（ｘ，ｙ，／３）　（４）　ｒｒ（ｘ，ｙ，／３）＝一ｘｓｉｎ／３＋ｙｃｏｓ／３　（５）　Ｙｒ（ｘ，Ｙ，ＪＢ）：ｘｃｏｓ／３＋ｙｓｉｎ／３　（６）　

６）＝（　，／Ｒ　ａ丽ｂ　＋　＋　

Ｐ　（ＪＢ，ａ，６））　（ｎ）　（７）　上述公式中　（　，ｙ，　）表示待重建的图　像；　表示在投影角度／３下的投影数据；　（ａ）　Ｄ　表示滤波函数；＿　兰Ｌ＿　表示预加权因子；　√　‘

＋口‘＋ｂ　

Ｐ　是经过预加权和滤波后的数据；（　，Ｙ，ｚ）表　示笛卡尔坐标系中的坐标；（ａ，ｂ）表示虚拟探　测器平面坐标系下的坐标。从上述公式可以看　出，ＦＤＫ算法的时间复杂度为Ｏ（　），随着图　像尺寸增大，该算法将消耗大量时间。但是对　于不同待重建点来说，其重建过程相对独立，这　为该过程的并行化提供了可能。因此，本文把　待重建点作为线程划分的基本依据并结合目前　最新型的ＣＵＤＡ构架，提出了一种基于现代　ＧＰＵ实现的实时ＦＤＫ加速算法。　２基于现代ＧＰＵ的ＦＤＫ算法　与传统ＧＰＵ加速算法不同，ＣＵＤＡ类加速　算法　可以更合理地划分ＦＤＫ算法的实现步　骤；并且对于该重建过程中的不同阶段采取与　之相适应的并行加速策略。但是文献［４］对　ＦＤＫ算法本身计算复杂度的简化程度不够，并　且与ＣＵＤＡ构架的结合程度不高。本文针对　上述两点不足做了充分的改进，提出了更高效　的ＣＵＤＡ—ＦＤＫ加速算法。该算法的实现可以　简单分成预处理（数据传输，加权和滤波）和反　投影两个步骤。　２．１　ＧＰＵ加速预处理　三维ＣＴ重建通常需要大量的投影数据。　而当今主流ＧＰＵ的显存容量并不足以用来同　时保存待重建图像和投影数据，并且ＧＰＵ设备　与主机间的数据传输的速度较慢。因此如何有　效利用ＧＰＵ设备内存成为首先要解决的内容。　本文采用了基于待重建图像的数据分配方式，　即首先在ＧＰＵ的全局内存中开辟待重建图像　的存储空间；然后将投影数据批量传输到ＧＰＵ　上，再根据已有投影数据对目标对象进行部分　重建；待所有重建操作完毕后再将重建好的图　像一并传输回主机或直接通过ＯＰＥＮＧＬ等图　形ＡＰＩ在显示器上显示。这种数据存储方式可　以与后文提出的并行策略有效地结合，从而获　得更高的加速性能。　在预加权和滤波过程中，本文改进了文献　［４］中所提出的两行同时滤波方案，将投影数　据中的两行实数数据合并成一行复数数据，利　用ＣＵＦＦ１＂库　进行傅里叶变换与反变换操作。　在ＧＰＵ设备中，一次数据访问往往需要上百个　时钟周期，这远大于流处理器进行普通运算所　花费的时间。但是在ＣＵＤＡ构架中，上述两种　操作是可以并发执行的。因此本文并没有单独　处理预加权过程，而是提出在数据合并的同时　对投影数据进行预加权处理。利用每个合并过　程中所需要的三次访存操作来掩盖预加权操作　的时间消耗。　２．２　ＧＰＵ加速反投影　本文采取基于体素点的反投影方案，并以　体素点的ｚ轴相关性作为线程划分的依据。　在这种方案下，并发的一组线程需要对经过预　处理的数据进行多次非线性拾取，这会造成　ＣＵＤＡ构架中访问冲突现象的发生，增加不必　要的访存等待时间。不过在另一方面，反投影　操作只会对上述数据进行只读访问。正是基于　这点，本文首先将经过预处理的数据绑定到二　１　２６１　维纹理内存上，利用纹理内存良好的全局随机　访问性能和硬件优化过的二维寻址系统加速对　数据的拾取和插值操作。　与文献［４］中的并行加速策略不同，本文　提出的并行策略可以更充分地利用体素点之间　的Ｚ轴相关性简化运算；并且与ＣＵＤＡ构架结　合得更加紧密，同时便于拓展。图２给出了该　策略所采用的线程划分方式。　

图２本文提出的线程划分方式　本文按照图２所示，根据栅格化后的ｘ．Ｙ　平面将重建任务划分成多个线程块。线程块中　的线程各自完成对应Ｚ轴上一系列点的重建任　务。通过公式（３）可以得出同一线程内各点之　间的ｚ轴相关性的数学表达式：　６（　ｏ，Ｙｏ，　１，　）＝６（　０，Ｙｏ，　０，　）　

．Ｌ　墨　＾　。　（　０，Ｙｏ，ＪＢ）　

其中　＝Ｚ０＋ａｚ　（８）　上文中给出的这种并行加速策略不仅简化　了ＦＤＫ算法本身的计算复杂度，同时有效地结　合了ＣＵＤＡ构架的相关特性，具有如下４个特　点：　（１）单一线程内的待重建点拥有相同的　，　Ｙ和ａ值，并且只需利用公式（３）计算一次体素　点对应的ｂ值，然后利用公式（８）简化线程内　其他体素点ｂ值计算，同时线程之间相互独立　不需要额外通信消耗。　（２）虽然线程块内线程对预处理后数据的　拾取和插值仍然是非线性的，但是其访问具有　方向一致性。这种一致性可以保证对纹理缓存　的充分利用，加快访问速度。　（３）线程内的连续访存操作可以与索引计　算并发执行，利用访存时间来掩盖大量索引计　算时间。　（４）在ＣＵＤＡ构架中，算法的性能会受到　１　２６２　

Ｂｌｏｃｋ（ｏ，ｏ】　多方面因素影响，单一线程对寄存器的消耗量　就是其中的重要指标之一。通过最大化线程内　的公共变量以及应用公式（８），本文提出的算　法可以最小限度地使用寄存器，从而获得更高　的运行效率。　２．３　ＧＰＵ加速反投影的进一步优化　由于上述加速策略对索引计算的有效掩　盖，使得内存访问变成反投影过程中消耗时间　最多的操作。因此，如何通过减少访存次数以　获得更高加速比成为当前要解决的最大问题。　基于这点，本文提出在一次反投影过程中同时　将多组预处理后的数据绑定到２Ｄ纹理内存　中，从多个角度下同时重建物体。并将所得到　的中间结果存放到ＧＰＵ的片上资源中，通过这　种方法减少访存次数和ｋｅｒｎｅｌ的启动次数。然　而单一地增加线程内的计算复杂度势必会导致　算法对ＧＰＵ资源过度消耗，从而限制算法性能　的发挥。因此，本文利用三角关系进一步优化　资源使用率，并简化索引运算。　在ＧＰＵ中，一个三角函数的运算需要３２　个时钟周期，而一个乘加运算只需要４个时钟　周期。所以可以通过公式（９）来简化三角函数　运算。　ｃｏｓ（　）＝ｃｏｓ（　０）ｃｏｓ（　１）一ｓｉｎ（￣ｏ）ｓｉｎ（ｆｌ】）　ｓｉｎ（　）＝ｓｉｎ（　）ｏｏ￣（ｐ　）＋ｃ０ｓ（　）ｓ　（　）　（９）