一种在有限视角断层成像中消除截断伪影的方法

1.本发明涉及图像处理领域，具体涉及一种在有限视角断层成像中消除截断伪影的方法。


背景技术：

2.传统的ct(computed tomography，计算机断层扫描)成像技术利用锥形x射线发射器绕目标组织180
°
，每隔一定的角度发射一次x射线，并通过x射线接收器在目标组织的另一侧接收信号，从而得到若干组投影数据，根据目标组织的不同部分对x射线的吸收或衰减系数的不同这一特性来重构目标组织的内部结构。
3.在实际使用时，由于目标组织结构的差异和应用环境的制约，通常难以对目标组织实现全角度范围的扫描，因此在实际应用时常采用有限视角ct成像技术。与需要采集180
°
完整数据集的传统ct技术不同，有限视角ct成像技术通常基于从15
°
到60
°
的有限角度范围内获得的投影数据重建图像。传统的有限角度ct重建方法，包括迭代算法和非迭代算法。非迭代算法包括bp(back projection，反投影)算法，和fbp(filtered back projection，滤波反投影)算法。bp算法的重建速度快，但是重建后的图像质量差。fbp算法在bp算法的基础上使用了滤波器用于滤除bp重建过程中的噪声信息，可以在计算速度快的同时提供令人满意的图像质量，但是该方法受成像的几何形状的限制，比如对于某些利用非常窄的投影角范围、非均匀角度采样方案和固定探测器的系统中，fbp的计算更复杂，计算误差也会降低图像的质量。迭代类算法，比如sart(simultaneous algebraic reconstruction technique，联合代数重建算法)，和基于统计模型的算法比如ml-convex(maximum likelihood convex，最大似然凸函数)算法，可以适应有限视角的断层扫描而不局限于几何形状，而且能有效增强目标组织图像的边缘和对比度。随着计算机硬件的升级和并行计算效率的提高，迭代算法的重建速度也有了一定的提高。
4.但是由于探测器的尺寸有限，且光源为点光源，所以部分大角度投影时接收器无法接收完整的投影数据，对应的pvs(projection views，投影视图)发生截断，更新发生截断的pvs后，fov(field of view，视场)外的体素值没有更新，导致边缘体素值不连续，如果不进行校正，重建后的图像将会出现tpa(truncated projection artifacts，截断投影伪影)。tpa显然影响重建图像的质量。


技术实现要素：

5.发明目的：本发明所需要解决的技术问题是针对有限视角ct成像结果在边缘处出现不连续的条状截断伪影的问题，提供一种在有限视角断层成像中消除截断伪影的方法。
6.为实现上述目的，本发明具体步骤如下：
7.步骤1，在目标组织上、下分别放置x射线发射器和接收器；
8.步骤2，在有限角度内使用x射线发射器绕组织移动并发射x射线，同时接收器采集原始投影数据；
9.步骤3，对原始投影数据进行重建，获得重建图像；
10.步骤4，对重建图像进行投影变换，获得完整投影数据；对所述完整投影数据取梯度信息；
11.步骤5，利用所述梯度信息对原始投影数据进行补全，获得补全后的投影数据；
12.步骤6，对补全后的投影数据进行重建，获得消除截断伪影的有限角度ct重建图像。
13.本发明中，优选地，步骤4中对sart重建图像进行投影变换，获得完整投影数据；所述完整投影数据未发生截断。对所述完整投影数据取梯度信息包括：
14.根据完整投影数据和原始投影数据的相对位置，将完整投影数据和原始投影数据中心对齐，对完整投影数据提取所需要的部分，该部分由左、右两块构成，分别为原始投影数据相对完整投影数据的截断处至完整投影数据两边的非零处，记为左侧重建投影数据f
l
和右侧重建投影数据fr；
15.对左侧重建投影数据f
l
和右侧重建投影数据fr取x方向梯度信息，左侧重建投影数据f
l
中，像素点(x
l
,yn)在x方向的梯度值g
l
(x
l
,yn)为：
16.g
l
(x
l
,yn)＝f
l
(x
l 1
,yn)-f
l
(x
l
,yn)
17.对右侧重建投影数据fr，像素点(xr,yn)在x方向的梯度值为：
18.gr(xr,yn)＝fr(xr,yn)-fr(x
r-1
,yn)
19.其中，l表示左侧重建投影数据f
l
中x方向的位置索引，r表示右侧重建投影数据fr中x方向的位置索引，n表示y方向的位置索引。
20.通过对重建图像进行投影变换，可以人为扩大接收器的面积，使投影数据不再发生截断，从而从根本上避免了伪影的产生。并且用此方法得到的投影图与原始投影数据具有更多的共性，可以更好地对原始投影数据进行补全。
21.本发明中，优选地，所述步骤5中，根据左侧重建投影数据f
l
和右侧重建投影数据fr的梯度信息，从内侧向两边对原始投影数据的缺失部分进行补全，获得补全后的投影数据。
22.本发明中，优选地，所述步骤5包括：
23.记补全后的投影数据为f0，补全后的投影数据f0中，中间为原始投影数据；由中间向两边对原始投影数据进行补全，对两边需要补全的部分分别进行如下迭代：
24.对左侧区域第yn行数据，有：
25.f0(x
l 1
,yn)＝f0(x
l
,yn) g
l
(x
l
,yn)
26.其中，记l＝l0处为原始投影数据，根据上述公式得到补全值再根据上述公式得到的值，以此类推，直到迭代结束；
27.对右侧区域第yn行数据，有：
28.f0(x
r-1
,yn)＝f0(xr,yn)-gr(xr,yn)
29.其中，记r＝r0处为原始投影数据，根据上述公式得到补全值再根据上述公式得到的值，以此类推，直到迭代结束；
30.对每行做上述左侧区域和右侧区域的迭代，完成对原始投影数据的补全。
31.在步骤5中，采用梯度信息补全投影数据，避免了投影数据截断处的不连续我们采
用的数据补全方法和现有的投影数据扩展的方法相比，可以更好地恢复出截断的投影数据，恢复出的数据也更接近真实值，也使得最终可以重建出良好的图像。在消除伪影的作用上，和现有的使用像素补偿的方法相比，本技术重建出来的图像更加的真实，可以完全的保留被伪影覆盖的信息，同时不会在图像的边缘造成锯齿状的伪影。
32.本发明中，优选地，还包括步骤7，对步骤4到步骤6进行一次以上的循环，获得消除截断伪影的有限角度ct重建图像，在循环过程中步骤4输入的重建图像为上一循环中步骤6获得的消除截断伪影的有限角度ct重建图像。当关注位置较低的切片时，步骤1至步骤6就能够有效地消除图像的边缘处条状截断伪影；当关注位置较高的切片时，执行步骤7对步骤4到步骤6进行多次循环，能够有效地减少高层图像的边缘处条状截断伪影。
33.本发明中，优选地，所述步骤1中x射线发射器应位于接收器的上方，使x射线束尽可能地覆盖目标组织。
34.本发明中，优选地，所述步骤2中，以接收器平面内的某一直线为旋转轴，x射线发射器绕此轴在垂直于接收器平面的运动平面内做有限角度的圆弧形移动，角度范围由发射器左右移动的边界点确定的圆心角表示，其中左右边界点关于轨迹的中心对称。采集过程中，x射线发射器从初始端点向另一运动端点移动，每移动一定的角度向接收器发射x射线，接收器接收后视为采集到一组投影数据。扫描过程即步骤2结束，采集到若干组原始投影数据。
35.本发明中，优选地，所述步骤3中使用sart重建算法对重建区域进行重建，使最终结果满足最优化条件，获得第一次有限视角ct成像结果，此时的成像结果在边缘处存在条状截断伪影。
36.本发明中，优选地，所述步骤6中，对补全后的投影数据使用sart算法重建，获得消除截断伪影的有限角度ct重建图像。
37.有益效果：
38.本发明提出的一种在有限视角断层成像中消除截断伪影的方法，以传统的sart重建算法为基础，利用重建图像投影值的梯度信息，对有截断的原始投影数据进行补全，对补全后的投影数据进行sart重建，得到消除截断伪影的有限视角ct重建结果，解决了有限视角下ct图像重建因接收器无法接收到所有x射线而造成真实投影数据缺失，进而导致的重建图像边缘处存在条状截断伪影的问题。
附图说明
39.下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。
40.图1是本发明方法的流程图。
41.图2是本发明的原始投影数据采集过程的示意图。
42.图3是本技术实施例中原始投影数据。
43.图4是本技术实施例中对原始投影数据做初次sart算法重建得到的图。
44.图5是本技术实施例中初次重建图的投影数据图。
45.图6是本技术实施例中补全过程的示意图。
46.图7是本技术实施例中补全后的投影数据图。
47.图8是本技术实施例中对补全后投影数据做sart算法重建得到的图。
具体实施方式
48.下面将结合附图，对本发明的实施例进行描述。
49.如图1所示，本发明提出的一种在有限视角断层成像中消除截断伪影的方法，包括以下步骤：
50.步骤1，在目标组织上、下分别放置x射线发射器和接收器；
51.步骤2，在有限角度内使用x射线发射器绕组织移动并发射x射线，同时接收器采集获得原始投影数据；
52.步骤3，对原始投影数据进行重建，获得重建图像；
53.步骤4，对重建图像进行投影变换，获得完整投影数据；对所述完整投影数据取梯度信息；
54.步骤5，利用所述梯度信息对原始投影数据进行补全，获得补全后的投影数据；
55.步骤6，对补全后的投影数据进行重建，获得消除截断伪影的有限角度ct重建图像。
56.本实施例中，步骤1，x射线接收器长23.04cm，宽19.20cm，每个像素的大小为一个边长为0.1mm的正方形，因此每次接收到的投影数据为一张分辨率为1920
×
2304的平面图。将目标组织放置于x射线发射器和接收器之间靠近接收器平面的位置，物体底面与接收器平面之间距离为2cm，x射线发射器到旋转支点的距离为64cm。x射线发射器、接收器与物体的摆放位置如图2所示。
57.本实施例中，步骤2，如图2所示，圆弧形运动轨迹中心位于点p，旋转支点p
′
为点p在旋转轴上的投影，x射线发射器初始位置位于运动轨迹的左端点l，运动终点为右端点r。在运动平面内，x射线发射器以点p
′
为旋转中心，绕该点作弧形运动，该弧的圆心角为60
°
，左端点l和右端点r以pp
′
为轴对称，与中轴各成30
°
。x射线发射器从左端点l出发向右端点r做圆弧形运动，每隔3
°
向x射线接收器发射一次x射线，接收器接收到的数据作为该角度的一组投影值。当发射器到达点r，扫描过程结束，共采集到21组投影数据，本实施例中的原始投影数据如图3所示。
58.本实施例中，步骤3，使用sart重建算法对采集到的原始投影数据进行初次重建。sart算法通过迭代求这样一类方程的解x:
59.ax＝b
ꢀꢀꢀ
(1)
60.公式(1)中，a为m
×
n矩阵，在图像重建过程中是投影变换的系统矩阵，x是n维列向量，由需要重建的体素的值构成，代表目标组织中每个体素对x射线的衰减系数，n为重建体素的数目，b是m维列向量，代表投影值，m是采集到的投影像素的总数目。在有限角度ct成像时，一般m的值远小于n，因此该方程有无数组解。sart算法首先对x给定一个初始的估计值x
′
，然后对这个估计值进行投影得到投影向量b的一个估计值b
′
，再用第一个角度下测量得到的投影值b1与b
′
之差b
1-b
′
作为修正值修正x
′
。对所有的21组投影值完成一次对x
′
的修正作为一次迭代，最终x
′
收敛到使最小的x，即最小二乘解。根据sart算法理论，加权矩阵w的引入可以有效减少伪影。
61.所以需要求使如下损失函数l(x)最小的x的值，如公式(2)所示：
[0062][0063]
利用梯度下降法求得使l(x)最小的x即得到sart算法的迭代公式：
[0064][0065]
其中，表示第k 1次迭代后第j个体素的值，i为所取的投影值的序号，λ为松弛因子，通过改变λ影响收敛速度和重建结果。j表示每条x射线穿过的体素个数，1≤j≤j；m表示此次迭代过程中第j个体素被更新的次数，mi表示此次迭代过程中第j个体素被更新的总次数，1≤m≤mi，a
mj,i
表示投影值序号为i时，第j个体素在本次迭代过程中第m次更新时的权重。为保证收敛，本例中λ取0.1，迭代5次。为了加快运算，将原始投影数据等比例缩小分辨率为480
×
576。
[0066]
如图4所示，得到初次重建结果，本例中重建结果的尺寸为480
×
576
×
100，重建图像是三维图像，将三维图像切割成100个二维图像处理。图4可以看到在边缘处存在明显的截断投影伪影。
[0067]
本实施例中，步骤4，对初次重建结果做投影变换，得到完整的投影数据图，如图5所示。对初次重建结果做投影变换使用本领域技术人员公知的现有技术计算得到，可以由计算机完成，所以对比步骤2中的扫描阶段，没有了接收器大小的限制，得到完整投影数据，因此原始投影数据与重建图像的投影数据的相对位置是已知的。将新的投影图和原始投影数据中心对齐，对重建图像的投影数据提取所需要的部分，该部分由左、右两块构成，为原始投影图相对重建投影图的截断处至重建投影图两边的非零处，对该部分取x方向梯度信息。如图6所示。梯度的计算方法为，某一像素点的梯度值表示为该点的像素值减去其右方一点的像素值。左侧重建投影数据f
l
中，点(x
l
,yn)在x方向的梯度值为：
[0068]gl
(x
l
,yn)＝f
l
(x
l 1
,yn)-f
l
(x
l
,yn)
ꢀꢀꢀ
(4)
[0069]
公式(4)中g
l
表示左区域中像素的梯度信息，是一个标量，下标l表示在左侧区域。
[0070]
对于右侧重建投影数据fr，点(xr,yn)在x方向的梯度值为：
[0071]gr
(xr,yn)＝fr(xr,yn)-fr(x
r-1
,yn)
ꢀꢀꢀ
(5)
[0072]
公式(5)中gr表示右区域中像素的梯度信息，是一个标量，下标r表示在右侧区域。
[0073]
本实施例中，步骤5，得到梯度信息后，由中间向两边对原始投影数据进行补全。重建区域f0中，中央为原始投影数据，对两边需要补全的部分分别进行如下迭代：
[0074]
对左侧区域第yn行数据，有：
[0075]
f0(x
l 1
,yn)＝f0(x
l
,yn) g
l
(x
l
,yn)
ꢀꢀꢀ
(6)
[0076]
其中,l＝l0处为原始的投影数据，根据公式(6)得到补全值再根据公式(6)得到的值，以此类推，直到迭代结束。
[0077]
对右侧区域第yn行数据，有：
[0078]
f0(x
r-1
,yn)＝f0(xr,yn)-gr(xr,yn)
ꢀꢀꢀ
(7)
[0079]
其中,r＝r0处为原始的投影数据，根据公式(7)得到补全值再根据公
式(7)得到的值，以此类推，直到迭代结束。
[0080]
重建区域f0的每行做上述迭代，完成对原始投影数据的补全。本例中，补全后的投影图如图7所示。
[0081]
本实施例中，步骤6，对补全后的投影数据做sart算法重建，得到新的重建图像。
[0082]
本实施例中，还包括步骤7，对上述步骤4到步骤6进行多次循环，以减少高层部分图片的截断伪影，使内部组织的轮廓更清晰。在循环过程中步骤4输入的sart重建图像为上一循环中步骤6获得的消除截断伪影的有限角度ct重建图像。根据截断伪影的残留情况，当关注位置较高的切片时，本实施例循环3次，最终得到消除边缘处截断投影伪影的有限视角ct重建结果，如图8所示。
[0083]
本发明提供了一种在有限视角断层成像中消除截断伪影的方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的具体实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。