CN101567090A - 一种自适应的锥束ct三维图像快速重建方法 - Google Patents

Abstract

本发明公开了一种自适应的锥束CT三维图像快速重建方法，对试件扫描采集一组投影图像；将该组所有投影图像按该正方形区域裁剪为边长为E个象素的正方形图像，并将裁剪后的投影图像进行对数运算；对该组投影图像进行滤波处理；设置需要重建的切片图像起止层号作为重建任务；对重建任务进行自适应分割，得到若干个相同大小的子重建任务；根据子重建任务的层数分配一个部分重建空间；判断是否还有未重建的子重建任务，若是，则在部分重建空间中进行下一个子重建任务的重建计算并继续判断，若否，则释放部分重建空间所占内存，完成重建任务。本发明解决了Z线优先重建算法的自适应处理问题，增强了适用性，减少了读取数据量，节省重建时间。

Description

一种自适应的锥束CT三维图像快速重建方法

技术领域

本发明属于CT系统图像重建领域，涉及一种对CT系统中三维图像自适应快速重建的方法。

背景技术

CT(Computed Tomography)技术利用射线(一般为X射线)透射被检测物体，并在探测器上获取一组投影图像，再配合相应的重建算法得到物体的切片图像。目前研究和应用的CT可分为二维CT和三维CT两大类，锥束CT属于三维CT。与传统二维CT相比，锥束CT扫描速度快、射线利用率高，能获得均匀、高精度的空间分辨率。锥束CT的重建是三维图像重建，即需要重建得到连续的很多层切片图像，涉及巨大的计算量。在锥束CT实际应用中，在保证良好的重建质量的基础上，提高重建速度成为必须突破的瓶颈。

目前在商用锥束CT系统中，应用最广泛的重建算法是FDK滤波反投影算法，该算法是一种近似重建，计算量相对较小，在小锥角的情况下可以获得质量良好的重建图像。即便如此，FDK算法的计算量依然较大，主要集中在反投影过程，计算复杂度为O(N⁴)，其中N为投影数据尺寸。近年来，CT系统硬件性能得到了较快发展，在平板探测器方面，Varian Medical System的PaxScan 4343R可以产生的最大投影图像为3072²，如果用360幅3072²投影图像重建3072³图像，重建算法的内部循环次数将达到10¹³次，即使重建少数层，运算量也非常大，这在普通PC机上计算是非常耗时的。

在现有的文献中所提到的提高图像重建速度的方法，主要可以分为以下三大类：(1)重建算法的改进，如几何变量分离法、增量法、图像细分法等；(2)采用特殊的硬件来加速，如采用专用的DSP模块，对反投影计算过程进行硬件化实施；(3)并行计算技术，如采用阵列处理器或分布式计算机对图像重建过程进行并行化。

毛海鹏、张定华、梁亮等人在《系统仿真学报》(2004，16(11)：2486-2489)的文章“一种基于PC的快速三维图像重建方法”中提出了一种Z线优先重建算法，它是对传统FDK算法进行改进并通过数据并行计算来实现快速三维图像重建的算法，其基本原理是：根据FDK反投影过程中反投影点的p值与z无关，重建图像空间的每条Z线可映射为滤波投影平面上的一条u向亚像素线，标记为L，L可由两条相邻的u线进行线性插值得到，最后根据L线与Z线的对应关系得到重建图像；在此基础上通过有效地组织和划分重建数据，使得对重建数据的内存访问非常连续，采用单指令多数据(SingleInstruction Multiple Data，SIMD)技术进行数据并行处理。该方法的计算效率很高，但需要对整个重建立方体空间进行内存整体分配，如果重建空间过大，就会因为内存不足而导致算法无法实施，而且实际应用中往往无需重建出整个立方体空间，只需重建出其中若干感兴趣层切片即可，这时就需要对Z线优先算法进行相应的改进。

另外，重建过程中制约重建速度的不仅仅是算法的效率，大量数据在硬盘与内存间的读写也是一个很重要的影响因素，特别是随着高分辨率探测器的使用，投影图像文件的尺寸也随之增大，这种影响更加明显。以现在的常规重建算法进行重建，约一半的时间花费在了硬盘与内存之前的数据读写。因此，提高数据读写速度或减少数据读写量也成为重建过程中重要的一个环节。

发明内容

为了克服现有技术需要对整个重建立方体空间进行内存整体分配，而且重建速度不高的不足，本发明提供一种自适应的锥束CT三维图像快速重建方法，使得Z线优先重建算法也能适用于局部重建的问题，同时通过适当的算法设计，减少硬盘与内存之间交互的数据量，从而缩短重建时间。

本发明解决其技术问题所采用的技术方案包括以下步骤：

(1)对试件进行圆周锥束CT扫描，采集一组用于重建的投影图像；

(2)在投影图像中选取一个边长为E个象素的正方形区域，将该组所有投影图像按该正方形区域裁剪为边长为E个象素的正方形图像，并将所有裁剪后的投影图像进行公知的对数运算，此时投影图像的象素灰度为单精度浮点型；

(3)对该组投影图像进行FDK算法中的滤波处理，滤波函数采用公知的S-L滤波器、R-L滤波器或SL-W滤波器；

(4)设置需要重建的切片图像起止层号，这段连续的切片图像即为重建任务；

(5)根据计算机当前的可用内存数对重建任务进行自适应分割，得到若干个相同大小的子重建任务；

(6)根据子重建任务的层数m，从计算机中分配一个E×E×m大小的单精度浮点型内存空间，此内存空间称为部分重建空间；

(7)判断：是否还有未重建的子重建任务，若是，则继续执行下一步，若否，则转第9步；

(8)在部分重建空间中进行下一个子重建任务的重建计算，完成后转第7步；

(9)释放部分重建空间所占内存，完成重建任务。

在上述方法第5步中，根据Z线优先重建算法的基本思想，总重建空间的每条Z线可映射为滤波投影平面上的一条u向亚像素L线，由此可以得出Z线上的一段Z’，必然可以对应于L线上的一段L’，因此在Z’所对应的部分重建空间中同样可以采用Z线优先重建的基本方法进行计算。换言之，Z线并非只能完整重建，而是可以根据需要分段重建，这是本发明方法的基础(图2)。根据计算机当前的可用内存数对重建任务进行自适应分割的具体步骤如下：

1)获取计算机物理内存当前的可用数，考虑到重建过程中可能会有其它程序的内存需求以及操作系统运行的内存需求，为确保计算过程的稳定，将该可用内存数减去50～150M，作为实际可用的最大内存数；

2)由重建任务的切片图像起止层号得到重建任务的总层数为N，设实际可用的最大内存数为T，重建任务所需的内存数R＝投影图像大小×N，则需要将重建任务分割为S个子重建任务，S＝int(R/T)+1；

3)考虑到子重建任务计算量的均衡性，应当使各子重建任务重建的切片层数相同，则子重建任务的层数m＝int(N/S+0.5)，考虑到后续的计算中会采用单指令多数据(SIMD)技术进行加速，因此需要将m取为m′，m′为大于或等于m的整数中能被4整除的最小数；

4)计算各子重建任务的起止层号，设重建任务的起始层号为L_s，则第i个子重建任务的起始层号为L_s+(i-1)m，终止层号为L_s+i·m-1，其中i的取值为1～S。

在上述方法第8步中，进行子重建任务的重建计算时，无需读取整幅滤波处理后的投影图像，只需读取子重建任务所需的滤波投影图像中的部分数据，以有效减少读取的数据量，从而节省总的重建时间，具体步骤如下：

1)选取部分重建空间的内切圆柱为感兴趣重建空间ROI，即每一层切片的重建区域只局限于切片图像的内切圆中，这样可将重建体素的数目减少至原始数目的π/4；

2)计算子重建任务的起止层切片图像在滤波投影图像中对应位置(图3)，分以下两种情况：

对于起始层切片图像，按以下公式计算：

对于终止层切片图像，按以下公式计算：

其中，h_s为起始层切片到中心切片的层数，h_e为终止层切片到中心切片的层数，r为子重建空间内切圆柱半径，dso为射线源焦点到工作台旋转中心的距离，dsd为射线源焦点到探测器的距离。将OZ_s取整得到P_s，OZ_e取整得到P_e，并设中心切片层号为P₀，并考虑一点冗余，则子重建任务所需的重建数据位于滤波投影图像中的第P₀±P_s-1行至第P₀±P_e+1行，当所计算层位于中心层之下时取-号，反之取+号；

3)判断：是否还有未读入的滤波投影图像，若是，则继续进行下一步，若否，则转第6步；

4)读入下一幅滤波投影图像中与子重建任务相对应的部分数据；

5)按Z线优先重建算法并采用单指令多数据(SIMD)技术进行加速计算，完成后转第3步；

6)将子重建空间的重建结果存储为Z向的序列切片图像，结束。

本发明的有益效果是：在保持Z线优先重建算法原有的重建图像质量基础上，分析出Z线并非只能完整重建，而是可以根据需要分段重建的性质，解决了Z线优先重建算法对各种分辨率三维图像重建的自适应处理问题，增强了算法的适用性。另外，子重建任务的重建计算只需读取其所需的滤波投影图像中的部分数据，有效减少了读取的数据量，从而节省总的重建时间。

下面结合附图和实施例对本发明进一步说明。

附图说明

图1是本发明技术流程图；

图2是三维图像重建中的Z-L线映射关系；

图3是重建空间与滤波投影图像的数据对应关系。

具体实施方式

在Intel Core II 2.33GHz处理器、2G内存的计算机上，实施本发明方法的步骤如下：

(1)对一零件进行圆周锥束CT扫描，采集一组360幅1536×1920的投影图像用于重建，其中dso＝780mm，dsd＝880mm，r＝80mm；

(2)在投影图像中选取一个边长为1024个象素的正方形区域，该正方形区域的中心位于整幅投影图像的中心，将该组所有投影图像按该正方形区域裁剪为边长为1024个象素的正方形图像，并将所有裁剪后的投影图像进行公知的对数运算，此时投影图像的象素灰度为单精度浮点型；

(3)对该组投影图像进行FDK算法中的滤波处理，滤波函数采用公知的S-L滤波器；

(4)设置需要重建的切片图像起止层号为第301-600层，这段连续的切片图像即为重建任务

(5)根据计算机当前的可用内存数对指定重建任务进行自适应分割，得到若干个相同大小的子重建任务，具体步骤如下：

1)获取计算机物理内存当前的可用数为1145M，考虑到重建过程中可能会有其它程序的内存需求以及操作系统运行的内存需求，为确保计算过程的稳定，将该可用内存数减去145M，即实际可用的最大内存数为1000M；

2)由重建任务的切片图像起止层号得到重建任务的总层数为N＝300，设实际可用的最大内存数为T，重建任务所需的内存数R＝4M×300＝1200M，则需要将重建任务分割为S个子重建任务，S＝int(R/T)+1＝int(1200/1000)+1＝2；

3)考虑到子重建任务计算量的均衡性，应当使各子重建任务重建的切片层数相同，则子重建任务的层数m＝int(N/S+0.5)＝int(300/2+0.5)＝150，考虑到后续的计算中会采用单指令多数据(SIMD)技术进行加速，因此需要将m取为m′＝152，m′为大于或等于m的整数中能被4整除的最小数；

4)计算各子重建任务的起止层号，得到第1个子重建任务为第301-452层，第2个子重建任务为第453-604层。

(6)根据子重建任务的层数m，从计算机中分配一个1024×1024×152×4/1024/1024＝608M大小的单精度浮点型内存空间，此内存空间称为部分重建空间；

(7)判断：是否还有未重建的子重建任务，若是，则继续执行下一步，若否，则转第9步；

(8)在部分重建空间中进行下一个子重建任务的重建计算，完成后转第7步，具体步骤如下：

1)选取部分重建空间的内切圆柱为感兴趣重建空间ROI，即每一层切片的重建区域只局限于切片图像的内切圆中，这样可将重建体素的数目减少至原始数目的π/4；

2)计算子重建任务的起止层切片图像在滤波投影图像中对应位置(图3)，分以下两种情况：

对于起始层切片图像，按以下公式计算：

对于终止层切片图像，按以下公式计算：

其中，h_s为起始层切片到中心切片的层数，h_e为终止层切片到中心切片的层数，r为子重建空间内切圆柱半径，dso为射线源焦点到工作台旋转中心的距离，dsd为射线源焦点到探测器的距离。将OZ_s取整得到P_s，OZ_e取整得到P_e，并设中心切片层号为P₀，并考虑一点冗余，则子重建任务所需的重建数据位于滤波投影图像中的第P₀±P_s-1行至第P₀±P_e+1行，当所计算层位于中心层之下时取-号，反之取+号。本例中P₀＝512，对于第1个子重建任务，P_s＝265，P_e＝75，则子重建任务所需的重建数据位于滤波投影图像中的第246行至第438行；

3)判断：是否还有未读入的滤波投影图像，若是，则继续进行下一步，若否，则转第6步；

4)读入下一幅滤波投影图像中与子重建任务相对应的部分数据；

5)按Z线优先重建算法并采用单指令多数据(SIMD)技术进行加速计算，完成后转第3步；

6)将子重建空间的重建结果存储为Z向的序列切片图像，结束。

(9)释放部分重建空间所占内存，完成重建任务。

表1给出了FDK算法、Z线优先方法和本发明方法在锥束CT三维图像重建的反投影部分的计算速度比较，可见本发明方法很好地实现了自适应三维图像快速重建计算，而且由于没有改变Z线优先算法本身的计算策略，所以重建图像质量与Z线优先算法保持一致。

表1二维图像重建的反投影部分的计算速度比较

	FDK算法	Z线优先方法	本发明方法
				时间(s)	24867.3	无法实施	451.34
加速比	1	/	55.1

Claims (3)

1、一种自适应的锥束CT三维图像快速重建方法，其特征在于包括下述步骤：

(1)对试件进行圆周锥束CT扫描，采集一组用于重建的投影图像；

(2)在投影图像中选取一个边长为E个象素的正方形区域，将该组所有投影图像按该正方形区域裁剪为边长为E个象素的正方形图像，并将所有裁剪后的投影图像进行对数运算，此时投影图像的象素灰度为单精度浮点型；

(3)对该组投影图像进行FDK算法中的滤波处理，滤波函数采用S-L滤波器、R-L滤波器或SL-W滤波器；

(4)设置需要重建的切片图像起止层号，这段连续的切片图像即为重建任务；

(5)根据计算机当前的可用内存数对重建任务进行自适应分割，得到若干个相同大小的子重建任务；

(6)根据子重建任务的层数m，从计算机中分配一个E×E×m大小的单精度浮点型内存空间，此内存空间称为部分重建空间；

(7)判断：是否还有未重建的子重建任务，若是，则继续执行下一步，若否，则转第(9)步；

(8)在部分重建空间中进行下一个子重建任务的重建计算，完成后转第7步；

(9)释放部分重建空间所占内存，完成重建任务。

2、根据权利要求1所述的一种自适应的锥束CT三维图像快速重建方法，其特征在于：所述的第(5)步中，根据计算机当前的可用内存数对重建任务进行自适应分割的具体步骤如下：

1)获取计算机物理内存当前的可用数，将该可用内存数减去50～150M，作为实际可用的最大内存数；

2)由重建任务的切片图像起止层号得到重建任务的总层数为N，设实际可用的最大内存数为T，重建任务所需的内存数R＝投影图像大小×N，则需要将重建任务分割为S个子重建任务，S＝int(R/T)+1；

3)使各子重建任务重建的切片层数相同，子重建任务的层数m＝int(N/S+0.5)，将m取为m′，m′为大于或等于m的整数中能被4整除的最小数；

4)计算各子重建任务的起止层号，设重建任务的起始层号为L_s，则第i个子重建任务的起始层号为L_s+(i-1)m，终止层号为L_s+i·m-1，其中i的取值为1～S。

3、根据权利要求1所述的一种自适应的锥束CT三维图像快速重建方法，其特征在于：所述的第(8)步具体步骤如下：

1)选取部分重建空间的内切圆柱为感兴趣重建空间ROI，即每一层切片的重建区域只局限于切片图像的内切圆中；

2)计算子重建任务的起止层切片图像在滤波投影图像中对应位置，分以下两种情况：

对于起始层切片图像，按以下公式计算：

对于终止层切片图像，按以下公式计算：

其中，h_s为起始层切片到中心切片的层数，h_e为终止层切片到中心切片的层数，r为子重建空间内切圆柱半径，dso为射线源焦点到工作台旋转中心的距离，dsd为射线源焦点到探测器的距离；将OZ_s取整得到P_s，OZ_e取整得到P_e，并设中心切片层号为P₀，子重建任务所需的重建数据位于滤波投影图像中的第P₀±P_s-1行至第P₀±P_e+1行，当所计算层位于中心层之下时取-号，反之取+号；

3)判断：是否还有未读入的滤波投影图像，若是，则继续进行下一步，若否，则转第6步；

4)读入下一幅滤波投影图像中与子重建任务相对应的部分数据；

5)按Z线优先重建算法并采用单指令多数据技术进行加速计算，完成后转第3)步；

6)将子重建空间的重建结果存储为Z向的序列切片图像，结束。

Images