CN106725565B - 一种稀疏投影下的锥束xct成像质量评估方法 - Google Patents

Abstract

本发明公开一种稀疏投影下的锥束XCT成像质量评估方法，该方法在得到重建出仿体的三维图像后，进行成像系统质量评估，包括：选取多个三维重建结果切片，沿z轴方向叠加平均三维重建结果切片得到二维图像，并得到极坐标系下图像；进行阈值分割，识别出目标体区域，然后通过质心法确定目标体区域的中心。设定过采样像素大小，将极坐标系下图像中各个同心圆环范围内的像素点进行像素值平均，而后按照离开圆心的距离排序，由此得到过采样边缘扩散函数曲线ESF。求调制传递函数；计算锥束XCT成像系统的二维噪声功率谱。得到一维噪声功率谱，计算噪声等效量子数NEQ。

Description

一种稀疏投影下的锥束XCT成像质量评估方法

技术领域

本发明属于生物医学工程及计算机领域，涉及一种锥束XCT成像质量评估方法。

背景技术

近年来，全球经济发展速度极快，人均寿命也在不断提升，人们对于医疗的依赖程度也在不断地加大，同时人们对于各种无创医疗手段也更加追求。X线设备作为医疗影像学中重要的仪器设备之一，自然成为了人们所关注的焦点。在提升X线设备成像质量的同时，降低X线照射剂量则是广大患者和医疗研究人员追求的目标之一。如何用最小的X线照射剂量，产生质量清晰的影像，是X线设备研究人员关注的重点。目前，由于锥束XCT具有成像速度快，射线剂量相对较低等优点，在临床应用中发展迅速，而低剂量锥束XCT研究的关注也愈发广泛。降低锥束XCT射线剂量主要有两种方式：一是改变管电流和管电压，这样的结果是图像噪声较大，图像的信噪比低；二是减少投影次数，但这种情况下采用传统的重建方法会造成重建时数据不一致性提高，重建质量下降。

针对上述问题，在近年来的医学影像领域，发展稀疏投影下的锥束XCT成像方法以保证重建图像质量的研究得到了广泛关注，这种稀疏投影为前提的XCT重建方法，由于减少了投影成像数量，不仅减少辐射剂量，减少扫描时间，也避免了因为长时间扫描中病人的移动而造成的运动伪影，因此对于新一代锥束XCT系统的发展具有重要的意义。但到目前为止，这类方法尚未得到临床的尝试和应用，其中一个重要的原因在于缺乏合理可靠的成像质量客观定量评估方法，目前相关的文献中只是给出稀疏投影锥束XCT重建与传统方法的定性比较，缺乏定量化标准。

对于稀疏投影锥束XCT图像重建而言，需要在三维成像空间内进行成像系统和成像方法的评估，比传统的投影成像系统，例如X射线透视这种二维成像系统的性能评估更加复杂。到目前为止，尚无被广泛认可的锥束XCT成像系统性能评估的标准，更没有针对稀疏投影下锥束XCT成像质量的评估。

当前国际上公认的成像系统性能评价参数包括调制传递函数(ModulationTransfer Function,MTF)，噪声功率谱(Noise Power of Spectra,NPS)，和噪声等效量子数(Noise Equivalent Quanta,NEQ)等，但仅给出了针对X射线投影成像系统的标准，且对应的标定仿体也仅适合于二维投影系统。而在三维的锥束XCT成像系统下，成像系统及方法的差别造成二维投影系统下采用的标定仿体和方法不再适用，更不适用于稀疏投影下的锥束XCT成像质量评估。目前有部分学者提出采用直径小于像素值大小的钨丝来实现锥束XCT的调制传递函数测量，但是这种方法中，由于钨丝直径小，成像结果质量差，且这种方法不易实现过采样，导致所得到的MTF结果缺失信息甚至失真。

因此，为了实现稀疏投影下的锥束XCT成像质量评估，设计新型仿体，发展出一种有效的评估方法，实现成像质量评估的标准化及定量化，已经成为稀疏投影下锥束XCT成像质量评估的研究重点，这对于成像系统的发展应用具有十分重要的意义。

发明内容

本发明的主旨是提供一种稀疏投影下的锥束XCT成像质量评估方法，以此解决锥束XCT成像中稀疏投影三维重建的成像质量评估的基本问题：实现稀疏投影下的锥束XCT成像结果的调制传递函数，噪声功率谱及噪声等效量子数的准确测量，为评估欠定条件下稀疏投影锥束XCT成像方法的性能提供有力条件。技术方案如下：

一种稀疏投影下的锥束XCT成像质量评估方法，该方法将仿体置于锥束XCT成像系统的点光源下方的成像腔，成像腔在水平平面上绕着成像腔中轴旋转，每旋转一定角度，进行一次曝光成像，利用TV正则化方法重建出仿体的三维图像；基于此三维图像进行成像系统质量评估，分为两个部分，

第一部分，计算稀疏投影下的锥束XCT成像系统的调制传递函数MTF

①选取多个三维重建结果切片，沿z轴方向叠加平均三维重建结果切片得到二维图像p(x,y)，(x,y)为图像像素在直角坐标系下位置，并将其从直角坐标系映射到极坐标系，得到极坐标系下图像p(r,θ)。

②对图像p(x,y)采用阈值分割将灰度图像p(x,y)映射为二值图像p^norm(x,y)，这里为图像中空气部分的像素平均值，为图像中成像质量评估仿体部分的像素平均值，将仿体和空气部分分割，识别出目标体区域，然后通过质心法确定目标体区域的中心。

③设定过采样像素大小Δr，由此形成半径分别为Δr，2Δr，3Δr，4Δr……的一系列同心圆环，将图像p(r,θ)中各个同心圆环范围内的像素点进行像素值平均，而后按照离开圆心的距离排序，由此得到过采样边缘扩散函数曲线ESF。

④对ESF求微分以后得到线扩散函数LSF，对LSF加hann窗后进行傅里叶变换，所得结果取模并进行零频率归一化后即为所求的调制传递函数MTF(u),这里的u为频率；

第二部分，求取稀疏投影下的锥束XCT成像系统的噪声功率谱NPS：

⑤对于所选取的三维重建结果切片，分别在每一个切片的圆形区域内的同一区域截取方形ROI，计算该方形ROI的平均像素值，而后将原ROI图像减去该平均像素值，以此获得各个切片的噪声图像。

⑥将每一幅噪声图像都分成互不重叠的几个子块ROI，而后分别对所有子块ROI计算噪声功率谱，计算噪声功率谱时采用的傅里叶变换的点数与计算MTF时采用的点数一致，再将所得的噪声功率谱叠加平均，得到锥束XCT成像系统的二维噪声功率谱。

⑦从二维噪声功率谱图中心沿横轴向右至终点，取横轴上下各7行功率谱值，而后将这14行功率谱值叠加平均，得到锥束XCT成像系统的一维噪声功率谱NPS_d(u)。

⑧根据公式计算得到噪声等效量子数NEQ，即得到成像质量空间分辨率和噪声水平的综合评价结果。

本发明的稀疏投影下的锥束XCT成像质量评估方法，该发明针对目前稀疏投影下锥束XCT成像缺乏定量评估方法的问题，利用圆柱体仿体，并设计相应的极坐标系下的过采样调制传递函数测量方法以及噪声功率谱测量方法，进而获取噪声等效量子数，以实现稀疏投影三维重建的成像质量的准确定量评估。该发明为准确获取稀疏投影锥束XCT三维重建的成像系统的噪声等效量子数，有效评估其成像性能奠定了方法学基础，从而为稀疏投影锥束XCT三维重建方法在影像学临床实践的应用推广提供有力支持。

附图说明

图1成像系统评估仿体设计

图2边缘扩散函数(ESF)形成示意图

图3噪声功率谱(NPS)计算中ROI选取示意图

图4锥束XCT成像系统示意图

图5(a)(b)(c)分别为.不同投影角度数情况下，滤波反投影算法和TV正则化算法的MTF评估结果对比

图6(a)(b)(c)分别为不同投影角度数情况下，滤波反投影算法和TV正则化算法的NPS评估结果对比

图7.(a)(b)(c)分别为不同投影角度数情况下，滤波反投影算法和TV正则化算法的NEQ评估结果对比

具体实施方式

针对目前在稀疏投影锥束XCT三维重建方法定量评估方面的问题，本发明提出一种稀疏投影下的锥束XCT成像质量评估方法，可以准确获得与成像质量相关的过采样调制传递函数，噪声功率谱以及噪声等效量子数。最后得到的技术方案如下：

1.设置稀疏投影下旋转步进角度(比如每隔2度旋转物体一次，则对应180个投影)，将设计的成像质量评估仿体安置于成像腔中，而后在物体旋转到设定角度时，进行曝光成像；

2.利用TV正则化方法重建出仿体的三维图像；

3.按照如下步骤计算稀疏投影下成像系统的调制传递函数MTF：

①.沿z轴方向叠加平均三维重建结果切片(选取6个切片叠加)得到二维图像p(x,y)，(x,y)为图像像素在直角坐标系下位置，并将其从直角坐标系映射到极坐标系，得到极坐标系下图像p(r,θ)。

②.对图像p(x,y)采用阈值分割将图像像素值映射到[0，1]，即将灰度图像p(x,y)映射为二值图像p^norm(x,y)，这里为图像中空气部分的像素平均值，为图像中成像质量评估仿体(圆形目标体)部分的像素平均值，将物体和空气部分分割，识别出目标体区域，然后通过质心法确定圆形目标体区域的中心。

③.设定过采样像素大小Δr，由此形成半径分别为Δr，2Δr，3Δr，4Δr……的一系列同心圆环，将图像p(r,θ)中各个同心圆环范围内的像素点进行像素值平均，而后按照离开圆心的距离排序，由此即可得到过采样边缘扩散函数曲线(ESF)。

④.对ESF求微分以后得到LSF，对LSF加hann窗后进行傅里叶变换，所得结果取模并进行零频率归一化后即为所求的MTF(u),这里的u为频率。

4.按照如下步骤计算稀疏投影下成像系统的噪声功率谱NPS：

①.选择计算MTF时采用的6个三维重建结果切片，而后分别在每一个切片的圆形区域内的同一区域截取方形ROI，计算该方形ROI的平均像素值，而后将原ROI图像减去该平均像素值，以此获得6幅噪声图像。

②.将每一幅噪声图像都分成互不重叠的4个方形ROI，而后分别对所有子块ROI计算噪声功率谱，计算噪声功率谱时采用的傅里叶变换的点数与计算MTF时采用的点数一致。再将所得的噪声功率谱叠加平均，可以得到锥束XCT成像系统的二维噪声功率谱。

③.从二维噪声功率谱图中心沿横轴向右至终点，取横轴上下各7行功率谱值，而后将这14行功率谱值叠加平均，得到锥束XCT成像系统的一维噪声功率谱NPS_d(u)。

5.根据公式计算得到噪声等效量子数NEQ，即得到成像质量空间分辨率和噪声水平的综合评价结果。

下面对本发明进行详细说明。

1稀疏投影下的锥束XCT成像方法

对于锥束CT成像系统，成像模型可以用下式来表示：

其中I是到达探测器的光子数，I₀是入射的光子数，l为从X射线源到探测器的积分直线，即X射线的路径，μ为待求目标区域吸收系数分布。实际采样和重建过程均是离散化模型，第i条射线穿过物体之后，到达探测器的光子数为：

其中，投影矩阵M的元素m_i,j表示第i条射线穿过图像第j个像素点的衰减系数。

得到的第i条射线的投影数据p_i可以表示为：

其中p_i为第i条射线对应的投影数据。那么将上式表示为矩阵形式，数字锥束XCT系统可以用如下模型描述：

P＝Mμ (4)

其中P为图像投影，M为CT系统投影矩阵。

对于稀疏投影下的锥束XCT投影数据重建图像问题，从数学角度讲是求解欠定线性方程组的问题，即用投影数据P和系统矩阵M求解原始图像μ。通过引入了待重建图像全变分(total variation，TV)先验信息，将上述问题等价于求解以下优化问题：

其中，第一项是保真项，为了满足成像模型，即公式(4)。第二项是基于对图像分段常数或者分段连续的假设，通过TV正则化方法引入分段常数的约束条件，锐化重建图像的边缘。||TV(·)||₁为TV范数，实际上是图像梯度的L₁-范数，τ为非负的正则化参数，用来平衡前述两项。另外，由于吸收系数μ是非负的，还需要在公式(3)满足的同时加入非负约束，即μ≥0。

2成像系统评估仿体设计

为了满足稀疏投影下的锥束XCT成像质量的评估，设计一种适合用于成像系统的调制传递函数(Modulation Transfer Function,MTF)及噪声功率谱(Noise PowerSpectrum)测量的仿体。仿体材料选择接近于人体软组织的聚丙烯酸酯，圆柱体边界平滑，内部材料均匀，所设计的仿体为直径20厘米，厚度4厘米的圆柱体，如图1所示。

通过改变投影数量，进行不同投影数量下的锥束XCT重建，以此实现基于稀疏投影的锥束XCT成像质量的准确评估。

3成像质量的调制传递函数(MTF)评估方法

利用所设计的仿体来测量MTF，采用的原理是求得径向过采样边缘扩散函数(EdgeSpread Function,ESF)，然后求微分得到线扩散函数(Line Spread Function,LSF)，对LSF进行傅里叶变换并取模求得MTF。基本的步骤如下：

1)沿z轴方向叠加平均三维重建结果切片(选取6个切片叠加)得到二维图像p(x,y)，此时的二维图像可以认为是系统的点扩散函数从三维降维到二维

这里(x,y)为图像像素在直角坐标系下位置，Δz是z轴叠加平均的切片数量。这里设定轴向分辨率是各向同性的，即沿x轴和y轴的空间分辨率相同。将二维图像p(x,y)从直角坐标系映射到极坐标系，得到极坐标系下图像p(r,θ)，其中代表着轴向的空间变量，是对应的极坐标下的极角。

2)用阈值分割的方法识别目标体区域，然后通过质心法确定圆形目标体区域的中心：阈值分割时将图像值映射到[0，1]内，0代表空气，1代表仿体区域的平均值，如下式：

即将灰度图像p(x,y)映射为二值图像p^norm(x,y)，这里为图像中空气部分的像素平均值，为图像中成像质量评估仿体(圆形目标体)部分的像素平均值。

3)设定过采样像素大小Δr(比如为原像素大小的1/10，则对应过采样倍数为10倍)，由此形成半径分别为Δr，2Δr，3Δr，4Δr……的一系列同心圆环，如图2所示给出了半径从Δr～10Δr的10个同心圆环，用不同颜色表示。将图像p(r,θ)中各个同心圆环范围内的像素点进行像素值平均，而后按照离开圆心的距离排序，由此即可得到图2右下角对应的过采样边缘扩散函数曲线(ESF)。

4)对ESF求微分以后得到线扩散函数(LSF)，对LSF加hann窗后进行傅里叶变换，所得结果取模并进行零频率归一化后即为所求的MTF(u),这里的u为频率。

4成像质量的噪声功率谱(NPS)评估方法

噪声功率谱被定义为空间频率域上随机信号在每一个频率位置的方差。根据维纳－辛钦理论，通过计算自相关函数的傅立叶变换就可以获得噪声功率谱。但在实际计算中，噪声功率谱通常通过对噪声密度图像的二维傅立叶变换取模来直接获取，具体计算公式如下：

这里FT[f(x,y)]是噪声图像f(x,y)的傅立叶变换，X和Y分别是图像在x方向和y方向的最大尺寸。

在数字成像系统中，上述的噪声功率谱的连续表达形式通常需要转换为离散形式进行计算。因此，需要重新计算连续噪声图像函数f(x,y)在一系列离散位置点的像素值，这些位置点是x＝mΔx,y＝nΔy,m＝0…M,n＝0…N，这里X＝MΔx,Y＝NΔy。这样就可以把连续图像表示为离散形式f(m,n)。相应的，空间频率变量u和v的函数也需要在离散频率位置重新计算，这些位置包括u＝kΔu,v＝lΔv,k＝0,±1,±2…±K,l＝0,±1,±2…±L。最大空间频率即奈奎斯特采样频率在两个频率轴方向分别为U＝1/2Δx及V＝1/2Δy。

为了估计数字成像系统地噪声功率谱，上述计算功率谱的连续表达式也需要转换为离散形式。根据上述定义，可以将公式(8)改写为：

这里FT_d[f(m,n)]是数字图像f(m,n)的二维傅立叶变换。

仍然采用图1所示的仿体，实现噪声功率谱测量的基本步骤如下：

1.选择计算MTF时采用的6个三维重建结果切片，而后分别在每一个切片的圆形区域内的同一区域截取方形ROI(如图3所示)，计算该方形ROI的平均像素值，而后将原ROI图像减去该平均像素值，以此获得6幅噪声图像。

2.将每一幅噪声图像都分成互不重叠的4个方形ROI，而后分别对所有子块ROI按照公式(9)计算噪声功率谱，计算噪声功率谱时采用的傅里叶变换的点数与计算MTF时采用的点数一致。再将所得的噪声功率谱叠加平均，可以得到锥束XCT成像系统的二维噪声功率谱。

3.从二维噪声功率谱图中心沿横轴向右至终点，取横轴上下各7行功率谱值，而后将这14行功率谱值叠加平均，得到锥束XCT成像系统的一维噪声功率谱NPS_d(u)。

5成像质量的噪声等效量子数(NEQ)评估方法

对于二维图像，其二维NEQ可以计算为：

因为具有空间频率依赖性的NEQ参数用量子来表示(或更精确地，每单位面积的量子)，该指标能够客观地比较在不同成像系统获取的图像的质量，对于在实际医学影像仪器和工业应用中可以方便地比较配置如散射光栅、探测器探元、光源X射线束、剂量等的不同以及重建算法的不同而带来的图像的差异。

在一维情况下，NEQ可以表示为如下形式

本发明采用美国BIOPTICS公司生产的Pixarray 100小动物数字放射成像系统构建一套锥束XCT成像系统，采用如图4采用的成像方式，所设计仿体安置在圆柱形物体承载容器位置。

我们将分别进行了180个角度、90个角度以及60个角度下对成像质量评估仿体进行投影。将上述三种情况下，分别采用目前锥束XCT采用的传统的滤波反投影算法以及TV正则化方法进行三维重建，而后依据本发明的技术流程依次对重建结果的中心6个切片求解MTF、NPS以及NEQ。

对不同投影角度数情况下的不同算法MTF结果之间比较如图5所示，从图中可以发现，在三种不同投影角度数情况下，TV正则化方法的MTF值均高于滤波反投影算法，滤波反投影方法得到的MTF在空间频率为0.3mm^-1时已经达到极限空间分辨率(10％MTF)，而对于TV正则化方法，三种情况下的MTF曲线中，极限空间分辨率都不低于0.7mm^-1，远高于传统的滤波反投影方法。

对不同投影角度数情况下的不同算法NPS结果之间比较如图6所示，从图中可以发现，在三种不同投影角度数情况下，TV正则化方法的NPS均远低于传统的滤波反投影算法。

图7显示在不同投影角度情况下，传统滤波反投影方法与TV正则化方法对应的系统NEQ结果之间的比较。从图中可以明显看到，TV正则化方法获得了远高于传统滤波反投影方法的NEQ值，即表明了TV正则化方法在获取的图像的空间分辨率及噪声方面具有更优的性能。

最终结果表明，在进行数字X射线成像系统的调制传递函数测量时，通过本发明的基于X射线成像系统探测器特性的调制传递函数测量新方法，充分利用了探测器ESF曲线的单调性和凹凸性，可以获得较传统刀口测量方法更为准确的MTF曲线。该方法的应用，将为准确测量系统调制传递函数，有效评估放射成像系统性能，深入开展放射影像学临床实践和研究提供有力支持。

Claims (1)

1.一种稀疏投影下的锥束XCT成像质量评估方法，该方法将仿体置于锥束XCT成像系统的点光源下方的成像腔，所述的仿体为由聚丙烯酸酯制成的边界平滑，内部材料均匀的圆柱体；成像腔在水平平面上绕着成像腔中轴旋转，每旋转一定角度，进行一次曝光成像，利用TV正则化方法重建出仿体的三维图像；基于此三维图像进行成像系统质量评估，分为两个部分，其中，

第一部分：计算稀疏投影下的锥束XCT成像系统的调制传递函数MTF

①选取多个三维重建结果切片，沿z轴方向叠加平均三维重建结果切片得到二维图像p(x,y)，(x,y)为图像像素在直角坐标系下位置，并将其从直角坐标系映射到极坐标系，得到极坐标系下图像p(r,θ)；

②对图像p(x,y)采用阈值分割将灰度图像p(x,y)映射为二值图像p^norm(x,y)，这里为图像中空气部分的像素平均值，为图像中成像质量评估仿体部分的像素平均值，将仿体和空气部分分割，识别出目标体区域，然后通过质心法确定目标体区域的中心；

③设定过采样像素大小Δr，由此形成半径分别为Δr，2Δr，3Δr，4Δr……的一系列同心圆环，将图像p(r,θ)中各个同心圆环范围内的像素点进行像素值平均，而后按照离开圆心的距离排序，由此得到过采样边缘扩散函数曲线ESF；

④对ESF求微分以后得到线扩散函数LSF，对LSF加hann窗后进行傅里叶变换，所得结果取模并进行零频率归一化后即为所求的调制传递函数MTF(u),这里的u为频率；

第二部分：求取稀疏投影下的锥束XCT成像系统的噪声功率谱NPS：

⑤对于所选取的三维重建结果切片，分别在每一个切片的圆形区域内的同一区域截取方形ROI，计算该方形ROI的平均像素值，而后将原ROI图像减去该平均像素值，以此获得各个切片的噪声图像；

⑥将每一幅噪声图像都分成互不重叠的几个子块ROI，而后分别对所有子块ROI计算噪声功率谱，计算噪声功率谱时采用的傅里叶变换的点数与计算MTF时采用的傅里叶变换的点数一致，再将所得的噪声功率谱叠加平均，得到锥束XCT成像系统的二维噪声功率谱；

⑦从二维噪声功率谱图中心沿横轴向右至终点，取横轴上下各7行功率谱值，而后将这14行功率谱值叠加平均，得到锥束XCT成像系统的一维噪声功率谱NPS_d(u)；

⑧根据公式计算得到噪声等效量子数NEQ，即得到成像质量空间分辨率和噪声水平的综合评价结果。