专利名称：一种ct并行重建系统及成像方法
技术领域：
本发明属于无损4企测技术领域，具体地说，本发明涉及一种CT并行重建成像系统及并行重建成像方法。
背景技术：
计算机断层成像术(CT)广泛应用于医学影像及工业无损探伤领域，其基本原理是通过采集不同角度下的被测物对射线的投影衰减信息，通过相应的滤波反投影算法(FBP),重建出被测物的二维或三维断层信息。在图像重建过程中，FBP算法需要先对原始数据进行滤波，然后计算相应重建坐标点对应的探测器像素单元，将对应角度下的探测器数据赋值给相应的坐标点，即为反投影。反投影过程涉及大量的三角及地址搜索运算，耗费大量时间。如对于用1800个投影角度信息重建1024 x 1024像素的断层图像，需要进行超过10,欠三角运算，直接对其进行图像重建计算，耗时超过l分钟。CT应用于过程工程领域的测量，4企测对象通常为运动的物体，只有毫秒级的准实时扫描及重建才能提供被测对象的真实信息。加速图像重建的一个方法是，事先将涉及三角运算的地址计算好存为地址文件供重建时查询，该方法可有效缩短重建所需时间，其代价是极大的增加了内存开销，如对于上述的1024 x1024 xl800重建，若地址为float型(4个字节)，其全部地址占用内存约为7.5G。地址查表法的另一个缺点是，当CT几何参凄丈发生变化时，需要全部重新计算投影地址。
另一方面，基于图形处理器(GPU)的高性能凝:值运算是近年来的一个研究热点。CPU通常使用半数以上的空间(如Inte^Xeoi^处理器60。/。的芯片空间)用于緩存及管理，GPU则将绝大多数的晶体管用于数值计算。除此之外，GPU的晶体管总数也远大于CPU，如AMD⑧的Athlon 64FX双核CPU处理器共有2. 27亿个晶体管，四核心Xeor^DPCPU处理器Harpertovra(45nm，代号Penryn)拥有8.2亿个晶体管，NVIDIA GeForce 9800GX2 ( G92核心)GPU图形处理器的晶体管数是15. 08亿。正因为GPU专注于数值运算，因此其计算能力远大于CPU。目前，峰值计算能力达500GFlops的GPU已经面市，约为16个处理器的Cray⑧C90 超级计算机峰值性能的100倍。(L Mueller, F. Xu， and N.Neophytou Why do GPUs work so well for acceleration of CT SPIEElectronic Imaging '07 San Jose, January 2007 )因此，GPU特别适用于并行度高数值运算量大的领域。
发明内容
本发明的目的是针对现有的CT图像重建系统的不足，提供一种基于图形处理器(GPU)集群的多进程多线程CT图像并行重建系统，并提供针对GPU集群特点设计的多进程多线程并行重建方法，使得图像重建时间大大缩短，可达到毫秒级准实时显示重建图像。
为实现上述发明目的，本发明提供的CT并行重建系统包括与CT探测器列阵连接的前端采样器，与所述前端采样器连接的中心节点，以及与所述中心节点连接的多个子节点，每个所述子节点均安装有图形处理器，所述中心节点用于将CT重建图像划分为多个block区域(即块区域)，并根据所划分的block区域对每个子节点的计算任务进行分配，所述子节点用于计算所分配到的block区域中各像素的重建值。
上述技术方案中，所述中心节点也安装有图形处理器。为实现上述发明目的，本发明提供的基于所述的CT并行重建系统的成像方法，包括如下步骤
1) 将重建图像划分为多个block区域；
2) 中心节点接收前端采样器所采集的各角度下的原始投影数据；
3) 中心节点为各节点分配计算任务，并将原始投影数据传输给分配了计算任务的各子节点；每个分配了计算任务的子节点计算一个或多个block区域的重建值；
4) 各子节点完成所分配的block区域内各像素的重建计算，将计算结果返回给中心节点；
5) 中心节点根据各节点的计算结果组合出完整的重建图像。上述技术方案中，所述步骤l)中，每个所述block区域含有的像素点
的数目是64的倍数。
上述技术方案中，所述步骤l)中，每个所述block区域的像素点组成方阵。
上述技术方案中，所述步骤l)中，每个所述block区域的像素点数目为256,所述block区i或为16*16方阵。
上述技术方案中，所述步骤2)中，所述前端采样器依次采集CT探测过程中各角度下的原始投影数据，所述中心节点分批次依次接收各角度下
的原始投影数据， 一个角度下的原始投影数据为一个批次；所述中心节点在接收到一个角度下的各探测器单元的原始投影数据后，根据当前各节点的计算负荷情况将计算任务动态分配给当前计算负荷较小的节点。
上述技术方案中，所述步骤3)中，当中心节点接收到一个角度下的原始投影数据时，将该角度下的重建图像的所有block区域的计算任务分配给同一节点进行计算；所述步骤4)中，各节点将各角度下的重建图像返回给中心节点；所述步骤5)中，中心节点将各角度下的重建图像组合成最终的重建图像。
上述技术方案中，所述步骤3)中，中心节点每接收到一个角度下的原始投影数据时，将该角度下的重建图像的各block区域的计算任务分配给不同的节点进行计算；所述步骤4)中，各节点将各自负责的block区域的重建图像返回给中心节点，所述步骤5)中，首先根据各节点返回的各block区域的重建图像组合出每个角度下的重建图像；然后再得出最终的重建图像。
上述技术方案中，所述步骤3)中，每个子节点负责该一个block区域的重建计算，每个子节点的图形处理器的一个线程负责一个像素的重建值计算。
与现有技术相比，本发明提供一种基于图形处理器(GPU)集群的多进程多线程CT图像并行重建系统，并且提供针对GPU集群硬件特性开发的多进程多线程的CT图像并行重建方法，采取划分重建图像子区域和线程块的方法，充分挖掘CT扫描及重建的并行特性，在数据采集的同时进行GPU的并行重建，本发明可将重建速度提高超过1000倍，重建时间由现有技术的分钟级降低为毫秒级，达到准实时的被测对象的断层图像重建显示。
以下，结合附图来详细说明本发明的实施例，其中图l是本发明的图像重建系统结构示意图，其中l为射线源，2为被测物，3为转台，4为^果测器。图2是本发明中的GPU工作原理示意图。
图3是本发明提供的基于GPU并行重建系统对三代CT扫描数据的FBP重建图，测量对象为有机玻璃管和楔形模型。
图4是单CPU标准FBP对三代CT扫描数据的重建图，测量对象为有机玻璃管和楔形模型。
图5是本发明提供的基于GPU并行重建系统对五代CT扫描数据的重建图，测量对象为有机玻璃圆盘模型。
图6是单CPU标准FBP对五代CT扫描数据的重建图，测量对象为有机玻璃圆盘模型。
具体实施例方式
实施例l
本实施例提供一种基于GPU集群的并行重建系统，包括一台前端数据釆集计算机，与CT系统的探测器连接；N台后端数据处理计算机(N不做特殊限定)，用Linux操作系统组成局域网，每台计算机上配有NVIDIA⑧公司生产的TeslaTM图形处理器；该系统还包括射线源、探测器阵列、转台等CT扫描所必需硬件。中心节点连接存储设备和图像输出设备。
本实施例中，每个GPU包括L个多处理器单元(multiprocessor )，每个多处理器单元同时可发起M个线程投入计算(不同型号的GPU的M数量不同)，因此，理论上，设计好针对GPU特点的CT图像重建计算方法，N个GPU可将计算速度提高L x N x M倍。
本实施例还提供了一种基于GPU集群的并行重建方法，包括以下步骤
1) 确定重建图像尺寸，并在各节点上进行相应的内存分配，分配的内存大小均为重建尺寸*数据类型，如重建图像尺寸为1024*1024像素，数据类型为float型(4个字节)，则各节点需划分的内存空间为1024*1024*4个字节，用于緩存各节点重建的图像。
2) 对重建计算任务进行预分配，每个节点分配一个投影角度数据的重建计算，各个节点的重建图像进一步划分为m*m个子图像区域，每个子图像区域由11*11个像素组成；将每个节点上的GPU分为n^m个线程块(block)，每个子图^象区域对应着一个GPU线程块，每个线程块划分为n*n个线程，每个线程计算一个像素点。
3) 将转台置于初始角度，转台有两种工作模式， 一种是射线源-探测器不动，转台带动被测物旋转，另一种是被测物不动，转台带动射线源-探测器绕被测物旋转，二者对图像重建不产生实质影响；
4) 测得当前角度下的投影原始数据，将投影原始数据传送至中心节点；中心节点侦测各子节点的计算负荷情况并对各子节点的计算任务进行分配；并将当前角度下的投影原始数据传送至被分配到计算任务的x号子节点；当前角度下的重建图像的计算任务可以分配给一个子节点，也可以分配给多个子节点共同完成，当中心节点也具有GPU时，中心节点本身也可以被分配计算任务；为了提高计算效率，本步骤中一个子节点的GPU负责计算一个block区域，GPU的每个线程计算该block区域的一个像素点的重建值。但值得注意的是以上任务分配方式并不是唯一的，比如每个线程也可以计算二个甚至多个像素点，这是本领域技术人员容易理解。
5) x号子节点将接收到的当前角度下的投影原始数据上传至本节点的GPU緩存，进行傅利叶变换对原始投影数据滤波；然后根据步骤2中划分的线程块区域，在GPU緩存中划分相应的子图像存储空间，调用编写好的GPU,利用与该子节点连接的图形处理器并行计算所指定block区域的重建数据值；所述图形处理器中的每个处理器单元负责所指定的block区域中的一个线程区域；
6) 调用驱动GPU工作的核函数，根据步骤2中划分好的线程块，并行发起重建计算线程，每个线程对应计算重建图像的一个像素点；
7) 侦测所有线程计算完成，将子图像数据合成，回传至中心节点，供最终重建图像的合成；
8) 判断是否遍历所有角度，如果判断为否，则按一定步长将转台置于下一个角度，返回步骤4);如果判断为是，将中心节点接收到的各节点处理后的重建图像进行合成，将其以图像形式输出至显示设备以及存储到硬盘、光盘、闪存盘等指定的存储设备上。
需说明的是，每个节点完成计算任务所需时间不同，因此在步骤4中，为了进一步提高计算效率，中心节点可以根据各子节点的负荷情况实时动态分配计算任务，选择当前时刻下计算负荷最小的x节点来进行计算，从而确保本并行重建系统达到计算速度的最优化。
本发明中，每个block内线程数分配不当，会导致GPU内存地址沖突。在采用Nvidia公司的Tesla C870 GPU时，当寄存器每个block内的线程数为64的倍数时，可最大程度的避免寄存器内存沖突。block内的线程数过少，线程进行读写及同步操作时，处理器中会有一部分寄存器处于空闲状态，计算性能不高，但随着block内的线程数的增加，每个线程可获得的寄存器数减少，其性能也会下降，因此，每个block内的线程数需有个最佳值，对于本实施例的重建系统，该最佳值为256,结合图像重建实际情况，将256个线程编为16*16方阵，分别对应重建图像的一部分。Block数过少，会导致部分multiprocessor(处理器单元)处于空闲状态， 一般block数至少为multiprocessor数的两倍，使得每个block进行读写及同步操作时，可以将multiprocessor切换至另 一个block进行计算。才艮据CT重建的特点(重建图像一般为正方形，其边长B通常为64、 128、 256、 512、 1024、2048个像素)，由于block内的线程一般编为16*16矩阵，对应的block划分也为n^m矩阵，m=B/16,对于1024*1024像素的重建图像，其block数为64*64。
本发明中对重建图像区域的划分是提高计算效率的重要步骤。每个block的线程数是64的倍数，以256与512是较佳值。同时，为了方便进行目标重建(英文为targeted reconstruction,是指针对重建图像某一特定部分进行，在实际应用中会经常用到这一点)，重建图像的每个block区域均设为方阵，因此以16*16的256线程的block为最佳。当然，block也可设为1*256，或其它的非方阵排列，但这样设置不容易进行目标重建。
本发明通过适当的划分，可实现用一个线程计算一个像素点，充分发挥了并行计算SIMD (单指令多数据)架构的优点，提高了重建运算的并行性。
实施例2
本实施实例是用三代单源等距扇束CT对有机玻璃管和楔形模型进行测量，探测器共有1536个像素点，像素点宽O. 4mm，步进电机驱动转台绕被测物进行360度全周扫描，共采集3600个角度投影数据。因此原始数据为1536*3600二维矩阵，以int类型存取，每个原始数据文件约为22M。重建图像image尺寸为1024"024。本发明的图像并行重建系统共有30个节点，每个节点安装有一块GPU图形处理器。
如图1所示，前端采样机控制转台的转动，以及射线源与探测器的同步工作，射线源1发射出的射线穿过被测物2后，由探测器3采集到衰减后的射线信号，传递至前端采样机，前端采样机经过模/数转换，将模拟电压线号转换成数字信号后，再实时地将数据传递至O号节点，O号节点(即中心节点)负责原始数据向其它29个节点的分发及重建后的数据的采集工作，最终重建图像素的合成也由零号节点完成，并将最终的图像数据以二进制编码进行存
储。根据本发明硬件系统的特点以及重建图像尺寸，将重建图像分为64*64个blocks,每个block由1646个线程组成，如图2所示。基于GPU的并行重建方法如下
1) 初始化图^象image，初始化后的image为零^i矩阵，前端采样才几读取第一个角度射线对板测物的衰减信息后(1536*1行向量，以P!表示)，将Pi通过以千兆以太网传至后端数据处理机群的O号节点；
2) O号节点侦测到数据传输完毕后，侦测当时各节点计算负荷情况，将P,发送至负荷最小的节点x。侦测方法是中心节点发送新计算任务前，先发出指令，统计每个节点完成计算的像素数，将新任务发到完成像素数最多的节点，即当前计算负荷最小的节点。
3 ) 节点x调用GPU使用FFT对Pi进行相应的滤波运算，此处滤波器可选标准的ramp滤波器，也可以选Ram-Lak、 cosine、 Hamming等滤波器； '
4) 节点x初始化该节点上的GPU，将PJi传至GPU的内存，调用数据处理核代码，每个GPU线程并行的计算其分辖的重建区域(即重建图像划分的blocks所对应的区域，每个线程对应一个像素点的计算)，其计算公式如下，d为重建图像中点(x, y)对应的投影探测器像素点，x、 y为图像的坐标，D为源到探测器的距离，本实例中为1000mm，a为像素点宽，本例中为O. 4mm, ^为投影角度，本步骤中为第一个角度，"=0;
」 Z) x (x x cos yff +少x sin々)
rf =-;-
(Z) + x x sin / — j; x cos / ) x a
5 ) 将步骤4中计算出的探测器像素点的值P!, d赋给重建图像上点
(x， y);
6 )对该GPU上所有线程进行同步操作，确保所有线程完成步骤4的计算
和步骤5的赋值，得102O1024矩阵ima {1};7) 节点x将第一个角度的重建矩阵imaU)传至O号节点，0号节点侦测到数据接收完毕，将ima {1}与image进行对应加和运算；8 ) 驱动步进电机旋转O. l度，按照步骤1-步骤7,采集下一个角度的投
影衰减数据，并对当前角度下投影数据进行重建；当完成所有3600个角度的数据采集及重建时，0号节点也完成imag^g/扁w运算，
得到重建图像的image矩阵；
9 ) 最后对重建图像进行平滑滤波、去除负值等后处理操作，将重建图
像的image矩阵发至图像显示机，显示最终重建图像。图3是使用本发明重建的二维图像，图4为使用单CPU+标准FBP算法的重建
结果。由于二者理论基础完全一致，因此成像质量是相同的。
使用本发明提供的重建系统和算法，每个角度的重建时间约为O. 5ms,完
成全部角度重建所需时间约为40ms。采用CPU单机重建约需160秒，采用30个
节点的CPU集群重建约需8秒，本发明可将重建效率分别提高约4000倍(相对
单CPU)和200倍(相对30个节点的CPU集群)。
实施例3
本实施实例是用五代等距扇束CT对有机玻璃圆盘模型进行测量，CT系统共有18组源-探测器，因此同时可采集18个投影角度的原始数据。每个探测器有640个像素，像素点宽O. 4隱，源与探测器成对安装于转台之上，数据通过滑环传至前端釆样机，驱动系统带动转台以667rpm的转速绕被测物进行360度全周扫描，共采集900个角度投影数据。因此原始数据为640*900二维矩阵。重建图像image尺寸为1024d024。本发明的图像并行重建系统共有30个节点，每个节点安装有一块GPU图形处理器，重建图像分为6"64个blocks,每个block由1646个线程(每个线程对应一个像素)组成，如图2所示。
基于GPU的并行重建方法如下
10) 初始化图^象image矩阵为零值矩阵，前端采样机同时读取第l、 51、101、…、851、 900个角度共18组投影数据，(640*1行向量,以P"b:l、
51、 101..... 851、 900}表示)，将Pb通过以千兆以太网传至后端
数据处理^L群的0号节点；
11) 0号节点侦测到凄t据传输完毕后，O号节点侦测当前各节点计算负荷情况，将Pb发送至18个负荷最小的子节点Gn;
12 ) 节点G。调用GPU使用FFT对Pb进行相应的滤波运算，此处滤波器可选才示准的ramp滤》皮器，也可以选Ram—Lak、 cosine、 Hamming等滤波器；
13) 初始化节点GJi的GPU，将PJi传至GPU的内存，调用数据处理核代
码，每个线程并行的计算其分辖的计算区域(即重建图像划分的
blocks中的线程所对应的区域)，其计算公式如下，d为重建图像中
'点(x， y)对应的投影探测器像素点，D为源到探测器的距离，本实例
中为1200腿，^为投影角度，"为像素点宽，本例中为O. 4mm;
j — Z) x (x x cos / + y x sin(Z) + x x sin / — x cosx "
14 )将步骤13中计算出的探测器像素点的值赋给点(x， y),即为Pb. d;
15) 对该GPU上所有线程进行同步操作，确保所有线程完成步骤14的计
算和步骤15的赋值，得18组102"1024矩阵Ib(b-l、 51、 101.....
851、 900};
16) 完成计算的子节点将重建矩阵Ib传至O号节点，O号节点侦测到数据接收完毕，将Ib与image进行对应加和运算；
17) 转台驱动系统带动源-探测器旋转O. 4度，按照步骤ll-步骤17，采集下一个角度下的18组投影衰减数据并对这些数据进行重建，在完成50个角度共9QQ组投影数据采集及重建后，0号节点也完成
image= 玄/ 运算，得到重建图像image矩阵；
18) 对重建图像进行平滑滤波、去除负值等后处理操作，将图像image发至图像显示机，显示最终重建图像。
图5是使用本发明重建的二维图像，图6为使用单CPU+标准FBP算法的重建结果。
由于该CT系统配有18组源-探测器，因此，转台仅需旋转1/18周即可完成数据采集，每个角度的投影采样时间约为O. lms，完成全部采样约需5ms。使用本基于GPU集群的图像并行重建系统，完成一个角度的重建约需O. 2ms,完全所有角度的图像重建约需10ms。由于数据采集与重建非串行完成，而是完成一个角度的18组投影采集后，在进行第个角度数据采集的同时进行图像重建，因此采集与重建是并行进行的，完成一个断面的扫描及重建，共耗时12ms
最后，以上所述的具体实施例，对本发明的目的、 一支术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实 施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的 任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
权利要求
1.一种CT并行重建系统，包括与CT探测器列阵连接的前端采样器，与所述前端采样器连接的中心节点，以及与所述中心节点连接的多个子节点，每个所述子节点均安装有图形处理器，所述中心节点用于将CT重建图像划分为多个块区域，并根据所划分的块区域对每个子节点的计算任务进行分配，所述子节点用于计算所分配到的块区域中各像素的重建值。
2. 根据权利要求1所述的CT并行重建系统，其特征在于，所述中心 节点也安装有图形处理器。
3. 基于权利要求1所述的CT并行重建系统的成像方法，包括如下步骤1) 将重建图像划分为多个块区域；2) 中心节点接收前端采样器所采集的各角度下的原始投影数据；3) 中心节点为各节点分配计算任务，并将原始投影数据传输给分配 了计算任务的各子节点；每个分配了计算任务的子节点计算一个或多个块 区域的重建值；4) 各子节点完成所分配的块区域内各像素的重建计算，将计算结果 返回给中心节点；5 )中心节点根据各节点的计算结果组合出完整的重建图像。
4. 根据权利要求3所述的成像方法，其特征在于，所述步骤l)中， 每个所述块区域含有的像素点的数目是64的倍数。
5. 根据权利要求3所述的成像方法，其特征在于，所述步骤l)中， 每个所述块区域的像素点组成方阵。
6. 根据权利要求4所述的成像方法，其特征在于，所述步骤l)中， 每个所述块区域的像素点数目为256，所述块区域为16*16方阵。
7. 根据权利要求3所述的成像方法，其特征在于，所述步骤2)中， 所述前端采样器依次采集CT探测过程中各角度下的原始投影数据，所述 中心节点分批次依次接收各角度下的原始投影数据， 一个角度下的原始投 影数据为一个批次；所述中心节点在接收到一个角度下的各探测器单元的 原始投影数据后，根据当前各节点的计算负荷情况将计算任务分配给当前计算负荷较小的节点。
8. 根据权利要求7所述的成像方法，其特征在于，所述步骤3)中，当中心节点接收到 一个角度下的原始投影数据时，将该角度下的重建图像的所有块区域的计算任务分配给同一节点进行计算；所述步骤4)中，各节 点将各角度下的重建图像返回给中心节点；所述步骤5)中，中心节点将 各角度下的重建图像组合成最终的重建图像。
9. 根据权利要求7所述的成像方法，其特征在于，所述步骤3)中， 中心节点每接收到 一个角度下的原始投影数据时，将该角度下的重建图像 的各块区域的计算任务分配给不同的节点进行计算；所述步骤4)中，各节 点将各自负责的块区域的重建图像返回给中心节点，所述步骤5)中，首先 根据各节点返回的各块区域的重建图像组合出每个角度下的重建图像；然 后再得出最终的重建图像。
10. 根据权利要求9所述的成像方法，其特征在于，所述步骤3)中， 每个子节点负责该一个块区域的重建计算，每个子节点的图形处理器的一 个线程负责一个像素的重建值计算。
全文摘要
本发明提供一种CT并行重建系统及其成像方法，所述系统包括前端采样器、中心节点、以及与所述中心节点连接的多个子节点，每个所述子节点均安装有图形处理器。所述成像方法包括如下步骤1)将重建图像划分为多个block区域；2)中心节点接收前端采样器所采集的各角度下的原始投影数据；3)中心节点为各节点分配计算任务，每个分配了计算任务的子节点计算一个或多个block区域的重建值；4)各子节点完成所分配的重建计算；5)中心节点组合出完整的重建图像。本发明采取划分重建图像子区域的方法充分挖掘CT扫描及重建的并行特性，在数据采集的同时进行GPU的并行重建，重建时间由现有技术的分钟级降低为毫秒级，达到准实时的被测对象的断层图像重建显示。
文档编号G06T5/00GK101596113SQ200810114478
公开日2009年12月9日 申请日期2008年6月6日 优先权日2008年6月6日
发明者孟凡勇, 李静海, 维 王, 陈飞国 申请人:中国科学院过程工程研究所