FPGA加速三维CT图像重建*

邓靖飞，李建新，李 磊，闫 镔

（解放军信息工程大学，河南 郑州 450002）

FPGA加速三维CT图像重建*

邓靖飞，李建新，李 磊，闫 镔

（解放军信息工程大学，河南 郑州 450002）

针对三维图像重建的经典算法（FDK算法）在FPGA上的加速，提出了并行无等待流水线的实现方法。实验结果表明，该方法获得了较高的加速比。

FPGA；CT重建；FDK；流水线

计算机断层成像技术CT（Computed Tomography）作为一种新型的成像方式已经被广泛应用于医学、工业等领域。三维CT相对于传统二维CT有空间分辨率高，各向同性的优势[1]。但是由于三维图像重建运算量大，重建时间长的问题已成为制约其走向实用的瓶颈。

目前，重建加速研究主要集中在通过改进算法的软件加速及利用GPU、FPGA进行的硬件加速。其中，FPGA由于具有极佳的并行计算能力及可重构可定制的特点[2]，利用FPGA实现CT重建加速正逐渐引起研究人员的注意。

2002年 Miriam Leeser[3]首次利用 FPGA对二维 CT重建进行了加速，重建规模为 512^2时需要3.6 s。2003年 Iain Goddard[4]首次对三维 CT重建 FDK算法中的反投影过程用FPGA实现加速，重建规模为512^3时，反投影过程需要 38.7 s；2008年 Benno Heigl[5]用 9块 FPGA协调配合完成了FDK算法中滤波及反投影部分的加速，重建规模为512^3时，该过程共需要9 s。2009年Nikhil Subramanian[6]利用FPGA作为协处理器用Impulse c语言开发实现了二维CT重建过程的加速，重建规模为512^2时，反投影过程需要38.4 ms。

在FPGA内实现硬件加速是通过全数据流的形式处理，脱离了指令的操作。为了充分利用FPGA的片内资源以获得更高的加速效果，本文设计了一种并行无等待流水线的处理结构，同时对核心算法电路进行资源优化，在保持高度并行性的同时保证了较高的资源利用率。

1 FDK算法

FDK算法作为一种经典的三维CT重建算法是由Feldkamp等人于 1984年针对锥束成像圆形扫描轨迹提出的。该算法作为近似算法的代表具有良好的重建效果以及极佳的并行性，适用于FPGA并行加速的原理。FDK投影几何意义如图1所示。S表示射线源，f（x，y，z）表示重建物体的某一点，P（θ，u，v）表示在物体旋转 θ角度时探测器采集到的投影数据。

FDK算法可以分为以下步骤[7]：

（1）对投影数据进行加权，修正由于锥束光源角度引起的误差。

其中D表示光源到探测器的距离。

（2）对加权后的投影数据按行进行滤波。

其中g（a）为Ramp滤波函数。

（3）将滤波后的投影数据反投影到重建物体上。

在FDK算法中，反投影的计算复杂度与时间消耗都是最大的，是制约重建速度的瓶颈所在，因此本文研究在FPGA内实现反投影部分的加速方法。

2 FPGA反投影加速实现

反投影过程需要对重建物体旋转一周所采集到的数据进行处理。实际情况中采集数据的过程是离散的，且每一个分度下投影数据的处理过程不相关。基于这种可并行性，用FPGA加速的思路是并行计算反投影过程，并且在保证每一个反投影单元速度最快时并行尽可能多的反投影单元。

2.1 无等待流水线的设计

重建物体体素的反投影流程有3个步骤。首先根据图1的几何关系定位出重建体素在探测器上的位置；然后从存储器中读取相应数据；最后对所读取数据进行双线性插值。

通过流水线设计，双线性插值部分可以在每一个时钟更新一个数据，但是每更新一个数据需要从数据存储空间读取4个数据来计算，如果花费4个时钟周期来读取这4个数据，就会造成前级数据读取时间大于后级双线性插值时间，双线性插值处理单元会产生空泡。空泡的产生，不但制约了后级处理单元的计算速度，也造成FPGA内资源利用率的降低。

为了解决上述问题，提高处理速度以及资源利用率，分析双线性插值过程所读取4个数据之间的关系，如图2所示。

X、Y分别代表探测器垂直的两个方向，F（X0，Y0）为重建目标点在探测器上投影位置，A（Xn，Yn），B（Xn+1，Yn），C（Xn，Yn+1），D（Xn+1，Yn+1）为探测器上与之相邻最近的 4个有效区域，有效区域的大小即为探测器分辨率的大小，双线性插值过程为：

其中x=X0-Xn，y=Y0-Yn+1，A、B、C、D四个点是探测器上相邻两行上相邻两列的数据。针对这种相关性，在FPGA内部设计一种双组奇偶双口RAM结构，如图3所示。

图3中两组共四个存储空间均为FPGA片内设计的双口RAM。通过数据预取技术以及两组RAM之间的乒乓操作避免了双线性插值时随机访问外部存储器带来的延迟。在反投影处理中，第一组RAM处于数据写入过程时，分别向该组两块RAM写入探测器上奇数行数据与偶数行数据。同时另一组RAM中的数据进行双线性插值。计算完成后，两组RAM进行读写状态的互换，完成一次乒乓循环。在投影寻址单元中，计算出A（Xn，Yn）的纵坐标Yn，对其奇偶性进行判断，当其为奇数时，从两块RAM中所取数据与地址分别为：

当Yn为偶数时，从两块RAM中所取数据与地址分别为：

其中：width代表探测器的宽度，addrA、addrB、addrC、addrD分别代表A（Xn，Yn）、B（Xn+1，Yn）、C（Xn，Yn+1）、D（Xn+1，Yn+1）四个数据在相应 RAM内的存储地址。

通过这种双组奇偶双口RAM缓存结构的设计在不增加片内存储资源消耗的前提下一个时钟周期内取出四个数据，消除了空泡，实现了无等待流水线，在一个时钟周期内可更新一个重建点的结果。

单条反投影流水线设计原理如图4所示。图中总体逻辑与时序控制模块通过状态机实现对流水线的控制；缓存写入控制单元对投影数据的写入操作进行判断与控制；重建点生成器产生重建点坐标，并根据此坐标由读地址生成器计算双线性插值数据的地址，同时通过查表找出空间系数sinθ与cosθ；循环累加控制器完成对各分度下反投影结果的归约过程。

2.2并行流水线处理结构

基于各个分度下反投影过程的不相关性，在FPGA内设计一种基于分度的并行流水线处理结构，如图5所示。数据控制接口以及时序状态控制模块完成对输入数据的分配调度，通过多条反投影流水线并行计算后，由循环归约单元完成反投影结果的归约。

为了尽量增加N的值，对单条流水线做优化来降低其所占的资源量。在XC5VLX330内实现一个32 bit的除法要占用约4%的逻辑资源，因此需要尽量避免除法运算。在反投影过程中有3次除法，分别为：

其中：υ为重建区域单元尺寸，τ为探测器单元尺寸。这三个除法运算在时间上是顺序执行的，前两个运算用于投影点位置的计算，第三个运算用于双线性插值。为了降低这三个运算在FPGA内的资源占有量，提取出三个式子的公因式：

在FPGA内设计时，预先算出该公因式的值，然后通过移位寄存的方法进行延迟同步，使之在相应的节拍打入到指定的计算单元。虽然这样增加了乘法运算，但是将三次除法运算优化为一次，节省了大量的资源。表1给出了优化前后资源占有情况对比。

表1优化前后资源比较

Block Ram/kbits乘法计算量除法计算量Number of slices LUTs Number of Flip Flops Number of occupied slices可并行流水线数量优化前1 620 15 3 11 468 19 806 5 560 7优化后1 656 19 1 7 516 10 347 2 985 14对比增加2.2%增加26.7%减少66.7%减少34.5%减少47.8%减少46.3%增加100%

通过优化，在XC5VLX330 FPGA内部可以并行14条反投影流水线，相比优化以前，流水线数增加1倍。

3 实验结果与分析

通过编写Verilog语言程序，在XC5VLX330 FPGA上进行综合仿真，采取14条并行反投影流水线对Shepp-Logan标准体模进行重建，得到图6所示结果。

图6（a）为CPU重建图像的一个切片，图6（b）为FPGA加速重建的图像切片，图 6（c），图 6（d）分别为图 6（a），图 6（b）切片图像中心横线的统计图，从图中可以看出FPGA加速重建结果的精度达到要求。表2给出FPGA与CPU对不同规模数据进行重建时速度对比，计算机CPU为Inter Xeon E5430 2.66 GB，内存为16 GB。

可以看出，通过FPGA加速FDK算法中的反投影过程加速比达到了115，具有显著效果。获得加速比的原因由FPGA的体系结构所决定，其基于数据流的处理方式脱离了指令的操作，保证每一个时钟周期都是用来计算。当采取多条流水线并行处理时，加速比进一步提高，加速比与在FPGA内并行流水线数目成正比。

表2重建不同规模物体时FPGA与CPU速度对比

重建规模128^3×360 256^3×360 512^3×360 CPU 41 s 330 s 2 656 s FPGA 0.36 s 2.87 s 23.1 s加速比114 115 115

本文在FPGA上实现了对三维CT重建过程的加速。针对FDK算法中计算复杂度最高的反投影过程，通过双组奇偶双口RAM的缓存结构实现了无等待流水线，达到每一个时钟周期可更新一个重建点的速度。另外通过优化电路设计降低了单条流水线的资源占有率，在XC5VLX330上实现了14条流水线的并行处理，在保证重建图像质量的同时，取得了115倍的加速比。在利用FPGA实现CT重建加速时，影响加速效果的主要因素是FPGA内部资源利用率以及数据传输效率，随着未来FPGA以及存储器技术的发展，FPGA可以实现更快的CT重建速度。

[1]包尚联.现代医学影像物理学[M].北京：北京大学医学出版社，2003.

[2]刘佳，焦斌亮.FPGA的发展趋势及其新应用[J].电子技术，2008（4）.

[3]LEESER M.Parallel-beam backprojection：an FPGA implementation optimized for medical imaging[J].Proc of the Tenth Int.Symposium on FPGA.2002（2）：217-226.

[4]GODDARD I.High-speed cone-beam reconstruction：an embedded systems approach[J].Proceedings of SPIE，2003：483-491.

[5]Benno Heig.High-speed reconstruction for C-arm computed tomography[C].Proceedings of the 9th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine.2008：25-28.

[6]Nikhil Subramanian.A C-to-FPGA solution for accelerating tomographic reconstruction[D].University of Washington，2009.

[7]张剑，陈志强.三维锥形束 CT成像 FDK重建算法发展综述[J].中国体视学与图像分析，2005（2）：116-121.

FPGA accelerate 3D CT reconstruction

DENG Jing Fei，LI Jian Xin，LI Lei，YAN Bin
（The PLA Information Engineering University,Zhengzhou 450002，China）

：In this paper,the acceleration of the classic three-dimensional reconstruction FDK algorithm based on FPGA is researched.The method of parallel no-wait pipelines is proposed.The experimental results show that the method obtains a higher speedup.

FPGA；CT reconstruction；FDK；pipeline

TP391

A

0258-7998（2010）09-0050-04

国家 863计划重点项目（2009AA012201）

2010-03-15）

邓靖飞，男，1985年生，硕士生，主要研究方向：CT成像技术与应用。

闫镔，男，1976年生，副教授，主要研究方向：X射线成像及三维图像处理。