基于NSCT稀疏编码的医学图像超分辨率重建算法

Super-Resolution Reconstruction of Medical Images Based on Sparse Representation in the NSCT Domain

摘要: 为了满足临床对医学图像分辨率的进一步需求，提出了一种非下采样Contourlet变换和稀疏表示结合的单帧医学图像超分辨率重建算法。该方法在改进字典学习步骤中，首先通过非下采样Contourlet变换中的方向滤波器提取训练样本集图像的高频特征；其次，采用K-SVD奇异值分解法和OMP正交匹配追踪算法来求取高、低分辨率字典和稀疏表示系数，并重建出高分辨率图像。最后，将提出的方法与Yang、Zeyde提出的稀疏表示重建算法都应用于医学图像超分辨率重建上，通过实验对比分析，提出算法重建结果的PSNR值与SSIM值分别达到45.389和0.907，充分展现出了提出算法的优势。

Abstract: In order to meet the additional demands of image resolution in clinical medicine, a super-resolution reconstruction algorithm based on sparse representation in the nonsubsampled contourlet transform (NSCT) domain is proposed. In this method, in the improved dictionary learning step, firstly, the high-frequency features of the training sample set images are extracted by the direction filter in the non-downsampling Contourlet transform; secondly, the K-SVD singular value decomposition method and the OMP orthogonal matching tracking algorithm are used to find the high- and low-resolution dictionaries and sparse representation coefficients, and reconstruct the high-resolution images. Finally, the proposed method and the sparse representation reconstruction algorithm proposed by Yang and Zeyde are both applied to the super-resolution reconstruction of medical images. Through experimental comparison and analysis, the PSNR and SSIM values of the reconstruction results of the proposed algorithm reach 45.389 and 0.907 respectively, fully demonstrating the advantages of the proposed algorithm.

1. 引言

随着医学图像的发展，对医学图像分辨率的研究具有重要的现实意义。然而，医学图像的成像过程受到医疗设备成像分辨率、图像采集时间、患者能承受的辐射剂量等因素的影响，获得图像的分辨率不高。目前，提高图像分辨率的方法包括增加传感器的像素数、增加传感器的感光区域、改变成像设备的扫描方式以及超分辨率重建技术。通过硬件的改进来提高图像的分辨率会带来辐射量大、医疗设备的尺寸大和成本高等问题，因此超分辨率重建算法成为首选。该方法通过算法提高了物理设备获取图像的分辨率，避免了硬件开发的限制，提高了图像质量，已广泛应用于军事、医学、公安、计算机视觉等领域。

超分辨率重建是指以低分辨率图像作为研究对象，使用算法来提高其分辨率，得到一幅高分辨率图像的过程。目前，在医学领域的图像超分辨率重建算法已经发展到一定程度，根据用于重建高分辨率图像的低分辨率图像数量，可分为基于单帧的超分辨率重建算法和基于序列的超分辨率重建算法；根据算法类型，分为传统类型算法和新型算法。

传统类型的超分辨率重建算法包括插值法、频域变换法、迭代反投影法(Iterative Back Projection, IBP) [1] 、凸集投影法(Projection onto Convex Sets, POCS) [2] 、最大后验概率法(Maximum A Posteriori, MAP) [3] 以及POCS-MAP混合的方法等等。插值法将周围像素代入插值函数中得到的值即为待插点的像素值，以此来增加图像的像素点，常用的插值函数有最近邻插值、双线性插值、三次样条插值和局部自适应插值等 [4] ，如Mahmoudzadeh等 [4] 就通过对脑部磁共振(Magnetic Resonance, MR)图像进行插值技术重建合成了三维脑部MR图像。频域变换法即在频域内对图像信息进行间接处理，即将对图像空间域中像素的操作转换为对变换域中频域系数的处理。常用的频域重建方法有离散余弦变换 [5] 、小波变换 [6] 、Contourlet变换 [7] 、非下采样Contourlet变换(Non-Subsampled Contourlet Transform, NSCT) [8] 等。Pak等 [9] 利用NSCT具有平移不变性、各向异性和多方向性的特点，对乳腺图像进行超分辨率重建，提高了图像质量。Zheng等 [10] 使用IBP方法改善了脑部MR图像的边缘；Xu等 [11] 提出的基于先验知识的改进POCS算法较好地恢复原始图像的高频；Jonghye提出了一种基于舌头MRI序列图像的重建方法，结合Markov随机场，将问题转化为最大后验概率的最优化问题 [12] 。

总体来看，传统类型的超分辨率重建算法存在一些不足。插值法针对单帧图像，虽然步骤简单，但具有受噪声影响较大、重建图像具有人工伪影等缺点，大多应用于对图像的分辨率要求不是特别高的场合，更广泛的作用是作为一种辅助算法应用到其他的新颖算法中。由于缺乏先验知识，频域方法只能处理序列图像间的全局或整体运动，对其他情况带来的图像分辨率问题解决效果不佳。IBP算法的解集不唯一、POCS依赖于初始估计图像，且收敛性慢 [13] 、MAP算法依赖先验信息和计算量大。传统类型的超分辨率重建算法局限于图像自身的关系进行处理，低分辨率图像本身包含的图像信息就很有限，因此，新型的超分辨率重建算法得到发展，该类算法能够突破低分辨率图像信息有限的局限性，充分利用高、低分辨率图像之间的对应关系。

新型超分辨率重建算法基于学习的理论，通过学习高、低图像样本之间的映射关系并将其引入到单帧低分辨率图像的像素重建中。最早的学习类方法是邻域嵌入法(Neighbor Embedding, NE) [14] ，除此之外，还有卷积神经网络法、稀疏字典学习法。

NE方法的基本思想是：针对待重建低分辨率图像块，寻找K个最相似样本块以及对应高分辨率样本块，然后，将这K个低分辨率样本块线性组合成当前低分辨率块，并使得重建误差最小，再将这种线性关系用于重建对应的高分辨率样本块，最后合成高分辨率图像。但该方法依赖于图像样本，且图像样本较少，难以保证寻求到的图像块是最佳的。

卷积神经网络(Super-Resolution Using Convolutional Neural Network, SRCNN)的超分辨率重建算法首先利用Cubic插值算法将图像放大，然后将训练好的卷积神经网络用于重建超分辨率图像上。Li等 [15] 针对一些图像存在较大噪声的问题，提出了一种基于抗噪声的SRCNN重建算法。Oktay等 [16] 提出了一种通用的训练策略，将新的正则化模型先验知识整合到神经网络中。

Yang等 [17] 提出了基于稀疏表达的方法，该方法源于信号的稀疏性求解得到广泛的关注，信号可由字典和稀疏表示系数重建得到，且在确定字典下，求得的系数越稀疏，重建得到的信号质量越好，即从字典中选择尽量少还能尽可能逼近原始信号结构的字典元素。稀疏字典学习法应用到图像重建中，包括字典训练和稀疏重建两个重要步骤。假设输入的低分辨率图像块可以被训练样本集稀疏线性表达，字典训练的过程中通过训练和学习大量与目标相似的图像样本来获得高、低分辨率字典，采用联合字典训练法，字典下的高、低分辨率图像块会具有相同的稀疏表示系数，因此在重建时，低分辨率字典下的稀疏系数结合高分辨率字典就可得到高分辨率重建图像。

稀疏字典重建法求解通常分为两步：① 字典学习算法有：最优方向法(Method of Optimal Direction, MOD) [18] 、快速奇异值分解法(K-singular Value Decomposition, K-SVD) [19] 、特征标志搜索方法(Feature Sign Search, FSS) [20] 、在线字典更新方法(Online Dictionary Learning, ODL) [21] 、主成分分析法等。② 稀疏分解优化算法求解在这个基字典下的稀疏表示系数，常用的稀疏表示系数求解算法又分为两大类，分别为贪婪算法、凸松弛算法。贪婪算法又可细分为基追踪法、匹配追踪法、正交匹配追踪法(Orthogonal Matching Pursuit, OMP)算法等。

其中，字典学习步骤中K-SVD算法明显提高了字典求解速度。在稀疏字典学习的应用上，Ota等 [22] 研究了基于稀疏字典学习的稀疏编码超分辨率重建方法，极大地提高了重建后电子计算机断层扫描(Computer Tomography, CT)图像质量。Asif等 [23] 也利用联合学习字典从单个低分辨率图像中提取高分辨率图像，提高了含噪声CT图像的质量。为了规避Yang的算法效率较低的特点，Zeyde等人 [24] 于2010年提出了使用K-奇异值分解(K-Singular Value Decomposition, K-SVD)算法训练字典的方法。

在新型超分辨率重建算法中，NE算法不能保证低分辨率图像块和寻找到的样本块具有较好的相似性；SRCNN算法重建效果好，但依赖于模型参数，且计算量大；Yang提出的稀疏表示重建算法充分利用图像的内容和结构，在重建过程中对大量的外部图像库进行训练，引入了大量的图像先验知识，是目前超分辨率重建领域较好的重建算法，但也存在重建时间长、算法效率低的缺点。

针对新型算法中Yang提出的稀疏字典学习算法存在的不足以及Zeyde提出的改进算法，本文提出了一种基于NSCT稀疏编码的医学图像超分辨率重建算法。该算法通过NSCT中的方向滤波器提取图像的高频信息，代替Zeyde方法中一、二阶梯度特征，得到一个表示图像块的特征向量，与原始方法相比，NSCT提取图像特征能够有效保留图像各个方向的高频信息；然后在稀疏编码阶段采用K-SVD奇异值分解法求解低分辨率字典和高分辨率字典，在重建阶段采用OMP算法计算稀疏表示系数；最后重建出高分辨率的医学图像。NSCT提取图像特征加上K-SVD求解稀疏字典的方法结合，提高了字典的有效性，使得重建的图像效果更好。

2. 理论

2.1. 图像的稀疏表示重建

稀疏编码的原理是：用较少的基本信号的线性组合去尽可能地表达原始信号。这些基本信号从字典中选择，进行线性表示的系数即为稀疏系数。通俗地说，字典中存在众多信号，有些与原始信号相关性较大，稀疏编码的目的就是用尽可能少的字典中信号去描述原始信号。将图像看作一种特殊的信号，图像的稀疏表示就是利用字典中尽量少的图像原子的线性组合去表示原图像信号。

在Yang等人提出的稀疏表示超分重建算法中，有如下假设：本质上图像是具有相似结构表示基元的，并且结构是高频信息，例如图像边界信息、纹理结构信息等。稀疏编码应用于图像超分辨率重建中，可以简单描述为：字典由多个图像特征原子组成，图像的稀疏字典重建即从字典中选择尽量少的能近似表现原始图像特征的图像特征原子。假设图像的特征信号 $y\in R^N$ 能够用字典 $D\in R^{N\times M}$ 中较少的原子来进行线性表示，模型表示为图1和公式(1) [25] ：

Figure 1. Sparse represents the model

图1. 稀疏表示模型

$y={\displaystyle {\sum }_{i=1}^L{d}_i{\alpha }_i}=D\alpha$ (1)

其中，y表示行数为N，列数为1的原始图像信号； $D=\left[d_1,d_2,\cdots ,d_M\right]$ 是大小为N*M的字典，D的每一列 $d_i$ 表示第i个图像特征； $\alpha ={\left[{\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_M\right]}^{\text{T}}$ 是行数为M的稀疏表示系数，每个系数对应的是一个值，且大部分系数的值都是0。此时，原始图像信号y可以被稀疏表示。

稀疏表示系数求解方式设定如公式(2)。

$\arg {\min }_{D,\alpha }{\Vert y-D\alpha \Vert }_2^2,\text{s}\text{.t}.{\Vert \alpha \Vert }_0<T$ (2)

即在 $\alpha$ 中非0系数个数不超过T的约束条件下，优化公式(2)中的目标方程，得到一个最优字典D和一个稀疏表示系数矩阵 $\alpha$ 。

稀疏表示用于图像超分辨率重建中，需要解决两个问题：一是如何依据给定的字典找出稀疏性最好的系数；二是如何利用样本图像集训练出最有用的字典。假设X表示原始高分辨率图像，x则表示对X进行分块后的高分辨率图像块，Y表示待重建的低分辨率图像，y表示低分辨率图像块， $\alpha$ 是稀疏系数，D_H和D_L分别表示高、低分辨率字典。

高分辨率图像块x可由高分辨率字典D_H和稀疏表示系数 ${\alpha }^{\prime }$ 表示，如公式(3)。

$x=D_H\ast {\alpha }^{\prime }{\Vert {\alpha }^{\prime }\Vert }_0\ll k$ (3)

其中k是稀疏表示系数 ${\alpha }^{\prime }$ 的行数。在稀疏表示系数中，值为非0的系数个数要远远小于k。

同样地，低分辨率图像块y可通过稀疏表示求出，对应的稀疏表示系数用 ${\alpha }^{″}$ 表示，如公式(4)。

$y=D_L\ast {\alpha }^{″}$ (4)

通过对对高分辨率图像块特征和低分辨率图像块特征进行联合字典训练，迭代优化得到具有相同稀疏表示系数的联合字典，即 $\alpha ={\alpha }^{\prime }={\alpha }^{″}$ 。此时，对于待重建的低分辨率图像特征，只要求得对应低分辨率字典D_L下的稀疏系数，再联合高分辨率字典D_H，便可重建出高分辨率图像特征。重建的高分辨率图像特征块 $\tilde{x}$ 可用公式(5)表示。

$\tilde{x}=D_H\ast \alpha$ (5)

2.2. 基于NSCT稀疏编码的医学图像超分辨率重建

基于NSCT特征提取的稀疏字典超分辨率重建算法的具体步骤包括NSCT特征提取、联合字典训练和图像稀疏表示重建三大步骤。

2.2.1. 算法流程

对图像进行稀疏编码重建的方法包含字典训练、图像重建两个阶段。其中，图像特征提取采用NSCT变换方法；低分辨率图像的构建通过对高分辨率图像进行下采样得到；稀疏字典求解采用K-SVD算法；稀疏表示系数则通过OMP算法得到。具体的重建流程如图2所示。

图像重建流程(图2中算法左边流程)：

1) 输入低分辨率图像LR、加载联合字典训练得到的低分辨率图像块字典D_L和高分辨率图像块字典D_H；

2) 结合LR和低分辨率字典D_L，采用OMP正交匹配追踪算法，求解出对应高分辨率字典D_H下的稀疏表示系数 $\alpha$ ；

3) 根据步骤(2)中求解出的稀疏表示系数 $\alpha$ 、高分辨率字典D_H，重建得到高分辨率医学图像 $\tilde{X}$ ；

字典训练流程(图2中虚线框内流程)：

1) 读入训练样本集，对其下采样，形成对应的低分辨率图像样本集；

2) 对低分辨率图像采用NSCT变换，得到各个频段子带图像；

3) 对训练样本中的每个低分辨率图像进行分块，每个块对应的中高频子带上系数即构成低分辨率图像块的样本特征；高分辨率图像块的特征直接由图像像素值减去均值图像得到；

4) 采用K-SVD奇异值分解法对高分辨率图像块特征和低分辨率图像块特征进行联合字典训练，迭代优化得到具有相同稀疏表示系数的高、低分辨率联合字典D_H和D_L。

Figure 2. The flow of algorithm

图2. 算法流程

2.2.2. NSCT特征提取

NSCT的每级分解会先采用Nonsubsampled Pyramid Filter Banks (NSPFB)进行塔型分解，得到低频、中高频子带；接着，中高频子带会进一步通过方向滤波器组(Directional Filter Bank, DFB)变为各个方向子带。与Contourlet变换相比，NSCT去掉了下采样和插值过程，有效减少了图像缩放造成的信息丢失。

采用NSCT对一幅图像进行特征提取，将该图像分解为对应的多频带图像，包含了低频、中频和高频等几个部分。NSCT变换后分别得到 个低频子带和 ${\displaystyle {\sum }_{l=1}^N{2}^M}$ 个中高频子带。通常用 $f^{\left(0\right)}\left(i,j\right)$ 表示低频子带坐标 $\left(i,j\right)$ 处的系数，用 $f^{\left(l,k\right)}\left(i,j\right)$ 表示l尺度下的第k个中高频子带坐标 $\left(i,j\right)$ 处的系数。其中l表示NSP分解尺度，N表示最大NSP分解尺度；k表示NSDFB方向子带，M表示最大NSDFB方向分解尺度。为减少计算量，本文提出算法只采用一级分解，即分解尺度 $M=1$ ，方向分解尺度 $M=2$ ，得到一个低频系数矩阵和4个中高频系数矩阵。进行稀疏编码时，高频信息更具有典型性，因此，4个中高频的子带才是本文提出的特征提取方法处理的对象。所用方法的NSCT分解方式示意图如图3所示 [26] 。

Figure 3. Structural diagram of the NSCT decomposition

图3. NSCT分解结构图

本文算法对图像特征提取应用在图像样本训练集准备和低分辨率图像重建两个过程中。

在图像训练样本集准备阶段，首先对收集到的HR图像样本进行下采样得到低分辨率样本对。针对低分辨率图像，进行NSCT变换后，每张图像的4个中高频子带图像划分为k × k的图像块，则低分辨率图像块的特征直接用每个中高频子带图像表示。其中，第i幅低分辨率图像的中高频系数，即低分辨率图像块的特征F_L，用公式(6)表示。

$F_L=\left[F_i^{\left(1,1\right)},F_i^{\left(1,2\right)},F_i^{\left(1,3\right)},F_i^{\left(1,4\right)}\right]$ (6)

针对高分辨率图像，由于超分辨率的本质是要找回图像降质过程丢失的高频细节信息，只需要将图像块内每个像素减去块内所有像素值的平均值，即可得到高分辨率图像块特征F_H，如公式(7)所示。

$F_H=\left[HR-mean\left(HR\right)\right]$ (7)

综上所述，F_L和F_H构成样本训练集图像特征对。

在低分辨率图像重建过程中，对待重建的低分辨率图像进行NSCT变换后选取4个中高频子带图像，将其划分同样大小的k × k图像块。待重建图像的特征块用LR表示为公式(8)。

$LR=\left[L{R}^{\left(1,1\right)},L{R}^{\left(1,2\right)},L{R}^{\left(1,3\right)},L{R}^{\left(1,4\right)}\right]$ (8)

2.2.3. 联合字典训练

在获得训练样本集特征F_L和F_H的基础上，可选择合适的字典学习算法对其进行训练得到高、低分辨率字典。

在保证字典训练有效性的同时，K-SVD方法采用奇异值分解法逐个更新字典中的原子，减少了MOD方法中矩阵求逆的步骤，同时更新相关表示系数和字典基元的策略也降低了时间复杂度。基于K-SVD方法求解字典具有的上述优点，本文采用该算法求解字典。

设 $Y=\left[y_1,y_2,\cdots ,y_n\right]\in R^{N\times n}$ 为训练样本集构成的矩阵， $D=\left[d_1,d_2,\cdots ,d_m\right]\in R^{N\times m}$ 是稀疏字典， $\alpha =\left[{\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_n\right]\in R^{m\times n}$ 为稀疏表示系数。训练样本中的第i个图像块特征 $y_i$ 可用字典D中的特征原子线性表示，线性系数为 ${\alpha }_i$ ，求解的目标函数为公式(9)。

$\left\{\tilde{D},\tilde{\alpha }\right\}=\arg \underset{\left\{D,\alpha \right\}}{\min }\left\{{\Vert Y-D\alpha \Vert }_F^2+\mu {\Vert \alpha \Vert }_1\right\}$ (9)

其中， $\tilde{D},\tilde{\alpha }$ 为字典D和稀疏表示系数 $\alpha$ 的近似值。

结合2.2.2节中的图像特征块， $F_L=\left[y_{L1},y_{L2},\cdots ,y_{Ln}\right]$ 和 $F_H=\left[y_{H1},y_{H2},\cdots ,y_{Hn}\right]$ 分别对应为样本集中的n个高、低分辨率样本特征块，则高、低分辨率字典的求解模型如公式(10)和公式(11)。

$\left\{{\tilde{D}}_H,\tilde{\alpha }\right\}=\arg \underset{\left\{D_H,\alpha \right\}}{\min }\left\{{\Vert F_H-D_H\alpha \Vert }_F^2+\mu {\Vert \alpha \Vert }_1\right\}$ (10)

$\left\{{\tilde{D}}_L,\tilde{\alpha }\right\}=\arg \underset{\left\{D_L,\alpha \right\}}{\min }\left\{{\Vert F_L-D_L\alpha \Vert }_F^2+\mu {\Vert \alpha \Vert }_1\right\}$ (11)

K-SVD算法求解字典的简要步骤如下：

1) 先初始化字典D为D₀，并设置迭代次数k；

2) 针对样本特征F_L、F_H以及字典D₀，采用OMP算法计算稀疏表示系数 $\alpha$ ；

3) 使用K-SVD算法逐列去更新字典D，当更新某一列时，字典矩阵的其他列固定，并得到最新的稀疏表示系数；

4) 判断D的所有列是否都被更新过，若未更新完，则重新回到步骤(3)；

5) 判断是否达到迭代次数k，若迭代次数小于k，则回到步骤(2)；

6) 得到最新的字典D和稀疏表示系数 $\alpha$ 。

K-SVD算法求解字典的伪代码描述如下：

2.2.4. 图像稀疏表示重建

经过2.2.3节字典训练过程，得到高、低分辨率字典对D_H和D_L，该字典对是基于学习的重建方法中常提到的映射关系，将其用于指导低分辨率图像的超分辨率重建过程即可。在图像稀疏表示重建中，首先对输入的单帧低分辨率图像采用相同的特征提取方法，并分成相同大小的图像块；然后结合低分辨率字典D_L，使用OMP优化算法求出对应的稀疏表示系数；最后，将稀疏表示系数与高分辨率字典D_H联合，重建出高分辨率图像块，进一步组合成高分辨率图像。

重建步骤为：

1) 输入待重建的低分辨率医学图像Y，采用NSCT提取图像的特征，并将其划分为大小是 $k\ast k$ 的图像块y；

2) 结合低分辨率字典D_L，根据公式(2)采用OMP算法求解图像块对应的稀疏表示系数 $\alpha$ ，本文算法求解用更具体的公式(12)进行表示。

$\arg {\min }_{\alpha }{\Vert y-D_L\alpha \Vert }_2^2,\text{s}\text{.t}.{\Vert \alpha \Vert }_0<T_0$ (12)

3) 联合字典训练的方法保证了高、低分辨率字典具有相同的系数表示系数，故联合高分辨率字典D_H，便可重建得到高分辨率图像块 $\tilde{x}$ ，如公式(13)所示。

$\tilde{x}=D_H\ast \alpha$ (13)

步骤(2)中OMP算法求解的过程实际上就是对字典中原子(图像特征)向量进行施密特正交化的过程，正交化的目的在于已经选择过的原子不会再被重复选择。将图像特征减去在正交化后的原子上各自的分量即可得残差，并通过最小二乘法不断更新残差最后实现准确重建稀疏向量。算法求解步骤可更具体地展开描述 [27] ：

1) 用LR表示图像特征，初始化残差 $e_0=LR$ ，初始化迭代次数 $m=1$ ；

2) 从字典中找到与 $e_0$ 的内积绝对值最大的原子，表示为 ${\phi }_1$ ；

3) 将选择的原子作为列组成矩阵 ${\Phi }_{mt}$ ；

4) 求最小二乘解，更新残差：

$x=\underset{e_m}{\arg \min }\Vert y-{\Phi }_{mt}x\Vert$

$e_m=LR-{\Phi }_{mt}x$

其中， ${\Phi }_{mt}$ 为第m次迭代过程中所有被选择过的原子组成的矩阵，随着m值的递增迭代更新， $e_m$ 为第m次迭代计算得到的残差。

5) 当达到最大迭代次数或者误差到达限制范围时，停止迭代并输出稀疏表示系数x。

本文使用OMP算法对低分辨率图像求解稀疏表示系数的伪代码如下：

3. 实验结果

3.1. 实验设定

3.1.1. 实验数据集的建立

为了验证本文提出的超分辨率重建算法有效性，对医学图像的训练集和测试样本进行了选择。用于训练字典的图像集由100张牙科自然光图像、30张CT图像、35张MRI图像和60张X-ray图像组成；用于验证和统计超分辨率重建算法实现效果的测试图像包括牙科图像200张，CT图像30张，MRI图像30张，X-ray图像40张，共计300张医学测试图像。

3.1.2. 实验平台的建立

本文在Windows7操作系统下验证算法效果，实验平台为：64位Intel Core i5处理器、主机频率2.5GHZ、8G运行内存。

3.1.3. 算法评价体系

超分辨率重建图像的质量评价方法主要有主观评价和客观评价。主观评价是通过人来观察图像，对图像的优劣做出主观评定。这种方法评价结果与观察条件和观察者的状态等诸多因素有关。因此人们提出了多种客观评价方法，最常用的是逼真度测量，即测量被评价图像偏离参考图像的程度。均方误差(Mean Squared Error, MSE)和峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)等误差度量方法基于图像像素灰度值进行统计和平均计算，是常用的衡量信号失真的指标，如公式(14)和公式(15)。另一种衡量两幅图像相似度的指标，即为结构相似性(Structural Similarity Index, SSIM)，定义为公式(16)。

$\text{MSE}=\frac{{\displaystyle {\sum }_{i=1}^M{\displaystyle {\sum }_{j=1}^N{\left(f\left(i,j\right)-\hat{f}\left(i,j\right)\right)}^2}}}{M\ast N}$ (14)

$\text{PSNR}=10{\log }_{10}\left[\frac{M\ast N\ast {255}^2}{{\displaystyle {\sum }_{i=1}^M{\displaystyle {\sum }_{j=1}^N{\left(f\left(i,j\right)-\hat{f}\left(i,j\right)\right)}^2}}}\right]$ (15)

$\text{SSIM}\left(I_1,I_2\right)=\frac{1}{M\ast N}{\displaystyle {\sum }_{x,y}\frac{\left(2{\mu }_1{\mu }_2+k_1^2{P}^2\right)\left(2{\sigma }_{12}+k_2^2{P}^2\right)}{\left({\mu }_1^2+{\mu }_2^2+k_1^2{P}^2\right)\left({\sigma }_1^2+{\sigma }_2^2+k_2^2{P}^2\right)}}$ (16)

其中，M和N分别为图像的高度和宽度，f和 $\hat{f}$ 对应源图像和重建图像。RMSE值越小表示两幅图像之间的差别越小，重建质量越好。相反，PSNR值越大表示两幅图像之间的差别越小。 ${\mu }_1$ ， ${\mu }_2$ 分别为滑动窗口内像素的均值； ${\sigma }_1$ ， ${\sigma }_2$ 分别为滑窗内像素的方差； ${\sigma }_{12}$ 为两幅图像对应滑窗内像素的协方差； $k_1$ ， $k_2$ 为常数，一般设定 $k_1=0.01$ ， $k_2=0.03$ 。SSIM取值区间为 $\left[-1,+1\right]$ ，其值越大表示图像质量越好。

3.2. 实验结果

3.2.1. 典型实验结果

图4展示了本文重建算法用于对牙齿光滑面、咬合面和牙齿根部三种牙齿图像进行超分辨率重建的效果：图4(a)~(c)分别为三种原始的高分辨率图像；图4(d)~(f)为对应降质后的低分辨率图像，与高分辨率图像相比，降质后的低分辨率图像在细节和边缘上都出现不同程度的迷糊现象，图像分辨率较低；图4(g)~(i)为超分辨率重建的效果，与低分辨率图像相比，图像具有较好的重建效果，极大地改善了图像的质量，以原高分辨率图像为参考，本文提出的算法较好地恢复了图像的分辨率。

Figure 4. Dental image reconstruction effect

图4. 牙科图像重建效果

为了定量分析三种图像超分辨率重建的结果，表1中列出了各图像对应的PSNR和SSIM两个重建指标值。从表中可以看出，本文提出算法对三种牙齿图像进行重建得到的各项指标都具有较好的效果。其中，PSNR值都超过了40 db，结构相似性指标SSIM也都保持在高于0.9的范围。

Table 1. Figure 4 Reconstruction effects

表1. 图4重建效果

为了进一步验证提出的算法对其他医学图像进行重建的能力，图5中显示了CT图像、MR图像和X-ray图像在视觉上达到的分辨率效果。图5(a)为原始的CT高分辨率图像，图5(d)为对应降质后的低分辨率图像，图像出现了视觉上的模糊，图5(g)是重建的CT超分辨率图像，可以看到，低分辨率重建图像有较好的效果；图5(b)是高分辨率脑部MR图像，图5(e)为降质后严重受损的低分辨率图像，图5(h)得到的图像完好地展示了本文提出算法的重建效果；图5(c)、图5(f)和图5(i)分别为X-ray图像的高分辨率效果、低分辨率效果和重建效果。

Figure 5. Other medical image reconstruction effects

图5. 其他医学图像重建效果

表2计算了上述三种测试图像重建效果的PSNR和SSIM指标，根据表中数据可知，本文提出的算法对MR图像和X-ray图像具有更好的重建效果，PSNR和SSIM得到了较大的值。

Table 2. Figure 5 Reconstruction effects

表2. 图5重建效果

3.2.2. 对比实验结果

为了充分验证本文算法的有效性与优越性，以如图6中牙齿图像、CT图像、MR图像和X-ray图像对应的低分辨率图像为研究对象，比较了4种超分辨率重建算法的效果：① 基于插值领域的Cubic插值方法；② Yang提出的稀疏编码超分辨率重建算法；③ Zeyde提出的稀疏表示重建算法；④ 本文提出的利用NSCT提取图像特征的稀疏编码重建算法。

将各种算法用于对该低分辨率图像进行超分辨率重建，不同算法的超分辨率重建结果如图6所示。图(a)、(e)、(i)和图(m)分别为对各低分辨率图像进行经典Cubic插值得到的重建图像，与低分辨率图像相比，通过引用原始图像周围像素点对图像进行插值放大，其像素增加，但图像的分辨率未得到根本的改善；图(b)、(f)、(g)、(n)为Yang提出的经典稀疏编码重建算法得到的结果，图像的分辨率有所改善，但与其他算法相比，重建效果比较一般；Zeyde提出的稀疏编码方法重建效果如图(c)、(g)、(k)、(o)所示，图像的清晰度效果优于前两种方法；图(d)、(h)、(l)、(p)是本文算法重建的高分辨率图像，与其他几种方法相比，图像在分辨率较好。

各种重建算法处理图6中图像得到的各指标定量分析如表3所示。从表3和表4可以看出，本文提出的算法对图中的低分辨率进行超分辨率重建在PSNR和SSIM指标上都达到了最好的效果。

Figure 6. Reconstruction effects of different reconstruction methods

图6. 不同重建方法重建效果

Table 3. Reconstruction effects

表3. 不同算法重建PSNR值

Table 4. SSIM values by different algorithms

表4. 不同算法重建SSIM值

3.2.3. 基于测试数据集的统计结果

为了使本文算法的重建效果更具有说服力，对图像库中300幅低分辨率图像的重建效果进行统计。其中，重建算法为上述4种算法，测试图像集包含了内窥镜牙齿图像、CT图像、MR图像和X-ray图像，并统计了重建的PSNR值和SSIM值的均值，如表5所示。由表可知，提出算法的平均精度与其他3种方法相比，均达到了最好的效果，本文提出重建算法的重建效果表明，PSNR为45.398，SSIM为0.907。其次，Cubic插值和Yang的方法得到的重建效果都相对较差，Zeyde的方法和本文的方法具有较好的效果，且本文的算法稍微具有优势。

Table 5. Average accuracy of the different super-resolution reconstruction algorithms

表5. 不同超分辨率重建算法的平均精度

为了更直观地展示各种重建算法的重建结果，首先将SSIM值乘以100，将其转换到与PSNR值相同数量级上，同时增加超分辨率图像结果的均方误差指标MSE，不同重建方法对应的3个指标在图7中体现。图中误差棒图统计了重建效果的指标平均值和方差。图中四种颜色表示了4种重建算法，横坐标上标注了PSNR、MSE和SSIM，所在位置对应的纵坐标展示了重建的统计效果。针对每个方差图，“*”所在的高度标识了该方法对应指标的平均值，位置越高，说明该方法重建效果具有较好的表现，其次，每个方差线的长度体现了算法的误差大小，长度越短，说明误差越小，即该算法具有较好的稳定性，对各种图像都有较好的重建效果。

Figure 7. Bad plots of reconstruction effect errors by different reconstruction methods

图7. 不同重建方法重建效果误差棒图

对误差棒图分析可以看出，从PSNR指标来看，本文提出的算法和Zeyde算法的重建效果高于另外两种算法，且本文提出算法以略微的优势达到最好的重建效果；从MSE指标来看，MSE指标与PSNR相反，值越小，重建效果越好，从误差图来看，从SSIM指标来看，Yang的方法效果最差，本文提出的算法在图像结构上具有更好的效果和更好的稳定性。

不论是从表5还是图7都可以直观的看出，本文提出算法具有有效性和一定的优越性。

4. 讨论

4.1. 本文算法的优点

本文提出的基于NSCT和稀疏表示的医学图像超分辨率重建算法提高了图像超重建的效果，无论在视觉上还是在客观指标评价上都具有一定的优势。

① NSCT变换相对于简单的一阶、二阶梯度算子来提取低分辨率图像的高频信息，具有明显的优势，NSCT变换采用方向滤波器，具有平移不变性，能够有效保留图像的各个方向的细节信息，使得训练出来的字典更具有有效性和典型性。

② K-SVD奇异值分解法来求取稀疏编码字典，相对于原始Yang采用的有效稀疏编码算法，效率更高；同等特征下，训练出来的字典效果相当。因此，不但能够保证最终的重建质量，而且也能整体提高重建效率。

4.2. 所提出算法的局限性

尽管提出的算法在对图像库中图像的超分辨率重建时，统计结果显示提出的算法具有一定的优势，但从某些特殊图像的重建效果以及算法的复杂度和算法运行时间上来看，仍然需要进行改进。

① 算法的复杂度和算法运行时间仍有进一步提升的空间。NSCT变换虽然能够更有效提取图像的高频信息，但会加大特征提取时的计算量，从而提高了本文算法的整体运行时间。

② K-SVD奇异值分解法，会涉及到迭代的过程，可能会出现陷入局部最小值的情况，因此能够选择一种更为有效的稀疏编码算法至关重要。

③ 某些医学图像的重建效果并未取得较好的效果。从测试数据库的统计结果来看，CT图像、MR图像和X-ray图像的重建效果与自然光牙齿图像的重建效果相比，具有较大差异，提出算法的重建效果指标仍然存在数值较小的情况。

5. 结论

本文提出了一种基于NSCT特征提取的稀疏编码超分辨率重建算法，该算法首先利用NSCT提取训练图像集中图像的特征，得到用于训练的特征集；然后，使用K-SVD算法和OMP算法对特征集进行训练得到高、低分辨率字典对和稀疏表示系数；最后，针对输入的低分辨率图像，通过提取其特征、使用上述训练得到的字典以及OMP算法等步骤，重建得到对应的高分辨率图像。实验结果表明，在主观的视觉效果上，以及峰值信噪比、结构相似性和均方误差等重建结果的客观重建指标上，本文的算法能够有效改善医学图像的超分辨率重建效果。

基金项目

上海市教委2016年青年教师培养计划项目(ZZs115012)资助；上海理工大学光电学院教师创新能力建设项目(1000302006)资助。