CN101520899B - 一种锥束ct三维图像的并行重建方法 - Google Patents
一种锥束ct三维图像的并行重建方法 Download PDFInfo
- Publication number
- CN101520899B CN101520899B CN2009100219130A CN200910021913A CN101520899B CN 101520899 B CN101520899 B CN 101520899B CN 2009100219130 A CN2009100219130 A CN 2009100219130A CN 200910021913 A CN200910021913 A CN 200910021913A CN 101520899 B CN101520899 B CN 101520899B
- Authority
- CN
- China
- Prior art keywords
- reconstruction
- space
- rebuild
- image
- son
- 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
- Image Processing (AREA)
Abstract
本发明公开了一种锥束CT三维图像的并行重建方法,对试件进行圆周锥束CT扫描,采集一组用于重建的投影图像;将该组所有投影图像按该正方形区域裁剪为边长为E个象素的正方形图像,并将投影图像进行对数运算;对该组投影图像进行滤波处理;从计算机中一次分配E3大小的单精度浮点型内存空间作为重建空间;对重建空间分割得若干个大小基本相同的子重建空间;对各子重建空间采用SIMD加速的Z线优先重建算法进行重建;根据需要将重建空间存储为公知的X向、Y向或Z向的序列切片图像。本发明解决了Z线优先重建算法不能直接在多核计算机中实行并行计算的问题。
Description
技术领域
本发明属于CT系统图像重建领域,涉及锥束CT系统并行化三维图像重建的方法。
背景技术
CT(Computed Tomography)技术利用射线(一般为X射线)透射被检测物体,并在探测器上获取一组投影图像,再配合相应的重建算法得到物体的切片图像。目前研究和应用的CT可分为二维CT和三维CT两大类,锥束CT属于三维CT。与传统二维CT相比,锥束CT扫描速度快、射线利用率高,能获得均匀、高精度的空间分辨率。锥束CT的重建是三维图像重建,即需要重建得到连续的很多层切片图像,涉及巨大的计算量。在锥束CT实际应用中,在保证良好的重建质量的基础上,提高重建速度成为必须突破的瓶颈。
目前在商用锥束CT系统中,应用最广泛的重建算法是FDK滤波反投影算法,该算法是一种近似重建,计算量相对较小,在小锥角的情况下可以获得质量良好的重建图像。即便如此,FDK算法的计算量依然较大,主要集中在反投影过程,计算复杂度为O(N4),其中N为投影数据的尺寸。例如,用360幅5122投影图像重建5123图像,在普通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线优先重建算法本身已经进行了SIMD并行处理,而一般并行计算的线程通信与数据保护技术很难应用于SIMD计算过程。
随着CPU技术的发展,目前双核CPU和四核CPU已经普及,而六核CPU、八核CPU也将很快得到广泛应用。因此,研究并实现Z线优先重建算法在多核计算机中的并行计算方法,是一个具有应用基础和实际价值的项目。
发明内容
为了克服现有技术不能通过简单增加线程的方法进行并行计算的不足,本发明提供一种锥束CT三维图像的并行重建方法,使得Z线优先重建算法能够直接在多核计算机中实行并行计算。
本发明解决其技术问题所采用的技术方案包括以下步骤:
(1)对试件进行圆周锥束CT扫描,采集一组用于重建的投影图像;
(2)在投影图像中选取一个边长为E个象素的正方形区域,将该组所有投影图像按该正方形区域裁剪为边长为E个象素的正方形图像,并将所有裁剪后的投影图像进行公知的对数运算,此时投影图像的象素灰度为单精度浮点型;
(3)对该组投影图像进行FDK算法中的滤波处理,滤波函数采用公知的S-L滤波器、R-L滤波器或SL-W滤波器;
(4)从计算机中一次分配E3大小的单精度浮点型内存空间,此内存空间称为重建空间,是一个边长为E个体素的立方体;
(5)对重建空间进行分割,得到若干个大小基本相同的子重建空间;
(6)对各个子重建空间同时采用SIMD加速的Z线优先重建算法进行重建,当所有子重建空间完成其重建计算时,重建结束;
(7)根据需要将重建空间存储为公知的X向、Y向或Z向的序列切片图像。
在上述方法第5步中,为减少不必要的重建计算,同时达到避免线程通信和实现负载平衡的目的,进行重建空间分割的具体步骤如下(图2):
1)选取立方体重建空间的内切圆柱为感兴趣重建空间ROI,即每一层切片的重建区域只局限于切片图像的内切圆中,这样可将重建体素的数目减少至原始数目的π/4;
2)设置:2≤子重建空间的个数≤计算机中所有CPU的总核心数;
3)对感兴趣重建空间ROI按子重建空间个数进行扇形等分,对位于等分线上的重建空间体素,将其归入包含该体素较多部分的子重建空间,这样得到的各子重建空间的大小基本相同,且各子重建空间之间无交叉重叠部分。
在上述方法第6步中,所有子重建空间同时采用SIMD加速的Z线优先重建算法进行重建的具体步骤如下:
1)分别为每一个子重建空间设置一个整数标志位,0表示空闲或已完成重建计算,1表示正在进行重建计算,并设置所有标志位的初始值为0;
2)判断:是否还有未读入的滤波处理后的投影图像,若是,则继续进行下一步,若否,则全部重建结束;
3)读入下一幅滤波处理后的投影图像;
4)设置所有标志位的值为1,同时从该滤波图像中读取所需的数据进行各子重建空间的重建,由于对滤波图像仅进行读操作,因此不会出现子重建空间重建任务间的排斥和滤波图像数据被改写的情况。另外,由于操作系统的原因,实际中不可能所有子重建空间同时完成各自的重建计算,因此当某个子重建空间完成其重建计算时,将其标志位设置为0;
5)判断:直到所有标志位的值均为0,转第2步。
本发明的有益效果是:在保持Z线优先重建算法原有的重建图像质量基础上,通过对重建空间进行简便实用的分割,同时达到避免线程通信和实现负载平衡的目的,解决了Z线优先重建算法不能直接在多核计算机中实行并行计算的问题。
下面结合附图和实施例对本发明进一步说明。
附图说明
图1是本发明方法流程图。
图2是三维图像重建中的Z-L线映射关系及重建空间按线程分割方式。
具体实施方式
在Intel Xeon 2.5GHz四核处理器、4G内存的计算机上,实施本发明方法的步骤如下:
(1)对一阶梯铝柱零件进行圆周锥束CT扫描,采集一组360幅768×960的投影图像用于重建;
(2)在投影图像中选取一个边长为512个象素的正方形区域,该正方形区域的中心位于整幅投影图像的中心,将该组所有投影图像按该正方形区域裁剪为边长为512个象素的正方形图像,并将所有裁剪后的投影图像进行公知的对数运算,此时投影图像的象素灰度为单精度浮点型;
(3)对该组投影图像进行FDK算法中的滤波处理,滤波函数采用公知的S-L滤波器;
(4)从计算机中一次分配5123大小的单精度浮点型内存空间,此内存空间称为重建空间,是一个边长为512个体素的立方体;
(5)对重建空间进行分割,得到若干个大小基本相同的子重建空间,为减少不必要的重建计算,同时达到避免线程通信和实现负载平衡的目的,进行重建空间分割的具体步骤为:
1)选取立方体重建空间的内切圆柱为感兴趣重建空间ROI,即每一层切片的重建区域只局限于切片图像的内切圆中,这样可将重建体素的数目减少至原始数目的π/4;
2)设置子重建空间的个数=计算机中所有CPU的总核心数=4;
3)对感兴趣重建空间ROI按子重建空间个数进行扇形4等分,对位于等分线上的重建空间体素,将其归入包含该体素较多部分的子重建空间,这样得到的各子重建空间的大小基本相同,且各子重建空间之间无交叉重叠部分。
(6)对各个子重建空间同时采用SIMD加速的Z线优先重建算法进行重建,当所有子重建空间完成其重建计算时,重建结束,具体步骤如下:
1)分别为每一个子重建空间设置一个整数标志位,0表示空闲或已完成重建计算,1表示正在进行重建计算,并设置所有标志位的初始值为0;
2)判断:是否还有未读入的滤波处理后的投影图像,若是,则继续进行下一步,若否,则全部重建结束;
3)读入下一幅滤波处理后的投影图像;
4)设置所有标志位的值为1,同时从该滤波图像中读取所需的数据进行各子重建空间的重建,由于对滤波图像仅进行读操作,因此不会出现子重建空间重建任务间的排斥和滤波图像数据被改写的情况。另外,由于操作系统的原因,实际中不可能所有子重建空间同时完成各自的重建计算,因此当某个子重建空间完成其重建计算时,将其标志位设置为0;
5)判断:直到所有标志位的值均为0,转第2步。
(7)根据需要将重建空间存储为公知的Z向序列切片图像。
针对上述裁剪后的投影图像,表1给出了FDK算法、Z线优先算法和本发明方法在锥束CT三维图像重建的反投影部分的计算速度比较,可见本发明方法很好地实现了Z线优先算法的并行加速计算,而且由于没有改变Z线优先算法本身的计算策略,所以重建图像质量与Z线优先算法保持一致。
表1三维图像重建的反投影部分的计算速度比较
FDK算法 | Z线优先算法 | 本发明方法(四线程) | |
时间(s) | 6636.219 | 135.274 | 37.42 |
加速比 | 1 | 49.06 | 177.34 |
Claims (1)
1.一种锥束CT三维图像的并行重建方法,其特征在于包括下述步骤:
(1)对试件进行圆周锥束CT扫描,采集一组用于重建的投影图像;
(2)在投影图像中选取一个边长为E个象素的正方形区域,将该组所有投影图像按该正方形区域裁剪为边长为E个象素的正方形图像,并将所有裁剪后的投影图像进行对数运算,此时投影图像的象素灰度为单精度浮点型;
(3)对该组投影图像进行FDK算法中的滤波处理,滤波函数采用S-L滤波器、R-L滤波器或SL-W滤波器;
(4)从计算机中一次分配E3大小的单精度浮点型内存空间,此内存空间称为重建空间,是一个边长为E个体素的立方体;
(5)对重建空间进行分割,得到若干个大小基本相同的子重建空间,具体包括以下步骤:
1)选取立方体重建空间的内切圆柱为感兴趣重建空间感兴趣区域,即每一层切片的重建区域只局限于切片图像的内切圆中;
2)设置:2≤子重建空间的个数≤计算机中所有CPU的总核心数;
3)对感兴趣重建空间感兴趣区域按子重建空间个数进行扇形等分,对位于等分线上的重建空间体素,将其归入包含该体素较多部分的子重建空间;
(6)对各个子重建空间同时采用单指令多数据技术加速的Z线优先重建算法进行重建,当所有子重建空间完成其重建计算时,重建结束;重建步骤如下:
1)分别为每一个子重建空间设置一个整数标志位,0表示空闲或已完成重建计算,1表示正在进行重建计算,并设置所有标志位的初始值为0;
2)判断:是否还有未读入的滤波处理后的投影图像,若是,则继续进行下一步,若否,则全部重建结束;
3)读入下一幅滤波处理后的投影图像;
4)设置所有标志位的值为1,同时从该滤波图像中读取所需的数据进行各子重建空间的重建,当某个子重建空间完成其重建计算时,将其标志位设置为0;
5)判断:直到所有标志位的值均为0,转第(6)步的第2)步;
(7)根据需要将重建空间存储为公知的X向、Y向或Z向的序列切片图像。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN2009100219130A CN101520899B (zh) | 2009-04-08 | 2009-04-08 | 一种锥束ct三维图像的并行重建方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN2009100219130A CN101520899B (zh) | 2009-04-08 | 2009-04-08 | 一种锥束ct三维图像的并行重建方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN101520899A CN101520899A (zh) | 2009-09-02 |
CN101520899B true CN101520899B (zh) | 2011-11-16 |
Family
ID=41081472
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN2009100219130A Expired - Fee Related CN101520899B (zh) | 2009-04-08 | 2009-04-08 | 一种锥束ct三维图像的并行重建方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN101520899B (zh) |
Families Citing this family (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101819679B (zh) * | 2010-04-19 | 2011-12-07 | 李楚雅 | 三维医学图像分割方法 |
CN101882319B (zh) * | 2010-06-24 | 2011-11-30 | 西北工业大学 | 基于最小三维凸包的锥束ct快速重建方法 |
CN102013116B (zh) * | 2010-11-30 | 2012-11-07 | 华润万东医疗装备股份有限公司 | 基于脑血管旋转造影术进行三维重建的方法 |
US20120159124A1 (en) * | 2010-12-15 | 2012-06-21 | Chevron U.S.A. Inc. | Method and system for computational acceleration of seismic data processing |
EP3242586A1 (en) * | 2015-01-07 | 2017-11-15 | C/o Canon Kabushiki Kaisha | Photoacoustic apparatus, image display method, and program |
CN106228601B (zh) * | 2016-07-21 | 2019-08-06 | 山东大学 | 基于小波变换的多尺度锥束ct图像快速三维重建方法 |
CN112085811B (zh) * | 2020-09-23 | 2021-04-23 | 赛诺威盛科技(北京)有限公司 | Ct局部重建的方法及装置 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6014419A (en) * | 1997-11-07 | 2000-01-11 | Hu; Hui | CT cone beam scanner with fast and complete data acquistion and accurate and efficient regional reconstruction |
CN1739455A (zh) * | 2005-09-16 | 2006-03-01 | 北京大学 | 三维锥束ct图像重建的处理系统及处理方法 |
WO2006131872A3 (en) * | 2005-06-07 | 2008-03-13 | Philips Intellectual Property | Fast reconstruction algorithm for cone-beam ct |
-
2009
- 2009-04-08 CN CN2009100219130A patent/CN101520899B/zh not_active Expired - Fee Related
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6014419A (en) * | 1997-11-07 | 2000-01-11 | Hu; Hui | CT cone beam scanner with fast and complete data acquistion and accurate and efficient regional reconstruction |
WO2006131872A3 (en) * | 2005-06-07 | 2008-03-13 | Philips Intellectual Property | Fast reconstruction algorithm for cone-beam ct |
CN1739455A (zh) * | 2005-09-16 | 2006-03-01 | 北京大学 | 三维锥束ct图像重建的处理系统及处理方法 |
CN100464707C (zh) * | 2005-09-16 | 2009-03-04 | 北京大学 | 三维锥束ct图像重建的处理系统 |
Non-Patent Citations (1)
Title |
---|
毛海鹏等.一种基于 PC 的快速三维图像重建方法.《系统 仿 真 学 报》.2004,第16卷(第11期), * |
Also Published As
Publication number | Publication date |
---|---|
CN101520899A (zh) | 2009-09-02 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN101520899B (zh) | 一种锥束ct三维图像的并行重建方法 | |
DE102013204062B4 (de) | Voll-Parallel-am-Platz-Konstruktion von 3D-Beschleunigungs-Strukturen in einer Grafik-Verarbeitungs-Einheit | |
Pratx et al. | Fast, accurate and shift-varying line projections for iterative reconstruction using the GPU | |
Hsu | Segmented ray casting for data parallel volume rendering | |
CN109754361A (zh) | 3d各向异性的混合网络:将来自2d图像的卷积特征传递到3d各向异性体积 | |
AU2005275463B2 (en) | System and method for object characterization of toboggan-based clusters | |
US10127710B2 (en) | Processor and method for accelerating ray casting | |
DE102018114286A1 (de) | Durchführen einer Traversierungs-Stack-Komprimierung | |
DE102020107080A1 (de) | Grafiksysteme und Verfahren zum Beschleunigen von Synchronisation mittels feinkörniger Abhängigkeitsprüfung und Planungsoptimierungen basierend auf verfügbarem gemeinsam genutztem Speicherplatz | |
Kainz et al. | Ray casting of multiple volumetric datasets with polyhedral boundaries on manycore GPUs | |
DE102020132272A1 (de) | Verfahren und vorrichtung zum codieren basierend auf schattierungsraten | |
Sudbø et al. | New algorithms based on the Voronoi diagram applied in a pilot study on normal mucosa and carcinomas | |
DE102007020879A1 (de) | Verfahren und Vorrichtung für die äußerst schnelle Symmetrie- und SIMD- gestützte Projektion/Rückprojektion für die 3D-PET-Bildrekonstruktion | |
Jiang et al. | Acceleration of CT reconstruction for wheat tiller inspection based on adaptive minimum enclosing rectangle | |
Buurlage et al. | A geometric partitioning method for distributed tomographic reconstruction | |
CN101882319B (zh) | 基于最小三维凸包的锥束ct快速重建方法 | |
CN101561937A (zh) | 基于普通微机的大数据量医学影像三维交互方法 | |
CN102314688B (zh) | 图像分割方法和图像集分割方法 | |
DE102021114013A1 (de) | Techniken zum effizienten abtasten eines bildes | |
DE102013021044A1 (de) | Erzeugung fehlerbefreiter Voxel-Daten | |
US10964094B1 (en) | Visualization system that transforms 2D images of objects slices into 3D point clouds | |
Shin et al. | Acceleration techniques for cubic interpolation MIP volume rendering | |
CN103745495A (zh) | 基于医学体数据的体绘制方法 | |
CN105844690A (zh) | 基于gpu的快速三维ct迭代重建系统 | |
CN101887591B (zh) | 基于矩形包围盒的锥束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: 20111116 Termination date: 20150408 |
|
EXPY | Termination of patent right or utility model |