CN101283913A - Ct图像重建的gpu加速方法 - Google Patents

Abstract

一种基于GPU硬件加速的CT图像重建方法，属于X射线CT技术领域。本发明软件部分通常包括基于GPU的CT数据预处理模块、基于GPU的CT数据滤波模块、基于GPU的CT正投影模块和基于GPU的CT图像重建反投影模块。本发明利用GPU硬件实现了对CT图像重建算法的加速，重建部分和数据处理部分均在GPU上完成。针对目前GPU显存不足以及内存到显存数据传输速度慢的瓶颈问题，提出了对大数据处理采用分块处理方法，与已有方法不同，本分块方法对每块待重建体重建时只需要用到部分投影数据，减少了数据传输，从而可提高整个体重建速度。

Description

CT图像重建的GPU加速方法

技术领域

本发明属于X射线CT技术领域，是一种基于计算机图形处理器(GraphicsProcessing Unit，缩写为GPU)的CT图像重建加速方法。

背景技术

如何提高CT图像的重建速度一直都是至关重要的问题，尤其针对三维锥束CT的体重建问题。锥束CT的数据量大，计算复杂度高，因此重建速度较慢，它是目前制约三维锥束CT应用的主要因素之一。在CPU上即使是优化的重建方法也要需要较长的时间，如从360个512×512的投影数据重建512×512×512的体数据大约也需要5分钟。欲重建高分辨的CT图像，需要更大规模的投影数据，如美国VARIAN公司PaxScan2520面板探测器图像矩阵1920*1536，探测器单元0.127*0.127平方毫米。为了使直径243.84mm的视野范围内分辨率约达到0.127～0.213mm，需要的投影个数(角度数)≥3600，相应的数据量将大于19GB.即使投影个数为720，相应的数据量也约4GB。

商业CT机中多采用专用集成电路(ASIC)和可编程逻辑阵列(FPGA)来实现CT算法加速，但其昂贵的价格一定程度上限制了在研究和医学等方面的广泛应用，更主要的缺点是ASIC或FPGA缺乏灵活性和可更改性，不利于升级或更新。

GPU的单指令多数据流(Single Instruction Multiple Data，缩写为SIMD)的处理模式可并行地对大规模数据进行同样的操作。由于计算机游戏和工程设计的巨大市场驱动，GPU的发展速度大大超过了CPU的发展速度，图形流水线的高速度和高带宽极大地提高了图形处理能力，近年来发展起来的可编程功能为图形处理以外的通用计算提供了高性价比的运行平台，使得基于GPU的通用计算成为近年来的研究热点之一。

CT数据形成的平行束投影和锥束投影的概念分别与计算机图形学中的正交投影和透视投影的概念非常相似，CT图像重建正投影和反投影过程与计算机图形学的渲染过程也很相似，CT理论中的X光源和探测器分别对应于计算机图形学中的摄像机和渲染目标，CT理论中的投影地址的概念也对应于图形学中纹理坐标的概念。由于GPU利用专用的纹理映射硬件计算纹理坐标(投影地址)，并能够并行处理多个被投影的数据，利用GPU可以将投影操作的速度提升一到两个数量级(与用CPU执行投影操作相比较)。反投影操作类似幻灯机将幻灯片投影到屏幕的过程，利用GPU也可以方便实现。因此利用GPU加速CT重建受到人们的广泛关注。

经对现有技术的文献检索发现：利用GPU加速CT重建速度的文章主要分为两类：一类是利用GPU加速CT重建的反投影过程，其中所用到的数据是在CPU上事先滤波好的，仅在GPU内完成反投影过程；另一类是利用GPU进行迭代法的加速，其中加权数据需要在CPU上事先计算好，迭代时需将加权数据表用纹理传入，因此占用了有限的显存空间，难以从大数据重建CT图像。美国纽约大学Klaus Mueller博士是国际上最早研究GPU加速CT重建的学者之一，2007年Klaus Mueller，Fang Xu和Neophytos Neophytou在SPIE Electronic Imaging 2007，Computational Imaging V Keynote发表《Why do Commodity Graphics HardwareBoards(GPUs)work so well for acceleration of Computed Tomography？》。他们认为GPU不适合进行FFT运算，因此采用CPU下的FFT进行滤波，投影数据预处理也在CPU上完成，没有充分利用GPU图形流水线的高速度和高带宽特点；2006年Klaus Mueller和Fang Xu在生物学和医学图像(IEEE Symp.Biomedical Imaging)发表的文章《PRACTICAL CONSIDERATIONS FORGPU-ACCELERATED CT》中为了提高重建速度进行了体数据分块，尽管在一定程度上提高了重建速度，但是由于重建每块体数据均用到所有的投影数据，因此投影数据的反复传输造成了CPU到GPU带宽的较大浪费。

发明内容

本发明的目的在于提供一种CT图像重建的GPU加速方法，针对目前GPU显存不足，以及内存到显存数据传输速度慢的瓶颈问题，提出了对大数据处理采用分块处理方法，与已有方法不同，本分块方法可以大幅度减少数据传输，从而可提高整个体重建速度。本发明的技术方案如下：

一种CT图像重建的GPU加速方法，软件部分通常包括基于GPU的CT数据预处理模块、基于GPU的CT数据滤波模块、基于GPU的CT正投影模块和基于GPU的CT图像重建反投影模块，通过如下步骤实现：

(1)从扫描的X射线强度数据获得投影数据并用基于GPU的CT数据预处理模块进行数据的去噪和去除环状伪影的处理；

(2)对投影数据用基于GPU的CT数据滤波模块和基于GPU的CT反投影模块采用滤波反投影型算法或反投影滤波型算法获得重建体，或者用基于GPU的CT正投影模块和基于GPU的CT反投影模块采用迭代算法获得重建体。

为了保证重建CT图像的质量，在上述方案中扫描数据、暗场数据、亮场数据和重建体数据均存储成32位浮点格式纹理，以保证重建图像的精度。利用GPU的渲染到纹理的技术，以提高处理数据的效率，采用多纹理技术，以减少GPU的渲染次数。所述的基于GPU的CT数据预处理模块用于数据的去噪、数据的提取、环状伪影的去除，以及从扫描的X射线强度数据获得投影数据等，所述的基于GPU的CT数据滤波模块用于对投影数据的各种滤波，所述的基于GPU的CT正投影模块用于模拟CT的扫描过程，用于获得扫描投影数据或用于迭代法求解，所述的基于GPU的CT反投影模块用于获得反投影图像或者反投影重建体，也用于迭代法求解。

基于GPU的CT数据滤波模块，将数据通过GPU上FFT变换到频域，频域滤波后通过GPU上的逆FFT获得滤波后的数据。基于GPU的CT数据滤波模块在GPU上实现FFT变换的方法为：FFT排序部分在GPU上完成，将预算好的FFT的系数和蝶式索引存成纹理，通过多次渲染实现蝶式运算，通过拷贝到纹理技术将处理后的数据替换原投影数据。

上述方案中，基于GPU的CT数据滤波模块和基于GPU的CT反投影模块组合可以实现多种解析类CT重建算法，包括滤波反投影类重建算法和反投影滤波类重建算法。

本发明的技术方案中，基于GPU的CT正投影模块，采用GPU的可编程管线实现多种扫描模式的模拟，其中基于GPU的CT正投影模块，采用了正投影权值实时计算的并行算法。基于GPU的CT反投影模块根据系统参数来设置变换矩阵，包括投影变换矩阵和模式变换矩阵。进一步讲，所述的基于GPU的CT反投影模块中将邻近四角度投影数据存储在一张纹理的四个通道，分别计算投影矩阵，所述的基于GPU的CT反投影模块中将邻近4片切层体数据存储在一张纹理的四个通道。另外，基于GPU的CT反投影模块根据计算机显存配置情况选择适宜方法对重建体进行分块重建。

上述技术方案中，相比于通过“依赖纹理”传递权值的方法，不仅提高了权值的计算速度，避免了权值在内存和显存之间的传递，且节省了显存用于保存权值所用的空间。所述的基于GPU的CT正投影模块，采用了GPU的可编程管线，算法灵活，可以实现多种扫描模式的模拟，包括圆轨迹断层扫描，圆轨迹锥束扫描，单层螺旋轨迹扫描，螺旋锥束扫描，断层大视野扫描，断层多次超视野扫描，锥束大视野扫描，锥束多次超视野扫描等扫描模式。所述的基于GPU的CT正投影模块，在一次渲染中即可计算出物体7片切片的正投影，而基于固定管线的正投影方法经过两次渲染才能计算出一片切片的正投影，本发明的算法所需计算步骤少，计算更简洁。所述的基于GPU的CT反投影模块中没有采用传统的四对称的投影数据存储方法，而是将邻近四角度投影数据存储在一张纹理的四个通道，分别计算投影矩阵。此方法直接重建出待重建切层，不需要额外在CPU上进行体数据的累加。所述的基于GPU的CT反投影模块中也没有采用传统的四对称的体数据存储方法，而是将邻近4片切层体数据存储在一张纹理的四个通道，分别计算投影矩阵。此方法可以一次渲染直接重建出4片待重建切层，减少了GPU的渲染次数。所述的基于GPU的CT反投影模块根据计算机显存配置对重建体进行分块重建，克服了显存大小对重建体的影响。

上述方案中，基于GPU的CT正投影模块和基于GPU的CT反投影模块组合可以实现多种迭代重建算法，包括代数迭代算法、EM迭代算法、子集排序的代数迭代算法和子集排序EM迭代算法等。

本发明充分利用GPU可编程管线(programmable pipeline)的特点，完全在GPU内实现了多种CT重建算法的加速。本发明包括基于GPU的CT数据预处理模块、基于GPU的CT数据滤波模块、基于GPU的CT的正投影模块(也称投影数据生成模块)和基于GPU的CT图像重建反投影模块。提出了合理的利用GPU实现数据滤波方案，数据处理部分和重建部分均在GPU内实现。针对目前GPU显存不足，以及内存到显存数据传输速度慢的瓶颈问题，提出了对大数据分块处理方法，与已有方法不同，本分块方法可以大幅度减少数据传输，从而可提高整个体重建速度。

下面将按模块更为详细的描述本发明的技术方案：

(1)基于GPU的CT数据预处理模块

①将相邻4个角度的数据分别存储到一张四通道矩形纹理的4个通道中，纹理格式GL_RGBA，设投影个数为NumAngle，每张纹理用tex[IAngle]表示，IAngle＝0，…，NumAngle/4-1；暗场数据存储成一张单通道矩形纹理，记作texNoise，纹理格式GL_INTENSITY；亮场数据存储成一张矩形纹理，记作texNormal，纹理格式GL_INTENSITY。纹理的长和宽分别均记作Wtex和Htex，纹理用32位浮点格式，以保证计算的精度。

②将数据纹理tex[IAngle]，暗场数据纹理texNoise和亮场数据纹理texNormal，利用多纹理技术贴到一长为Wtex，宽为Htex的矩形框上，投影矩阵为正交投影矩阵。

③顶点变换用固定管线的变换，编写元处理器，使得三片纹理按 $\mathrm{gl}\_\mathrm{FragColor}=\log \frac{\mathrm{texNormal}-\mathrm{texNoise}}{\mathrm{tex}\left[\mathrm{IAngle}\right]-\mathrm{texNoise}}$ 形式混合，gl_FragColor为最后的混合结果。

④利用拷贝到纹理技术，将混合后的结果拷贝到数据纹理tex[IAngle]。

不难理解，同样原理和方法，修改片元shader程序，即可实现对数据的其它方式的预处理。

(2)基于GPU的CT数据滤波模块

①利用GPU上的FFT将数据变换到频域。由于本发明的目的是对大量数据进行一维频域滤波，因此采用FFT系数和蝶式索引预先计算的方法，而没有采用通常FFT系数和蝶式索引临时计算的方法。GPU上数据的FFT变换分以下几步实现：

第一步，对数据纹理tex[IAngle]补零，使其长度为2的整次幂记作WtexNew＝2^M，宽度不变为Htex，补零的方法是将纹理tex[IAngle]贴在长为WtexNew，宽为Htex的矩形框上，投影矩阵为正交投影矩阵。利用渲染到纹理技术，渲染到一新的纹理texNew，格式为GL_RGBA，长为WtexNew，宽为Htex；

第二步，在CPU上计算出FFT系数和蝶式索引以及滤波器的频域数据，将这些数据存储成一张纹理，记作texCoff，纹理格式GL_RGBA.长为WtexNew，宽为M+1；

第三步，顶点变换用固定管线的变换，编写片元处理器，在片元shader中访问纹理texCoff，利用它所存储的FFT系数和蝶式索引经过M次渲染实现蝶式算法。获得数据的频域表示。

②将频域数据与滤波函数的频域数据点乘，实现频域滤波。由于纹理texCoff已存储了滤波函数的频域数据，利用多纹理混合可容易实现点乘运算。

③利用GPU上FFT的逆变换将数据变换到实域。众所周知，傅里叶正变换与傅里叶逆变换具有对称关系，记f(t)原函数，F(ω)为f(t)的傅里叶变换

则有

因此很容易利用GPU上FFT实现GPU上FFT逆变换。

④将纹理texNew贴在长为Wtex，宽为Htex的矩形框上，投影矩阵为正交投影矩阵，利用拷贝到纹理技术，将渲染后的结果拷贝到数据纹理tex[IAngle]。

(3)基于GPU的CT正投影模块

本发明的正投影模块完全在GPU的可编程管线中实现，不仅支持32bit浮点运算精度，而且算法需要的计算步骤少，实现更加灵活；采用GPU的正投影权值实时并行算法，相比于通过“依赖纹理”传递权值的方法，不仅提高了权值的计算速度，避免了权值在内存和显存之间的传递，且节省了显存用于保存权值所用的空间。

①将被测物体按竖切片刨切成NumRow片切片，每片切片和与它有对称关系的切片的数据存储在一张四通道的矩形纹理中，纹理格式GL_RGBA，长为NumRow，宽为NumSlice。

②根据被测物体的切片模型，将切片纹理映射到所对应得几何体切片上，纹理坐标的映射关系通过可编程管线计算获得，利用GPU渲染到纹理技术，将多个切片渲染到平板探测器对应的纹理上，经过多次渲染，将所有切片渲染到平板探测器对应的纹理上获得没有校正的投影数据。

③根据正投影的几何参数，利用GPU可编程管线计算正投影权值，对投影数据进行校正。

④采用了OpenGL快速数据传输扩展PBO(Pixel Buffer Object)进行切片纹理数据和投影数据的传输，该数据传输扩展不仅能加快数据在内存和显存之间传输速度，而且支持数据传输和GPU计算同步进行。

(4)基于GPU的CT反投影模块

①投影数据的存储，投影数据尺寸为NumWidth×NumHeight×NumAngle，存储成NumAngle/4片四通道矩形纹理，纹理格式GL_RGBA，长为NumWidth，宽为NumHeight，每4个相邻角度的投影数据存储在一张纹理的四个通道中。

②重建体数据的存储，假设重建体为NumRow×NumRow×NumSlice，横切片的存储方法是将重建体存储成NumSlice/4片GL_RGBA格式的矩形纹理，长为NumRow，宽为NumRow.每4个相邻横层切片的重建体数据存储在一张纹理的四个通道中，竖切片的存储方法是将重建体存储成NumRow/4片GL_RGBA格式的矩形纹理，长为NumRow，宽为NumSlice，每4个相邻竖层切片的重建体数据存储在一张纹理的四个通道中。

③根据CT系统的参数设定透视投影矩阵，CT系统的参数包括：探测器单元的尺寸，面阵探测器的行单元个数，面阵探测器的列单元个数，射线源到转台的距离，射线源到探测器的距离等。CT系统中射线源与图形学中视点概念相对应，探测器与图形学中投影平面概念相对应，CT锥束扫描与透视投影模式相对应，因此根据CT系统参数容易计算透视投影矩阵。

④利用投影纹理技术实现反投影，将投影数据纹理贴到待重建切层(横切层或竖切层)上，纹理坐标利用待重建切层上的重建点在投影平面上的投影坐标得到。

在实际应用中有两种方法减少显存消耗，一种是将所有投影角度下的投影数据纹理全部反投重建一片重建体纹理后，将重建体纹理传出显存，因此显存中除放置所有投影数据纹理外，只需要放置一片重建体纹理即可。另一种方法是将一片投影数据纹理或几片投影数据纹理对所有的重建切层进行反投影，然后释放这几片投影数据纹理，从内存重新读入新的投影数据纹理，继续对重建切层进行反投影，因此显存中除了放置重建切层纹理外，只需要放置几片投影数据纹理即可。

上述两种方法的选择需要根据实际情况中重建体大小和投影数据大小来定，如果重建体数据较大，则选择前一种将重建切层分次重建的方法，如果投影数据较大，则选择后一种将投影数据纹理分次读入的办法。

⑤大数据量下的重建问题，尽管④中的两种方法都可降低显存消耗，但是都无法解决投影数据和重建体数据均超出显存限制的问题。当投影数据和重建体数据均超出显存限制时，本发明采用对投影数据和重建体数据同时分块的方法，重建体的纹理存储采用横切片的纹理存储方式，然后根据显存大小计算可以容纳的切层数量，依据重建体的横切层分块情况对每个投影角度下的投影数据进行分块，也就是计算重建每一块时所需要的每个投影下的投影数据，将它分次传入显存。之所以选择横切片的体数据存储方式是因为每片横切片在每个投影角度小只需要少量几行投影数据，因此可以达到解决大数据量下的图像重建问题。

(5)模块组合实现多种GPU上的CT重建算法

①圆轨迹FDK重建算法

首先利用基于GPU的CT数据预处理模块对射线强度数据进行处理，在处理的过程中投影数据的获得和投影数据预校正同时进行；然后利用基于GPU的CT数据滤波模块选择斜坡滤波器沿探测器行方向对投影数据进行一维滤波；利用基于GPU的CT反投影模块进行加权反投影，权重采用实时计算的方法，模型变换矩阵设置为绕探测器的竖直方向旋转的旋转变换矩阵，透视投影矩阵根据系统参数设置。

②锥束多次扫描数据的BPF型重建算法

首先利用基于GPU的CT数据预处理模块对射线强度数据进行处理，在处理的过程中投影数据的获得和投影数据沿探测器行方向的求导同时进行；然后利用基于GPU的CT反投影模块生成投影数据导数的三维反投影图像(也称三维DBP图像)，其中设置模型变换矩阵为先沿探测器的水平方向平移，再绕探测器的竖直方向旋转的平移旋转变换矩阵，根据系统参数设置透视投影矩阵；最后利用基于GPU的CT数据滤波模块对三维DBP图像沿行方向进行一维滤波，滤波核函数为Hilbert逆变换函数。

③代数迭代重建算法

首先利用基于GPU的CT数据预处理模块将探测器获取的射线强度数据转换成投影数据；然后给待重建的体数据付初始值，利用基于GPU的CT正投影模块获得投影数据的估计值；接着利用基于GPU的CT数据预处理模块计算估计值和实际投影值的残差；再利用基于GPU的CT反投影模块将残差值反投影给初始的待重建体数据。这样就完成了一次迭代，将重建结果作为下次迭代的初始值实现多步迭代过程。

附图说明

图1为实采数据利用本发明重建的效果图；

图2为图1的局部体显示图；

图3为用于测试的模型的照片；

图4为基于GPU的大视野重建的DBP效果图，图4(a)为主视方向的中心切层图，图4(b)为俯视方向的中心切层图，图4(c)为侧视的中心切层图；

图5为基于GPU的大视野重建的效果图，图5(a)为重建体的体渲染图，图5(b)为正视图的中心切层图，图5(c)为俯视图的中心切层图，图5(d)为侧视的中心切层图。

具体实施方式

下面将结合实施例对本发明作进一步的说明。利用本发明的锥束多次扫描数据的BPF型重建算法对一次偏置大视野的CT数据进行重建，扫描模型为Shepp-logan模型，GPU使用Nvidia的Quadro FX 4600芯片。

(1)投影数据大小为512×512×360，待重建体的大小为512×512×512，重建体数据大于投影数据，因此选取数据纹理一次存储成纹理，重建切层分次重建的反投影方法；

(2)将投影数据存成90片长为512，宽为512，格式为GL_RGBA的四通道32位浮点矩形纹理，每4个相邻角度的投影数据存储在一张纹理的四个通道中；

(3)利用基于GPU的CT数据预处理模块对投影数据纹理沿行方向进行求导，利用拷贝到纹理技术将求导后的投影数据替换原投影数据；

(4)CT系统参数为：探测器单元尺寸为0.127mm，转台中心到探测器的距离为240mm，转台中心到主射线方向的距离为50.5mm，射线源到探测器的距离为2150mm。根据上述参数设定变换矩阵，包括透视投影矩阵和模型视图矩阵；

(5)利用基于GPU的CT反投影模块对重建体进行反投影，获取导数的反投影三维图像，在反投影中需要根据投影地址计算符号函数和权重。重建效果如图4所示，投影数据导数的反投影图像也称DBP图像。

(6)利用基于GPU的CT数据滤波模块对重建的三维图像沿行方向进行滤波，滤波核函数为Hilbert逆变换函数。

整个重建过程从硬盘数据的读入到重建数据的硬盘保存共耗时30秒，重建效果如图5所示。

需要说明的是：虽然上述的实施例阐述了本发明的模块组合能够实现的几种重建算法，但是本领域的研究人员可以理解，本发明可以将滤波反投影类重建算法、反投影滤波类重建算法、迭代类重建算法通过本发明的模块移植到GPU上实现，进而来提高它们的重建速度。

Claims (10)

1.一种CT图像重建的GPU加速方法，软件部分通常包括基于GPU的CT数据预处理模块、基于GPU的CT数据滤波模块、基于GPU的CT正投影模块和基于GPU的CT图像重建反投影模块，通过如下步骤实现：

(1)从扫描的X射线强度数据获得投影数据并用基于GPU的CT数据预处理模块进行数据的去噪和去除环状伪影的处理；

(2)对投影数据用基于GPU的CT数据滤波模块和基于GPU的CT反投影模块采用滤波反投影型算法或反投影滤波型算法获得重建体，或者用基于GPU的CT正投影模块和基于GPU的CT反投影模块采用迭代算法获得重建体。

2.如权利要求1所述的CT图像重建的GPU加速方法，其特征在于扫描数据、暗场数据、亮场数据和重建体数据均存储成32位浮点格式纹理。

3.如权利要求2所述的CT图像重建的GPU加速方法，其特征在于对投影数据用基于GPU的CT数据滤波模块和基于GPU的CT反投影模块采用滤波反投影型算法或反投影滤波型算法获得重建体，基于GPU的CT数据滤波模块，将数据通过GPU上FFT变换到频域，频域滤波后通过GPU上的逆FFT获得滤波后的数据。

4.如权利要求3所述的CT图像重建的GPU加速方法，其特征在于所述的基于GPU的CT数据滤波模块中在GPU上实现FFT，其中FFT排序部分在GPU上完成，将预算好的FFT的系数和蝶式索引存成纹理，通过多次渲染实现蝶式运算，通过拷贝到纹理技术将处理后的数据替换原投影数据。

5.如权利要求2所述的CT图像重建的GPU加速方法，其特征在于用基于GPU的CT正投影模块和基于GPU的CT反投影模块采用迭代算法获得重建体，基于GPU的CT正投影模块，采用了GPU的可编程管线，实现多种扫描模式的模拟。

6.如权利要求5所述的CT图像重建的GPU加速方法，其特征在于所述的基于GPU的CT正投影模块，采用正投影权值实时计算的并行算法。

7.如权利要求3或5所述的CT图像重建的GPU加速方法，其特征在于所述的基于GPU的CT反投影模块根据系统参数来设置变换矩阵，包括投影变换矩阵和模式变换矩阵。

8.如权利要求7所述的CT图像重建的GPU加速方法，其特征在于所述的基于GPU的CT反投影模块中将邻近四角度投影数据存储在一张纹理的四个通道，分别计算投影矩阵。

9.如权利要求7所述的CT图像重建的GPU加速方法，其特征在于所述的基于GPU的CT反投影模块中将邻近4片切层体数据存储在一张纹理的四个通道。

10.如权利要求8或9所述的CT图像重建的GPU加速方法，其特征在于所述的基于GPU的CT反投影模块根据计算机显存配置对重建体进行分块重建。

Images