Iterative CBCT reconstruction using Hessian penalty

TV正则作为惩罚项具有明显的阶梯效应，会使图像不自然，这篇文章为了解决这个问题，提出了利用图像Hessian矩阵的Frobenius范数的二阶导数作为惩罚项进行CBCT重建，可以有效抑制TV惩罚项的阶梯效应。
TV惩罚项优先最小化一阶导数，因此往往具有分段函数结果。在这次研究当中，我们提出了使用hessian惩罚来进行CBCT重建。这就会涉及到图像的hessian矩阵的F范数来抑制TV惩罚中的阶梯效应。其中，Hessian惩罚是二阶导数惩罚。这背后的动机是二阶导数对于相邻像素之间的绝对差有较弱的惩罚，并且允许分段地平滑重建结果。
对于之前的重建方法，大多数的惩罚项都是二次式的，因此相关目标函数的优化可以直接进行，例如高斯-赛德尔和共轭梯度法。
但是不幸的是hessian惩罚不是二次式的，这就会让目标函数的最小化具有挑战性。作者发展了一种有效的算法来最小化目标函数—MM算法。
图像去噪的能量泛函：

第一项式加权最小二乘元作为（WLS）或数据保真项。对角线矩阵Σ的元素在WLS损失函数中起着加权的角色，决定着每项的贡献。第二项是先验约束或一个惩罚项。图像重建任务是通过最小化具有正约束的目标函数来找到一个衰减系数图像μ
hessian惩罚
由于上式中的先验约束强制了平滑约束，图像的TV最小化在CBCT重建中在噪声抑制和边缘保持中显示了良好的作用。让μ是连续可微的三维图像。TV惩罚可以定义为：

TV的阶梯效应是由于TV惩罚的一阶导数会倾向于最小化相邻元素之间的差异，由此导致了分段常数图像。尝试在惩罚项中选择二阶导数，提出的Hessian惩罚定义为：

矩阵A的Frobenius范数定义为矩阵A各项元素的绝对值平方的总和，即

这种Hessian惩罚满足凸性（产生于线性算子范数的积分）、均匀性、旋转和平移不变性。
这里控制（每个特征方向）步长的，有两个东西：原来的一阶梯度和对应的Hessian矩阵特征值。
Hessian惩罚满足：

其中，“u”和“v”是以球面坐标表示的单位范数向量：

其中

Du1和Dv1分别表示由u和v指定方向的一阶导数。
足迹分离向前投影
之前，我们将矩阵A的元素Aij定义为从X射线源到图像像素μj的第i个检测器的投影射线交点长度。但是这种线交模型对于我们像素级更新策略是不精确的，计算也复杂。而足迹分离映射算法的主要思想是真正的像素足迹函数可以用二维可分离函数近似。SF-TR模型在横向使用梯形函数，纵向使用矩形函数，投影矩阵A的元素可以表示为：

其中s，t分别是探测器平面的横向和纵向坐标，r是投影的极角，(x,y,z)是像素的坐标，F1是横向的梯形函数，F2是纵向的矩形函数。足迹函数F1和F2近似与真实足迹的形状相同。
目标函数最小化
衍生滤波器
在离散情况中，TV惩罚可以被写为：

提出的Hessian惩罚可以写为：

上式的精确形式可以由被表示为：

其中Li是在像素点(x,y,z)产生了所有可能的二阶导数组成的滤波器。
主要化-最小化（MM）方法
这个方法的主要思想是得到一系列二次函数，作为原目标函数上界，并收敛到最佳状态。
对于最小化问题：

不是直接最小化Ф(u)，而是线用一系列Q(μ；μ(t))来作为它的上界，给定初始值为μ(0)，通过方程：

得到μ(t+1)，然后利用不等式：
得到基于集合{μ(0),μ(1),μ(2)…μ(t)…}的单调不增的Ф(u)。在条件Q(μ；μ’)对于μ和μ‘是连续的下，Ф(u(t))单调收敛到一个固定的点μ对应的Ф(u)上。因此，我们可以得到下列不等关系：

这种不平等构成了主要化-最小化（MM）方法的基础。
惩罚函数的二次主项
为了表示出TV和Hessian的惩罚项，我们利用不等式：

当且仅当g(x)=g(y)时取等号。
由于

所以可以得到TV惩罚的主项：

同理也可以得到Hessian惩罚的主项：

由于数据保真项为：

所以得到完整的二次主项：

让上式对μ求导为0，主项的最小化等效于求下面的线性方程组：

上述方程可以通过高斯-赛德尔策略有效解决。注意到对于一个固定的像素μj，交叉项μjμi源于衍生操作，并且仅仅在对于TV的主项当μj与μi相距不超过一个像素的时候存在，对于Hessian的主项当μj与μi相距不超过两个像素的时候存在。因此，对于TV的主项和Hessian的主项可以等同表示为：

材料和评价
比较算法
除了TV惩罚，还将Hessian惩罚与Huber惩罚进行了比较，后者被广泛用于CT重建中。Huber惩罚形式如下：
当相邻像素之间的差异的绝对值大于阈值时采用线性惩罚，对于出现在过渡区域的小差异采用二次惩罚。选择合适的阈值就可以对CT重建有较好的效果。
评价指标
为了量化重建图像的质量，我们选择了峰值信噪比（PSNR）、信噪比（ISNR）、对比噪声比（CNR）和结构相似度（SSIM）。在我们的实验中，FDK重建图像被用作迭代重建算法的初始猜测。
PSNR：

其中MSE是重建图像与参考图像之间的均方误差，而μmax是图像的最大可能值。较高的PSNR意味着重建图像与指定图像之间的误差较小。
ISNR：

其中MSEin和MSEout是在FDK重建图像和指定图像之间的均方误差，和迭代重建图像和指定图像之间的均方误差。ISNR越高意味着迭代重建图像比FDK重建图像质量上提高越多。
CNR：

其中μROI和μRef是选择的ROI和指定均匀背景区域的平均值，σROI和σRef是ROI和指定区域标准差。CNR越高表示能更好的保护信号不受噪声的影响。
SSIM：

其中a和b是两幅图像中尺寸为8×8像素的局部窗口。这两个窗口位置相同，其中μa和σa，μb和σb是他们的均值和方差。σab是两个窗口之间的协方差。C1和C2是避免不稳定的两个常数。在此可以选择
当两个窗口在重建图像和指定图像中移动时，我妈得到一个SSIM映射，一般，我妈使用均值SSIM（MSSIM）来评价整个图像。