CT图像重建反投影过程中的线性插值方法与流程

本发明涉及一种ct图像重建反投影过程中的线性插值方法。

背景技术：

反投影(backprojection)是ct图像重建算法中的主要过程之一，现有的反投影方法大体上分为两类，一类是基于射线的反投影方法(ray-drivenbackprojection)，一类是基于像素的反投影方法(pixel-drivenbackprojection)，但实践中，从简单性和便于实现等方面考虑，多使用基于像素的反投影方法。

传统的基于像素的反投影方法是假定给定的像素点没有大小尺寸，只是纯粹的一个点，考虑x-光源点和给定像素点p的连线到探测器的交点通常可能会介于两个探测器单元之间，为了获得该像素点对应于探测器上的数据，通常需要使用插值方法，其中最常用的两种插值方法是最近邻(nearestneighbour)插值法和线性(linear)插值法，其中最近邻插值法以与该交点最近邻的探测器单元γi的值g(γi)作为该像素点的值(其计算原理可参见图3)，线性插值法则以与该交点相邻的两个探测器单元γi和γi+1的值g(γi)和g(γi+1)的加权平均值作为该像素点p的值gp。这两种插值方法最重要的一个特点是不考虑像素点的大小，其结果是光源和像素点的连线最多只能落在两个探测器单元之间，噪声特性无法满足人们日益提高的要求。

为了改进反投影技术，美国专利7,227,982提出了一种新的基于距离的反投影技术。相对于过去的基于像素的反投影技术，这种技术不再认为像素点没有大小，而是假定像素点是一个边长为δ的方块。假定需要重建的ct图像的视野大小为fov(mm)，这个视野中水平(x-方向)和竖直(y-方向)两个方向各有n个像素，则方块边长δ的计算方法为：

公式(1)表明方块的长度实际上就是图像空间一个网格的边长。这种技术将像素边界和探测器单元投影到一个公共轴上，使用这两个投影的重叠部分作为反投影的权重，这种反投影法得到的ct图像信噪比(signal-to-noiseratio)较高，即图像的噪声较低，但该方法实现比较复杂，实施起来有较大的难度，影响实践中的推广使用。

同样为了改进反投影技术，美国专利8,116,426中提出了另一种方法。和美国专利7,227,982相似的想法，这种方法中也是考虑了像素点的大小，对于有n×n个像素点、视野为fov(mm)的图像，将每个像素点看作是半径为δs的圆盘，δs用下式计算：

由于该方法与本发明密切相关，在此详细描述此方法(其计算原理可参见图2)。

图2中，对于给定的像素点p有坐标(x，y)，将其视为半径为δs的圆盘，p点所在的位置为该圆盘的圆心。圆盘半径δs为像素大小的一半，由公式(2)计算得到。在反投影过程中，该方法考虑到这个半径为δs的圆盘投影到探测器上所占的通道数目。

首先需要使用下式计算出p点在探测器上对应的通道号γ：

c＝int(γ)(4)

公式(3)中r是光源所在位置s到旋转中心o的距离，即r＝|so|，β是线段so与y-轴所成的角度。公式(4)是将公式(3)得到的浮点通道值γ取为整数。要注意的是公式(3)是以so所对应的探测器通道号为0(图2中线段so延长线与探测器弧交点处编号为0)来进行编号的，所以通道编号有正的编号和负的编号之分。如果采用其它的编号方法需要对公式(3)进行相应的调整则应依据编号规则与公式(3)体现的编号规则的对应关系对公式(3)进行调整，例如，通常可以在公式(3)的右侧加上一个起始号γ0，γ0为so所对应的探测器通道号，或者在不加起始号的情形下，将γ视为相对于线段so延长线与探测器弧交点处编号的相对通道号，即实际通道号减去γ0的相对通道号作为公式(3)的通道编号。

其次考虑半径为δs的圆盘投影到探测器上所占通道的数目。

使用γi和γj来表示圆盘投影到探测器上的边界通道值，位于γi和γj之间的通道均属于该圆盘的投影区域。使用下式计算γi和γj

公式(5)和(6)中的l表示光源s到p点的距离，即l＝|sp|，其计算公式为

定义

δc＝γj-γi(8)

δc表示半径为δs的圆盘投影到探测器上所占通道的数目。

为了得到该圆盘所对应的探测器数据值gp，使用g(·)来表示某个探测器通道所对应的值，采用如下方法进行插值得到gp：

ci＝int(γi)(9)

cj＝int(γj)(10)

wi＝ci+1-γi(11)

wj＝γj-cj(12)

如果ci＝c＝cj，则

gp＝g(c)(13)

否则

公式(9)和(10)是将公式(5)和(6)计算出的浮点通道号γi和γj转化成整数通道号值。公式(11)和(12)计算出的是边界通道γi和γj在加权平均中的权重。

相对于其他现有技术而言，美国专利8,116,426的图像重建方法将像素点看成是半径为网格大小一半的圆盘，实现过程比较简单，易于使用，但是该像素模型弹性不够。首先是该方法将像素点看成是半径为网格大小一半的圆盘，网格大小一旦固定，圆盘的半径也就固定了(参见公式(2))，缺乏弹性，因此在视野较大而像素数目相对不变的情况下就表现为像素尺寸太大，由于像素尺寸太大，投影到探测器上表现为覆盖的通道数目太多，从而会造成重建出来的图像虽然噪声低，但分辨率太低，细节丢失；其次是该方法在网格尺寸变小时会退变成最近邻插值法，这一点可以从极限情况δs＝0看出，当δs＝0时，公式(5)和(6)给出γi＝γj＝γ，从而ci＝c＝cj，而从公式(4)和(15)看出这实际上就是最近邻插值法，因此在像素尺寸变小的情况下使用这种方法得到的图像会存在较大的噪声。而像素尺寸变小最起码有两种情况，从公式(2)可以看到，第一种情况是如果图像的像素点数目n不变，则视野大小fov变小时会造成像素尺寸变小；第二种情况是如果视野大小fov不变，而像素点数目n变大，则像素尺寸也会变小，而第二种情况目前正成为ct图像重建中的一种趋势，即使用大矩阵(如1k或2k矩阵)重建ct图像，因此，美国专利8,116,426难以适应于ct图像重建技术的发展要求。

技术实现要素：

为解决上述技术问题，本发明提供了一种ct图像重建反投影过程中的线性插值方法，该方法简便易行，且具有良好的弹性，在像素尺寸较小的情况下不会导致过高的噪音。

本发明的技术方案是：一种ct图像重建反投影过程中的线性插值方法，其对于像素点p，定义一个以像素点p为圆心的半径可变的圆盘，所述圆盘的半径δs′依据下列公式计算：

其中

ρ≥0，为设定的半径系数，fov为视野范围，n为fov范围内的像素点数量，

计算并获得该圆盘投影到探测器后覆盖的探测器通道边界，依据投影所覆盖的各探测器通道所对应的值g(·)，以下列公式计算像素点p的值gp，

如果c′i＝c′j，则

否则(即c′i+1≤c′j)

其中

c′i＝int(γ′i)(22)

c′j＝int(γ′j)(23)

w′i＝γ′i-c′i(24)

w′j＝γ′j-c′j(25)

其中，

x和y分别为像素点p的x-轴y-轴坐标；

r是光源所在位置s到旋转中心o的距离，即r＝|so|；

β是线段so与y-轴所成的角度；

l为光源s到p点的距离，即l＝|sp|；

γ为像素点p对应的相对通道号计算值，相对通道号计算值为γ的通道称为通道γ；

γi和γj为圆盘投影到探测器后覆盖的边界相对通道号计算值，介于γi和γj之间的通道(如果有的话)均为圆盘投影到探测器后覆盖的通道，相对通道号计算值为γi和γj的通道分别称为通道γi和通道γj；

c′i和c′j分别为对γi和γj取整数；

wi和wj分别为边界通道γi和γj的权重系数，位于γi和γj之间的通道(如果有的话)为全部覆盖的通道，权重系数均为1。

可以依据图像视野和像素点数目设定或改变半径系数ρ。

通常，可以设定ρ＞0。

可以设定0＜ρ＜1，相应的圆盘半径在0到图像网格尺寸(边长)一半之间。

可以设定ρ＞1，相应的圆盘半径超过图像网格尺寸的一半，不同像素点的圆盘存在重叠区域。

所称相对通道号为相对于线段so延长线与探测器弧交点处通道编号的通道号，等于以线段so延长线与探测器弧交点处通道编号为零情形下的通道号或者以实际通道编号减去线段so延长线与探测器弧交点处的通道编号。

本发明的有益效果是：将像素点看成是-个半径可变的圆盘，依据图像视野和像素点数目选择半径系数ρ，使其符合分辨率要求及数据处理能力，由此避免了将像素点视为矩形导致的复杂处理过程，同时使图像重建具有很好的信噪比和分辨率，对不同情形具有良好的弹性和适应性，当像素尺寸变小时，插值方法退化成为近似于线性插值法，而不会像美国专利8,116,426那样退化成近似于最近邻插值，保证了在使用大矩阵(如1k或2k矩阵)重建ct图像时具有更好的信噪比，能够更好地适用于ct图像重建技术的发展。

附图说明

图1是反投影模型和计算原理示意图，该图所示的圆盘半径小于图像网格尺寸(边长)一半；

图2是反投影模型和计算原理示意图，该图所示的圆盘半径等于图像网格尺寸(边长)的一半；

图3是反投影模型和计算原理示意图，该图对应于圆盘半径为零时的极限情形，即不考虑像素点大小的情形，可以将该模型视为圆盘半径为零(无限趋近于零)的极限情形；

图4是在相同情形下反投影过程中不考虑像素尺寸大小(现有技术)和考虑像素尺寸大小(本发明)的图像质量比较，其中左侧图像为不考虑像素尺寸大小的图像，右侧图像为考虑像素尺寸大小的图像；

图5是在相同情形下反投影过程中选择不同半径系数的图像质量比较，其中左侧图像为半径系数为1(即圆盘半径等于图像网格边长一半)的图像，右侧图像为半径系数为0.5(即圆盘半径等于图像网格边长四分之一)的图像。

具体实施方式

参见图1-3，为了增加像素模型的弹性，本发明在现有技术的基础上，将像素点看成是一个半径可变的圆盘，其半径为：

其中ρ≥0。

当ρ＝0时，该像素模型是将像素点看成一个没有大小的点；当ρ＝1时，像素点被看成半径为图像网格大小的一半的圆盘；当0＜ρ＜1时，圆盘的半径在0到图像网格尺寸一半之间；当ρ＞1时，则表示圆盘的半径大于网格尺寸的一半，即相邻像素点之间有重叠。图2和图3中描述的像素模型为本发明在两种情形下的具体实施方式，分别对应于现有技术下像素点被看成半径为图像网格大小的一半的圆盘和将像素点视为一个点的情形。而图1所示的具体实施方式则为0＜ρ＜1时的情形，圆盘的半径在0到图像网格尺寸一半之间。

设计了这样的像素模型后，计算该圆盘投影到探测器后覆盖的通道边界，具体公式为：

由上可知，γi≤γj，γi和γj分别代表圆盘圆盘投影到探测器上的边界通道值，对应于圆盘在探测器投影的边界位置，

为了得到p点的值(探测器数据值或图像像素值)gp，采用下列插值公式：

c′i＝int(γ′i)

w′i＝γ′i-c′i

w′j＝γ′j-c′j

如

c′i＝c′j

则

否则，即

除特别说明或另有其他明显含义，下标i、j分别表示相应参数为两个边界通道的变量。

公式(16)中，虽然有c′i＝c′j的条件，但是通常情况下wi≠wj，故而分母(1-w′i)+w′j≠1，公式(16)其实是一种线性插值的变体。公式(17)是公式(16)的推广。

现在考虑一种极限情况δs′＝0，即ρ＝0的情况，计算得到yi＝yj和c′i＝c′j，进而w′i＝w′j，在此情形下，公式(16)退化为线性插值法的计算公式：

gp＝(1-w′i)g(c′i)+w′ig(c′i+1)

通过上面的分析可以看出，当像素尺寸变得很小时，本发明的插值方法会趋近于线性插值法，而不同于美国专利8,116,426中的插值法那样趋近于最近邻插值法，因此，在像素尺寸变得很小时，本发明重建图像的噪音会明显低于美国专利8,116,426重建图像的噪音。

实验证实，本发明在反投影中通过考虑像素尺寸大小可以得到信噪比较高的图像。图4中左侧图像的反投影过程中没有考虑像素大小，即在反投影时认为像素是没有大小的物理点，右侧图像则采用考虑了像素大小的模型，是使用本发明插值方法得到的图像，右侧图像的噪声特性相比左侧图像得到了改善，表明反投影中依据本发明模型考虑像素模型可以提高信噪比。

实验还证实，本发明的像素模型具有较强的弹性，在反投影中通过调节像素尺寸大小可以得到噪声和分辨率的平衡，以获得更好的检测结果。图5显示了反投影过程中选用不同像素尺寸大小图像质量的比较，左侧图使用的像素模型是像素圆盘的半径为图像网格大小的一半(ρ＝1)，右侧图像使用的像素模型是像素圆盘的半径为图像网格大小的四分之一(ρ＝0.5)，结果显示左侧图像和右侧图像箭头所指处的分辨率是不一样的，左侧箭头处的黑色小点状结构被模糊掉了，而右侧图像则可以看到黑色小点状结构，即分辨率得到了改善。