基于广义全变分的X射线有限角CT图像重建方法和装置

基于广义全变分的x射线有限角ct图像重建方法和装置
技术领域
1.本发明涉及x射线ct图像重建技术领域，特别是关于一种基于广义全变分的x射线有限角ct图像重建方法和装置。


背景技术：

2.在工业无损检测中，常遇到电路板、机翼等板状物的检测问题，由于待测物品在某个方向具有较长的边界，可能出现以下情况：1、为了能够重建高分辨率的图像，要求射线源到物品之间的距离较小，为了避免射线源与物品之间的碰撞，扫描过程中转台无法旋转整圈；2、即使转台能够旋转整圈，但部分x射线由于穿过物体的较长边界而衰减严重，造成对应角度的探测数据信噪比差而无法使用。这些情况均导致仅能在一定的角度范围内采集得到有效的投影数据。
3.使用一定角度范围内采集的数据重建得到图像即为有限角ct图像重建问题，隶属于不完全数据图像重建问题。采用传统的图像重建算法(如fdk、sart等)得到的图像，存在部分边界严重模糊的问题，影响ct图像的可用性。
4.因此，希望有一种技术方案来克服或至少减轻现有技术的上述缺陷中的至少一个。


技术实现要素：

5.本发明的目的在于提供一种基于广义全变分的x射线有限角ct图像重建方法来克服或至少减轻现有技术的上述缺陷中的至少一个，提高有限角ct图像的质量和可用性。
6.为实现上述目的，本发明提供一种基于广义全变分的x射线有限角ct图像重建方法，其包括：
7.步骤1，利用有限角ct扫描数据集p，通过与扫描几何参数集g相关的图像重建算子rg更新估计图像u
(k)
，获得中间图像u
(k+1/2)
，k为迭代次数；
8.步骤2，按照式(1)，利用梯度的稀疏性约束中间图像u
(k+1/2)
，得到估计图像u
(k+1)
；
9.u
(k+1)
＝p(u
(k+1/2)
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
10.式(1)中，p表示约束图像梯度稀疏性的算子，由如式(2)提供的最优化问题求解得到：
[0011][0012]
式(2)中：
[0013]
为正则化项；
[0014]
为数据保真项；
[0015]
分别表示输入图像u的水平方向x和竖直方向y的梯度；
[0016]
u表示输入图像；
[0017]
diag(
·
)表示权重矩阵，为每个像素点的水平方向梯度和竖直方向梯度赋予惩罚权重；
[0018]
||
·
||2表示l2范数；
[0019]
表示l2范数的平方；
[0020]
||
·
||1表示l1范数；
[0021]
||
·
||
2,1
表示梯度模长的l1范数；
[0022]
|
·
|表示逐像素做绝对值操作；
[0023]
κ表示用于调节惩罚权重的参数；
[0024]gσ
表示高斯卷积核；
[0025]
*表示卷积运算；
[0026]
t表示转置；
[0027]
λ表示用于控制正则化项和数据保真项权重的参数；
[0028]
步骤3，判断是否达到迭代要求，若是，则终止迭代，输出估计图像u
(k+1)
，反之，则转至步骤1。
[0029]
进一步地，像素点水平方向的惩罚权重的方法包括：
[0030]
根据像素点水平方向的梯度，通过式(2)中权重矩阵diag(
·
)中的第一部分计算，得到该像素点水平方向的惩罚权重。
[0031]
进一步地，像素点竖直方向的惩罚权重的方法包括：
[0032]
根据像素点竖直方向的梯度，利用式(2)中权重矩阵diag(
·
)中的第二部分计算，得到该像素点竖直方向的惩罚权重。
[0033]
进一步地，估计图像u
(k)
的初始值u
(0)
的所有像素值设置为0。
[0034]
进一步地，扫描几何参数集g包括射线源到探测器中心的距离、射线源到转台中心的距离、探测器单元个数、探测器单元尺寸、扫描角度数和角度采样间隔。
[0035]
本发明还提供一种基于广义全变分的x射线有限角ct图像重建装置，其包括：
[0036]
迭代处理估计模块，其用于利用有限角ct扫描数据集p，通过与扫描几何参数集g相关的图像重建算子rg更新估计图像u
(k)
，获得中间图像u
(k+1/2)
，k为迭代次数；
[0037]
梯度稀疏性约束模块，其用于按照式(1)，利用梯度的稀疏性约束中间图像u
(k+1/2)
，得到估计图像u
(k+1)
；
[0038]u(k+1)
＝p(u
(k+1/2)
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0039]
式(1)中，p表示约束图像梯度稀疏性的算子，由如式(2)提供的最优化问题求解得到：
[0040][0041]
式(2)中：
[0042]
为正则化项；
[0043]
为数据保真项；
[0044]
分别表示输入图像u的水平方向x和竖直方向y的梯度；
[0045]
u表示输入图像；
[0046]
diag(
·
)表示权重矩阵，为每个像素点的水平方向梯度和竖直方向梯度赋予惩罚权重；
[0047]
||
·
||2表示l2范数；
[0048]
表示l2范数的平方；
[0049]
||
·
||1表示l1范数；
[0050]
||
·
||
2,1
表示梯度模长的l1范数；
[0051]
|
·
|表示逐像素做绝对值操作；
[0052]
κ表示用于调节惩罚权重的参数；
[0053]gσ
表示高斯卷积核；
[0054]
*表示卷积运算；
[0055]
t表示转置；
[0056]
λ表示用于控制正则化项和数据保真项权重的参数；
[0057]
判断模块，其用于判断是否达到迭代要求，若是，则终止迭代，输出估计图像u
(k+1)
，反之，则由迭代处理估计模块继续更新估计图像u
(k+1)
。
[0058]
进一步地，像素点水平方向的惩罚权重的方法包括：
[0059]
根据像素点水平方向的梯度，通过式(2)中权重矩阵diag(
·
)中的第一部分计算，得到该像素点水平方向的惩罚权重。
[0060]
进一步地，像素点竖直方向的惩罚权重的方法包括：
[0061]
根据像素点竖直方向的梯度，利用式(2)中权重矩阵diag(
·
)中的第二部分计算，得到该像素点竖直方向的惩罚权重。
[0062]
进一步地，估计图像u
(k)
的初始值u
(0)
的所有像素值设置为0。
[0063]
进一步地，扫描几何参数集g包括射线源到探测器中心的距离、射线源到转台中心的距离、探测器单元个数、探测器单元尺寸、扫描角度数和角度采样间隔。
[0064]
本发明由于采取以上技术方案，其具有以下优点：
[0065]
本发明提供的方法，通过对图像u理想放置状态(水平方向)以及非理想放置状态下(非水平方向)的梯度自适应地赋予不同的惩罚权重，利用图像梯度的稀疏性，逐步恢复图像的边界信息，能够有效地抑制图像的模糊和伪影，提高有限角ct图像的重建质量，进而克服了现有的有限角ct图像重建算法中边界模糊的问题。
附图说明
[0066]
图1为本发明实施例的流程图；
[0067]
图2为本发明实施例中扫描的模体图像；
[0068]
图3为本发明实施例中sart的重建图像；
[0069]
图4为利用本发明实施例获得的重建图像。
具体实施方式
[0070]
下面结合附图和实施例对本发明进行详细的描述。
[0071]
如图1所示，本发明实施例提供的基于广义全变分的x射线有限角ct图像重建方法包括：
[0072]
步骤1，利用有限角ct扫描数据集p，通过与扫描几何参数集g相关的图像重建算子rg更新估计图像u
(k)
，获得中间图像u
(k+1/2)
，k为迭代次数，k＝0,1,
…
,n。其中，更新方法可以是art、sart或fbp等现有的更新方法。ct扫描几何参数集g包括射线源到探测器中心的距离(sdd)、射线源到转台中心的距离(sod)、探测器单元个数、探测器单元尺寸、扫描角度数和角度采样间隔。rg选择迭代类重建算法或解析类重建算法获得，所述迭代类重建算法为art、sart、em；所述解析类重建算法为fbp、bpf。
[0073]
步骤2，按照式(1)，利用梯度的稀疏性约束中间图像u
(k+1/2)
，得到估计图像u
(k+1)
；
[0074]u(k+1)
＝p(u
(k+1/2)
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0075]
式(1)中，p表示约束图像梯度稀疏性的算子，由如式(2)提供的最优化问题求解得到，该优化问题可以通过admm方法求解：
[0076][0077]
式(2)中：
[0078]
为正则化项；
[0079]
为数据保真项；
[0080]
分别表示输入图像u的水平方向x和竖直方向y的梯度，以输入图像u为参考，通常可以理解为该图像的水平方向，即左右方向，视为本发明中的水平方向，而该图像的竖直方向，即上下方向，视为本发明中的竖直方向，因此，的大小为当前像素值与其相邻的右侧的像素值的差值，的大小为当前像素值与其相邻的下方的像素值的差值，而遇到边界时，的大小为当前像素值与第一行像素值的差值，的大小为当前像素值与第一列像素值的差值；
[0081]
u表示输入图像，要与u
(k+1/2)
充分接近；
[0082]
diag(
·
)表示权重矩阵，依赖于迭代过程中输入图像u的局部信息，为每个像素点的水平方向梯度和竖直方向梯度赋予不同的惩罚权重，惩罚权重的具体确定方法可以包括：根据像素点水平方向的梯度，通过式(2)中权重矩阵diag(
·
)中的第一部分计算，得到该像素点水平方向的惩罚权重；与此类似的方法，根据像素点竖直方向的梯度，利用式(2)中权重矩阵diag(
·
)中的第二部分计算得到该像素点竖直方向的惩罚权重。因此，当同一个像素点水平和竖直方向的梯度不同时，diag(
·
)中的两个数不同，那这个像素点水平和竖直方向的惩罚权重就不同。对于不同像素点，不同像素点之间水平和竖直方向的梯度的值不一样，那么不同像素点的惩罚权重也不同；
[0083]
||
·
||2表示l2范数，其表示各个参数的平方的和的开方值，也就是欧氏距离，l2范数越小，可以使得参数的每个元素都很小，接近于0，正则化项处处可导，例如：u表示m
×
m大小的输入图像，u
i,j
表示u的第i行第j列的像素，||u||2＝
[0084]
表示l2范数的平方；
[0085]
||
·
||1表示l1范数，其表示各参数绝对值的和，可以防止过拟合，l1范数不会使得参数等于0而是接近于0，有参数平滑的作用，优化后的参数向量往往比较稀疏；还例如：u＝(u1,u2)，u1、u2分别表示m
×
m大小的矩阵，
[0086]
||
·
||
2,1
表示梯度模长的l1范数，梯度模长定义为的平方和再开根号获得；
[0087]
|
·
|表示逐像素做绝对值操作；
[0088]
κ表示用于调节惩罚权重的参数，该参数能够更好地区分图像的光滑区域与边界区域，具体数值通过实验得到；
[0089]gσ
表示高斯卷积核；
[0090]
*表示卷积运算；
[0091]
t表示转置；
[0092]
λ表示用于控制正则化项和数据保真项权重的参数，具体数值通过实验得到。
[0093]
通过式(2)，可以解决现有有限角ct图像重建算法重建图像沿特定方向模糊，并存在有限角伪影的问题。
[0094]
当然，也可以通过现有的施加保边界扩散和保边界光滑正则项、施加曲率约束正则项等方法获得约束图像梯度稀疏性的算子p。
[0095]
步骤3，判断是否达到迭代要求，若是，则终止迭代，输出估计图像u
(k+1)
，反之，则转至步骤1。
[0096]
本实施例提供的方法，通过对图2示意的输入图像u进行水平和竖直方向的梯度自适应地赋予不同的惩罚权重，利用图像梯度的稀疏性，逐步恢复图像的边界信息，从而能够有效地抑制图像的模糊和伪影，提高有限角ct图像的重建质量，进而克服了现有的有限角ct图像重建算法中边界模糊的问题。
[0097]
在一个实施例中，估计图像u
(k)
的初始值u
(0)
的所有像素值设置为0，当然，也可以为其他值。
[0098]
在一个实施例中，还可以设置迭代终止阈值ε以及迭代次数上限n，那么，步骤3判断是否达到迭代要求时，既可以判断相邻两次迭代图像间的差别是否小于给定阈值，即||u
(k+1)-u
(k)
||≤ε，也可以判断迭代次数否达到迭代次数上限n。
[0099]
本发明实施例还提供一种基于广义全变分的x射线有限角ct图像重建装置，其包括迭代处理估计模块、梯度稀疏性约束模块和判断模块，其中：
[0100]
迭代处理估计模块用于利用有限角ct扫描数据集p，通过与扫描几何参数集g相关
的图像重建算子rg更新估计图像u
(k)
，获得中间图像u
(k+1/2)
，k为迭代次数。
[0101]
梯度稀疏性约束模块用于按照式(1)，利用梯度的稀疏性约束中间图像u
(k+1/2)
，得到估计图像u
(k+1)
。其中使用的约束图像梯度稀疏性的算子p可以由如式(2)提供的最优化问题求解。
[0102]
判断模块用于判断是否达到迭代要求，若是，则终止迭代，输出估计图像u
(k+1)
，反之，则由迭代处理估计模块继续更新估计图像u
(k+1)
。
[0103]
为了更好地体现本发明提供的有限角ct图像重建方法在重建效果方面的优势，下面结合一具体实施例将本发明所述的方法与已存在的典型算法sart进行比较。
[0104]
本实施例数据采集自工业ct系统，扫描样品为pcb模体，pcb板非理想放置于扫描台，其图像如图2所示。实验电压与电流分别为140kv以及100ua，在120度范围内每隔0.5度采集一次投影数据。
[0105]
分别采用sart方法和本发明方法对pcb板的扫描数据进行重建，重建结果分别如图3、图4所示。其中：图3是sart方法的重建结果；图4是本发明方法的重建结果。可以看出，sart方法重建的图像因扫描数据不完全存在严重的边界模糊和条状伪影，pcb模体中的边界信息几乎完全丢失，从图4可以看到，本发明实施例获得的重建图像边界清晰，准确地重建了具有较高重建难度的长条状边界。
[0106]
最后需要指出的是：以上实施例仅用以说明本发明的技术方案，而非对其限制。本领域的普通技术人员应当理解：可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。