首页> 中国专利> 伪影修正辅助的CBCT迭代重建方法

伪影修正辅助的CBCT迭代重建方法

页面导航

摘要
著录项
法律信息
说明书
相似文献

摘要

本发明公开了一种伪影修正辅助的CBCT迭代重建方法，该方法首先使用病人或者模体的CBCT投影数据重建出初始图像；其次引入偏置场项，将伪影修正的概念引入到迭代重建的理论框架中去；偏置场为现实中受到伪影污染的重建图像和理想情况下没有伪影污染，满足分段常数性质的图像之间的差距；然后使用轮寻的思想求解引入偏置场项的迭代重建框架的目标函数，得到待重建图像以及偏置场：最后将求得的偏置场叠加到待重建图像上，最终实现图像域阴影伪影修正。本发明在原有的迭代重建理论框架下引入伪影修正的基本概念，构造偏置场项来平衡大面积阴影伪影对重建算法的不良影响，实现了没有阴影伪影的精确、稳定、快速的迭代重建。

著录项

公开/公告号CN105631909A

专利类型发明专利
公开/公告日2016-06-01

原文格式PDF
申请/专利权人浙江大学;
展开▼

申请/专利号CN201510980575.9
发明设计人牛田野;吴蓬威;龚书涛;
展开▼

申请日2015-12-23
分类号G06T11/00(20060101);G06T5/00(20060101);
代理机构33200 杭州求是专利事务所有限公司;
代理人邱启旺
地址 310027 浙江省杭州市西湖区浙大路38号
入库时间 2023-12-18 15:29:29

法律信息

法律状态公告日

法律状态信息

法律状态
2018-05-29

授权

授权
2016-06-29

实质审查的生效 IPC(主分类):G06T11/00 申请日:20151223

实质审查的生效
2016-06-01

公开

公开

说明书

技术领域

本发明属于医学成像技术领域，尤其涉及一种伪影修正辅助的Cone‐beamCT (CBCT)迭代重建方法。

背景技术

采用压缩感知技术的CBCT迭代重建算法具有相对完善的理论框架，相比于传统的滤波反投影重建，它可以在很大程度上减少图像重建所需要的投影数量，从而减少病人在成像过程中摄入的放射剂量，保证病人的生命安全，因而它在医学影像领域有着重要的应用和潜在的商业价值。整个迭代重建算法的基础是：重建图像的CT值具有分段常数或者近似分段常数的兴致。但是在CBCT成像过程中，由于锥角和接收X光部分的面积较大，散射污染和射术硬化等不良因素会导致重建图像上具有严重的阴影伪影。这些伪影并不满足分段常数的性质，在一定程度上破坏了迭代重建的理论基础，影响了整个重建算法的稳定性和速度，同时重建结果的CT数精度也在很大程度上被阴影伪影所破坏。

发明内容

本发明的目的在于针对现有技术的不足，提供一种伪影修正辅助的CBCT迭代重建方法。

本发明的目的是通过以下技术方案来实现的：一种伪影修正辅助的CBCT迭代重建方法，包括以下步骤：

(1)使用病人或者模体CBCT投影数据重建出初始图像f₀；

(2)建立含有偏置场项的迭代重建框架：引入偏置场项，将伪影修正的概念引入到迭代重建的理论框架中去；将现实中受到伪影污染的重建图像和理想情况下没有伪影污染，满足分段常数性质的图像之间的差距称为偏置场；在引入偏置场后，重建算法的目标函数如下所示：

$(f, f_{b i a s}) = argmin [\frac{1}{2} | | M f - b | |_{2}^{2} + λ | | f + f_{b i a s} | |_{T V}], - - - (1)$

其中f是需要求取的重建图像，f_bias是需要求取的偏置场。M是写成矩阵形式的投影，也被称为正投矩阵，Mf代表对重建图像进行正投操作。b代表原始投影转换为线积分之后的数据，λ为正则化项因子。‖■‖₂代表求取二范数，‖■‖_TV代表求取全变分，为置信项。

(3)求解迭代重建框架的目标函数，得到待重建图像f以及偏置场f_bias：使用轮寻的思想来求解目标函数，在每次迭代的过程中，首先计算偏置场，然后固定偏置场，将目标函数转变为单变量优化问题(2)来进行求解：

$f = argmin [\frac{1}{2} | | M f - b | |_{2}^{2} + λ | | f + f_{b i a s} | |_{T V}], - - - (2)$

(4)将求得的偏置场f_bias叠加到待重建图像f上，最终实现图像域阴影伪影修正。

进一步地，所述步骤1中使用滤波反投影技术重建出初始图像f₀。

进一步地，所述步骤3具体包括以下子步骤：

(3.1)通过公式(3)计算偏置场：

f_bias＝H(f_seg-f),(3)

其中，f_seg是在待重建图像上进行图像分割而获得的模板图像，模板图像的灰度值被填充为不同组织CT数的标准值，通过图像分割将高对比度物质和软组织区分开来；H是一个低通且连续的滤波器，在保持偏置场强度的情况下将偏置场从残差图像中提取出来；残差图像指模板图像与待重建图像的差距。

(3.2)使用GP‐BB方法进行目标函数最小化，具体包括以下步骤：

(3.2.1)按照如下公式进行目标函数梯度g的求取：

$g = ▿ | | f + f_{b i a s} | |_{T V} + M^{T} (M f - b), - - - (4)$

公式中T是用来进行矩阵转置的算子，使用如下公式来进行全变分的求导：

$\begin{matrix} ▿ | | F | |_{T V} = \frac{3 F_{k, m, n} - F_{k - 1, m, n} - F_{k, m - 1, n} - F_{k, m, n - 1}}{u (m, n)} - \frac{F_{k + 1, m, n} - F_{k, m, n}}{u (k + 1, m, n)} - \frac{F_{k, m + 1, n} - F_{k, m, n}}{u (k, m + 1, n)} \\ - \frac{F_{k, m, n + 1} - F_{k, m, n}}{u (k, m, n + 1)} \\ = \frac{3 F_{k, m, n} - F_{k - 1, m, n} - F_{k, m - 1, n} - F_{k, m, n - 1}}{\sqrt{δ + {(F_{k - 1, m, n} - F_{k, m, n})}^{2} + {(F_{k, m - 1, n} - F_{k, m, n})}^{2} + {(F_{k, m, n - 1} - F_{k, m, n})}^{2}}} \\ - \frac{F_{k + 1, m, n} - F_{k, m, n}}{\sqrt{δ + {(F_{k, m, n} - F_{k + 1, m, n})}^{2} + {(F_{k + 1, m - 1, n} - F_{k + 1, m, n})}^{2} + {(F_{k + 1, m, n - 1} - F_{k + 1, m, n})}^{2}}} \\ - \frac{F_{k, m + 1, n} - F_{k, m, n}}{\sqrt{δ + {(F_{k - 1, m + 1, n} - F_{k, m + 1, n})}^{2} + {(F_{k, m, n} - F_{k, m + 1, n})}^{2} + {(F_{k, m + 1, n - 1} - F_{k, m + 1, n})}^{2}}} \\ - \frac{F_{k, m, n + 1} - F_{k, m, n}}{\sqrt{δ + {(F_{k - 1, m, n + 1} - F_{k, m, n + 1})}^{2} + {(F_{k, m - 1, n + 1} - F_{k, m, n + 1})}^{2} + {(F_{k, m, n} - F_{k, m, n + 1})}^{2}}}, \end{matrix} - - - (5)$

公式中的δ是一个小的正数。GP算法使用如下公式来对重建图像进行更新：

f_n+1＝max(f_n-a_np_n,0),(6)

其中，α是每次迭代中的步长，投影后的梯度被记作p_n，计算公式如下：

$p_{n} (l) = (\begin{matrix} g_{n} (l), & i f >g_{n}(l)\leq0,or>f_{n}(l)>00, & o t h e r w i s e,---(7) \end{matrix})$

其中l是体素点的位置坐标。

(3.2.2)使用BB算法来解析的计算每一次迭代步长的α，在每次迭代中计算出两个步长，公式如下：

$α_{n}^{B B 1} = \frac{{(f_{n} - f_{n - 1})}^{T} (f_{n} - f_{n - 1})}{{(f_{n} - f_{n - 1})}^{T} (p_{n} - p_{n - 1})}, - - - (8)$

$α_{n}^{B B 2} = \frac{{(f_{n} - f_{n - 1})}^{T} (p_{n} - p_{n - 1})}{{(p_{n} - p_{n - 1})}^{T} (p_{n} - p_{n - 1})}, - - - (9)$

下标n代表本次迭代，而下标n-1代表前一次迭代。

在两个步长使用公式(10)选择一个步长：

$α_{n} = (\begin{matrix} α_{n}^{B B 1}, & i f >α_{n}^{B B 1}α_{n}^{B B 2}<κ \\ α_{n}^{B B 2}, & o t h e r w i s e \end{matrix}), - - - (10)$

其中，κ是一个小于1的正数。

(3.3)停止判据：通过判断全变分的作用和置信项的作用是否达到平衡来判断迭代算法是否停止。确定这两种作用的公式如下：

$d_{T V} = d i a g (f_{i n d i c a t o r}) (▿ | | f + f_{b i a s} | |_{T V}), - - - (11)$

$d_{d a t a} = d i a g (f_{i n d i c a t o r}) (▿ | | M f - b | |_{2}^{2}), - - - (12)$

其中，diag(x)是用来产生对角阵的函数，x上的元素会被填充到对角线上。 f_indicator是一个每当f不等于0就取1的指示函数，等于0则取0，确定了这两种作用后，按照下述公式计算停止判据c_α：

$c_{α} = \frac{d_{T V} \cdot d_{d a t a}}{| d_{T V} | | d_{d a t a} |}, - - - (13)$

在c_α小于设定阈值并保持一段时间后终止算法。

进一步地，所述步骤3.1中，使用阈值化算法和二相水平集算法结合的方法来进行图像分割。

进一步地，所述二相水平集算法公式如下所示：

$\begin{matrix} F (φ, c, b) = \int (Σ_{i = 1}^{2} \int K (y - x) {| I (x) - b (y) c_{i} |}^{2} M_{i} (φ (x)) d x) d y + v \int | ▿ H (φ) | d x \\ + u \int p | ▿ φ | d x, \end{matrix} - - - (14)$

式中，φ是水平集函数，该函数通过其函数值的符号进行两个不同区域的分割，I 是需要进行分割的图像，c_i是需要填充到第i个区域的函数值，y代表图像上每一个点，而x 代表y的邻域内的每个点。b是用来对图像本身的不均匀度进行补偿的补偿项。K(y-x)是一个非负的窗函数，该函数在不属于y的邻域的部分取0。M_i是对第i种组织用符号函数定义的成员函数，在属于这个组织的部分该函数会取为1。u，v是用来进行效果调节的参量，通过对这两个参量进行调节来实现置信项和平滑项之间的平衡。为了将公式(14)最小化，依次基于φ,c,b来执行梯度降算法，公式(4)在被最小化之后就可以提供分割区域，将特定组织的标准CT值填入到分割区域中，从而得到模板图像。

进一步地，所述步骤3.1中，所述滤波器为二维的Savitzky‐Golay滤波器。

进一步地，所述步骤3中，并没有在每次迭代中对偏置场进行更新，而是每间隔固定次数的迭代后再对偏置场进行更新。如果在某一次迭代中偏置场没有被更新，那么算法将沿用上一次迭代中的偏置场。

本发明的有益效果是：本发明在原有的迭代重建理论框架下引入伪影修正的基本概念，构造偏置场项来平衡大面积阴影伪影对重建算法的不良影响，实现了没有阴影伪影的精确、稳定、快速的迭代重建。

附图说明

图1为本发明在600模体上的实施结果，(a):使用655个投影进行传统FBP 重建的结果，(b):首先使用655个投影进行传统FBP重建，然后使用计划CT图像进行散射修正的结果，(c):在迭代算法结束后的最终偏置场，(d):使用传统迭代重建算法基于92个投影得到的重建结果，(e):使用本发明中的重建算法基于92个投影得到的重建结果，显示窗 (除(c))为[‐250250]HU。

图2为本发明在600模体低对比度层上的实施结果，(a):使用655个投影进行传统FBP重建的结果，(b):首先使用655个投影进行传统FBP重建，然后使用计划CT图像进行散射修正的结果，(c):使用传统迭代重建算法基于92个投影得到的重建结果，(d):使用本发明中的重建算法基于92个投影得到的重建结果。在(b)中使用黑色虚线框标准的部分在(z.a)，(z.b)，(z.c)，(z.d)中被放大显示。(a)，(b)，(c)中的显示窗为[‐2500]HU，(d) 中的显示窗为[‐100150]HU。

图3为本发明在600模体高对比度线对层上的实施结果，(a):使用655个投影进行传统FBP重建的结果，(b):首先使用655个投影进行传统FBP重建，然后使用计划CT 图像进行散射修正的结果，(c):使用传统迭代重建算法基于92个投影得到的重建结果， (d):使用本发明中的重建算法基于92个投影得到的重建结果。在(b)中使用黑色虚线框标注的部分在(z.a)，(z.b)，(z.c)，(z.d)中被放大显示，(a)，(c)，(d)中的显示窗为[‐250 250]HU，(b)中的显示窗为[‐50450]HU。

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细说明。

本发明提出的一种伪影修正辅助的CBCT迭代重建方法，包括以下步骤：

(1)使用滤波反投影技术产生初始图像

将实际测得的病人或者模体CBCT投影数据使用滤波反投影技术重建出初始图像 f₀，注意，在减少剂量的情况下，采集到的投影数据相对较少，滤波反投影算法重建得到的图像具有较多条纹和噪声污染。

(2)建立含有偏置场项的迭代重建框架

通过引入偏置场项，将伪影修正的概念引入到迭代重建的理论框架中去。将现实中受到伪影污染的重建图像和理想情况下没有伪影污染，满足分段常数性质的图像之间的差距称为偏置场。在引入偏置场后，本发明中重建算法的目标函数如下所示：

$(f, f_{b i a s}) = argmin [\frac{1}{2} | | M f - b | |_{2}^{2} + λ | | f + f_{b i a s} | |_{T V}], - - - (1)$

其中f是需要求取的重建图像，f_bias是需要求取的偏置场。M是写成矩阵形式的投影，也被称为正投矩阵，Mf代表对重建图像进行正投操作，该操作为线性因而可以被写成矩阵相乘的形式。b代表原始投影转换为线积分之后的数据，λ为手工选取的正则化项因子，用来控制重建图像的平滑程度。‖■‖₂代表求取二范数，‖■‖_TV代表求取全变分，全变分在数学被定义为空间梯度图像的一范数，为置信项。

通过求解该优化问题，就可以同时得到待重建图像f，以及偏置场f_bias。随后简单的将偏置场叠加到待重建图像上就可以完成精确的图像域阴影伪影修正了。

(3)求解迭代重建框架的目标函数，得到待重建图像f以及偏置场f_bias

本发明使用的目标函数(1)由于有两个自变量的存在，很难直接进行最小化操作。因此本发明提出使用轮寻的思想来进行求解。也就是在每次迭代的过程中，首先使用信号与图像处理的手段来计算偏置场，偏置场在计算完成后即被固定为不变，从而将目标函数转变为单变量优化问题来进行求解。本发明计算偏置场的方法如下所示：

f_bias＝H(f_seg-f),(2)

其中，f_seg是在待重建图像上进行图像分割而获得的模板图像，模板图像的灰度值被填充为不同组织CT数的标准值，因而不包含伪影信息。f是待重建图像，H是一个低通且连续的滤波器，该滤波器可以将f_seg和f_p之间的插值图像进行滤波，从而得到估计的偏置场。

以上就是轮寻求解的第一步，也就是在每次迭代开始的时候首先对偏置场进行更新。在进行过偏置场的更新后，本发明将偏置场固定为不变，这样，只剩下一个变量的优化函数如公式(3)可以按照传统迭代重建的算法进行求解，这就是轮寻求解的第二步，下面本发明将分别对这两步进行介绍。

$f = argmin [\frac{1}{2} | | M f - b | |_{2}^{2} + λ | | f + f_{b i a s} | |_{T V}], - - - (3)$

在第一次迭代时，将待重建图像f置为初始图像f₀；

(3.1)偏置场求解

(3.1.1)图像分割算法

为了对公式(2)中的f_seg进行计算，本发明使用阈值化算法和二相水平集算法结合的方法来进行图像分割。本发明认为，骨头和空气等高对比度物质，由于其CT数和软组织 (包括脂肪，肌肉等)具有较大程度的区别，可以直接通过手工设定阈值的做法将其分割出来，也就是CT数大于某一个设定值的部分认定为骨骼，CT数小于某一个设定值的部分认定为空气。

在完成高对比物体的分割之后，本发明使用二相水平集算法进行软组织的分割，也就是进行脂肪和肌肉的分割。本发明中的二相水平集算法公式如下所示：

$\begin{matrix} F (φ, c, b) = \int (Σ_{i = 1}^{2} \int K (y - x) {| I (x) - b (y) c_{i} |}^{2} M_{i} (φ (x)) d x) d y + v \int | ▿ H (φ) | d x \\ + u \int p | ▿ φ | d x, \end{matrix} - - - (4)$

式中，φ是水平集函数，该函数通过其函数值的符号进行两个不同区域的分割，I 是需要进行分割的图像，c_i是需要填充到第i个区域的函数值，y代表图像上每一个点，而x 代表y的邻域内的每个点。b是用来对图像本身的不均匀度进行补偿的补偿项。K(y-x)是一个非负的窗函数，该函数在不属于y的邻域的部分取0。M_i是对第i种组织用符号函数定义的成员函数，在属于这个组织的部分该函数会取为1。u，v是用来进行效果调节的参量，本发明通过对这两个参量进行调节来实现置信项和平滑项之间的平衡。为了将公式(4)最小化，本发明依次的基于φ,c,b来执行梯度降算法，公式(4)在被最小化之后就可以提供分割区域，本发明将特定组织的标准CT值填入到分割区域中，从而得到模板图像。

(3.1.2)滤波器设计

本发明使用的滤波器是二维的Savitzky‐Golay滤波器(也就是公式(2)中的H)，该滤波器具有较好的轮廓保持能力，因而可以在保持偏置场强度的情况下将偏置场从残差图像中提取出来。残差图像指模板图像与待重建图像的差距。该滤波器连续的使用相邻点信息，基于低次多项式来对未知点进行最小二乘拟合。这种效果可以通过卷积来实现，因而具有较低的计算复杂度。

(3.2)使用GP‐BB方法进行目标函数最小化

求得偏置场后的下一个步骤是固定偏置场，来最小化公式(3)。本发明使用GP‐BB 方法来进行目标函数的最小化，具体如下：

首先，按照如下公式进行目标函数(公式(3))梯度g的求取：

$g = ▿ | | f + f_{b i a s} | |_{T V} + M^{T} (M f - b), - - - (5)$

公式中T是用来进行矩阵转置的算子，本发明中使用如下公式来进行全变分的求导：

公式中的δ是一个小的正数，用来防止除0的错误。GP算法使用如下公式来对重建图像进行更新：

f_n+1＝max(f_n-α_np_n,0),(7)

其中，α是每次迭代中的步长，投影后的梯度被记作p_n，计算公式如下：

$p_{n} (l) = (\begin{matrix} g_{n} (l), & i f >g_{n}(l)\leq0,or>f_{n}(l)>00, & o t h e r w i s e,---(8) \end{matrix})$

其中l是体素点的位置坐标。

本发明使用BB算法来解析的计算每一次迭代步长的α，在每次迭代中，本发明计算出两个步长，使用的公式如下：

$α_{n}^{B B 1} = \frac{{(f_{n} - f_{n - 1})}^{T} (f_{n} - f_{n - 1})}{{(f_{n} - f_{n - 1})}^{T} (p_{n} - p_{n - 1})}, - - - (9)$

$α_{n}^{B B 2} = \frac{{(f_{n} - f_{n - 1})}^{T} (p_{n} - p_{n - 1})}{{(p_{n} - p_{n - 1})}^{T} (p_{n} - p_{n - 1})}, - - - (10)$

下标n代表本次迭代，而下标n-1代表前一次迭代。本发明使用如下公式来确定在上述两个步长中选择哪一个：

$α_{n} = (\begin{matrix} α_{n}^{B B 1}, & i f >α_{n}^{B B 1}α_{n}^{B B 2}<κ \\ α_{n}^{B B 2}, & o t h e r w i s e \end{matrix}), - - - (11)$

其中，κ是一个小于1的正数。

(3.3)停止判据

本发明通过判断全变分的作用和置信项的作用是否达到平衡来判断迭代算法是否停止。确定这两种作用的公式如下：

$d_{T V} = d i a g (f_{i n d i c a t o r}) (▿ | | f + f_{b i a s} | |_{T V}), - - - (12)$

$d_{d a t a} = d i a g (f_{i n d i c a t o r}) (▿ | | M f - b | |_{2}^{2}), - - - (13)$

其中，diag(x)是用来产生对角阵的函数，x上的元素会被填充到对角线上。 f_indicator是一个每当f不等于0就取1的指示函数(等于0则取0)，确定了这两种作用后，按照下述公式计算停止判据c_α：

$c_{α} = \frac{d_{T V} \cdot d_{d a t a}}{| d_{T V} | | d_{d a t a} |}, - - - (14)$

在理想情况下，c_α＝-1的时候算法应该停止，但是由于噪声等非理想因素的存在，本发明在c_alpha小于某一设定阈值并保持一段时间后终止算法。

(4)将求得的偏置场叠加到待重建图像上，最终实现图像域阴影伪影修正。

由于在本发明中，轮寻更新偏置场所需要的迭代次数远小于轮寻更新目标函数所需要的迭代次数。为了减少计算负担，本发明并没有在每次迭代中对偏置场进行更新，而是每间隔固定次数的迭代后再对偏置场进行更新。如果在某一次迭代中偏置场没有被更新，那么算法将沿用上一次迭代中的偏置场。

偏置场的更新是基于公式(2)的，也就是说，模板图像或者重建图像的更新都会造成偏置场的更新。本发明中更新模板图像的频率小于更新偏置场的频率。但是由于重建图像在每次迭代中都会变化，所以在不更新模板图像的时候也可以更新偏置场。

实施例

1.在三个实施例中，选择的停止条件阈值为：c_alpha＜-0.8，并且保持100次迭代。模板图像的更新为每隔33次迭代进行一次，只在前100次迭代中进行，也就是在第[1,34,67, 100]次迭代的时候进行模板图像的更新。偏置场的更新为每隔10次迭代进行一次。

2.在600模体中，迭代重建的正则化项参数：λ＝0.2。

3.在三个实施例中，全变分求导公式(6)中使用的δ为10^-8。步长选择公式(11)中使用的κ为0.3。

4.上述技术方案中，使用传统迭代重建算法得到的重建图像如图(2)(3)中的(c)，以及图(1)中的(d)所示，可以发现使用传统的迭代重建算法无法保证足够的CT数精度，重建后的图像具有较大面积的阴影伪影。

5.上述技术方案中，所述的600模体，使用本发明中的重建算法得到的结果，相比使用传统迭代重建算法得到的结果，CT数精度从224HU降低到了6HU，空间均匀度从42.74％提升到了63.92％，算法完成所需要的迭代次数从320次减小到了160次。

以上所述的具体实施方式对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的最优选实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换等，均应包含在本发明的保护范围之内。

去获取专利，查看全文>

相似文献

专利
中文文献
外文文献

1. 伪影修正辅助的CBCT迭代重建方法 [P] . 中国专利： CN105631909B . 2018.05.29
2. 伪影修正辅助的CBCT迭代重建方法 [P] . 中国专利： CN105631909A . 2016-06-01
3. METHODS OF IMAGE RECONSTRUCTION TO REDUCE ARTIFACTS IN RAPID CBCT SCANS [P] . US2021264591A1 . 2021-08-26

机译：减少快速CBCT扫描中伪影的图像重建方法
4. Method for creating reconstruction image from previously recorded image data by iterative reconstruction method under consideration of auxiliary conditions, involves detecting position dependent standard deviation of image data [P] . 德国专利： DE102008027749A1 . 2009-12-24

机译：在考虑辅助条件的情况下通过迭代重建方法从先前记录的图像数据中创建重建图像的方法，包括检测图像数据的位置相关标准偏差
5. RECONSTRUCTION METHOD FOR MOTION COMPENSATED HIGH QUALITY 4D-CBCT IMAGE BASED ON BILATERAL FILTERING [P] . US2021248790A1 . 2021-08-12

机译：基于双侧滤波的运动补偿高质量4D-CBCT图像的重建方法