CN105631909B - 伪影修正辅助的cbct迭代重建方法 - Google Patents

伪影修正辅助的cbct迭代重建方法 Download PDF

Info

Publication number
CN105631909B
CN105631909B CN201510980575.9A CN201510980575A CN105631909B CN 105631909 B CN105631909 B CN 105631909B CN 201510980575 A CN201510980575 A CN 201510980575A CN 105631909 B CN105631909 B CN 105631909B
Authority
CN
China
Prior art keywords
mrow
msub
msup
mtd
bias
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201510980575.9A
Other languages
English (en)
Other versions
CN105631909A (zh
Inventor
牛田野
吴蓬威
龚书涛
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University ZJU
Original Assignee
Zhejiang University ZJU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang University ZJU filed Critical Zhejiang University ZJU
Priority to CN201510980575.9A priority Critical patent/CN105631909B/zh
Publication of CN105631909A publication Critical patent/CN105631909A/zh
Application granted granted Critical
Publication of CN105631909B publication Critical patent/CN105631909B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/008—Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20172—Image enhancement details
    • G06T2207/20192—Edge enhancement; Edge preservation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)
  • Medicines Containing Antibodies Or Antigens For Use As Internal Diagnostic Agents (AREA)

Abstract

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

Description

伪影修正辅助的CBCT迭代重建方法
技术领域
本发明属于医学成像技术领域，尤其涉及一种伪影修正辅助的Cone‐beam CT(CBCT)迭代重建方法。
背景技术
采用压缩感知技术的CBCT迭代重建算法具有相对完善的理论框架，相比于传统的滤波反投影重建，它可以在很大程度上减少图像重建所需要的投影数量，从而减少病人在成像过程中摄入的放射剂量，保证病人的生命安全，因而它在医学影像领域有着重要的应用和潜在的商业价值。整个迭代重建算法的基础是：重建图像的CT值具有分段常数或者近似分段常数的兴致。但是在CBCT成像过程中，由于锥角和接收X光部分的面积较大，散射污染和射术硬化等不良因素会导致重建图像上具有严重的阴影伪影。这些伪影并不满足分段常数的性质，在一定程度上破坏了迭代重建的理论基础，影响了整个重建算法的稳定性和速度，同时重建结果的CT数精度也在很大程度上被阴影伪影所破坏。
发明内容
本发明的目的在于针对现有技术的不足，提供一种伪影修正辅助的CBCT迭代重建方法。
本发明的目的是通过以下技术方案来实现的：一种伪影修正辅助的CBCT迭代重建方法，包括以下步骤：
(1)使用病人或者模体CBCT投影数据重建出初始图像f0；
(2)建立含有偏置场项的迭代重建框架：引入偏置场项，将伪影修正的概念引入到迭代重建的理论框架中去；将现实中受到伪影污染的重建图像和理想情况下没有伪影污染，满足分段常数性质的图像之间的差距称为偏置场；在引入偏置场后，重建算法的目标函数如下所示：
其中f是需要求取的重建图像，fbias是需要求取的偏置场。M是写成矩阵形式的投影，也被称为正投矩阵，Mf代表对重建图像进行正投操作。b代表原始投影转换为线积分之后的数据，λ为正则化项因子。‖■‖2代表求取二范数，‖■‖TV代表求取全变分，为置信项。
(3)求解迭代重建框架的目标函数，得到待重建图像f以及偏置场fbias：使用轮寻的思想来求解目标函数，在每次迭代的过程中，首先计算偏置场，然后固定偏置场，将目标函数转变为单变量优化问题(2)来进行求解：
(4)将求得的偏置场fbias叠加到待重建图像f上，最终实现图像域阴影伪影修正。
进一步地，所述步骤1中使用滤波反投影技术重建出初始图像f0。
进一步地，所述步骤3具体包括以下子步骤：
(3.1)通过公式(3)计算偏置场：
fbias＝H(fseg-f), (3)
其中，fseg是在待重建图像上进行图像分割而获得的模板图像，模板图像的灰度值被填充为不同组织CT数的标准值，通过图像分割将高对比度物质和软组织区分开来；H是一个低通且连续的滤波器，在保持偏置场强度的情况下将偏置场从残差图像中提取出来；残差图像指模板图像与待重建图像的差距。
(3.2)使用GP‐BB方法进行目标函数最小化，具体包括以下步骤：
(3.2.1)按照如下公式进行目标函数梯度g的求取：
公式中T是用来进行矩阵转置的算子，使用如下公式来进行全变分的求导：
公式中的δ是一个小的正数。GP算法使用如下公式来对重建图像进行更新：
fn+1＝max(fn-αnpn,0), (6)
其中，α是每次迭代中的步长，投影后的梯度被记作pn，计算公式如下：
其中l是体素点的位置坐标。
(3.2.2)使用BB算法来解析的计算每一次迭代步长的α，在每次迭代中计算出两个步长，公式如下：
下标n代表本次迭代，而下标n-1代表前一次迭代。
在两个步长使用公式(10)选择一个步长：
其中，κ是一个小于1的正数。
(3.3)停止判据：通过判断全变分的作用和置信项的作用是否达到平衡来判断迭代算法是否停止。确定这两种作用的公式如下：
其中，diag(x)是用来产生对角阵的函数，x上的元素会被填充到对角线上。findicator是一个每当f不等于0就取1的指示函数，等于0则取0，确定了这两种作用后，按照下述公式计算停止判据cα：
在cα小于设定阈值并保持一段时间后终止算法。
进一步地，所述步骤3.1中，使用阈值化算法和二相水平集算法结合的方法来进行图像分割。
进一步地，所述二相水平集算法公式如下所示：
式中，φ是水平集函数，该函数通过其函数值的符号进行两个不同区域的分割，I是需要进行分割的图像，ci是需要填充到第i个区域的函数值，y代表图像上每一个点，而x代表y的邻域内的每个点。b是用来对图像本身的不均匀度进行补偿的补偿项。K(y-x)是一个非负的窗函数，该函数在不属于y的邻域的部分取0。Mi是对第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个投影得到的重建结果，(f):在迭代算法结束后的最终模板图，显示窗(除(c))为[‐250 250]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)中的显示窗为[‐250 0]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)中的显示窗为[‐250250]HU，(b)中的显示窗为[‐50 450]HU。
具体实施方式
下面结合附图和具体实施例对本发明作进一步详细说明。
本发明提出的一种伪影修正辅助的CBCT迭代重建方法，包括以下步骤：
(1)使用滤波反投影技术产生初始图像
将实际测得的病人或者模体CBCT投影数据使用滤波反投影技术重建出初始图像f0，注意，在减少剂量的情况下，采集到的投影数据相对较少，滤波反投影算法重建得到的图像具有较多条纹和噪声污染。
(2)建立含有偏置场项的迭代重建框架
通过引入偏置场项，将伪影修正的概念引入到迭代重建的理论框架中去。将现实中受到伪影污染的重建图像和理想情况下没有伪影污染，满足分段常数性质的图像之间的差距称为偏置场。在引入偏置场后，本发明中重建算法的目标函数如下所示：
其中f是需要求取的重建图像，fbias是需要求取的偏置场。M是写成矩阵形式的投影，也被称为正投矩阵，Mf代表对重建图像进行正投操作，该操作为线性因而可以被写成矩阵相乘的形式。b代表原始投影转换为线积分之后的数据，λ为手工选取的正则化项因子，用来控制重建图像的平滑程度。‖■‖2代表求取二范数，‖■‖TV代表求取全变分，全变分在数学被定义为空间梯度图像的一范数，为置信项。
通过求解该优化问题，就可以同时得到待重建图像f，以及偏置场fbias。随后简单的将偏置场叠加到待重建图像上就可以完成精确的图像域阴影伪影修正了。
(3)求解迭代重建框架的目标函数，得到待重建图像f以及偏置场fbias
本发明使用的目标函数(1)由于有两个自变量的存在，很难直接进行最小化操作。因此本发明提出使用轮寻的思想来进行求解。也就是在每次迭代的过程中，首先使用信号与图像处理的手段来计算偏置场，偏置场在计算完成后即被固定为不变，从而将目标函数转变为单变量优化问题来进行求解。本发明计算偏置场的方法如下所示：
fbias＝H(fseg-f), (2)
其中，fseg是在待重建图像上进行图像分割而获得的模板图像，模板图像的灰度值被填充为不同组织CT数的标准值，因而不包含伪影信息。f是待重建图像，H是一个低通且连续的滤波器，该滤波器可以将fseg和fp之间的插值图像进行滤波，从而得到估计的偏置场。
以上就是轮寻求解的第一步，也就是在每次迭代开始的时候首先对偏置场进行更新。在进行过偏置场的更新后，本发明将偏置场固定为不变，这样，只剩下一个变量的优化函数如公式(3)可以按照传统迭代重建的算法进行求解，这就是轮寻求解的第二步，下面本发明将分别对这两步进行介绍。
在第一次迭代时，将待重建图像f置为初始图像f0；
(3.1)偏置场求解
(3.1.1)图像分割算法
为了对公式(2)中的fseg进行计算，本发明使用阈值化算法和二相水平集算法结合的方法来进行图像分割。本发明认为，骨头和空气等高对比度物质，由于其CT数和软组织(包括脂肪，肌肉等)具有较大程度的区别，可以直接通过手工设定阈值的做法将其分割出来，也就是CT数大于某一个设定值的部分认定为骨骼，CT数小于某一个设定值的部分认定为空气。
在完成高对比物体的分割之后，本发明使用二相水平集算法进行软组织的分割，也就是进行脂肪和肌肉的分割。本发明中的二相水平集算法公式如下所示：
式中，φ是水平集函数，该函数通过其函数值的符号进行两个不同区域的分割，I是需要进行分割的图像，ci是需要填充到第i个区域的函数值，y代表图像上每一个点，而x代表y的邻域内的每个点。b是用来对图像本身的不均匀度进行补偿的补偿项。K(y-x)是一个非负的窗函数，该函数在不属于y的邻域的部分取0。Mi是对第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的求取：
公式中T是用来进行矩阵转置的算子，本发明中使用如下公式来进行全变分的求导：
公式中的δ是一个小的正数，用来防止除0的错误。GP算法使用如下公式来对重建图像进行更新：
fn+1＝max(fn-αnpn,0), (7)
其中，α是每次迭代中的步长，投影后的梯度被记作pn，计算公式如下：
其中l是体素点的位置坐标。
本发明使用BB算法来解析的计算每一次迭代步长的α，在每次迭代中，本发明计算出两个步长，使用的公式如下：
下标n代表本次迭代，而下标n-1代表前一次迭代。本发明使用如下公式来确定在上述两个步长中选择哪一个：
其中，κ是一个小于1的正数。
(3.3)停止判据
本发明通过判断全变分的作用和置信项的作用是否达到平衡来判断迭代算法是否停止。确定这两种作用的公式如下：
其中，diag(x)是用来产生对角阵的函数，x上的元素会被填充到对角线上。findicator是一个每当f不等于0就取1的指示函数(等于0则取0)，确定了这两种作用后，按照下述公式计算停止判据cα：
在理想情况下，cα＝-1的时候算法应该停止，但是由于噪声等非理想因素的存在，本发明在calpha小于某一设定阈值并保持一段时间后终止算法。
(4)将求得的偏置场叠加到待重建图像上，最终实现图像域阴影伪影修正。
由于在本发明中，轮寻更新偏置场所需要的迭代次数远小于轮寻更新目标函数所需要的迭代次数。为了减少计算负担，本发明并没有在每次迭代中对偏置场进行更新，而是每间隔固定次数的迭代后再对偏置场进行更新。如果在某一次迭代中偏置场没有被更新，那么算法将沿用上一次迭代中的偏置场。
偏置场的更新是基于公式(2)的，也就是说，模板图像或者重建图像的更新都会造成偏置场的更新。本发明中更新模板图像的频率小于更新偏置场的频率。但是由于重建图像在每次迭代中都会变化，所以在不更新模板图像的时候也可以更新偏置场。
实施例
1.在三个实施例中，选择的停止条件阈值为：calpha＜-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次。
以上所述的具体实施方式对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的最优选实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种伪影修正辅助的CBCT迭代重建方法，其特征在于，包括以下步骤：
(1)使用病人或者模体CBCT投影数据重建出初始图像f0；
(2)建立含有偏置场项的迭代重建框架：引入偏置场项，将伪影修正的概念引入到迭代重建的理论框架中去；将现实中受到伪影污染的重建图像和理想情况下没有伪影污染，满足分段常数性质的图像之间的差距称为偏置场；在引入偏置场后，重建算法的目标函数如下所示：
<mrow> <mo>(</mo> <mi>f</mi> <mo>,</mo> <msub> <mi>f</mi> <mrow> <mi>b</mi> <mi>i</mi> <mi>a</mi> <mi>s</mi> </mrow> </msub> <mo>)</mo> <mo>=</mo> <mi>arg</mi> <mi>min</mi> <mo>&amp;lsqb;</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <mi>M</mi> <mi>f</mi> <mo>-</mo> <mi>b</mi> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>+</mo> <mi>&amp;lambda;</mi> <mo>|</mo> <mo>|</mo> <mi>f</mi> <mo>+</mo> <msub> <mi>f</mi> <mrow> <mi>b</mi> <mi>i</mi> <mi>a</mi> <mi>s</mi> </mrow> </msub> <mo>|</mo> <msub> <mo>|</mo> <mrow> <mi>T</mi> <mi>V</mi> </mrow> </msub> <mo>&amp;rsqb;</mo> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
其中f是待重建图像，fbias是需要求取的偏置场；M是写成矩阵形式的投影操作，也被称为正投矩阵，Mf代表对重建图像进行正投操作；b代表原始投影转换为线积分之后的数据，λ为正则化项因子；‖■‖2代表求取二范数，‖■‖TV代表求取全变分，为置信项；
(3)求解迭代重建框架的目标函数，得到待重建图像f以及偏置场fbias：使用轮寻的思想来求解目标函数，在每次迭代的过程中，首先计算偏置场，然后固定偏置场，将目标函数转变为单变量优化问题(2)来进行求解：
<mrow> <mi>f</mi> <mo>=</mo> <mi>arg</mi> <mi>min</mi> <mo>&amp;lsqb;</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <mi>M</mi> <mi>f</mi> <mo>-</mo> <mi>b</mi> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>+</mo> <mi>&amp;lambda;</mi> <mo>|</mo> <mo>|</mo> <mi>f</mi> <mo>+</mo> <msub> <mi>f</mi> <mrow> <mi>b</mi> <mi>i</mi> <mi>a</mi> <mi>s</mi> </mrow> </msub> <mo>|</mo> <msub> <mo>|</mo> <mrow> <mi>T</mi> <mi>V</mi> </mrow> </msub> <mo>&amp;rsqb;</mo> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
(4)将求得的偏置场fbias叠加到待重建图像f上，最终实现图像域阴影伪影修正。
2.根据权利要求1所述一种伪影修正辅助的CBCT迭代重建方法，其特征在于，所述步骤(1)中使用滤波反投影技术重建出初始图像f0。
3.根据权利要求1所述一种伪影修正辅助的CBCT迭代重建方法，其特征在于，所述步骤(3)具体包括以下子步骤：
(3.1)通过公式(3)计算偏置场：
fbias＝H(fseg-f), (3)
其中，fseg是在待重建图像上进行图像分割而获得的模板图像，模板图像的灰度值被填充为不同组织CT数的标准值，通过图像分割将高对比度物质和软组织区分开来；H是一个低通且连续的滤波器，在保持偏置场强度的情况下将偏置场从残差图像中提取出来；残差图像指模板图像与待重建图像的差距；
(3.2)使用GP‐BB方法进行目标函数最小化，具体包括以下步骤：
(3.2.1)按照如下公式进行目标函数梯度g的求取：
<mrow> <mi>g</mi> <mo>=</mo> <mo>&amp;dtri;</mo> <mo>|</mo> <mo>|</mo> <mi>f</mi> <mo>+</mo> <msub> <mi>f</mi> <mrow> <mi>b</mi> <mi>i</mi> <mi>a</mi> <mi>s</mi> </mrow> </msub> <mo>|</mo> <msub> <mo>|</mo> <mrow> <mi>T</mi> <mi>V</mi> </mrow> </msub> <mo>+</mo> <msup> <mi>M</mi> <mi>T</mi> </msup> <mrow> <mo>(</mo> <mi>M</mi> <mi>f</mi> <mo>-</mo> <mi>b</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mo>,</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
公式中T是用来进行矩阵转置的算子，使用如下公式来进行全变分的求导：
<mrow> <mtable> <mtr> <mtd> <mrow> <mo>&amp;dtri;</mo> <mo>|</mo> <mo>|</mo> <mi>F</mi> <mo>|</mo> <msub> <mo>|</mo> <mrow> <mi>T</mi> <mi>V</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <mn>3</mn> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> </mrow> <mrow> <mi>u</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>-</mo> <mfrac> <mrow> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> </mrow> <mrow> <mi>u</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>-</mo> <mfrac> <mrow> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> </mrow> <mrow> <mi>u</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>-</mo> <mfrac> <mrow> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> </mrow> <mrow> <mi>u</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>=</mo> <mfrac> <mrow> <mn>3</mn> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> </mrow> <msqrt> <mrow> <mi>&amp;delta;</mi> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>-</mo> <mfrac> <mrow> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> </mrow> <msqrt> <mrow> <mi>&amp;delta;</mi> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>m</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>-</mo> <mfrac> <mrow> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> </mrow> <msqrt> <mrow> <mi>&amp;delta;</mi> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>-</mo> <mfrac> <mrow> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> </mrow> <msqrt> <mrow> <mi>&amp;delta;</mi> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>F</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>
公式中的δ是一个小的正数；
GP算法使用如下公式来对重建图像进行更新：
fn+1＝max(fn-αnpn,0), (6)
其中，αn表示第n次迭代中的步长，投影后的梯度被记作pn，计算公式如下：
<mrow> <msub> <mi>p</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>g</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>i</mi> <mi>f</mi> <mi> </mi> <msub> <mi>g</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>&amp;le;</mo> <mn>0</mn> <mo>,</mo> <mi>o</mi> <mi>r</mi> <mi> </mi> <msub> <mi>f</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>&gt;</mo> <mn>0</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>0</mn> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>o</mi> <mi>t</mi> <mi>h</mi> <mi>e</mi> <mi>r</mi> <mi>w</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>
其中l是体素点的位置坐标；
(3.2.2)使用BB算法来解析计算每一次迭代步长的α，在每次迭代中计算出两个步长，公式如下：
<mrow> <msubsup> <mi>&amp;alpha;</mi> <mi>n</mi> <mrow> <mi>B</mi> <mi>B</mi> <mn>1</mn> </mrow> </msubsup> <mo>=</mo> <mfrac> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mi>n</mi> </msub> <mo>-</mo> <msub> <mi>f</mi> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mi>n</mi> </msub> <mo>-</mo> <msub> <mi>f</mi> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mi>n</mi> </msub> <mo>-</mo> <msub> <mi>f</mi> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mi>n</mi> </msub> <mo>-</mo> <msub> <mi>p</mi> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>
<mrow> <msubsup> <mi>&amp;alpha;</mi> <mi>n</mi> <mrow> <mi>B</mi> <mi>B</mi> <mn>2</mn> </mrow> </msubsup> <mo>=</mo> <mfrac> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mi>n</mi> </msub> <mo>-</mo> <msub> <mi>f</mi> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mi>n</mi> </msub> <mo>-</mo> <msub> <mi>p</mi> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mi>n</mi> </msub> <mo>-</mo> <msub> <mi>p</mi> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mi>n</mi> </msub> <mo>-</mo> <msub> <mi>p</mi> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>
下标n代表本次迭代，而下标n-1代表前一次迭代；
在两个步长使用公式(10)选择一个步长：
<mrow> <msub> <mi>&amp;alpha;</mi> <mi>n</mi> </msub> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msubsup> <mi>&amp;alpha;</mi> <mi>n</mi> <mrow> <mi>B</mi> <mi>B</mi> <mn>1</mn> </mrow> </msubsup> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <msubsup> <mi>if&amp;alpha;</mi> <mi>n</mi> <mrow> <mi>B</mi> <mi>B</mi> <mn>1</mn> </mrow> </msubsup> <msubsup> <mi>&amp;alpha;</mi> <mi>n</mi> <mrow> <mi>B</mi> <mi>B</mi> <mn>2</mn> </mrow> </msubsup> <mo>&lt;</mo> <mi>&amp;kappa;</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msubsup> <mi>&amp;alpha;</mi> <mi>n</mi> <mrow> <mi>B</mi> <mi>B</mi> <mn>2</mn> </mrow> </msubsup> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>o</mi> <mi>t</mi> <mi>h</mi> <mi>e</mi> <mi>r</mi> <mi>w</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow>
其中，κ是一个小于1的正数；
(3.3)停止判据：通过判断全变分的作用和置信项的作用是否达到平衡来判断迭代算法是否停止；确定这两种作用的公式如下：
<mrow> <msub> <mi>d</mi> <mrow> <mi>T</mi> <mi>V</mi> </mrow> </msub> <mo>=</mo> <mi>d</mi> <mi>i</mi> <mi>a</mi> <mi>g</mi> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mrow> <mi>i</mi> <mi>n</mi> <mi>d</mi> <mi>i</mi> <mi>c</mi> <mi>a</mi> <mi>t</mi> <mi>o</mi> <mi>r</mi> </mrow> </msub> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mo>&amp;dtri;</mo> <mo>|</mo> <mo>|</mo> <mi>f</mi> <mo>+</mo> <msub> <mi>f</mi> <mrow> <mi>b</mi> <mi>i</mi> <mi>a</mi> <mi>s</mi> </mrow> </msub> <mo>|</mo> <msub> <mo>|</mo> <mrow> <mi>T</mi> <mi>V</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow>
<mrow> <msub> <mi>d</mi> <mrow> <mi>d</mi> <mi>a</mi> <mi>t</mi> <mi>a</mi> </mrow> </msub> <mo>=</mo> <mi>d</mi> <mi>i</mi> <mi>a</mi> <mi>g</mi> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mrow> <mi>i</mi> <mi>n</mi> <mi>d</mi> <mi>i</mi> <mi>c</mi> <mi>a</mi> <mi>t</mi> <mi>o</mi> <mi>r</mi> </mrow> </msub> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mo>&amp;dtri;</mo> <mo>|</mo> <mo>|</mo> <mi>M</mi> <mi>f</mi> <mo>-</mo> <mi>b</mi> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow>
其中，diag(x)是用来产生对角阵的函数，x上的元素会被填充到对角线上；findicator是一个每当f不等于0就取1的指示函数，等于0则取0，确定了这两种作用后，按照下述公式计算停止判据cα：
<mrow> <msub> <mi>c</mi> <mi>&amp;alpha;</mi> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mi>d</mi> <mrow> <mi>T</mi> <mi>V</mi> </mrow> </msub> <mo>&amp;CenterDot;</mo> <msub> <mi>d</mi> <mrow> <mi>d</mi> <mi>a</mi> <mi>t</mi> <mi>a</mi> </mrow> </msub> </mrow> <mrow> <mo>|</mo> <msub> <mi>d</mi> <mrow> <mi>T</mi> <mi>V</mi> </mrow> </msub> <mo>|</mo> <mo>|</mo> <msub> <mi>d</mi> <mrow> <mi>d</mi> <mi>a</mi> <mi>t</mi> <mi>a</mi> </mrow> </msub> <mo>|</mo> </mrow> </mfrac> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow>
在cα小于设定阈值并保持一段时间后终止算法。
4.根据权利要求3所述一种伪影修正辅助的CBCT迭代重建方法，其特征在于，所述步骤(3.1)中，使用阈值化算法和二相水平集算法结合的方法来进行图像分割。
5.根据权利要求4所述一种伪影修正辅助的CBCT迭代重建方法，其特征在于，所述二相水平集算法公式如下所示：
<mrow> <mtable> <mtr> <mtd> <mrow> <mi>F</mi> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>,</mo> <mi>c</mi> <mo>,</mo> <mi>b</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>&amp;Integral;</mo> <mrow> <mo>(</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>2</mn> </munderover> <mo>&amp;Integral;</mo> <mi>K</mi> <mo>(</mo> <mi>y</mi> <mo>-</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>|</mo> <mi>I</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>b</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> <msub> <mi>c</mi> <mi>i</mi> </msub> <msup> <mo>|</mo> <mn>2</mn> </msup> <msub> <mi>M</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>(</mo> <mi>x</mi> <mo>)</mo> <mo>)</mo> </mrow> <mi>d</mi> <mi>x</mi> <mo>)</mo> <mi>d</mi> <mi>y</mi> <mo>+</mo> <mi>v</mi> <mo>&amp;Integral;</mo> <mo>|</mo> <mo>&amp;dtri;</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>)</mo> </mrow> <mo>|</mo> <mi>d</mi> <mi>x</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mi>u</mi> <mo>&amp;Integral;</mo> <mi>p</mi> <mo>|</mo> <mo>&amp;dtri;</mo> <mi>&amp;phi;</mi> <mo>|</mo> <mi>d</mi> <mi>x</mi> <mo>,</mo> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> </mrow>
式中，φ是水平集函数，该函数通过其函数值的符号进行两个不同区域的分割，I是需要进行分割的图像，ci是需要填充到第i个区域的函数值，y代表图像上每一个点，而x代表y的邻域内的每个点；b是用来对图像本身的不均匀度进行补偿的补偿项；K(y-x)是一个非负的窗函数，该函数在不属于y的邻域的部分取0；Mi是对第i种组织用符号函数定义的成员函数，在属于这个组织的部分该函数会取为1；u，v是用来进行效果调节的参量，通过对这两个参量进行调节来实现置信项和平滑项之间的平衡；为了将公式(14)最小化，依次基于φ,c,b来执行梯度降算法，公式(14)在被最小化之后就可以提供分割区域，将特定组织的标准CT值填入到分割区域中，从而得到模板图像。
6.根据权利要求5所述一种伪影修正辅助的CBCT迭代重建方法，其特征在于，所述步骤(3.1)中，所述滤波器为二维的Savitzky‐Golay滤波器。
7.根据权利要求1所述一种伪影修正辅助的CBCT迭代重建方法，其特征在于，所述步骤(3)中，并没有在每次迭代中对偏置场进行更新，而是每间隔固定次数的迭代后再对偏置场进行更新；如果在某一次迭代中偏置场没有被更新，那么算法将沿用上一次迭代中的偏置场。
CN201510980575.9A 2015-12-23 2015-12-23 伪影修正辅助的cbct迭代重建方法 Active CN105631909B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510980575.9A CN105631909B (zh) 2015-12-23 2015-12-23 伪影修正辅助的cbct迭代重建方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510980575.9A CN105631909B (zh) 2015-12-23 2015-12-23 伪影修正辅助的cbct迭代重建方法

Publications (2)

Publication Number Publication Date
CN105631909A CN105631909A (zh) 2016-06-01
CN105631909B true CN105631909B (zh) 2018-05-29

Family

ID=56046796

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510980575.9A Active CN105631909B (zh) 2015-12-23 2015-12-23 伪影修正辅助的cbct迭代重建方法

Country Status (1)

Country Link
CN (1) CN105631909B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106408543B (zh) * 2016-11-07 2019-04-23 深圳先进技术研究院 一种用于锥束ct图像散射修正的阻挡光栅优化方法及装置
WO2018126434A1 (zh) * 2017-01-06 2018-07-12 深圳先进技术研究院 Ct图像阴影校正方法、装置及电子设备
CN107871332A (zh) * 2017-11-09 2018-04-03 南京邮电大学 一种基于残差学习的ct稀疏重建伪影校正方法及系统
CN108109185B (zh) * 2017-12-18 2021-07-20 上海联影医疗科技股份有限公司 一种生成用于消除ct伪影的校正系数的方法，以及一种基于校正系数消除ct伪影的方法
US10922855B2 (en) 2017-11-30 2021-02-16 Shanghai United Imaging Healthcare Co., Ltd. Systems and methods for determining at least one artifact calibration coefficient
CN109118439B (zh) * 2018-07-03 2022-02-18 浙江大学 基于线积分的锥束ct去模糊方法
WO2020061814A1 (en) * 2018-09-26 2020-04-02 Shanghai United Imaging Healthcare Co., Ltd. Systems and methods for image generation
CN110097517B (zh) * 2019-04-28 2022-12-27 东软医疗系统股份有限公司 去除图像伪影的方法及装置
CN113256500B (zh) * 2021-07-02 2021-10-01 北京大学第三医院（北京大学第三临床医学院） 一种多模态图像合成的深度学习神经网络模型系统

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101777177A (zh) * 2009-12-29 2010-07-14 上海维宏电子科技有限公司 基于衰减滤波的ct图像去金属伪影混合重建法
CN104665862A (zh) * 2015-02-16 2015-06-03 清华大学 在cbct中消除几何伪影的方法以及使用该方法的cbct系统

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7747057B2 (en) * 2006-05-26 2010-06-29 General Electric Company Methods and apparatus for BIS correction

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101777177A (zh) * 2009-12-29 2010-07-14 上海维宏电子科技有限公司 基于衰减滤波的ct图像去金属伪影混合重建法
CN104665862A (zh) * 2015-02-16 2015-06-03 清华大学 在cbct中消除几何伪影的方法以及使用该方法的cbct系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Image-domain shading correction for cone-beam CT without prior patient information;Qiyong Fan et al;《JOURNAL OF APPLIED CLINICAL MEDICAL PHYSICS》;20151108;第16卷(第6期);第65-75页 *
一种去除医学图像中金属伪影的新方法;李艳丽 等;《电视技术》;20131231;第37卷(第17期);第59-61页 *

Also Published As

Publication number Publication date
CN105631909A (zh) 2016-06-01

Similar Documents

Publication Publication Date Title
CN105631909B (zh) 伪影修正辅助的cbct迭代重建方法
Xu et al. A practical cone-beam CT scatter correction method with optimized Monte Carlo simulations for image-guided radiation therapy
Li et al. Motion correction for improved target localization with on-board cone-beam computed tomography
Wei et al. X-ray CT high-density artefact suppression in the presence of bones
CN105188543B (zh) X射线ct装置、重构运算装置以及重构运算方法
Zhao et al. Patient-specific scatter correction for flat-panel detector-based cone-beam CT imaging
CN106960429B (zh) 一种ct图像金属伪影校正方法及装置
CN104751429B (zh) 一种基于字典学习的低剂量能谱ct图像处理方法
US10282872B2 (en) Noise reduction in tomograms
Riblett et al. Data‐driven respiratory motion compensation for four‐dimensional cone‐beam computed tomography (4D‐CBCT) using groupwise deformable registration
Marchant et al. Reduction of motion artefacts in on-board cone beam CT by warping of projection images
CN103617598B (zh) 一种基于轨迹的ct图像金属伪影去除方法
Lin et al. An efficient polyenergetic SART (pSART) reconstruction algorithm for quantitative myocardial CT perfusion
Qi et al. Performance studies of four-dimensional cone beam computed tomography
Zhang et al. PET image reconstruction using a cascading back-projection neural network
CN110335325A (zh) 一种ct图像重建方法及其系统
Hansen et al. Fast 4D cone-beam CT from 60 s acquisitions
Van Der Heyden et al. A Monte Carlo based scatter removal method for non-isocentric cone-beam CT acquisitions using a deep convolutional autoencoder
KR102493193B1 (ko) 다중 수집을 통해 콘트라스트를 정규화하는 방법 및 시스템
Chee et al. McSART: an iterative model-based, motion-compensated SART algorithm for CBCT reconstruction
Soimu et al. A novel approach for distortion correction for X-ray image intensifiers
Nuyts et al. Evaluation of maximum-likelihood based attenuation correction in positron emission tomography
Peterlik et al. Reducing residual‐motion artifacts in iterative 3D CBCT reconstruction in image‐guided radiation therapy
Shi et al. A Virtual Monochromatic Imaging Method for Spectral CT Based on Wasserstein Generative Adversarial Network With a Hybrid Loss.
Abadi et al. Emphysema quantifications with CT scan: Assessing the effects of acquisition protocols and imaging parameters using virtual imaging trials

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant