0 150化U。
[0049] 图3为本发明在Ca化han@600模体高对比度线对层上的实施结果，（a):使用655个 投影进行传统FBP重建的结果，（b):首先使用655个投影进行传统FBP重建，然后使用计划CT 图像进行散射修正的结果，（C):使用传统迭代重建算法基于92个投影得到的重建结果， (d):使用本发明中的重建算法基于92个投影得到的重建结果。在(b)中使用黑色虚线框标 注的部分在^.曰），^.6)，^.(3)，^.(1)中被放大显示，（曰），（(3)，（(1)中的显示窗为[-250 250化U，（b)中的显示窗为[-50 450化U。
【具体实施方式】
[0050] 下面结合附图和具体实施例对本发明作进一步详细说明。
[0051] 本发明提出的一种伪影修正辅助的CBCT迭代重建方法，包括W下步骤：
[0052] (1)使用滤波反投影技术产生初始图像
[0053] 将实际测得的病人或者模体CBCT投影数据使用滤波反投影技术重建出初始图像 fo,注意，在减少剂量的情况下，采集到的投影数据相对较少，滤波反投影算法重建得到的 图像具有较多条纹和噪声污染。
[0054] (2)建立含有偏置场项的迭代重建框架
[0055] 通过引入偏置场项，将伪影修正的概念引入到迭代重建的理论框架中去。将现实 中受到伪影污染的重建图像和理想情况下没有伪影污染，满足分段常数性质的图像之间的 差距称为偏置场。在引入偏置场后，本发明中重建算法的目标函数如下所示：
[0化6]
他
[0057] 其中f是需要求取的重建图像，fbias是需要求取的偏置场。Μ是写成矩阵形式的投 影，也被称为正投矩阵，Mf代表对重建图像进行正投操作，该操作为线性因而可W被写成矩 阵相乘的形式。b代表原始投影转换为线积分之后的数据，λ为手工选取的正则化项因子，用 来控制重建图像的平滑程度。II ·ΙΙ2代表求取二范数，II ·ΙΙτν代表求取全变分，全变分在数学 被定义为空间梯度图像的一范数，ΙΙΜ/ - fclli为置信项。
[005引通过求解该优化问题，就可W同时得到待重建图像f，W及偏置场fbias。随后简单 的将偏置场叠加到待重建图像上就可W完成精确的图像域阴影伪影修正了。
[0059] (3)求解迭代重建框架的目标函数，得到待重建图像fW及偏置场fbias
[0060] 本发明使用的目标函数（1)由于有两个自变量的存在，很难直接进行最小化操作。 因此本发明提出使用轮寻的思想来进行求解。也就是在每次迭代的过程中，首先使用信号 与图像处理的手段来计算偏置场，偏置场在计算完成后即被固定为不变，从而将目标函数 转变为单变量优化问题来进行求解。本发明计算偏置场的方法如下所示：
[0061] fbias = H(fseg-f), (2)
[0062] 其中，fseg是在待重建图像上进行图像分割而获得的模板图像，模板图像的灰度值 被填充为不同组织CT数的标准值，因而不包含伪影信息。f是待重建图像，Η是一个低通且连 续的滤波器，该滤波器可W将fseg和fp之间的插值图像进行滤波，从而得到估计的偏置场。
[0063] W上就是轮寻求解的第一步，也就是在每次迭代开始的时候首先对偏置场进行更 新。在进行过偏置场的更新后，本发明将偏置场固定为不变，运样，只剩下一个变量的优化 函数如公式(3)可W按照传统迭代重建的算法进行求解，运就是轮寻求解的第二步，下面本 发明将分别对运两步进行介绍。
[0064]
巧）
[0065] 在第一次迭代时，将待重建图像f置为初始图像fo;
[0066] (3.1)偏置场求解
[0067] (3.1.1)图像分割算法
[0068] 为了对公式(2)中的fseg进行计算，本发明使用阔值化算法和二相水平集算法结合 的方法来进行图像分割。本发明认为，骨头和空气等高对比度物质，由于其CT数和软组织 (包括脂肪，肌肉等)具有较大程度的区别，可W直接通过手工设定阔值的做法将其分割出 来，也就是CT数大于某一个设定值的部分认定为骨骼，CT数小于某一个设定值的部分认定 为空气。
[0069] 在完成高对比物体的分割之后，本发明使用二相水平集算法进行软组织的分割， 也就是进行脂肪和肌肉的分割。本发明中的二相水平集算法公式如下所示：
[0070]
[0071] 式中，Φ是水平集函数，该函数通过其函数值的符号进行两个不同区域的分割，I 是需要进行分割的图像，Cl是需要填充到第i个区域的函数值，y代表图像上每一个点，而X 代表y的邻域内的每个点。b是用来对图像本身的不均匀度进行补偿的补偿项。K(y-x)是一 个非负的窗函数，该函数在不属于y的邻域的部分取OdMi是对第i种组织用符号函数定义的 成员函数，在属于运个组织的部分该函数会取为1dU，v是用来进行效果调节的参量，本发明 通过对运两个参量进行调节来实现置信项和平滑项之间的平衡。为了将公式(4)最小化，本 发明依次的基于Φ，(3,6来执行梯度降算法，公式(4)在被最小化之后就可W提供分割区域， 本发明将特定组织的标准CT值填入到分割区域中，从而得到模板图像。
[0072] (3.1.2)滤波器设计
[0073] 本发明使用的滤波器是二维的Savitzky-Golay滤波器(也就是公式(2)中的H)，该 滤波器具有较好的轮廓保持能力，因而可W在保持偏置场强度的情况下将偏置场从残差图 像中提取出来。残差图像指模板图像与待重建图像的差距。该滤波器连续的使用相邻点信 息，基于低次多项式来对未知点进行最小二乘拟合。运种效果可W通过卷积来实现，因而具 有较低的计算复杂度。
[0074] (3.2)使用GP-BB方法进行目标函数最小化
[0075] 求得偏置场后的下一个步骤是固定偏置场，来最小化公式(3)。本发明使用GP-BB 方法来进行目标函数的最小化，具体如下：
[0076] 首先，按照如下公式进行目标函数(公式(3))梯度g的求取：
[0077]
,巧）
[0078] 公式中T是用来进行矩阵转置的算子，本发明中使用如下公式来进行全变分的求 导：
[0079]
[0080] 公式中的δ是一个小的正数，用来防止除0的错误。GP算法使用如下公式来对重建 图像进行更新：
[0081] fn+i=max(fn-anPn，〇)， (7)
[0082] 其中，α是每次迭代中的步长，投影后的梯度被记作pn，计算公式如下：
[0083]
货）
[0084] 其中1是体素点的位置坐标。
[0085] 本发明使用BB算法来解析的计算每一次迭代步长的α，在每次迭代中，本发明计算 出两个步长，使用的公式如下：
[0088] 下标η代表本次迭代，而下标n-1代表前一次迭代。本发明使用如下公式来确定在 上述两个步长中选择哪一个：
[0089]
(1： I
[0090] 其中，K是一个小于1的正数。
[0091] (3.3)停止判据
[0092] 本发明通过判断全变分的作用和置信项的作用是否达到平衡来判断迭代算法是 否停止。确定运两种作用的公式如下：
[0095] 其中，diag(x)是用来产生对角阵的函数，X上的元素会被填充到对角线上。 findwatDr是一个每当f不等于0就取1的指示函数(等于0则取0)，确定了运两种作用后，按照 下述公式计算停止判据ca:
[0096]
(14)
[0097] 在理想情况下，ca = -l的时候算法应该停止，但是由于噪声等非理想因素的存在， 本发明在Ca邮a小于某一设定阔值并保持一段时间后终止算法。
[009