在小动物计算机断层扫描(CT)实验中, 因需考虑小动物存活率以及实验的连续性等问题, 一般较少采用高剂量的X射线进行实验; 而低剂量的X射线会导致重建图像被噪声污染, 影响图像质量, 不利于后续实验分析。为解决此问题, 本文介绍了一种基于全局字典学习的降噪方法, 并将其应用于提升低剂量小动物CT重建图像质量的研究中。针对真实的小动物CT重建数据, 选择高剂量的小动物CT重建图像作为训练样本, 利用逐列更新的字典学习算法(K-SVD), 构建包含图像信息的全局字典; 利用正交匹配追踪算法(OMP)将低剂量重建图像利用全局字典进行稀疏分解, 分离噪声, 最后将重建图像复原, 达到降噪、提升图像质量、降低小动物CT实验的拍摄剂量、提高小动物存活率的目的。实验结果表明, 本文提出的方法能够有效减少低剂量动物CT图像的噪声, 并能够较好地保留图像细节。
引用本文： 李中源, 李光, 孙翌, 陈功, 罗守华. 一种基于全局字典学习的小动物低剂量计算机断层扫描降噪方法. 生物医学工程学杂志, 2016, 33(2): 279-286. doi: 10.7507/1001-5515.20160048 复制
引言
动物模型是现代生物医学研究中的重要组成部分。通过对动物接种肿瘤观察其生理变化、利用各类新型造影剂进行动物实验、研究药物对动物生理的影响等[1]，能够有效地帮助人类认识疾病发生规律，并采取相应的防治措施。近年来，各类医学影像技术广泛地应用于各类动物实验中，其中计算机断层扫描(computed tomography，CT)技术由于具有分辨率高、成像时间短等特点，成为动物医学影像研究的重要手段之一。
与传统的临床CT相比，小动物CT具有更高的分辨率，但却受到射线剂量方面的限制。小动物本身十分脆弱，高剂量的X射线辐射往往会对诸如小鼠之类的动物造成损伤，甚至死亡，不利于实验的连续性，因此低剂量小动物CT将是未来研究的一个发展趋势。但是低剂量的CT重建图像会被噪声干扰，影响图像质量，因而如何提高低剂量小动物CT的重建图像质量，成为了亟待解决的问题。
传统的低剂量CT重建问题的解决方案主要分为三类，分别应用于重建过程的三个部分：重建前处理(对CT投影图或者Sinogram图进行降噪处理)、重建(对CT重建过程中的算法进行改进以抑制噪声)以及重建后处理(对CT重建结果进行降噪处理)。在前处理过程中，由于投影图是通过X射线直接扫描得到的，因而噪声类型较好统计，如基于惩罚加权的最小二乘(penalized weighted least-squares，PWLS)算法通过统计方法获得投影图像噪声模型，针对噪声模型进行降噪处理[2]；此外也有对Sinogram图像进行噪声统计，并利用非线性滤波进行降噪处理的办法[3]。由于此类方法应用于重建算法之前，有可能会造成图像信息的丢失或者引发图像的畸变。另外，由于投影图的数据量非常庞大，统计过程耗时较长。而图像重建过程中，对重建算法进行改进的一种基本思路是通过在传统的滤波反投影过程中添加平滑窗函数来获取较为平滑的结果，比如信号处理领域中常见的汉宁窗、海明窗等。此外，也有采用基于统计迭代图像重建的算法，这类算法通过将CT重建图像数据与投影数据相关联，对投影数据的噪声模型建模，在迭代过程中引入权重，达到图像降噪的目的[4-5]。虽然迭代类算法对噪声的抑制作用比较明显，但其计算量庞大、耗时长，在临床上应用并不广泛。在重建图像的后处理中，一般将重建图像视为普通图像，采用传统的线性或非线性滤波方法，例如利用核函数的办法对图像进行降噪滤波，以及利用总变分最小化(total variation，TV)的办法对图像进行降噪等[6-7]；也有学者提出将图像变换到小波域进行降噪处理[8]，但是由于重建图像需要通过反投影等复杂计算才能得到，噪声模型不便统计，传统的图像处理办法很难取得一个好的结果。
近阶段兴起的信号稀疏表示与字典学习等方法为上述问题提供了新的解决方案。字典的概念最早由1993年Mallat等[9]提出。1998年Chen等[10]提出利用一范数评估稀疏性以及利用凸函数求解系数问题，推进了该理论的发展。2001年，Donoho等[11]验证了追踪算法的确定性，使得稀疏表示理论取得长足发展，并广泛应用于图像的降噪、压缩与恢复等各个领域。在图像降噪方面，基本的算法思路可以概括为将图像利用字典稀疏的表达，在给定噪声误差的前提下，利用稀疏化的参数和字典将图像稀疏重构，达到滤除噪声的目的。字典解决方案直接应用于重建图像，保证了图像结果的真实性，也省去了前处理庞大的计算量，十分具有应用前景。
本文利用字典降噪的基本思路，使用高剂量图像作为训练样本，利用逐列更新的字典学习算法(K-means singular value decomposition, K-SVD)获取合适的全局字典，并将其应用于低剂量小动物CT的图像降噪研究之中。
1 算法
在小鼠低剂量CT图像降噪过程中，主要会涉及全局字典的训练学习、图像降噪算法实现以及小鼠CT应用框架搭建等具体操作问题。
1.1 K-SVD字典训练
利用K-SVD算法进行字典训练的过程可以归结于求解如下目标函数[12]
| $ \mathop {\min }\limits_{D, \alpha } \left\{ {\left\| {\boldsymbol{Y}-\boldsymbol{D}\alpha } \right\|_F^2} \right\}, st.{\forall _i}, \left\| {{\alpha _i}} \right\| \le T $ |
其中，Y是用于训练的样本集，Y的每一个列向量都为一个训练样本；D是初始字典，一般为大小为n×K的一个矩阵，K表示字典总共包含K个原子，n表示样本的维度。常用的初始字典有离散余弦字典或者小波字典；α是对样本集Y进行稀疏分解后获得的稀疏编码系数，T为式(1)的稀疏度约束，表明表达式中每一次稀疏表示的过程使用的字典原子数量不能超过T个。由于K-SVD算法是对字典D逐个原子进行更新，稀疏编码计算也是针对单个训练样本进行计算，因此我们将式(1)细化到每一个信号，重新改写为：
| $ \mathop {\min }\limits_{{\alpha _I}} \left\{ {\left\| {{\boldsymbol{y}_i}-\boldsymbol{D}{\alpha _i}} \right\|_2^2} \right\}, st.\left\| {{\alpha _i}} \right\| \le T $ |
其中i代表样本序号，K-SVD算法的目的在于使得式(2)取得最小值，联合两式可以将目标函数继续修改进行求解：
| $ \left\| {\boldsymbol{Y}-\boldsymbol{D}\alpha } \right\|_F^2=\left\| {\boldsymbol{Y}-\sum\limits_{j=1}^K {{\boldsymbol{d}_j}\alpha _T^j} } \right\|_2^2=\left\| {\left({\boldsymbol{Y}-\sum\limits_{j \ne k} {{\boldsymbol{d}_j}\alpha _T^j} } \right)-{\boldsymbol{d}_k}\alpha _T^k} \right\|_F^2 $ |
其中dj表示第j个字典原子，为一个列向量。令Ek表示，则Ek代表了在字典中去除第k个原子后式(3)所产生的误差，对Ek进行奇异值分解(singular value decomposition，SVD)，求解最大的满秩列向量，并利用该向量对原字典中第k个原子进行更新；将字典中所有原子遍历一次，便可得到迭代一次后更新的字典，为了使字典更为平滑、有效地抑制噪声信息，一般需要迭代多次。
1.2 图像降噪框架
通过K-SVD算法获取合适的全局字典后，则可以利用全局字典对噪声图像进行降噪处理[13]，利用全局字典进行降噪的问题可以表述为：
| $ \left\{ {{{\hat \alpha }_{ij}},\boldsymbol{\hat X}} \right\} = \mathop {{\rm{argmin}}\lambda }\limits_{{\alpha _{ij}},X} \left\| {\boldsymbol{X} - \boldsymbol{Y}} \right\|_2^2 + \sum\limits_{ij} {{\mu _{ij}}} {\left\| {{\alpha _{ij}}} \right\|_0} + \sum\limits_{ij} {\left\| {{\boldsymbol{R}_{ij}}\boldsymbol{X} - \boldsymbol{D}{\alpha _{ij}}} \right\|_2^2} $ |
其中，式(4)右侧第一项衡量的是原始图像X与噪声图像Y之间的整体相似程度，λ为拉格朗日系数；第二项为字典降噪过程中的稀疏度约束，其中μij为约束项的权重系数，由于整图较大，所以需要将原图划分为一块块的小图像并重排为与全局字典的原子维度一致的列向量；第三项中的RijX表示了原图中的第ij幅子图像，Rij是用于提取子图的矩阵，ij为子图像序列下标，Dαij为对应子图像的稀疏表达，第三项用于衡量这两者之间的误差，它们的误差越小越好。
由于已经通过K-SVD算法获取到了合适的全局字典D，为了获取每张子图的稀疏表达系数，需要求解下式最优化问题：
| $ {\alpha _{ij}}=\mathop {{\mathop{\rm a}\nolimits} \hat rgmin{\mu _{ij}}}\limits_\alpha {\left\| \alpha \right\|_0} + \sum\limits_{ij} {\left\| {\boldsymbol{D}{\alpha _{ij}}-{\boldsymbol{R}_{ij}}\boldsymbol{X}} \right\|_2^2} $ |
式(5)一般通过贪婪算法来求解。本文采用正交匹配追踪(orthogonal matching pursuit，OMP)来求解式(5)的问题[14]，其基本思路是对每一次逼近的残差进行新一轮的逼近，直到满足设定的误差限或者稀疏度要求，得到近似解âij。如此，式(5)的求解更新即可转换为如下形式：
| $ \boldsymbol{\hat X}=\mathop {{\mathop{\rm argmin}\nolimits} \lambda }\limits_\boldsymbol{X} \left\| {\boldsymbol{X}-\boldsymbol{Y}} \right\|_2^2 + \sum\limits_{ij} {\left\| {\boldsymbol{D}{{\hat \alpha }_{ij}}-{\boldsymbol{R}_{ij}}\boldsymbol{X}} \right\|_2^2} $ |
这是一个简单的凸函数问题，对其进行求解得到降噪后的图像：
| $ \boldsymbol{\hat X}={\left({\lambda \boldsymbol{I} + \sum\limits_{ij} {\boldsymbol{R}_{ij}^T{\boldsymbol{R}_{ij}}} } \right)^{-1}}\left({\lambda \boldsymbol{Y} + \sum\limits_{ij} {\boldsymbol{R}_{ij}^T\boldsymbol{D}{{\hat \alpha }_{ij}}} } \right) $ |
其中，I为单位矩阵，为降噪后的结果图像。
1.3 基于全局字典降噪的小动物CT低剂量图像降噪框架
在Aharon等[12]的工作中，利用噪声图像自身作为训练集，经过K-SVD算法训练过后的字典在降噪方面比初始字典更有优势。这主要是因为经由噪声图像训练自身得到的字典原子包含了图像的结构信息，这使得图像在做稀疏表达时能够获得更好的稀疏度，从而取得更好的结果。但是在实际应用中，如果对每一幅噪声图像都进行字典训练，无疑大大增加了计算量，并不可取。考虑到同种类的小动物CT图像结构相对比较接近，同种类动物的不同个体在CT图像上的结构差异更多地体现于结构整体的不同，比如组织的轮廓或者位置等，但是其组织及其骨骼等构成成分都类似，因而在图像细节方面表现大多类似，比如组织的边缘、灰度级、骨骼内部细微结构等，因此通过提前对高剂量的同种类CT图像进行字典训练获取包含结构信息的全局字典，对同种类小动物的低剂量噪声图像进行降噪处理，既能达到噪声图像自身训练字典的降噪效果，也能省去单独每幅图像全局字典训练的时间，使得全局字典在低剂量小动物CT上的应用具有了可行性。
因此，基于全局字典降噪的低剂量小动物CT降噪方案可以归纳为如下几个步骤：
(1) 获取训练样本集：利用高剂量X射线拍摄小动物获得质量好的小动物的CT图像；
(2) 获取全局字典：对高剂量小动物的CT图像进行处理，利用K-SVD算法获取合适的全局字典D；
(3) 低剂量图像稀疏分解：将低剂量图像分割为图像子集，利用D以及OMP算法对低剂量图像进行稀疏分解，获取每一个子集的稀疏编码α；
(4) 低剂量图像重构：利用式(7)以及步骤(3)中获得结果对低剂量图像进行重构，获取降噪后的纯净图像。
2 实验与结果分析
本文实验主要分为两个部分，模拟实验以及小动物CT数据实验。为了验证全局字典降噪的有效性，本文采用小波软阈值、TV算法降噪进行对比实验。实验结果通过峰值信噪比(peak signal-noise ratio，PSNR)以及结构相似度(structure similarity index measurement，SSIM)两个评价指标进行衡量[15]。PSNR以及SSIM都是图像客观评价的标准，其中PSNR反映的是最大值信号与噪声之间的关系，单位为dB，PSNR值越大，噪声滤除效果越好；而SSIM则综合考量两幅图像之间亮度、对比度以及结构相似度作为整体评判指标，SSIM值越接近1，代表两幅图像越相似。
2.1 模拟数据实验结果
为了尽可能地接近小动物CT实验的真实情况，模拟实验采用人体胸腔CT图像作为模拟数据(数据来源于文献[16])，图像大小为512×512。CT图像噪声是由于到达探测器的光子数量不同(被扫描物体各部分对X光吸收不一致)并且有暗电流背景噪声等因素影响而形成。为了使模拟实验更加接近真实情况，对CT重建图像的噪声产生过程进行模拟，首先图像进行正投影后，投影图添加泊松噪声再进行重建[17]，重建后的CT图像利用全局字典、小波软阈值以及TV算法进行降噪处理，并计算处理后的结果与无噪声标准重建CT图像之间的PSNR值与SSIM值，观察实验结果。
其中，全局字典由无噪声标准重建图像训练而来，以8×8图像块大小对原图像素点进行采样，总共获取255 025个训练样本构成全局字典的训练集。全局字典训练参数稀疏度设置为3，迭代20次。字典大小为64×256，表示总共256个原子，每个原子代表了8×8大小的图像块。原图与全局字典如图 1所示。不同方法的降噪结果如图 2所示。
图1
模拟实验全局字典
Figure1.
The globe dictionary of simulation experiments
图2
模拟实验降噪结果(1 000 000光子)
Figure2.
The denoising results of simulation expreriments (with 1 000 000 photons)
如图 2所示，小波软阈值在滤除噪声的同时，使得胸腔边界模糊并且产生了伪影，TV算法对于重建产生的条纹噪声不能很好地抑制，而全局字典获得了相对平滑的效果。从局部放大图中可以看出，胸腔的边缘以及肺叶部分，小波软阈值处理结果模糊不清，TV算法处理结果仍然保留了许多放射条纹噪声，同时也导致了边缘一定程度的模糊，全局字典在去除噪声的同时很好地保留了边缘细节信息。不同量级噪声的处理结果如表 1所示。
表1
模拟实验结果
Table1.
Results of simulation experiments
由模拟数据的实验结果可以看出，全局字典降噪的效果，不仅在噪声滤除方面，在结构性保持方面也优于目前较为流行的一些降噪算法。
2.2 动物实验结果
本节中所有小动物的实际CT图像数据由苏州海斯菲德信息科技有限公司Hiscan M-1000型CT拍摄小鼠获得。图像大小为884×884，像素的空间分辨率为50 μm。CT扫描管电压40 kVp，由于CT的X射线剂量与管电流成线性正比关系，实验过程中高剂量的CT图像管电流为200 μA，低剂量CT图像的管电流依次为100 μA、50 μA以及20 μA。CT图像的数据位数为12位，显示窗宽为900，窗位为600。
全局字典参数设置与模拟实验一致，选取高剂量扫描图序列中结构较为丰富的一幅图像作为训练样本，同样以8×8大小的图像块对高剂量图像进行采样，构成包含777 899个训练样本的训练集；字典大小为64×256，训练时稀疏度为3，迭代20次。训练样本及全局字典如图 3所示。需要特别说明的是，在实际实验中，由于需要计算PSNR以及SSIM值对实验结果进行分析，这要求不同剂量组别的CT图像结构一致。小老鼠活体CT实验由于小鼠呼吸及运动等因素难以做到这一点，考虑CT图像反映的是不同组织对X射线的吸收率，组织凋亡与否对吸收率影响并不明显，因此包含数值计算验证实验可行性的小鼠实验是利用刚死亡的小鼠进行拍摄的。在一次实验过程中，对小鼠尸体分别进行剂量从高到低的连续四组拍摄，每组重建20层，对低剂量的三组CT图像利用全局字典、小波软阈值以及TV算法进行处理并与高剂量组进行对比；总共重复3次实验，并观察实验结果。实验结果如图 4及表 2所示。其中图 4展示了一次实验中一层图像的处理结果；表 2展示为三次试验不同剂量组每组各60幅图像的PSNR与SSIM数值的平均值。验证方法的可行性后，再对活体小鼠进行降噪方法的实验对比，保证实验的完整性。
图3
动物实验全局字典
Figure3.
The global dictionary of animal experiments
图4
动物实验结果(100 μA)
Figure4.
Results of animal experiments (100 μA)
表2
动物实验结果(平均)
Table2.
Results of animal experiments (averaged)
从实验结果中可以看出，小波软阈值法并不能很好地将噪声滤除，因此PSNR提升并不明显；TV降噪算法能够有效地将噪声滤除，并较为明显地提高PSNR，但在滤除噪声的过程中，图像的边缘变得平滑，部分有效信息被滤除，不能很好地保持图像的结构信息，因此在指标SSIM上的表现结果不尽如人意。全局字典降噪的效果相较于对比算法，在提高相同水平的PSNR的前提下，能够更好地保持结构信息，在PSRN以及SSIM两项衡量标准下都取得了较为满意的结果。如图 4所示的局部放大图中可以发现，原图噪声中微小的脏器间隙(圈内区域)以及骨骼内部的细小结构(方框内)在小波软阈值与TV算法的降噪过程中，细节信息被一并滤除了，全局字典算法则很好地将间隙及骨小梁结构保留了下来。活体小鼠实验的细节对比结果如图 5所示，可以看出全局字典算法在结构保持、边缘保留上均具有较好的结果，也说明了全局字典降噪算法在活体实验中的可行性。
图5
动物实验结果(活体)
Figure5.
Results of animal experiments(in vivo)
3 结语
本文提出了一种基于K-SVD字典学习的全局字典小动物CT图像降噪方法。通过对高剂量小动物CT图像进行训练获得全局字典，并将其用于低剂量小动物CT图像的稀疏分解中，最后利用稀疏分解获得的稀疏编码对低剂量图像进行重建。实验结果表明，与传统的主流算法相比，本文提出的算法在降噪的同时，能够较好地保留小动物CT图像中的结构信息，较为显著地提升低剂量小动物CT图像质量，有利于动物模型实验的顺利进行。
需要指出的是，在低剂量情况下，小动物CT图像不仅会被噪声干扰，也会由于X射线光子不足而产生放射状伪影，如何利用全局字典算法将伪影去除，进一步提升低剂量小动物CT图像质量，将是课题组今后的研究方向。
引言
动物模型是现代生物医学研究中的重要组成部分。通过对动物接种肿瘤观察其生理变化、利用各类新型造影剂进行动物实验、研究药物对动物生理的影响等[1]，能够有效地帮助人类认识疾病发生规律，并采取相应的防治措施。近年来，各类医学影像技术广泛地应用于各类动物实验中，其中计算机断层扫描(computed tomography，CT)技术由于具有分辨率高、成像时间短等特点，成为动物医学影像研究的重要手段之一。
与传统的临床CT相比，小动物CT具有更高的分辨率，但却受到射线剂量方面的限制。小动物本身十分脆弱，高剂量的X射线辐射往往会对诸如小鼠之类的动物造成损伤，甚至死亡，不利于实验的连续性，因此低剂量小动物CT将是未来研究的一个发展趋势。但是低剂量的CT重建图像会被噪声干扰，影响图像质量，因而如何提高低剂量小动物CT的重建图像质量，成为了亟待解决的问题。
传统的低剂量CT重建问题的解决方案主要分为三类，分别应用于重建过程的三个部分：重建前处理(对CT投影图或者Sinogram图进行降噪处理)、重建(对CT重建过程中的算法进行改进以抑制噪声)以及重建后处理(对CT重建结果进行降噪处理)。在前处理过程中，由于投影图是通过X射线直接扫描得到的，因而噪声类型较好统计，如基于惩罚加权的最小二乘(penalized weighted least-squares，PWLS)算法通过统计方法获得投影图像噪声模型，针对噪声模型进行降噪处理[2]；此外也有对Sinogram图像进行噪声统计，并利用非线性滤波进行降噪处理的办法[3]。由于此类方法应用于重建算法之前，有可能会造成图像信息的丢失或者引发图像的畸变。另外，由于投影图的数据量非常庞大，统计过程耗时较长。而图像重建过程中，对重建算法进行改进的一种基本思路是通过在传统的滤波反投影过程中添加平滑窗函数来获取较为平滑的结果，比如信号处理领域中常见的汉宁窗、海明窗等。此外，也有采用基于统计迭代图像重建的算法，这类算法通过将CT重建图像数据与投影数据相关联，对投影数据的噪声模型建模，在迭代过程中引入权重，达到图像降噪的目的[4-5]。虽然迭代类算法对噪声的抑制作用比较明显，但其计算量庞大、耗时长，在临床上应用并不广泛。在重建图像的后处理中，一般将重建图像视为普通图像，采用传统的线性或非线性滤波方法，例如利用核函数的办法对图像进行降噪滤波，以及利用总变分最小化(total variation，TV)的办法对图像进行降噪等[6-7]；也有学者提出将图像变换到小波域进行降噪处理[8]，但是由于重建图像需要通过反投影等复杂计算才能得到，噪声模型不便统计，传统的图像处理办法很难取得一个好的结果。
近阶段兴起的信号稀疏表示与字典学习等方法为上述问题提供了新的解决方案。字典的概念最早由1993年Mallat等[9]提出。1998年Chen等[10]提出利用一范数评估稀疏性以及利用凸函数求解系数问题，推进了该理论的发展。2001年，Donoho等[11]验证了追踪算法的确定性，使得稀疏表示理论取得长足发展，并广泛应用于图像的降噪、压缩与恢复等各个领域。在图像降噪方面，基本的算法思路可以概括为将图像利用字典稀疏的表达，在给定噪声误差的前提下，利用稀疏化的参数和字典将图像稀疏重构，达到滤除噪声的目的。字典解决方案直接应用于重建图像，保证了图像结果的真实性，也省去了前处理庞大的计算量，十分具有应用前景。
本文利用字典降噪的基本思路，使用高剂量图像作为训练样本，利用逐列更新的字典学习算法(K-means singular value decomposition, K-SVD)获取合适的全局字典，并将其应用于低剂量小动物CT的图像降噪研究之中。
1 算法
在小鼠低剂量CT图像降噪过程中，主要会涉及全局字典的训练学习、图像降噪算法实现以及小鼠CT应用框架搭建等具体操作问题。
1.1 K-SVD字典训练
利用K-SVD算法进行字典训练的过程可以归结于求解如下目标函数[12]
| $ \mathop {\min }\limits_{D, \alpha } \left\{ {\left\| {\boldsymbol{Y}-\boldsymbol{D}\alpha } \right\|_F^2} \right\}, st.{\forall _i}, \left\| {{\alpha _i}} \right\| \le T $ |
其中，Y是用于训练的样本集，Y的每一个列向量都为一个训练样本；D是初始字典，一般为大小为n×K的一个矩阵，K表示字典总共包含K个原子，n表示样本的维度。常用的初始字典有离散余弦字典或者小波字典；α是对样本集Y进行稀疏分解后获得的稀疏编码系数，T为式(1)的稀疏度约束，表明表达式中每一次稀疏表示的过程使用的字典原子数量不能超过T个。由于K-SVD算法是对字典D逐个原子进行更新，稀疏编码计算也是针对单个训练样本进行计算，因此我们将式(1)细化到每一个信号，重新改写为：
| $ \mathop {\min }\limits_{{\alpha _I}} \left\{ {\left\| {{\boldsymbol{y}_i}-\boldsymbol{D}{\alpha _i}} \right\|_2^2} \right\}, st.\left\| {{\alpha _i}} \right\| \le T $ |
其中i代表样本序号，K-SVD算法的目的在于使得式(2)取得最小值，联合两式可以将目标函数继续修改进行求解：
| $ \left\| {\boldsymbol{Y}-\boldsymbol{D}\alpha } \right\|_F^2=\left\| {\boldsymbol{Y}-\sum\limits_{j=1}^K {{\boldsymbol{d}_j}\alpha _T^j} } \right\|_2^2=\left\| {\left({\boldsymbol{Y}-\sum\limits_{j \ne k} {{\boldsymbol{d}_j}\alpha _T^j} } \right)-{\boldsymbol{d}_k}\alpha _T^k} \right\|_F^2 $ |
其中dj表示第j个字典原子，为一个列向量。令Ek表示，则Ek代表了在字典中去除第k个原子后式(3)所产生的误差，对Ek进行奇异值分解(singular value decomposition，SVD)，求解最大的满秩列向量，并利用该向量对原字典中第k个原子进行更新；将字典中所有原子遍历一次，便可得到迭代一次后更新的字典，为了使字典更为平滑、有效地抑制噪声信息，一般需要迭代多次。
1.2 图像降噪框架
通过K-SVD算法获取合适的全局字典后，则可以利用全局字典对噪声图像进行降噪处理[13]，利用全局字典进行降噪的问题可以表述为：
| $ \left\{ {{{\hat \alpha }_{ij}},\boldsymbol{\hat X}} \right\} = \mathop {{\rm{argmin}}\lambda }\limits_{{\alpha _{ij}},X} \left\| {\boldsymbol{X} - \boldsymbol{Y}} \right\|_2^2 + \sum\limits_{ij} {{\mu _{ij}}} {\left\| {{\alpha _{ij}}} \right\|_0} + \sum\limits_{ij} {\left\| {{\boldsymbol{R}_{ij}}\boldsymbol{X} - \boldsymbol{D}{\alpha _{ij}}} \right\|_2^2} $ |
其中，式(4)右侧第一项衡量的是原始图像X与噪声图像Y之间的整体相似程度，λ为拉格朗日系数；第二项为字典降噪过程中的稀疏度约束，其中μij为约束项的权重系数，由于整图较大，所以需要将原图划分为一块块的小图像并重排为与全局字典的原子维度一致的列向量；第三项中的RijX表示了原图中的第ij幅子图像，Rij是用于提取子图的矩阵，ij为子图像序列下标，Dαij为对应子图像的稀疏表达，第三项用于衡量这两者之间的误差，它们的误差越小越好。
由于已经通过K-SVD算法获取到了合适的全局字典D，为了获取每张子图的稀疏表达系数，需要求解下式最优化问题：
| $ {\alpha _{ij}}=\mathop {{\mathop{\rm a}\nolimits} \hat rgmin{\mu _{ij}}}\limits_\alpha {\left\| \alpha \right\|_0} + \sum\limits_{ij} {\left\| {\boldsymbol{D}{\alpha _{ij}}-{\boldsymbol{R}_{ij}}\boldsymbol{X}} \right\|_2^2} $ |
式(5)一般通过贪婪算法来求解。本文采用正交匹配追踪(orthogonal matching pursuit，OMP)来求解式(5)的问题[14]，其基本思路是对每一次逼近的残差进行新一轮的逼近，直到满足设定的误差限或者稀疏度要求，得到近似解âij。如此，式(5)的求解更新即可转换为如下形式：
| $ \boldsymbol{\hat X}=\mathop {{\mathop{\rm argmin}\nolimits} \lambda }\limits_\boldsymbol{X} \left\| {\boldsymbol{X}-\boldsymbol{Y}} \right\|_2^2 + \sum\limits_{ij} {\left\| {\boldsymbol{D}{{\hat \alpha }_{ij}}-{\boldsymbol{R}_{ij}}\boldsymbol{X}} \right\|_2^2} $ |
这是一个简单的凸函数问题，对其进行求解得到降噪后的图像：
| $ \boldsymbol{\hat X}={\left({\lambda \boldsymbol{I} + \sum\limits_{ij} {\boldsymbol{R}_{ij}^T{\boldsymbol{R}_{ij}}} } \right)^{-1}}\left({\lambda \boldsymbol{Y} + \sum\limits_{ij} {\boldsymbol{R}_{ij}^T\boldsymbol{D}{{\hat \alpha }_{ij}}} } \right) $ |
其中，I为单位矩阵，为降噪后的结果图像。
1.3 基于全局字典降噪的小动物CT低剂量图像降噪框架
在Aharon等[12]的工作中，利用噪声图像自身作为训练集，经过K-SVD算法训练过后的字典在降噪方面比初始字典更有优势。这主要是因为经由噪声图像训练自身得到的字典原子包含了图像的结构信息，这使得图像在做稀疏表达时能够获得更好的稀疏度，从而取得更好的结果。但是在实际应用中，如果对每一幅噪声图像都进行字典训练，无疑大大增加了计算量，并不可取。考虑到同种类的小动物CT图像结构相对比较接近，同种类动物的不同个体在CT图像上的结构差异更多地体现于结构整体的不同，比如组织的轮廓或者位置等，但是其组织及其骨骼等构成成分都类似，因而在图像细节方面表现大多类似，比如组织的边缘、灰度级、骨骼内部细微结构等，因此通过提前对高剂量的同种类CT图像进行字典训练获取包含结构信息的全局字典，对同种类小动物的低剂量噪声图像进行降噪处理，既能达到噪声图像自身训练字典的降噪效果，也能省去单独每幅图像全局字典训练的时间，使得全局字典在低剂量小动物CT上的应用具有了可行性。
因此，基于全局字典降噪的低剂量小动物CT降噪方案可以归纳为如下几个步骤：
(1) 获取训练样本集：利用高剂量X射线拍摄小动物获得质量好的小动物的CT图像；
(2) 获取全局字典：对高剂量小动物的CT图像进行处理，利用K-SVD算法获取合适的全局字典D；
(3) 低剂量图像稀疏分解：将低剂量图像分割为图像子集，利用D以及OMP算法对低剂量图像进行稀疏分解，获取每一个子集的稀疏编码α；
(4) 低剂量图像重构：利用式(7)以及步骤(3)中获得结果对低剂量图像进行重构，获取降噪后的纯净图像。
2 实验与结果分析
本文实验主要分为两个部分，模拟实验以及小动物CT数据实验。为了验证全局字典降噪的有效性，本文采用小波软阈值、TV算法降噪进行对比实验。实验结果通过峰值信噪比(peak signal-noise ratio，PSNR)以及结构相似度(structure similarity index measurement，SSIM)两个评价指标进行衡量[15]。PSNR以及SSIM都是图像客观评价的标准，其中PSNR反映的是最大值信号与噪声之间的关系，单位为dB，PSNR值越大，噪声滤除效果越好；而SSIM则综合考量两幅图像之间亮度、对比度以及结构相似度作为整体评判指标，SSIM值越接近1，代表两幅图像越相似。
2.1 模拟数据实验结果
为了尽可能地接近小动物CT实验的真实情况，模拟实验采用人体胸腔CT图像作为模拟数据(数据来源于文献[16])，图像大小为512×512。CT图像噪声是由于到达探测器的光子数量不同(被扫描物体各部分对X光吸收不一致)并且有暗电流背景噪声等因素影响而形成。为了使模拟实验更加接近真实情况，对CT重建图像的噪声产生过程进行模拟，首先图像进行正投影后，投影图添加泊松噪声再进行重建[17]，重建后的CT图像利用全局字典、小波软阈值以及TV算法进行降噪处理，并计算处理后的结果与无噪声标准重建CT图像之间的PSNR值与SSIM值，观察实验结果。
其中，全局字典由无噪声标准重建图像训练而来，以8×8图像块大小对原图像素点进行采样，总共获取255 025个训练样本构成全局字典的训练集。全局字典训练参数稀疏度设置为3，迭代20次。字典大小为64×256，表示总共256个原子，每个原子代表了8×8大小的图像块。原图与全局字典如图 1所示。不同方法的降噪结果如图 2所示。
图1
模拟实验全局字典
Figure1.
The globe dictionary of simulation experiments
图2
模拟实验降噪结果(1 000 000光子)
Figure2.
The denoising results of simulation expreriments (with 1 000 000 photons)
如图 2所示，小波软阈值在滤除噪声的同时，使得胸腔边界模糊并且产生了伪影，TV算法对于重建产生的条纹噪声不能很好地抑制，而全局字典获得了相对平滑的效果。从局部放大图中可以看出，胸腔的边缘以及肺叶部分，小波软阈值处理结果模糊不清，TV算法处理结果仍然保留了许多放射条纹噪声，同时也导致了边缘一定程度的模糊，全局字典在去除噪声的同时很好地保留了边缘细节信息。不同量级噪声的处理结果如表 1所示。
表1
模拟实验结果
Table1.
Results of simulation experiments
由模拟数据的实验结果可以看出，全局字典降噪的效果，不仅在噪声滤除方面，在结构性保持方面也优于目前较为流行的一些降噪算法。
2.2 动物实验结果
本节中所有小动物的实际CT图像数据由苏州海斯菲德信息科技有限公司Hiscan M-1000型CT拍摄小鼠获得。图像大小为884×884，像素的空间分辨率为50 μm。CT扫描管电压40 kVp，由于CT的X射线剂量与管电流成线性正比关系，实验过程中高剂量的CT图像管电流为200 μA，低剂量CT图像的管电流依次为100 μA、50 μA以及20 μA。CT图像的数据位数为12位，显示窗宽为900，窗位为600。
全局字典参数设置与模拟实验一致，选取高剂量扫描图序列中结构较为丰富的一幅图像作为训练样本，同样以8×8大小的图像块对高剂量图像进行采样，构成包含777 899个训练样本的训练集；字典大小为64×256，训练时稀疏度为3，迭代20次。训练样本及全局字典如图 3所示。需要特别说明的是，在实际实验中，由于需要计算PSNR以及SSIM值对实验结果进行分析，这要求不同剂量组别的CT图像结构一致。小老鼠活体CT实验由于小鼠呼吸及运动等因素难以做到这一点，考虑CT图像反映的是不同组织对X射线的吸收率，组织凋亡与否对吸收率影响并不明显，因此包含数值计算验证实验可行性的小鼠实验是利用刚死亡的小鼠进行拍摄的。在一次实验过程中，对小鼠尸体分别进行剂量从高到低的连续四组拍摄，每组重建20层，对低剂量的三组CT图像利用全局字典、小波软阈值以及TV算法进行处理并与高剂量组进行对比；总共重复3次实验，并观察实验结果。实验结果如图 4及表 2所示。其中图 4展示了一次实验中一层图像的处理结果；表 2展示为三次试验不同剂量组每组各60幅图像的PSNR与SSIM数值的平均值。验证方法的可行性后，再对活体小鼠进行降噪方法的实验对比，保证实验的完整性。
图3
动物实验全局字典
Figure3.
The global dictionary of animal experiments
图4
动物实验结果(100 μA)
Figure4.
Results of animal experiments (100 μA)
表2
动物实验结果(平均)
Table2.
Results of animal experiments (averaged)
从实验结果中可以看出，小波软阈值法并不能很好地将噪声滤除，因此PSNR提升并不明显；TV降噪算法能够有效地将噪声滤除，并较为明显地提高PSNR，但在滤除噪声的过程中，图像的边缘变得平滑，部分有效信息被滤除，不能很好地保持图像的结构信息，因此在指标SSIM上的表现结果不尽如人意。全局字典降噪的效果相较于对比算法，在提高相同水平的PSNR的前提下，能够更好地保持结构信息，在PSRN以及SSIM两项衡量标准下都取得了较为满意的结果。如图 4所示的局部放大图中可以发现，原图噪声中微小的脏器间隙(圈内区域)以及骨骼内部的细小结构(方框内)在小波软阈值与TV算法的降噪过程中，细节信息被一并滤除了，全局字典算法则很好地将间隙及骨小梁结构保留了下来。活体小鼠实验的细节对比结果如图 5所示，可以看出全局字典算法在结构保持、边缘保留上均具有较好的结果，也说明了全局字典降噪算法在活体实验中的可行性。
图5
动物实验结果(活体)
Figure5.
Results of animal experiments(in vivo)
3 结语
本文提出了一种基于K-SVD字典学习的全局字典小动物CT图像降噪方法。通过对高剂量小动物CT图像进行训练获得全局字典，并将其用于低剂量小动物CT图像的稀疏分解中，最后利用稀疏分解获得的稀疏编码对低剂量图像进行重建。实验结果表明，与传统的主流算法相比，本文提出的算法在降噪的同时，能够较好地保留小动物CT图像中的结构信息，较为显著地提升低剂量小动物CT图像质量，有利于动物模型实验的顺利进行。
需要指出的是，在低剂量情况下，小动物CT图像不仅会被噪声干扰，也会由于X射线光子不足而产生放射状伪影，如何利用全局字典算法将伪影去除，进一步提升低剂量小动物CT图像质量，将是课题组今后的研究方向。
