CN101887591B - 基于矩形包围盒的锥束ct快速重建方法 - Google Patents
基于矩形包围盒的锥束ct快速重建方法 Download PDFInfo
- Publication number
- CN101887591B CN101887591B CN2010102101102A CN201010210110A CN101887591B CN 101887591 B CN101887591 B CN 101887591B CN 2010102101102 A CN2010102101102 A CN 2010102101102A CN 201010210110 A CN201010210110 A CN 201010210110A CN 101887591 B CN101887591 B CN 101887591B
- Authority
- CN
- China
- Prior art keywords
- image
- reconstruction
- rectangular bounding
- cutpoint
- projection
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Landscapes
- Apparatus For Radiation Diagnosis (AREA)
Abstract
本发明公开了一种基于矩形包围盒的锥束CT快速重建方法，对试件进行锥束CT圆周扫描采集一组投影图像，裁剪为边长为E个象素的一组正方形投影图像；计算重建空间中所述试件的矩形包围盒大小，对正方形投影图像分别计算其对数图像；对对数图像进行FDK算法中的滤波处理；根据矩形包围盒的大小，从计算机中分配一个单精度浮点型内存空间，按Z线优先重建算法并采用单指令多数据技术进行锥束CT重建计算；将重建结果按其坐标系存储为X向、Y向或Z向的序列切片图像。本发明有效提高了锥束CT的重建速度，降低了重建算法对内存大小的要求。
Description
技术领域
本发明属于CT***图像重建领域，涉及对锥束CT***中进行CT图像快速重建的一整套解决方案。
背景技术
计算机断层成像技术(Computed Tomography，CT)是通过对物体不同角度的射线投影重建而获取被测物体内部断层图像信息的成像技术。锥束CT利用锥形束X射线源和面阵探测器采集被测物体的投影数据，与传统二维CT相比，锥束CT一次扫描即可重建出数百甚至上千个断层图像，具有射线利用率高、切片连续、切片内和切片间空间分辨率相同、精度高等特点。
锥束CT应用中的一个关键问题就是提高图像重建速度。目前在商业领域中，应用最广泛的是FDK滤波反投影重建算法，该算法相比于其它重建算法具有较高的运算效率，但其反投影过程的计算复杂度仍然达到O(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，ROI)法是图像重建中一种非常实用的降低计算量的方法。使用某些先验知识来生成图像重建的ROI边界范围，如圆柱域或球形域等，仅重建位于ROI内的体素，而忽略那些位于ROI之外的体素。采用内切圆柱域ROI可将重建体素的数目减少至原数目的π/4(Z线优先重建算法即采用了这种方法)，而采用内切球形域ROI可减少至原数目的π/6。张顺利、张定华、赵歆波等人在《计算机辅助设计与图形学学报》(2009，21(2)：160-164)的文章“基于最小区域的快速CT图像重建”针对线阵平行束CT重建方式提出了一种基于最小区域的快速CT重建方法，在二维CT重建中进一步减少了重建的计算量，但该方法不适用于锥束CT。
发明内容
为了克服现有技术仍然存在大量非物体部分体素重建计算的不足，本发明提供一种基于矩形包围盒的锥束CT快速重建方法，在确保物体完整重建的前提下，进一步减少重建体素的数量，达到提高锥束CT重建速度的目的。
本发明解决其技术问题所采用的技术方案包括以下步骤：
(1)对试件进行锥束CT圆周扫描，采集一组投影图像，将该组所有投影图像按相同位置和大小裁剪为边长为E个象素的一组正方形投影图像，并确保该组正方形投影图像的边长E比试件投影的边长至少大20个象素；
(2)根据步骤(1)所得的一组正方形投影图像计算重建空间中所述试件的矩形包围盒尺寸，设其长、宽、高分别为L、W、H，单位为象素；
(3)对步骤(1)所得的一组正方形投影图像分别计算其对数图像，得到一组象素灰度为单精度浮点型的对数图像；
(4)对上一步所得的对数图像进行FDK算法中的滤波处理，滤波函数采用公知的S-L滤波器、R-L滤波器或SL-W滤波器；
(5)根据矩形包围盒的大小，从计算机中分配一个L×W×H大小的单精度浮点型内存空间；
(6)在该内存空间中按Z线优先重建算法并采用单指令多数据(SIMD)技术进行锥束CT重建计算；
(7)将重建结果按其坐标系存储为X向、Y向或Z向的序列切片图像，并释放重建所占内存空间。
在上述步骤(2)中，根据步骤(1)所得的一组正方形投影图像计算重建空间中该试件的矩形包围盒尺寸的具体步骤如下：
1)设该组投影图像的列方向为扫描轴方向(即Z向)，分别对各幅投影图像按列叠加成一行图像；
2)分别对上一步获取的每一行图像，从其左端点往右取连续的10-20个象素(该部分是不含试件投影的背景区域)，计算这些象素的灰度均值，并以该均值的0.95-0.98倍为阈值，由该行图像的左端点向右查找，所得到的第一个比该阈值小的象素即为该行图像的左分割点；该行图像右分割点的确定方法与上述方法相似，仅是其中的方向相反；
3)分别计算每一行图像的左分割点与右分割点的间隔象素个数＝右分割点位置-左分割点位置+1；
4)查找每一行图像的间隔象素个数的最小值，该最小值即为矩形包围盒的W，设该最小值所对应的投影图像的投影方位角度为A；
5)沿扫描旋转方向查找投影方位角度为A+90°的投影图像，其叠加成一行图像的左分割点与右分割点间隔的象素个数即为矩形包围盒的L；
6)矩形包围盒的H默认为整个重建空间的高度，并可指定为重建空间中任意感兴趣的Z向第S1层至第S2层的高度(S2＞S1)，此时H＝S2-S1+1；
7)按扫描旋转方向将该组投影图像重新排序，A+90°投影方位的投影图像为重排后的第一幅投影图像，此时矩形包围盒的各面分别平行于重建空间立方体的各面，便于后续的锥束CT重建计算。
本发明的有益效果是：由于通过确定试件锥束CT重建时的矩形包围盒，并结合Z线优先重建算法，本发明在大大减少重建体素数量的同时，也减少了重建所需的内存消耗，从而有效提高了锥束CT的重建速度，降低了重建算法对内存大小的要求。
下面结合附图和实施例对本发明进一步说明。
附图说明
图1为本发明方法的实现流程。
具体实施方式
利用主要部件为Varian PaxScan 2520探测器、YXLON 450kV X射线源的锥束CT***，在Intel Core II 2.33GHz处理器、2G内存的计算机上，实施本发明方法的步骤如下：
(1)对一试件进行锥束CT圆周扫描，采集一组360幅1536×1920的投影图像，将该组所有投影图像按中心位置和大小裁剪为边长为512个象素的一组正方形投影图像，此时可确保该组正方形投影图像比试件投影至少大20个象素；
(2)根据步骤(1)所得的一组正方形投影图像计算重建空间中所述试件的矩形包围盒尺寸，设其长、宽、高分别为L、W、H，单位为象素，具体步骤如下：
1)设该组投影图像的列方向为扫描轴方向(即Z向)，分别对各幅投影图像按列叠加成一行图像；
2)分别对上一步获取的每一行图像，从其左端点往右取连续的15个象素(该部分是不含试件投影的背景区域)，计算这些象素的灰度均值，并以该均值的0.96倍为阈值，由该行图像的左端点向右查找，所得到的第一个比该阈值小的象素即为该行图像的左分割点；该行图像右分割点的确定方法与上述方法相似，仅是其中的方向相反；
3)分别计算每一行图像的左分割点与右分割点的间隔象素个数＝右分割点位置-左分割点位置+1；
4)查找每一行图像的间隔象素个数，得到其最小值为192，该最小值即为矩形包围盒的W，该最小值所对应的投影图像的投影方位角度为128°；
5)沿扫描旋转方向查找投影方位角度为128°+90°的投影图像，其叠加成一行图像的左分割点与右分割点间隔的象素个数为355，即为矩形包围盒的L；
6)矩形包围盒的H默认为整个重建空间的高度512，这里指定为重建空间中感兴趣的Z向第45层至第444层的高度，此时H＝444-45+1＝400；
7)按扫描旋转方向将该组投影图像重新排序，128°+90°投影方位的投影图像为重排后的第一幅投影图像，此时矩形包围盒的各面分别平行于重建空间立方体的各面，便于后续的锥束CT重建计算。
(3)对步骤(1)所得的一组正方形投影图像分别计算其公知的对数图像，得到一组象素灰度为单精度浮点型的对数图像；
(4)对上一步所得的对数图像进行FDK算法中的滤波处理，滤波函数采用公知的S-L滤波器；
(5)根据矩形包围盒的大小，从计算机中分配一个L×W×H大小的单精度浮点型内存空间；
(6)在该内存空间中按Z线优先重建算法并采用单指令多数据(SIMD)技术进行锥束CT重建计算；
(7)将重建结果按其坐标系存储为Z向的序列切片图像，并释放重建所占内存空间。
对于该试件，表1给出了采用Z线优先算法(重建512层)和本发明方法(重建400层)在锥束CT图像重建的反投影部分的计算速度比较，可见本发明方法获得了可观的加速比。表2给出了Z线优先算法完全重建和本发明方法重建所需内存的对比，可见本发明方法重建对内存的需求大大降低。由于本发明方法没有改变Z线优先算法本身的计算策略，所以重建图像质量与Z线优先算法保持一致。
表1锥束CT图像重建的反投影计算速度比较
| Z线优先算法(重建512层) | 本发明方法(重建400层) | |
| 时间(s) | 195.76 | 52.72 |
| 加速比 | 1 | 3.71 |
表2重建所需内存对比
| Z线优先算法(重建512层) | 本发明方法(重建400层) | |
| 内存大小(M) | 512 | 104 |
| 内存占用比 | 1 | 4.92 |
Claims (1)
1.基于矩形包围盒的锥束CT快速重建方法，其特征在于包括下述步骤：
(1)对试件进行锥束CT圆周扫描，采集一组投影图像，将该组所有投影图像按相同位置和大小裁剪为边长为E个象素的一组正方形投影图像，并确保该组正方形投影图像的边长E比试件投影的边长至少大20个象素；
(2)根据步骤(1)所得的一组正方形投影图像计算重建空间中所述试件的矩形包围盒大小，设其长、宽、高分别为L、W、H个象素，具体步骤如下：
1)设该组投影图像的列方向为扫描轴方向，分别对各幅投影图像按列叠加成一行图像；
2)分别对上一步获取的每一行图像，从其左端点往右取连续的10-20个象素，这些象素是不含试件投影的背景区域，计算这些象素的灰度均值，并以该均值的0.95-0.98倍为阈值，由该行图像的左端点向右查找，所得到的第一个比该阈值小的象素即为该行图像的左分割点；该行图像右分割点的确定方法与上述方法相似，仅是其中的方向相反；
3)分别计算每一行图像的左分割点与右分割点的间隔象素个数＝右分割点位置-左分割点位置+1；
4)查找每一行图像的间隔象素个数的最小值，该最小值即为矩形包围盒的W，设该最小值所对应的投影图像的投影方位角度为A；
5)沿扫描旋转方向查找投影方位角度为A+90°的投影图像，其叠加成一行图像的左分割点与右分割点间隔的象素个数即为矩形包围盒的L；
6)矩形包围盒的H默认为整个重建空间的高度，并可指定为重建空间中任意感兴趣的Z向第S1层至第S2层的高度，S2＞S1，此时H＝S2-S1+1；
7)按扫描旋转方向将该组投影图像重新排序，A+90°投影方位的投影图像为重排后的第一幅投影图像；
(3)对步骤(1)所得的一组正方形投影图像分别计算其对数图像，得到一组象素灰度为单精度浮点型的对数图像；
(4)对上一步所得的对数图像进行FDK算法中的滤波处理，滤波函数采用S-L滤波器、R-L滤波器或SL-W滤波器；
(5)根据矩形包围盒的大小，从计算机中分配一个L×W×H大小的单精度浮点型内存空间；
(6)在该内存空间中按Z线优先重建算法并采用单指令多数据技术进行锥束CT重建计算；
(7)将重建结果按其坐标系存储为X向、Y向或Z向的序列切片图像，并释放重建所占内存空间。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN2010102101102A CN101887591B (zh) | 2010-06-24 | 2010-06-24 | 基于矩形包围盒的锥束ct快速重建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN2010102101102A CN101887591B (zh) | 2010-06-24 | 2010-06-24 | 基于矩形包围盒的锥束ct快速重建方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN101887591A CN101887591A (zh) | 2010-11-17 |
| CN101887591B true CN101887591B (zh) | 2011-09-14 |
Family
ID=43073501
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN2010102101102A Expired - Fee Related CN101887591B (zh) | 2010-06-24 | 2010-06-24 | 基于矩形包围盒的锥束ct快速重建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN101887591B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102243766A (zh) * | 2011-06-08 | 2011-11-16 | 无锡引速得信息科技有限公司 | 一种ct图像重建加速的新方法 |
| CN110000812A (zh) * | 2019-03-21 | 2019-07-12 | 深圳点猫科技有限公司 | 一种用于机器人的物体碰撞检测方法及其*** |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP3950782B2 (ja) * | 2002-05-22 | 2007-08-01 | ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー | ３次元逆投影方法およびｘ線ｃｔ装置 |
| US7170966B2 (en) * | 2003-08-05 | 2007-01-30 | Gioietta Kuo-Petravic | Practical implementation of a CT cone beam algorithm for 3-D image reconstruction as applied to nondestructive inspection of baggage, live laboratory animal and any solid materials |
| JP2005137390A (ja) * | 2003-11-04 | 2005-06-02 | Ge Medical Systems Global Technology Co Llc | Ｃｔ画像生成方法およびｘ線ｃｔ装置 |
| CN100464707C (zh) * | 2005-09-16 | 2009-03-04 | 北京大学 | 三维锥束ct图像重建的处理*** |
| CA2643931A1 (en) * | 2006-02-27 | 2007-08-30 | University Of Rochester | Method and apparatus for cone beam ct dynamic imaging |
-
2010
- 2010-06-24 CN CN2010102101102A patent/CN101887591B/zh not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| CN101887591A (zh) | 2010-11-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Noël et al. | GPU-based cone beam computed tomography | |
| Pratx et al. | Fast, accurate and shift-varying line projections for iterative reconstruction using the GPU | |
| CN101596113B (zh) | 一种ct并行重建***及成像方法 | |
| Abella et al. | Software architecture for multi-bed FDK-based reconstruction in X-ray CT scanners | |
| CN101625767B (zh) | 非平行轮廓医学影像的点云式三维表面重建方法 | |
| Liu et al. | GPU-based branchless distance-driven projection and backprojection | |
| US9111381B2 (en) | Shift-varying line projection using graphics hardware | |
| CN105118039B (zh) | 实现锥束ct图像重建的方法及*** | |
| CN104899903A (zh) | 基于cuda架构的gpu并行加速双能谱ct重建方法 | |
| CN102779350A (zh) | 一种锥束ct迭代重建算法投影矩阵构建方法 | |
| CN101520899B (zh) | 一种锥束ct三维图像的并行重建方法 | |
| Ma et al. | An encoder-decoder network for direct image reconstruction on sinograms of a long axial field of view PET | |
| Zhao et al. | GPU based iterative cone-beam CT reconstruction using empty space skipping technique | |
| Park et al. | A fully GPU-based ray-driven backprojector via a ray-culling scheme with voxel-level parallelization for cone-beam CT reconstruction | |
| CN101882319B (zh) | 基于最小三维凸包的锥束ct快速重建方法 | |
| CN101567090B (zh) | 一种自适应的锥束ct三维图像快速重建方法 | |
| CN101887591B (zh) | 基于矩形包围盒的锥束ct快速重建方法 | |
| Nagai et al. | Tomographic surface reconstruction from point cloud | |
| US7272205B2 (en) | Methods, apparatus, and software to facilitate computing the elements of a forward projection matrix | |
| Wang et al. | Sparse-view cone-beam CT reconstruction by bar-by-bar neural FDK algorithm | |
| Yan et al. | Fast Katsevich algorithm based on GPU for helical cone-beam computed tomography | |
| Noël et al. | Clinical evaluation of GPU-based cone beam computed tomography | |
| Lin et al. | Calibration method of center of rotation under the displaced detector scanning for industrial CT | |
| Pratx et al. | 3-D tomographic image reconstruction from randomly ordered lines with CUDA | |
| Wang et al. | Fast system matrix calculation based on voxel projection in cone-beam CT |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| C14 | Grant of patent or utility model | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20110914 Termination date: 20180624 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |