专利名称：基于最小三维凸包的锥束ct快速重建方法
技术领域：
本发明属于CT系统图像重建领域，涉及对锥束CT系统中基于最小三维凸包进行 CT图像快速重建的方法。
背景技术：
计算机断层成像技术(Computed Tomography, CT)是通过对物体不同角度的射线 投影重建而获取被测物体内部断层图像信息的成像技术。锥束CT利用锥形束X射线源和 面阵探测器采集被测物体的投影数据，与传统二维CT相比，锥束CT 一次扫描即可重建出数 百甚至上千个断层图像，具有射线利用率高、切片连续、切片内和切片间空间分辨率相同、 精度高等特点。锥束CT应用中的一个关键问题就是提高图像重建速度。目前在商业领域中，应用 最广泛的是FDK滤波反投影重建算法，该算法相比于其它重建算法具有较高的运算效率， 但其反投影过程的计算复杂度仍然达到0(N4)，其中N为投影数据的尺寸，反投影运行时间 占整个重建过程的98%以上。因此，要想提高FDK算法的重建速度，关键在于反投影过程的 优化。从目前的文献资料来看，FDK快速重建主要分为以下两类一是FDK重建算法的改 进及代码优化，该类方法主要通过一些近似处理以减少反投影过程的运算量，可能引入一 些新的重建误差，如P-FDK、S-FDK、T-FDK等FDK衍生算法；二是采用并行计算技术，如采用 工作站配以专门的阵列处理器实现并行计算。毛海鹏、张定华、梁亮等人在《系统仿真学报》 (2004,16(11) =2486-2489)的文章“一种基于PC的快速三维图像重建方法”中对FDK重建 算法进行改进并结合数据并行处理提出了 Z线优先重建算法，在PC机上实现了三维图像快 速重建。传统的锥束CT图像重建方法是将重建空间限定在一个立方体中，针对重建空间 中的每个体素进行反投影计算，但有的体素并没有被检测物体所覆盖，从而导致不必要的 计算量。感兴趣区域(Region of Interest, R0I)法是图像重建中一种非常实用的降低计 算量的方法。使用某些先验知识来生成图像重建的ROI边界范围，如圆柱域或球形域等，仅 重建位于ROI内的体素，而忽略那些位于ROI之外的体素。采用内切圆柱域ROI可将重建 体素的数目减少至原数目的η /4 (Ζ线优先重建算法即采用了这种方法)，而采用内切球形 域ROI可减少至原数目的π/6。张顺利、张定华、赵歆波等人在《计算机辅助设计与图形学 学报》(2009，21⑵160-164)的文章“基于最小区域的快速CT图像重建”针对线阵平行束 CT重建方式提出了一种基于最小区域的快速CT重建方法，在二维CT重建中进一步减少了 重建的计算量，但该方法不适用于锥束CT。
发明内容
为了克服现有技术仍然存在大量非物体部分体素重建计算的不足，本发明提供一 种基于最小三维凸包的锥束CT快速重建方法，在确保物体完整重建的前提下，进一步减少重建体素的数量，达到提高锥束CT重建速度的目的。
本发明解决其技术问题所采用的技术方案包括以下步骤(1)对试件进行锥束CT圆周扫描，采集一组投影图像，将该组所有投影图像按相 同位置和大小裁剪为边长为E个象素的一组正方形投影图像，并确保该组正方形投影图像 的边长E比试件投影边长至少大20个象素；(2)根据步骤(1)所得的一组正方形投影图像计算重建空间中该试件的最小三维 凸包参数；(3)对步骤(1)所得的一组正方形投影图像分别计算其公知的对数图像，得到一 组象素灰度为单精度浮点型的对数图像；(4)对上一步所得的对数图像进行FDK算法中的滤波处理，滤波函数采用公知的 S-L滤波器、R-L滤波器或SL-W滤波器；(5)将最小三维凸包体素化，分配重建内存空间，再按Z线优先重建算法并采用单 指令多数据(SIMD)技术重建最小三维凸包内的体素；(6)将重建结果按其坐标系存储为X向、Y向或Z向的序列切片图像，并释放重建 所占内存空间。在上述步骤⑵中，将最小三维凸包定义为包含试件的最小柱状体，且其Z向截面 为相同的二维凸包。确定最小三维凸包，就是要确定其3个基本参数截面形状、高度和在 重建空间中的位置，具体步骤如下1)设该组投影图像的列方向为扫描轴方向(即Z向)，分别对各幅投影图像按列 叠加成一行图像(假设每幅图像的高为M,宽为N，H幅图像分别按列叠加，“一行图像”是形 成了 H个IXN的图像)；2)对上一步获取的每一行图像，分别计算其试件投影区域的左分割点P1*右分割 点P2，计算步骤如下①生成一个与叠加后的一行图像象素个数相等的单精度浮点型一维数组，并将该 数组全部置1，该数组元素序列与该行图像象素序列相对应，用于存储对应象素的标识值；②对该行图像除左右两端η个象素外的每个象素分别取其左右等量的η个象素， 设其左边象素灰度之和为IV，右边象素灰度之和为τκ，则该象素的标识值为ιγ/τκ并存入数 组中与该像素对应的位置；η取3 7 ；③在数组中从左向右查找第一个局部极大值礼，准则为礼大于1.01且同时大于 或等于其左右的各3 5个值，然后取左分割阈值& = (1+MJ/2，从礼向左查找得到的第 一个小于&的值所对应的象素就是左分割点P1 ；④在数组中从右向左查找第一个局部极小值Μκ，准则为Mr小于0. 99且同时小于 或等于其左右的各3 5个值，然后取右分割阈值Sk = (1+MK)/2，WMK向右查找得到的第 一个大于Sk的值所对应的象素就是右分割点P2 ；⑤分别将P1往左移3 5个象素，P2往右移3 5个象素，以增强算法可靠性。3)以锥束CT重建空间的旋转中心为原点、以0°扫描时的中心射束为Y轴、以 90°扫描时的中心射束为X轴建立直角坐标系，射线源和探测器绕原点旋转，设射线源到 旋转中心的距离为Dso，探测器到旋转中心的距离为Dod。以射线源在0°扫描位置成像为 例，射线源R的坐标为(0，Dso)，P1的坐标为(-0' P1, -Dod)，P2的坐标为(0' P2, -Dod)，分别计算射线RP1和RP2与重建空间Z向切片正方形内切圆的交点(共4个)，相应两个交 点连线即得到该位置处的两条包络线段。对于其余的每次成像(设其扫描角度为θ )，先将 其放在0°扫描位置并按上述方法计算出4个交点，然后将这4个交点分别绕原点旋转θ 角，再把相应两个交点连线即得到该扫描角度处的两条包络线段。最后由总数为投影图像 数量2倍的包络线段所围成的中间图形就是最小三维凸包的截面形状，是一个二维凸包；4)确定最小三维凸包的高度任取一幅投影图像，将其按行叠加成一列，从该列 的上端点往下取连续的10 20个象素(该部分是不含试件投影的背景区域)，计算这些象 素的灰度均值，并以该均值的0. 95 0. 98倍为分割阈值，由该列图像的上端点向下查找， 所得到的第一个比分割阈值小的象素位置即为最小三维凸包的上限位置Zmax ；由该列图像 的下端点向上查找，所得到的第一个比分割阈值小的象素位置即为最小三维凸包的下限位 置Zmin ；则(Zmax_Zmin+l)即为最小三维凸包的高度；为便于Z线优先的SIMD计算，需通过Zmin 或Zmax的增大或减小使最小三维凸包的高度规整为与原值最接近的4的整数倍；5)最小三维凸包在重建空间中的位置，由上述步骤3)中最小三维凸包的截面形 状在XOY平面中的位置和上述步骤4)中的最小三维凸包的上限位置Zmax以及下限位置Zmin 共同确定。在上述步骤(5)中，考虑到Z线优先重建算法的特点，不必对整个最小三维凸包进行体素化，只须对其截面的二维凸包进行象素化即可，具体步骤为1)生成一个与重建空间Z向切片图像大小相同的标识切片图像，所有象素初值置 0 ；2)在标识切片图像中采用Bresenham直线生成方法将上述步骤⑵中的包络线段 象素化，包络线段通过的象素值置1，包络线段所围成的中间图形就是最小三维凸包的截面 二维凸包；3)确定一个位于二维凸包内的种子点取0°方位投影图像的P1P2的中点横坐标值 为种子点的横坐标值，取90°方位投影图像的P1P2的中点纵坐标值为种子点的纵坐标值；4)采用扫描线种子填充算法在标识切片图像中的截面二维凸包内进行区域填充， 将填充后的象素值置2;5)查找标识切片图像中灰度值为2的象素的最小X坐标、最大X坐标、最小Y坐标 和最大Y坐标，由其确定二维凸包的外接矩形；6)分配重建内存空间，大小为外接矩形长X外接矩形宽X三维凸包高度，类型为 单精度浮点型，并全部置0;7)对标识切片图像中象素值为2位置处的Z线(其长度为最小三维凸包的高度)， 按Z线优先重建算法并采用单指令多数据(SIMD)技术重建其上的体素。本发明的有益效果是通过确定试件锥束CT重建时的最小三维凸包，并结合Z线 优先重建算法，在大大减少重建体素数量的同时，也减少了重建所需的内存消耗，从而有效 提高了锥束CT的重建速度，降低了重建算法对内存大小的要求。下面结合附图和实施例对本发明进一步说明。
图1是锥束CT扫描时试件V及其最小三维凸包V'的投影成像关系图2是最小三维凸包的XOY截面成像图。
具体实施例方式利用主要部件为Varian PaxScan 2520探测器、YXLON 450kV X射线源的锥束CT 系统，在Intel Core II 2. 33GHz处理器、2G内存的计算机上，实施本发明方法的步骤如 下(1)对一试件进行锥束CT圆周扫描，采集一组360幅1536X 1920的投影图像，将 该组所有投影图像按中心位置和大小裁剪为边长为512个象素的一组正方形投影图像，此 时可确保该组正方形投影图像比试件投影至少大20个象素；(2)根据步骤(1)所得的一组正方形投影图像计算重建空间中该试件的最小三维 凸包参数，具体步骤如下1)设该组投影图像的列方向为扫描轴方向(即Z向)，分别对各幅投影图像按列 叠加成一行图像；2)对上一步获取的每一行图像，分别计算其试件投影区域的左分割点P1*右分割 点P2，计算步骤如下①生成一个与叠加后的一行图像象素个数相等的单精度浮点型一维数组并全部 置1，该数组元素序列与该行图像象素序列相对应，用于存储对应象素的标识值；②对该行图像的每个象素(左右两端的5个象素除外)，分别取其左右等量的5个 象素，设左边象素灰度之和为IV，右边象素灰度之和为τκ，则该象素的标识值为ιγ/τκ并存入 数组对应位置；③在数组中从左向右查找第一个局部极大值Ml，准则为ML大于1.01且同时大于 或等于其左右的各3个值，然后取左分割阈值& = (Ι+Μ /2，从礼向左查找得到的第一个 小于&的值所对应的象素就是左分割点P1 ；④在数组中从右向左查找第一个局部极小值Μκ，准则为Mr小于0. 99且同时小于 或等于其左右的各3个值，然后取右分割阈值Sk = (1+MK)/2，WMK向右查找得到的第一个 大于Sk的值所对应的象素就是右分割点P2 ；⑤分别将P1往左移3个象素，P2往右移3个象素，以增强算法可靠性。3)以锥束CT重建空间的旋转中心为原点、以0°扫描时的中心射束为Y轴、以 90°扫描时的中心射束为X轴建立直角坐标系，射线源和探测器绕原点旋转，设射线源到 旋转中心的距离为Dso，探测器到旋转中心的距离为Dod。以射线源在0°扫描位置成像为 例，射线源R的坐标为(0，Dso)，P1的坐标为(-0' P1, -Dod)，P2的坐标为(0' P2, -Dod)， 分别计算射线RP1和RP2与重建空间Z向切片正方形内切圆的交点(共4个)，相应两个交 点连线即得到该位置处的两条包络线段。对于其余的每次成像(设其扫描角度为θ )，先将 其放在0°扫描位置并按上述方法计算出4个交点，然后将这4个交点分别绕原点旋转θ 角，再把相应两个交点连线即得到该扫描角度处的两条包络线段。最后由总数为投影图像 数量2倍的包络线段所围成的中间图形就是最小三维凸包的截面形状，是一个二维凸包；4)确定最小三维凸包的高度任取一幅投影图像，将其按行叠加成一列，从该列 的上端点往下取连续的15个象素(该部分是不含试件投影的背景区域)，计算这些象素 的灰度均值，并以该均值的0. 96倍为分割阈值，由该列图像的上端点向下查找，所得到的第一个比分割阈值小的象素位置即为最小三维凸包的上限位置Zmax = 462 ；由该列图像的 下端点向上查找，所得到的第一个比分割阈值小的象素位置即为最小三维凸包的下限位置 Zmin = 0(因试件下方有夹具支撑)；则(Zmax-Zmin+1)即为最小三维凸包的高度463，为便于 Z线优先的SIMD计算，将Zmax增大为463使其规整为与原值最接近的4的整数倍464 ；5)最小三维凸包在重建空间中的位置，由上述步骤3)中最小三维凸包的截面形 状在XOY平面中的位置和上述步骤4)中的最小三维凸包的上限位置Zmax以及下限位置Zmin 共同确定。
(3)对步骤(1)所得的一组正方形投影图像分别计算其公知的对数图像，得到一 组象素灰度为单精度浮点型的对数图像；(4)对上一步所得的对数图像进行FDK算法中的滤波处理，滤波函数采用公知的 S-L滤波器；(5)将最小三维凸包体素化，分配重建内存空间，再按Z线优先重建算法并采用单 指令多数据(SIMD)技术重建最小三维凸包内的体素，具体步骤如下1)生成一个与重建空间Z向切片图像大小相同的标识切片图像，所有象素初值置 0 ；2)在标识切片图像中采用Bresenham直线生成方法将上述步骤⑵中的包络线段 象素化，包络线段通过的象素值置1，包络线段所围成的中间图形就是最小三维凸包的截面 二维凸包；3)确定一个位于二维凸包内的种子点取0°方位投影图像的P1P2的中点横坐标 值为种子点的横坐标值，取90°方位投影图像的P1P2的中点纵坐标值为种子点的纵坐标 值；4)采用扫描线种子填充算法在标识切片图像中的二维凸包内进行区域填充，将填 充后的象素值置2;5)查找标识切片图像中灰度值为2的象素，得到其最小X坐标为90、最大X坐标 为407、最小Y坐标为142、最大Y坐标为425，由其确定二维凸包的外接矩形；6)分配重建内存空间，大小为外接矩形长(407-90+1 = 318) X外接矩 形宽(425-142+1 = 284) X三维凸包高度464，为其分配相应的重建内存空间为 318X 284X464X4/ (1024*1024) = 159. 8M ；类型为单精度浮点型，并全部置0 ；7)对标识切片图像中象素值为2位置处的Z线(其长度为最小三维凸包的高度)， 按Z线优先重建算法并采用单指令多数据(SIMD)技术重建其上的体素。(6)将重建结果按其坐标系存储为Z向的序列切片图像，并释放重建所占内存空 间。对于该试件，表1给出了采用Z线优先算法(重建512层)和本发明方法(重建 464层)在锥束CT图像重建的反投影部分的计算速度比较，可见本发明方法获得了可观的 加速比。表2给出了 Z线优先算法完全重建和本发明方法重建所需内存的对比，可见本发 明方法重建对内存的需求大大降低。由于本发明方法没有改变Z线优先算法本身的计算策 略，所以重建图像质量与Z线优先算法保持一致。表1锥束CT图像重建的反投影计算速度比较Z线优先算法 本发明方法 (重建512层)(重建464层)
时间(s) 195. 7656.41
加速比 3 47表2重建所需内存对比
Z线优先算法 本发明方法 (重建512层) (重建464层)
内存大小
512159.8
(M)
内存占用比 权利要求
基于最小三维凸包的锥束CT快速重建方法，其特征在于包括下述步骤(1)对试件进行锥束CT圆周扫描，采集一组投影图像，将该组所有投影图像按相同位置和大小裁剪为边长为E个象素的一组正方形投影图像，并确保该组正方形投影图像的边长E比试件投影边长至少大20个象素；(2)根据步骤(1)所得的一组正方形投影图像计算重建空间中该试件的最小三维凸包参数；(3)对步骤(1)所得的一组正方形投影图像分别计算其公知的对数图像，得到一组象素灰度为单精度浮点型的对数图像；(4)对上一步所得的对数图像进行FDK算法中的滤波处理，滤波函数采用公知的S-L滤波器、R-L滤波器或SL-W滤波器；(5)将最小三维凸包体素化，分配重建内存空间，再按Z线优先重建算法并采用单指令多数据技术重建最小三维凸包内的体素；(6)将重建结果按其坐标系存储为X向、Y向或Z向的序列切片图像，并释放重建所占内存空间。
2.根据权利要求1所述的基于最小三维凸包的锥束CT快速重建方法，其特征在于所 述的最小三维凸包为包含试件的最小柱状体，且其Z向截面为相同的二维凸包；确定最小 三维凸包的截面形状、高度和在重建空间中的位置具体步骤如下1)设该组投影图像的列方向为扫描轴方向，分别对各幅投影图像按列叠加成一行图像；2)对上一步获取的每一行图像，分别计算其试件投影区域的左分割点P1和右分割点 P2，计算步骤如下①生成一个与叠加后的一行图像象素个数相等的单精度浮点型一维数组，并将该数组 全部置1，该数组元素序列与该行图像象素序列相对应，用于存储对应象素的标识值；②对该行图像除左右两端η个象素外的每个象素分别取其左右等量的η个象素，设其 左边象素灰度之和为IV，右边象素灰度之和为Tk，则该象素的标识值为ιγ/τκ并存入数组中 与该像素对应的位置；η取3 7 ；③在数组中从左向右查找第一个局部极大值My准则为Ml大于1.01且同时大于或等 于其左右的各3 5个值，然后取左分割阈值& = (1+MJ/2，从礼向左查找得到的第一个 小于&的值所对应的象素就是左分割点P1 ；④在数组中从右向左查找第一个局部极小值Mr,准则为Mk小于0.99且同时小于或等 于其左右的各3 5个值，然后取右分割阈值Sk = (1+MK)/2，WMK向右查找得到的第一个 大于Sk的值所对应的象素就是右分割点P2 ；⑤分别将P1往左移3 5个象素，P2往右移3 5个象素；3)以锥束CT重建空间的旋转中心为原点、以0°扫描时的中心射束为Y轴、以90°扫 描时的中心射束为X轴建立直角坐标系，射线源和探测器绕原点旋转，设射线源到旋转中 心的距离为Dso，探测器到旋转中心的距离为Dod ；射线源在0°扫描位置成像时，射线源R 的坐标为(0，Dso)，P1的坐标为(-0' P1, -Dod)，P2的坐标为(0' P2, -Dod)，分别计算射线 RP1和RP2与重建空间Z向切片正方形内切圆的交点，相应两个交点连线即得到该位置处的 两条包络线段；对于扫描角度为θ的每次成像，先将其放在0°扫描位置并按上述方法计算出4个交点，然后将这4个交点分别绕原点旋转θ角，再把相应两个交点连线即得到该 扫描角度处的两条包络线段；最后由总数为投影图像数量2倍的包络线段所围成的中间图 形就是最小三维凸包的截面形状，是一个二维凸包； 4)确定最小三维凸包的高度；任取一幅投影图像，将其按行叠加成一列，从该列的上 端点往下取连续的10 20个象素的不含试件投影的背景区域，计算这些象素的灰度均 值，并以该均值的0. 95 0. 98倍为分割阈值，由该列图像的上端点向下查找，所得到的第 一个比分割阈值小的象素位置即为最小三维凸包的上限位置Zmax ；由该列图像的下端点向 上查找，所得到的第一个比分割阈值小的象素位置即为最小三维凸包的下限位置Zmin ；则 (Zmax_Zmin+l)即为最小三维凸包的高度；为便于Z线优先的SIMD计算，需通过Zmin或Zmax的 增大或减小使最小三维凸包的高度规整为与原值最接近的4的整数倍；5)最小三维凸包在重建空间中的位置，由上述步骤3)中最小三维凸包的截面形状在 XOY平面中的位置和上述步骤4)中的最小三维凸包的上限位置Zmax以及下限位置Zmin共同 确定。
3.根据权利要求1所述的基于最小三维凸包的锥束CT快速重建方法，其特征在于所 述的步骤(5)只须对最小三维凸包截面的二维凸包进行象素化即可，具体步骤为1)生成一个与重建空间Z向切片图像大小相同的标识切片图像，所有象素初值置0；2)在标识切片图像中采用Bresenham直线生成方法将上述步骤(2)中的包络线段象素 化，包络线段通过的象素值置1，包络线段所围成的中间图形就是最小三维凸包的截面二维 凸包；3)确定一个位于二维凸包内的种子点取0°方位投影图像的P1P2的中点横坐标值为 种子点的横坐标值，取90°方位投影图像的P1P2的中点纵坐标值为种子点的纵坐标值；4)采用扫描线种子填充算法在标识切片图像中的截面二维凸包内进行区域填充，将填 充后的象素值置2;5)查找标识切片图像中灰度值为2的象素的最小X坐标、最大X坐标、最小Y坐标和最 大Y坐标，由其确定二维凸包的外接矩形；6)分配重建内存空间，大小为外接矩形长X外接矩形宽X三维凸包高度，类型为单精 度浮点型，并全部置0;7)对标识切片图像中象素值为2位置处的Z线，按Z线优先重建算法并采用单指令多 数据技术重建其上的体素。
全文摘要
本发明公开了一种基于最小三维凸包的锥束CT快速重建方法，对试件进行锥束CT圆周扫描采集一组投影图像，将该组所有投影图像裁剪为边长为E个象素的一组正方形投影图像；计算重建空间中该试件的最小三维凸包参数；对正方形投影图像计算得到一组象素灰度为单精度浮点型的对数图像；对对数图像进行FDK算法中的滤波处理，将最小三维凸包体素化，分配重建内存空间，再按Z线优先重建算法并采用单指令多数据技术重建最小三维凸包内的体素；将重建结果按其坐标系存储为X向、Y向或Z向的序列切片图像，并释放重建所占内存空间。本发明有效提高了锥束CT的重建速度，降低了重建算法对内存大小的要求。
文档编号G06T11/00GK101882319SQ20101021030
公开日2010年11月10日 申请日期2010年6月24日 优先权日2010年6月24日
发明者卜昆, 张定华, 李明君, 程云勇, 黄魁东 申请人:西北工业大学