CN101596113B - 一种ct并行重建***及成像方法 - Google Patents

Abstract

本发明提供一种CT并行重建***及其成像方法，所述***包括前端采样器、中心节点、以及与所述中心节点连接的多个子节点，每个所述子节点均安装有图形处理器。所述成像方法包括如下步骤：1)将重建图像划分为多个block区域；2)中心节点接收前端采样器所采集的各角度下的原始投影数据；3)中心节点为各节点分配计算任务，每个分配了计算任务的子节点计算一个或多个block区域的重建值；4)各子节点完成所分配的重建计算；5)中心节点组合出完整的重建图像。本发明采取划分重建图像子区域的方法充分挖掘CT扫描及重建的并行特性，在数据采集的同时进行GPU的并行重建，重建时间由现有技术的分钟级降低为毫秒级，达到准实时的被测对象的断层图像重建显示。

Description

一种CT并行重建***及成像方法

技术领域

本发明属于无损检测技术领域，具体地说，本发明涉及一种CT并行重建成像***及并行重建成像方法。

背景技术

计算机断层成像术(CT)广泛应用于医学影像及工业无损探伤领域，其基本原理是通过采集不同角度下的被测物对射线的投影衰减信息，通过相应的滤波反投影算法(FBP)，重建出被测物的二维或三维断层信息。在图像重建过程中，FBP算法需要先对原始数据进行滤波，然后计算相应重建坐标点对应的探测器像素单元，将对应角度下的探测器数据赋值给相应的坐标点，即为反投影。反投影过程涉及大量的三角及地址搜索运算，耗费大量时间。如对于用1800个投影角度信息重建1024×1024像素的断层图像，需要进行超过10¹⁰次三角运算，直接对其进行图像重建计算，耗时超过1分钟。CT应用于过程工程领域的测量，检测对象通常为运动的物体，只有毫秒级的准实时扫描及重建才能提供被测对象的真实信息。加速图像重建的一个方法是，事先将涉及三角运算的地址计算好存为地址文件供重建时查询，该方法可有效缩短重建所需时间，其代价是极大的增加了内存开销，如对于上述的1024×1024×1800重建，若地址为float型(4个字节)，其全部地址占用内存约为7.5G。地址查表法的另一个缺点是，当CT几何参数发生变化时，需要全部重新计算投影地址。

另一方面，基于图形处理器(GPU)的高性能数值运算是近年来的一个研究热点。CPU通常使用半数以上的空间(如

处理器60％的芯片空间)用于缓存及管理，GPU则将绝大多数的晶体管用于数值计算。除此之外，GPU的晶体管总数也远大于CPU，如

的Athlon 64FX双核CPU处理器共有2.27亿个晶体管，四核心

DPCPU处理器Harpertown(45nm，代号Penryn)拥有8.2亿个晶体管，

GeForce 9800GX2(G92核心)GPU图形处理器的晶体管数是15.08亿。正因为GPU专注于数值运算，因此其计算能力远大于CPU。目前，峰值计算能力达500GFlops的GPU已经面市，约为16个处理器的C90超级计算机峰值性能的100倍。(K.Mueller，F.Xu，and N.Neophytou Why do GPUs work so well for acceleration of CT？SPIEElectronic Imaging′07San Jose，January 2007)因此，GPU特别适用于并行度高数值运算量大的领域。

发明内容

本发明的目的是针对现有的CT图像重建***的不足，提供一种基于图形处理器(GPU)集群的多进程多线程CT图像并行重建***，并提供针对GPU集群特点设计的多进程多线程并行重建方法，使得图像重建时间大大缩短，可达到毫秒级准实时显示重建图像。

为实现上述发明目的，本发明提供的CT并行重建***包括与CT探测器列阵连接的前端采样器，与所述前端采样器连接的中心节点，以及与所述中心节点连接的多个子节点，每个所述子节点均安装有图形处理器，所述中心节点用于将CT重建图像划分为多个block区域(即块区域)，并根据所划分的block区域对每个子节点的计算任务进行分配，所述子节点用于计算所分配到的block区域中各像素的重建值。

上述技术方案中，所述中心节点也安装有图形处理器。

为实现上述发明目的，本发明提供的基于所述的CT并行重建***的成像方法，包括如下步骤：

1)将重建图像划分为多个block区域；

2)中心节点接收前端采样器所采集的各角度下的原始投影数据；

3)中心节点为各节点分配计算任务，并将原始投影数据传输给分配了计算任务的各子节点；每个分配了计算任务的子节点计算一个或多个block区域的重建值；

4)各子节点完成所分配的block区域内各像素的重建计算，将计算结果返回给中心节点；

5)中心节点根据各节点的计算结果组合出完整的重建图像。

上述技术方案中，所述步骤1)中，每个所述block区域含有的像素点的数目是64的倍数。

上述技术方案中，所述步骤1)中，每个所述block区域的像素点组成方阵。

上述技术方案中，所述步骤1)中，每个所述block区域的像素点数目为256，所述block区域为16*16方阵。

上述技术方案中，所述步骤2)中，所述前端采样器依次采集CT探测过程中各角度下的原始投影数据，所述中心节点分批次依次接收各角度下的原始投影数据，一个角度下的原始投影数据为一个批次；所述中心节点在接收到一个角度下的各探测器单元的原始投影数据后，根据当前各节点的计算负荷情况将计算任务动态分配给当前计算负荷较小的节点。

上述技术方案中，所述步骤3)中，当中心节点接收到一个角度下的原始投影数据时，将该角度下的重建图像的所有block区域的计算任务分配给同一节点进行计算；所述步骤4)中，各节点将各角度下的重建图像返回给中心节点；所述步骤5)中，中心节点将各角度下的重建图像组合成最终的重建图像。

上述技术方案中，所述步骤3)中，中心节点每接收到一个角度下的原始投影数据时，将该角度下的重建图像的各block区域的计算任务分配给不同的节点进行计算；所述步骤4)中，各节点将各自负责的block区域的重建图像返回给中心节点，所述步骤5)中，首先根据各节点返回的各block区域的重建图像组合出每个角度下的重建图像；然后再得出最终的重建图像。

上述技术方案中，所述步骤3)中，每个子节点负责该一个block区域的重建计算，每个子节点的图形处理器的一个线程负责一个像素的重建值计算。

与现有技术相比，本发明提供一种基于图形处理器(GPU)集群的多进程多线程CT图像并行重建***，并且提供针对GPU集群硬件特性开发的多进程多线程的CT图像并行重建方法，采取划分重建图像子区域和线程块的方法，充分挖掘CT扫描及重建的并行特性，在数据采集的同时进行GPU的并行重建，本发明可将重建速度提高超过1000倍，重建时间由现有技术的分钟级降低为毫秒级，达到准实时的被测对象的断层图像重建显示。

附图说明

以下，结合附图来详细说明本发明的实施例，其中：

图1是本发明的图像重建***结构示意图，其中1为射线源，2为被测物，3为转台，4为探测器。

图2是本发明中的GPU工作原理示意图。

图3是本发明提供的基于GPU并行重建***对三代CT扫描数据的FBP重建图，测量对象为有机玻璃管和楔形模型。

图4是单CPU标准FBP对三代CT扫描数据的重建图，测量对象为有机玻璃管和楔形模型。

图5是本发明提供的基于GPU并行重建***对五代CT扫描数据的重建图，测量对象为有机玻璃圆盘模型。

图6是单CPU标准FBP对五代CT扫描数据的重建图，测量对象为有机玻璃圆盘模型。

具体实施方式

实施例1

本实施例提供一种基于GPU集群的并行重建***，包括一台前端数据采集计算机，与CT***的探测器连接；N台后端数据处理计算机(N不做特殊限定)，用Linux操作***组成局域网，每台计算机上配有

公司生产的Tesla^TM图形处理器；该***还包括射线源、探测器阵列、转台等CT扫描所必需硬件。中心节点连接存储设备和图像输出设备。

本实施例中，每个GPU包括L个多处理器单元(multiprocessor)，每个多处理器单元同时可发起M个线程投入计算(不同型号的GPU的M数量不同)，因此，理论上，设计好针对GPU特点的CT图像重建计算方法，N个GPU可将计算速度提高L×N×M倍。

本实施例还提供了一种基于GPU集群的并行重建方法，包括以下步骤：

1)确定重建图像尺寸，并在各节点上进行相应的内存分配，分配的内存大小均为重建尺寸^*数据类型，如重建图像尺寸为1024*1024像素，数据类型为float型(4个字节)，则各节点需划分的内存空间为1024*1024*4个字节，用于缓存各节点重建的图像。

2)对重建计算任务进行预分配，每个节点分配一个投影角度数据的重建计算，各个节点的重建图像进一步划分为m*m个子图像区域，每个子图像区域由n*n个像素组成；将每个节点上的GPU分为m*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划分也为m*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个像素点，像素点宽0.4mm，步进电机驱动转台绕被测物进行360度全周扫描，共采集3600个角度投影数据。因此原始数据为1536*3600二维矩阵，以int类型存取，每个原始数据文件约为22M。重建图像image尺寸为1024*1024。本发明的图像并行重建***共有30个节点，每个节点安装有一块GPU图形处理器。

如图1所示，前端采样机控制转台的转动，以及射线源与探测器的同步工作，射线源1发射出的射线穿过被测物2后，由探测器3采集到衰减后的射线信号，传递至前端采样机，前端采样机经过模/数转换，将模拟电压线号转换成数字信号后，再实时地将数据传递至0号节点，0号节点(即中心节点)负责原始数据向其它29个节点的分发及重建后的数据的采集工作，最终重建图像素的合成也由零号节点完成，并将最终的图像数据以二进制编码进行存储。根据本发明硬件***的特点以及重建图像尺寸，将重建图像分为64*64个blocks，每个block由16*16个线程组成，如图2所示。

基于GPU的并行重建方法如下：

1)初始化图像image，初始化后的image为零值矩阵，前端采样机读取第一个角度射线对被测物的衰减信息后(1536*1行向量，以P₁表示)，将P₁通过以千兆以太网传至后端数据处理机群的0号节点；

2)0号节点侦测到数据传输完毕后，侦测当时各节点计算负荷情况，将P₁发送至负荷最小的节点x。侦测方法是：中心节点发送新计算任务前，先发出指令，统计每个节点完成计算的像素数，将新任务发到完成像素数最多的节点，即当前计算负荷最小的节点。

3)节点x调用GPU使用FFT对P₁进行相应的滤波运算，此处滤波器可选标准的ramp滤波器，也可以选Ram-Lak、cosine、Hamming等滤波器；

4)节点x初始化该节点上的GPU，将P₁上传至GPU的内存，调用数据处理核代码，每个GPU线程并行的计算其分辖的重建区域(即重建图像划分的blocks所对应的区域，每个线程对应一个像素点的计算)，其计算公式如下，d为重建图像中点(x，y)对应的投影探测器像素点，x、y为图像的坐标，D为源到探测器的距离，本实例中为1000mm，a为像素点宽，本例中为0.4mm，β为投影角度，本步骤中为第一个角度，β＝0；

$d=\frac{D\×(x\×\cos \mathrm{\β}+y\×\sin \mathrm{\β})}{(D+x\×\sin \mathrm{\β}-y\×\cos \mathrm{\β})\×a}$

5)将步骤4中计算出的探测器像素点的值P_1，d赋给重建图像上点(x，y)；

6)对该GPU上所有线程进行同步操作，确保所有线程完成步骤4的计算和步骤5的赋值，得1024*1024矩阵ima{1}；

7)节点x将第一个角度的重建矩阵ima{1}传至0号节点，0号节点侦测到数据接收完毕，将ima{1}与image进行对应加和运算；

8)驱动步进电机旋转0.1度，按照步骤1-步骤7，采集下一个角度的投影衰减数据，并对当前角度下投影数据进行重建；当完成所有3600个角度的数据采集及重建时，0号节点也完成 $\mathrm{image}=\underset{n=1}{\overset{3600}{\mathrm{\Σ}}}\mathrm{ima}\left\{n\right\}$ 运算，得到重建图像的image矩阵；

9)最后对重建图像进行平滑滤波、去除负值等后处理操作，将重建图像的image矩阵发至图像显示机，显示最终重建图像。

图3是使用本发明重建的二维图像，图4为使用单CPU+标准FBP算法的重建结果。由于二者理论基础完全一致，因此成像质量是相同的。

使用本发明提供的重建***和算法，每个角度的重建时间约为0.5ms，完成全部角度重建所需时间约为40ms。采用CPU单机重建约需160秒，采用30个节点的CPU集群重建约需8秒，本发明可将重建效率分别提高约4000倍(相对单CPU)和200倍(相对30个节点的CPU集群)。

实施例3

本实施实例是用五代等距扇束CT对有机玻璃圆盘模型进行测量，CT***共有18组源-探测器，因此同时可采集18个投影角度的原始数据。每个探测器有640个像素，像素点宽0.4mm，源与探测器成对安装于转台之上，数据通过滑环传至前端采样机，驱动***带动转台以667rpm的转速绕被测物进行360度全周扫描，共采集900个角度投影数据。因此原始数据为640*900二维矩阵。重建图像image尺寸为1024*1024。本发明的图像并行重建***共有30个节点，每个节点安装有一块GPU图形处理器，重建图像分为64*64个blocks，每个block由16*16个线程(每个线程对应一个像素)组成，如图2所示。

基于GPU的并行重建方法如下：

10)初始化图像image矩阵为零值矩阵，前端采样机同时读取第1、51、101、...、851、900个角度共18组投影数据，(640*1行向量，以P_b{b＝1、51、101、...、851、900}表示)，将P_b通过以千兆以太网传至后端数据处理机群的0号节点；

11)0号节点侦测到数据传输完毕后，0号节点侦测当前各节点计算负荷情况，将P_b发送至18个负荷最小的子节点G_n；

12)节点G_n调用GPU使用FFT对P_b进行相应的滤波运算，此处滤波器可选标准的ramp滤波器，也可以选Ram-Lak、cosine、Hamming等滤波器；

13)初始化节点G_n上的GPU，将P₁上传至GPU的内存，调用数据处理核代码，每个线程并行的计算其分辖的计算区域(即重建图像划分的blocks中的线程所对应的区域)，其计算公式如下，d为重建图像中点(x，y)对应的投影探测器像素点，D为源到探测器的距离，本实例中为1200mm，β为投影角度，a为像素点宽，本例中为0.4mm；

$d=\frac{D\×(x\×\cos \mathrm{\β}+y\×\sin \mathrm{\β})}{(D+x\×\sin \mathrm{\β}-y\×\cos \mathrm{\β})\×a}$

14)将步骤13中计算出的探测器像素点的值赋给点(x，y)，即为P_b，d；

15)对该GPU上所有线程进行同步操作，确保所有线程完成步骤14的计算和步骤15的赋值，得18组1024*1024矩阵I_b{b＝1、51、101、...、851、900}；

16)完成计算的子节点将重建矩阵I_b传至0号节点，0号节点侦测到数据接收完毕，将I_b与image进行对应加和运算；

17)转台驱动***带动源-探测器旋转0.4度，按照步骤11-步骤17，采集下一个角度下的18组投影衰减数据并对这些数据进行重建，在完成50个角度共900组投影数据采集及重建后，0号节点也完成 $\mathrm{image}=\underset{n=1}{\overset{900}{\mathrm{\Σ}}}{I}_n$ 运算，得到重建图像image矩阵；

18)对重建图像进行平滑滤波、去除负值等后处理操作，将图像image发至图像显示机，显示最终重建图像。

图5是使用本发明重建的二维图像，图6为使用单CPU+标准FBP算法的重建结果。

由于该CT***配有18组源-探测器，因此，转台仅需旋转1/18周即可完成数据采集，每个角度的投影采样时间约为0.1ms，完成全部采样约需5ms。使用本基于GPU集群的图像并行重建***，完成一个角度的重建约需0.2ms，完全所有角度的图像重建约需10ms。由于数据采集与重建非串行完成，而是完成一个角度的18组投影采集后，在进行第个角度数据采集的同时进行图像重建，因此采集与重建是并行进行的，完成一个断面的扫描及重建，共耗时12ms。

最后，以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种CT并行重建***的成像方法，所述CT并行重建***包括与CT探测器列阵连接的前端采样器，与所述前端采样器连接的中心节点，以及与所述中心节点连接的多个子节点，每个所述子节点均安装有图形处理器，所述中心节点用于将CT重建图像划分为多个块区域，并根据所划分的块区域对每个子节点的计算任务进行分配，所述子节点用于计算所分配到的块区域中各像素的重建值，所述CT并行重建***的成像方法包括如下步骤：

1)将重建图像划分为多个块区域；

2)中心节点接收前端采样器所采集的各角度下的原始投影数据；

3)中心节点为各子节点分配计算任务，并将原始投影数据传输给分配了计算任务的各子节点；每个分配了计算任务的子节点计算一个或多个块区域的重建值；

4)各子节点完成所分配的块区域内各像素的重建计算，将计算结果返回给中心节点；

5)中心节点根据各子节点的计算结果组合出完整的重建图像；

所述步骤2)中，所述前端采样器依次采集CT探测过程中各角度下的原始投影数据，所述中心节点分批次依次接收各角度下的原始投影数据，一个角度下的原始投影数据为一个批次；

所述步骤3)中，所述中心节点在接收到一个角度下的各探测器单元的原始投影数据后，根据当前各子节点的计算负荷情况将计算任务分配给当前计算负荷较小的子节点。

2.根据权利要求1所述的成像方法，其特征在于，所述步骤1)中，每个所述块区域含有的像素点的数目是64的倍数。

3.根据权利要求1所述的成像方法，其特征在于，所述步骤1)中，每个所述块区域的像素点组成方阵。

4.根据权利要求2所述的成像方法，其特征在于，所述步骤1)中，每个所述块区域的像素点数目为256，所述块区域为16*16方阵。

5.根据权利要求4所述的成像方法，其特征在于，所述步骤3)中，当中心节点接收到一个角度下的原始投影数据时，将该角度下的重建图像的所有块区域的计算任务分配给同一子节点进行计算；所述步骤4)中，各 子节点将各角度下的重建图像返回给中心节点；所述步骤5)中，中心节点将各角度下的重建图像组合成最终的重建图像。

6.根据权利要求4所述的成像方法，其特征在于，所述步骤3)中，中心节点每接收到一个角度下的原始投影数据时，将该角度下的重建图像的各块区域的计算任务分配给不同的子节点进行计算；所述步骤4)中，各子节点将各自负责的块区域的重建图像返回给中心节点，所述步骤5)中，首先根据各子节点返回的各块区域的重建图像组合出每个角度下的重建图像；然后再得出最终的重建图像。

7.根据权利要求6所述的成像方法，其特征在于，所述步骤3)中，每个子节点负责该一个块区域的重建计算，每个子节点的图形处理器的一个线程负责一个像素的重建值计算。 

Images