一种基于三维加权滤波反投影和统计迭代的图像重建方法与流程

本发明属于图像处理技术领域，更为具体地讲，涉及一种基于三维加权滤波反投影和统计迭代的图像重建方法。

背景技术：

计算机断层成像技术(computedtomogrophy，ct)是一种通过对物体进行不同角度的射线投影测量以获取物体横截面信息的技术，被广泛应用于工业无损检测、医学诊断、地球物理以及农林业等领域。经过多年的发展,ct的扫描方式己经从二维平行束扫描发展到三维锥束扫描，三维锥束ct(conebeamct，cbct)。与二维ct相比，具有射线利用率高、纵向分辨率高、扫描时间短、部分容积效应低等优点,成为目前ct重建领域的研究热点。

重建算法是ct中十分重要的问题，它的本质是利用某种重建算法从投影数据解算出成像断层平面上各个像素点的衰减系数。通过几十年的发展，基于滤波反投影的解析法己逐步成熟，目前的ct重建技术可以对无噪声的完备的投影数据以较快的速度完成重建。在cbct中，由于圆轨迹的短扫描方式不满足完全精确重建条件，基于的fdk算法的短扫描重建算法存在锥束伪影。目前改善锥束伪影的方法主要有以下两种：采用“圆弧+圆弧”、“圆弧+线段”、“圆弧+螺旋”等扫描方式以满足完备重建条件，利用精确重建算法来进行图像重建，但这种方法要求轨迹精度高，增加了扫描的机械设计难度；或者采用改进的fdk算法，很多学者提出了基于fdk的改进算法，这些改进方法的中心思想是通过一些操作来减少锥角增大时不精确重建形成的伪影，这些算法对锥角较小情况有较高的重建质量，但当锥角大于5度时，改善效果有限。

另一种重建算法是迭代重建法，能够有效抑制噪声，迭代重建算法又分为代数迭代算法和统计迭代算法。在迭代中，根据计算得到的正投影值与实际投影数据之间的偏差更新图像。可以在采样不规则和数据缺失的情况下重建出图像。但是迭代法中需要进行多次的正投影和反投影的计算，成像速度较解析重建算法慢很多，限制了其在实际中的应用。

技术实现要素：

本发明的目的在于克服现有技术的不足，提供一种基于三维加权滤波反投影和统计迭代的图像重建方法，通过三维加权和统计迭代实现图像重建，可以在使用较少投影角度的情况下，有效解决锥束伪影问题，抑制噪声，而且提高重建效率。

为实现上述发明目的，本发明一种基于三维加权滤波反投影和统计迭代的图像重建方法，其特征在于，包括以下步骤：

(1)、获取投影数据并进行标记；

(1.1)、以待重建物体的中心为原点建立三维物体坐标系，将待重建物体分成一个个体素，单个体素表示为xi，对应坐标设为(xi,yi,zi)，i＝1,2,…,n表示待重建物体分成单个体素的总个数；

(1.2)、获取待重建物体的投影数据，计算第j个投影数据在投影时射线源相应于物体坐标系的旋转角度，即投影角度ηj、投影射线相对于当前投影角度中心射线的扇角γj、以及相对于中心平面的锥角αj，j＝1,2,…,m，m为投影数据总个数，最后将每一组投影数据标记为pj(ηj,γj,αj)。

(2)、对投影数据进行预加权处理：

p'j(ηj,γj,αj)＝cosγjpj(ηj,γj,αj)；

(3)、将加权的投影数据通过两个滤波核进行滤波处理，滤波后分别得到：

(4)、三维加权；

(4.1)、

记投影角度总数为l，则穿过体素xi的射线对应的投影数据一共l个，按照投影角度从小到大的方式排列，将对应的投影角、扇角、锥角记为(ηil,γil,αil)，l＝1,2,...l；其中ηi1＜ηi2＜…＜ηil；

(4.2)、获取共轭射线及其锥角

将任意一单体素xi与投影角度为ηil时源点s之间的连线sxi在扫描轨迹内的垂直投影的延长线与扫描轨迹的交点记为s'，则单体素xi与交点s'的连线s'xi为sxi的共轭射线；连线sxi相对于中心平面的锥角为αil，则共轭射线s'xi相对于中心平面的锥角记为α'il；

(4.3)、计算投影角度为ηil时xi处的三维加权因子为：

(5)、还原初始重建图像；

(5.1)、对单体素xi进行反投影，得到重建后的单体素f(xi)；

其中，r是扫描半径，l²(ηil,xi)是在ηil投影角度下，连接x射线源s到xi的线段sxi在中心平面垂直投影的长度；

(5.2)、将重建后的所有单体素按照三维物体坐标系还原，得到待重建物体的初始重建图像；

(6)、通过统计迭代方式修正重建图像；

(6.1)、对初始重建图像进行投影，得到投影估计值；

(6.1.1)、将重建后的所有单体素f(xi)排列成列向量f^(k)＝(f(x1)^k,f(x2)^k,…,f(xn)^k)^t，将所有的投影数据排列成列向量p＝(p1,p2,…,pm)^t，其中，k表示迭代次数；

(6.1.2)、利用siddon算法计算系统矩阵w＝(wji)m×n，其中，元素wji代表第i个重建体素f(xi)对第j个投影数据pj的贡献值；

对初始重建图像进行正向投影，得到投影估计值p'：

p'＝wf^(k)

(6.2)、计算实际投影值和投影估计值之间的相对偏差：

p^rel＝p./p'

(6.3)、根据相对偏差修正初始重建图像；

(6.3.1)、计算归一化：

(6.3.2)、将投影的相对偏差反投影到图像空间，得到：f^back＝w^tp^rel；

(6.3.3)、更新重建图像：

其中，是f^back的第i个元素；

(6.4)、重复步骤(6.1)-(6.3)，直到相邻两次迭代的相对偏差小于10^-3，则停止迭代，并将最后一次迭代后更新的图像作为满足条件的重建图像。

本发明的发明目的是这样实现的：

本发明基于三维加权滤波反投影和统计迭代的图像重建方法，先获取投影数据，获得每个投影数据对应射线的投影角度、锥角以及扇角，然后对投影数据进行预加权，再对预加权后的投影数据进行逐行滤波，然后进行三维加权，进行反投影得到初始的重建图像。再根据得到的重建图像，计算所有角度下的投影估计值，与实际的投影数据进行比较，计算误差，根据误差修正重建图像，然后再重复迭代此过程，直到满足条件的重建图像。

同时，本发明基于三维加权滤波反投影和统计迭代的图像重建方法还具有以下有益效果：

(1)、由于迭代法对数据完备性要求低，所以在三维加权滤波反投影中，所需的投影角度只要能满足迭代法的重建要求，比单独使用解析法所需的投影少，使得扫描时间更快，x光照射剂量降低。

(2)、在三维加权滤波反投影中，采用了一种改进的三维加权权重，相较之前的加权函数，除了考虑锥角的影响，还进一步考虑了扇角γ的影响，更有效的消除了重建的伪影问题。

(3)、将三维加权滤波反投影重建的图像作为迭代初值，减少了迭代所需次数，大大节省了重建时间，有效的抑制了噪声，结合了两种方法的优点，在重建质量与重建效率上都得到了令人满意的效果。

附图说明

图1是本发明基于三维加权滤波反投影和统计迭代的图像重建方法流程图；

图2是锥束圆轨迹扫描几何示意图；

图3是共轭射线几何示意图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

图1是本发明基于三维加权滤波反投影和统计迭代的图像重建方法流程图。

在本实施例中，如图1所示，本发明一种基于三维加权滤波反投影和统计迭代的图像重建方法，包括以下步骤：

s1、获取投影数据并进行标记；

s1.1、以待重建物体的中心为原点建立三维物体坐标系，将待重建物体分成一个个体素，单个体素表示为xi，对应坐标设为(xi,yi,zi)，i＝1,2,…,n表示待重建物体分成单个体素的总个数；

s1.2、在本实施例中，如图2所示，获取待重建物体的投影数据，计算第j个投影数据在投影时射线源相应于物体坐标系的旋转角度，即投影角度ηj、投影射线相对于当前投影角度中心射线的扇角γj、以及相对于中心平面的锥角αj，j＝1,2,…,m，m为投影数据总个数，

其中，投影角度ηj为已知量，通过扫描可以直接读出旋转角度，而扇角γj和锥角αj可以通过如下公式计算：

其中，r为扫描半径，xi、yi为xi的x、y轴坐标。

对于任意一个投影数据，其所对应的η，γ，α组合是惟一的，故每一组投影数据可以表示为pj(ηj,γj,αj)。

s2、对投影数据进行预加权处理，具体为投影数据乘以扇角的余弦值，获得加权后的数据：

p'j(ηj,γj,αj)＝cosγjpj(ηj,γj,αj)；

s3、将加权的投影数据通过两个滤波核进行滤波处理，此处两个滤波核为：

通过两个滤波核滤波后分别得到：

s4、三维加权；

s4.1、记投影角度总数为l，则穿过体素xi的射线对应的投影数据一共l个，按照投影角度从小到大的方式排列，将对应的投影角、扇角、锥角记为(ηil,γil,αil)，l＝1,2,...l；其中ηi1＜ηi2＜…＜ηil；

s4.2、获取共轭射线及其锥角

如图3所示，将任意一单体素xi与投影角度为ηil时源点s之间的连线sxi在扫描轨迹内的垂直投影的延长线与扫描轨迹的交点记为s'，则单体素xi与交点s'的连线s'xi为sxi的共轭射线；连线sxi相对于中心平面的锥角为αil，则共轭射线s'xi相对于中心平面的锥角记为α'il；

s4.3、计算投影角度为ηil时xi处的三维加权因子为：

s5、还原初始重建图像；

s5.1、对单体素xi进行反投影，得到重建后的单体素f(xi)；

其中，r是扫描半径，l²(ηil,xi)是在ηil投影角度下，连接x射线源s到xi的线段sxi在中心平面垂直投影的长度；

s5.2、将重建后的所有单体素按照三维物体坐标系还原，得到待重建物体的初始重建图像；

s6、通过统计迭代方式修正重建图像；

s6.1、对初始重建图像进行投影，得到投影估计值；

s6.1.1、将重建后的所有单体素f(xi)排列成列向量f^(k)＝(f(x1)^k,f(x2)^k,…,f(xn)^k)^t，将所有的投影数据排列成列向量p＝(p1,p2,…,pm)^t，其中，k表示迭代次数；

s6.1.2、利用siddon算法计算系统矩阵w＝(wji)m×n，其中，元素wji代表第i个重建体素f(xi)对第j个投影数据pj的贡献值；

对初始重建图像进行正向投影，得到投影估计值p'：

p'＝wf^(k)

s6.2、计算实际投影值和投影估计值之间的相对偏差：

p^rel＝p./p'

s6.3、根据相对偏差修正初始重建图像；

s6.3.1、计算归一化：

s6.3.2、将投影的相对偏差反投影到图像空间，得到：f^back＝w^tp^rel；

s6.3.3、更新重建图像：

其中，是f^back的第i个元素；

s6.4、重复步骤s6.1-s6.3，直到相邻两次迭代的相对偏差小于10^-3，则停止迭代，并将最后一次迭代后更新的图像作为满足条件的重建图像。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。