数字图像修复技术是图像处理和计算机视觉领域中的一个重要研究内容。图像修复是利用图像中未破损的信息对图像中信息完全丢失的区域进行自动修补的过程，其目的是恢复图像破损区域的信息，并使观察者不易察觉图像曾被修复过。图像修复的应用领域越来越广泛，主要在文物保护、多余目标物体移除、影视制作、虚拟现实等方面有着重要的应用价值。

近年来，有很多学者提出了不同的图像修复算法，主要分为3类：1)基于变分偏微分方程的图像修复；2)基于纹理合成的图像修复；3)混合方法。Bertalmio等人^[1]首次提出了基于偏微分方程的图像修复技术，该方法通过对待修补区域的边界进行不同方向的扩散将未破损区域的信息扩散到待修补区域的内部来对图像进行修复。该算法仅对照片中的细小划痕等小尺度破损具有较好的修复效果。Chan等人^[2]提出了TV(total variation)修复算法，并在此基础上提出了一些改进算法^[3-4]。随后，Criminisi等人^[5]提出了基于样本块的图像修复算法，其利用待修补区域的边界信息计算待修复块的优先级，然后在图像未破损区域中寻找与待修复块最匹配的样本块来填充破损区域。该算法对大面积的破损区域具有较好的修复效果，但是修复时间较长，大大降低了算法的修复效率。因此，针对Criminisi算法的一些不足，学者提出了不同的改进方法^[6-9]。文献[10]首先根据图像的颜色或者纹理特征对破损图像进行分块，生成一个区域图，以便对给定目标块的搜索区域做进一步的限定，提高整体的搜索效率。文献[11]通过引入调节因子来重新定义置信度项，得到新的优先权函数，确保待修复块信息的正确填充；并利用图像的颜色信息和梯度信息的相互结合来寻找最佳匹配块，该算法对结构简单的破损图像修复效果较好，对于背景结构复杂的图像则会出现结构不连续的现象。文献[12]利用图像的边缘结构特征和纹理特征引入差别因子，以改进优先权函数来提高算法填充顺序的准确性，并改进了传统算法的搜索方式，来提高算法的修复效率。文献[13]针对传统方法的时间复杂度较高及错误累积等缺点，通过引入结构张量来构建局部结构度量函数，以此来优化优先权函数，并且提出了一种新的匹配准则。该算法对图像的强边缘结构进行修复时有较好的保持，而在修复弱边缘结构时出现了结构断裂的现象。文献[14]通过结合图像的结构张量和梯度来更加准确地区分图像的显著结构，并定义新的优先级函数来使得图像的结构部分优先修复，该算法在修复复杂结构时，其效果不理想，而且运行时间也很长。文献[15]利用调节因子来决定每个块的优先权，并在匹配条件中加入归一化互相关函数(NCC)来衡量目标块和源块之间的相似性。该算法对简单的结构区域有较好的修复，但在修复纹理区域时，容易产生错误填充块。文献[16]提出了基于结构的图像修复算法，主要利用结构张量场来表示和分析图像的结构层，首先修复破损图像的结构层，然后利用修复后的结构层来约束后续图像的修复过程，该算法对图像的线性结构有较的保持。

然而上述改进算法在修复结构复杂的破损图像时，仍会出现结构不连续的现象，分析其原因，一是块的优先级函数影响待修复块的修复顺序；二是在寻找最佳匹配块时，只利用了图像像素点的颜色信息，而没有充分利用图像的局部特性，同时样本块大小采用固定形式。针对这些不足，本文基于结构张量理论，充分利用彩色图像的结构张量特征，在优先权函数、匹配准则中引入结构张量，以及采用自适应大小的样本块，提高彩色图像的修复效果。

Criminisi算法是一种趋于等照度线的图像采样过程^[5]。$\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}$为图像的待修复区域，${\bf{ \pmb{\mathsf{ δ}} }}{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}}$为填充边缘，$\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}$为已知区域，${\mathit{\boldsymbol{\psi }}_p}$是以像素$p$∈${\bf{ \pmb{\mathsf{ δ}} }}{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}}$为中心的大小为9×9像素的待修复像素块，如图 1所示。这样彩色破损图像的修复主要包含以下3个步骤^[5]：

1) 计算待修复块的优先权，具体公式为

$ P\left( p \right) = C\left( p \right) \cdot D\left( p \right) $

(1)

式中，$C(p)$和$D(p)$分别称为置信度项和数据项，分别定义为

$ C\left( p \right) = \sum\limits_{q \in {\mathit{\boldsymbol{\psi }}_p} \cap \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}} {C\left( q \right)/\left| {{\mathit{\boldsymbol{\psi }}_p}} \right|} $

(2)

$ D\left( p \right) = \left| {\nabla {\bf{I}}_p^ \bot \cdot {\mathit{\boldsymbol{n}}_p}} \right|/\alpha $

(3)

式中，$p$为填充边缘${\bf{ \pmb{\mathsf{ δ}} }}{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}}$上的像素点，$\left| {{\mathit{\boldsymbol{\psi }}_p}} \right|$为待修复块${{\mathit{\boldsymbol{\psi }}_p}}$的面积，$α$为归一化因子，${\mathit{\boldsymbol{n}}_p}$为$p$点的法向量，$\nabla {\bf{I}}_p^ \bot $代表$p$点的等照度线向量。初始化时$C(p)$的值为$C(p)$=0, $\forall p \in \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}$，$C(p)$=1, $\forall p \in \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}$。

2) 搜索最佳匹配块并填充。找到具有最高优先级的待修复块${\mathit{\boldsymbol{\psi }}_{\hat p}}$后，在已知区域中搜索${\mathit{\boldsymbol{\psi }}_{\hat p}}$的最佳匹配块${\mathit{\boldsymbol{\psi }}_q}$，使其${\mathit{\boldsymbol{\psi }}_{\hat p}}$与${\mathit{\boldsymbol{\psi }}_q}$的颜色差值的平方和(SSD)最小，匹配准则为

$ {\mathit{\boldsymbol{\psi }}_{\hat q}} = \arg \mathop {\min }\limits_{{\mathit{\boldsymbol{\psi }}_q} \in \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}} d\left( {{\mathit{\boldsymbol{\psi }}_{\hat p}},{\mathit{\boldsymbol{\psi }}_q}} \right) $

(4)

式中，$d$(${\mathit{\boldsymbol{\psi }}_{\hat p}}$, ${\mathit{\boldsymbol{\psi }}_q}$)表示待修复块${{\psi }_{\hat p}}$和最佳匹配块${\mathit{\boldsymbol{\psi }}_q}$中对应已知像素点的颜色差值的平方和。

$ d\left( {{\mathit{\boldsymbol{\psi }}_{\hat p}},{\mathit{\boldsymbol{\psi }}_q}} \right) = \sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {{{\left( {{p_{ij}} - {q_{ij}}} \right)}^2}} } $

(5)

式中，$m$, $n$分别为待修复块${\mathit{\boldsymbol{\psi }}_{\hat p}}$的长和宽，${p_{ij}}$和${q_{ij}}$分别为待修复块和最佳匹配块中已知的像素值。

3) 更新置信度。具体公式为

$ C\left( p \right) = C\left( {\hat p} \right),p \in {\mathit{\boldsymbol{\psi }}_{\hat p}} \cap \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} $

(6)

式中，${\hat p}$是具有最高优先权的像素点，$p$是待修复块${\mathit{\boldsymbol{\psi }}_{\hat p}}$的破损像素点。

重复步骤1)—3)，直至待修复区域的像素全部被填充。

结构张量通常被用来分析图像的局部几何结构^[17]，它不仅包含了局部区域的强度信息，同时也包含了特定像素邻域梯度的主要方向和这些方向的相干程度。

设彩色图像$\mathit{\boldsymbol{I}}$的梯度为^[5]

$ \nabla {\mathit{\boldsymbol{I}}_i} = {\left[ {{I_{xi}},{I_{yi}}} \right]^{\rm{T}}} $

(7)

式中，i=1, 2, 3分别代表彩色图像的R、G、B分量；$\nabla {\mathit{\boldsymbol{I}}_i}$代表彩色图像的每个通道的梯度，${I_{xi}}$, ${I_{yi}}$分别为$x$, $y$方向的偏导数，${I_{xi}} = \partial {\mathit{\boldsymbol{I}}_i}/\partial x, {I_{yi}} = \partial {\mathit{\boldsymbol{I}}_i}/\partial y$。

则彩色图像的结构张量${\mathit{\boldsymbol{J}}_\rho }$定义为^[17]

$ \mathit{\boldsymbol{J}} = \sum\limits_{i = 1}^n {\nabla {\mathit{\boldsymbol{I}}_i}\nabla \mathit{\boldsymbol{I}}_i^{\rm{T}}} = \left[ {\begin{array}{*{20}{c}} \begin{array}{l} I_x^2\\ {I_x}{I_y} \end{array}&\begin{array}{l} {I_x}{I_y}\\ I_y^2 \end{array} \end{array}} \right] $

(8)

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{J}}_\rho } = \left[ {\begin{array}{*{20}{c}} {{j_{11}}}&{{j_{12}}}\\ {{j_{21}}}&{{j_{22}}} \end{array}} \right] = {\mathit{\boldsymbol{G}}_\rho } * \mathit{\boldsymbol{J}} = }\\ {\left[ {\begin{array}{*{20}{c}} \begin{array}{l} {G_\rho } * I_x^2\\ {G_\rho } * {I_x}{I_y} \end{array}&\begin{array}{l} {G_\rho } * {I_x}{I_y}\\ {G_\rho } * I_y^2 \end{array} \end{array}} \right]} \end{array} $

(9)

式中，${G_\rho }$是均值为0、方差为$ρ$的2维高斯函数^[18]，即${G_\rho } = \frac{1}{{2{\rm{ \mathit{ π} }}{\rho ^2}}}{\rm{exp}}\left({-\frac{{{x^2} + {y^2}}}{{2{\rho ^2}}}} \right)$。由于矩阵${\mathit{\boldsymbol{J}}_\rho }$是对称且半正定的2维矩阵，因此有两个非负特征值，其大小分别为

$ {\lambda _{1,2}} = \frac{1}{2}\left( {{j_{11}} + {j_{22}} \pm \sqrt {{{\left( {{j_{11}} - {j_{22}}} \right)}^2} + 4j_{12}^2} } \right) $

(10)

${\lambda _1}$、${\lambda _2}$分别代表结构张量在像素点处的最大和最小特征值，表明了图像局部边缘的强度^[18]。当${\lambda _1}$和${\lambda _2}$两特征值都相对较小时，意味着在像素点邻域内灰度值变化较小，即像素点位于平坦区域；当${\lambda _1}$较大而${\lambda _2}$较小时，则表明在一个主要方向上，灰度值有很强的变化，该像素点处在图像的边缘区域；当两个特征值都较大时，则灰度值在相应的特征向量所指定的方向上存在很强的变化，像素点处在图像的角点区域。因此，根据特征值${\lambda _1}$、${\lambda _2}$的特点，本文从以下几个方面对Criminisi算法进行改进。

由于结构张量可以准确估计图像的边缘方向和该方向上的变化强度，其特征值之间的大小能够区分图像的边缘区域、纹理区域和平坦区域。因此，将结构张量特征值引入到数据项中，即

$ D\left( p \right) = 1 - \exp \left[ { - {{\left( {{\lambda _1} - {\lambda _2}} \right)}^2}} \right] $

(11)

从式(11)中可以看出，当像素点$p$处在结构边缘时，${\lambda _1}$-${\lambda _2}$>>0，则$D(p)$≈1；当$p$处在平坦区域时，${\lambda _1}$-${\lambda _2}$≈0，则$D(p)$≈0；当$p$处在纹理区域时，${\lambda _1}$－${\lambda _2}$>0，则0 < $D(p)$ < 1。这样，数据项$D(p)$能更加精准地反映图像结构特征。

在传统算法中，由于优先权函数$P(p)$是利用置信度项与数据项的乘积来确定(如式(1)所示)，但当数据项或置信度项在某点处为零时，其优先权也随着为零，从而影响图像的填充顺序，产生错误的修复结果。为此，本文则采用改进的数据项与置信度项的线性加权和，重新定义优先权函数，即

$ P\left( p \right) = \alpha \cdot C\left( p \right) + \beta \cdot D\left( p \right) $

(12)

式中，$α$、$β$为置信度项和数据项的权重，且$α$+$β$=1。为使图像的结构部分优先修复，应使数据项的权重大于置信度项的权重，即$α$ < $β$。

图 2分别给出了采用不同优先权函数对破损彩色图像Lena的修复效果图。从图 2对比可以看出，改进的优先权函数使得填充顺序更加准确，图像的纹理部分及结构部分都得到了较好的修复。

在Criminisi算法中，通常采用固定的样本块大小进行寻找最佳匹配块，而没有考虑图像的局部结构特征，从而导致在修复结构比较复杂的区域时，会出现部分结构不连续或延伸现象。为此，本文利用结构张量特征值的大小，将图像划分为边缘区域、纹理区域和平坦区域，并在不同区域，采用不同大小的样本块进行最佳匹配块的搜索。首先，参照文献[18]定义一个以像素点$p$为中心的待修复块${{\mathit{\boldsymbol{\psi }}_p}}$的平均相干因子$av{e_{moa}}$，其大小为

$ av{e_{moa}} = {\left( {\frac{{{{\bar \lambda }_1} - {{\bar \lambda }_2}}}{{{{\bar \lambda }_1} + {{\bar \lambda }_2}}}} \right)^2} $

(13)

式中，${{\bar \lambda }_1}$、${{\bar \lambda }_2}$为待修复块${{\mathit{\boldsymbol{\psi }}_p}}$中已知像素点的结构张量的平均特征值。表 1分别给出了图 3中各像素块的平均相干性$av{e_{moa}}$大小。在图 3中标出的像素块大小为9×9像素，其中像素块1、2、3、4、5处在边缘区域，像素块6、7、8、9处在纹理区域，像素块10、11、12处在平滑区域。从表 1中可以看出，像素块处在边缘区域时，$av{e_{moa}}$的值在0.90以上；当处在平滑区域时，$av{e_{moa}}$的值在0.70以下；而处在纹理区域时，$av{e_{moa}}$的值在0.80左右。因此，图像像素块的平均相干因子$av{e_{moa}}$的大小，能够较准确地区分图像的平滑区域、纹理区域和边缘区域。

为了在修复结构复杂的破损图像时，既可以对图像的平坦区域进行很好的修复，也可以保持边缘结构的连续性。因此，本文利用$av{e_{moa}}$的变化来自适应地选择样本块大小，其选择原则是，当待修复块的平均相干性较大时，则处于图像的边缘区域，应采用较小的样本块；当其平均相干性较小时，则处在图像的平坦区域，应采用较大的样本块。即样本块大小$size(p)$为

$ size\left( p \right) = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} 11 \times 11\;像素\\ 9 \times 9\;像素\\ 7 \times 7\;像素\\ 5 \times 5\;像素 \end{array}&\begin{array}{l} \left| {av{e_{moa}}} \right| \le 0.7\\ 0.7 < \left| {av{e_{moa}}} \right| < 0.8\\ 0.8 \le \left| {av{e_{moa}}} \right| < 0.9\\ \left| {av{e_{moa}}} \right| \ge 0.9 \end{array} \end{array}} \right. $

(14)

图 4分别给出了采用不同样本块大小时，对Lena破损图像(图 2(b))修复结果的对比图。通过对比得出，采用自适应样本块大小较好地克服了固定大小的样本块容易在结构处产生错误延伸的现象。

为了能更准确地找到最佳匹配块，避免Criminisi算法只根据待修复块和已知区域中样本块的像素颜色值来判断最佳匹配块的不足，本文结合结构张量特征值的特性，在式(5)的基础上，利用像素的结构张量相干因子${moa}$，将匹配准则重新定义为

$ \begin{array}{*{20}{c}} {d\left( {{\mathit{\boldsymbol{\psi }}_{\hat p}},{\mathit{\boldsymbol{\psi }}_q}} \right) = \sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {\left[ {{{\left( {{p_{ij}} - {q_{ij}}} \right)}^2} + } \right.} } }\\ {\left. {10 \times {{\left( {mo{a_{ij}} - mo{{a'}_{ij}}} \right)}^2}} \right]} \end{array} $

(15)

式中，$moa = {\left({\frac{{{\lambda _1}-{\lambda _2}}}{{{\lambda _1} + {\lambda _2}}}} \right)^2}$为像素的结构张量相干因子；$mo{a_{ij}}$、$moa{'_{ij}}$分别为待修复块和匹配块中已知像素点的结构张量的相干因子。

综上所述，本文算法的流程如图 5所示，具体实现步骤如下：

1) 确定彩色图像待修复区域的边界；

2) 根据式(12)，计算边界上像素点的优先权，并找出具有最大优先权的待修复块${\mathit{\boldsymbol{\psi }}_{\hat p}}$；

3) 根据式(14)，自适应确定待修复块${\mathit{\boldsymbol{\psi }}_{\hat p}}$的最佳大小；

4) 根据式(15)，分别在R、G、B的已知区域中寻找最佳匹配块；

5) 将最佳匹配块分别复制到R、G、B分量的待修复块中；

6) 根据式(6)更新置信度项；

7) 重复步骤1)— 6)，直到破损区域中的像素被填充完。

为了验证改进算法的有效性，在3.30 GHz处理器、2 GB内存的计算机上，利用Matlab R2011软件进行了仿真实验。下面给出参数为$ρ$=0.4、$α$=0.4、$β$=0.6的实验仿真结果，同时也给出了本文算法与Criminisi算法^[5]、文献[13]算法(LLZ)、文献[14]算法(SYM)和文献[15]算法(WLP)的对比结果，彩色图像的修复效果采用峰值信噪比(PSNR)、结构相似度(SSIM)[19]及运行时间来进行客观评判。其中两幅图像$\mathit{\boldsymbol{X}}$和$\mathit{\boldsymbol{Y}}$的结构相似度定义为^[19]

$ S = \frac{{\left( {2{u_X}{u_Y} + {c_1}} \right)\left( {2{\delta _{XY}} + {c_2}} \right)}}{{\left( {u_X^2 + u_Y^2 + {c_1}} \right)\left( {\delta _X^2 + \delta _Y^2 + {c_2}} \right)}} $

(16)

式中，${u_X}$、${u_Y}$分别为图像$\mathit{\boldsymbol{X}}$、$\mathit{\boldsymbol{Y}}$的均值，${\delta _X}$和${\delta _Y}$为其方差，${\delta _{XY}}$为$X$和$Y$的协方差，${c_1}$=6.502 5，${c_2}$=58.522 5。

图 6分别给出了Criminisi算法^[5]、文献[13]算法、文献[14]算法、文献[15]算法和本文算法对Lena彩色破损图像的修复结果。为了比较方便，并对修复结果作了局部放大。从图 6中看出，Criminisi算法的修复结果产生了一些额外的填充块，并在曲线结构处(如肩膀部位)未能保持结构的连续性；文献[13]和文献[14]算法对图像的纹理部分修复较好，但对曲线结构的修复不理想，出现了断裂或延伸现象；文献[15]算法对图像的结构边缘有较好的保持，但在纹理区域出现了模糊现象；从图 6(g)中可以看出，本文算法对不同的破损区域都有较好的修复结果。

图 7是5种算法对纹理比较丰富的Baboon彩色破损图像的修复结果。从图 7中看出，Criminisi算法^[5]和文献[14]算法在修复图像的边缘(如左脸颊部位)时产生了明显的结构延伸现象，并在纹理区域填充了错误的匹配块；文献[13]和文献[15]算法在修复图中的鼻子时，其结果不太符合人眼的视觉特性。而本文算法对线性结构和纹理区域都有较好的保持，克服了Criminisi算法所存在的不足，同时也说明，在寻找最佳匹配块时，采用自适应大小的样本块能对复杂结构有较好的修复。

图 8是对结构和纹理都比较复杂的Barb彩色划痕图像进行修复的结果。从图 8(c) (d) (e) (f)可以看出，4种算法都对图像的纹理部分修复较好，但Criminisi算法^[5]和文献[13]算法在胳膊与膝盖的交界处等结构边缘出现了明显的错误延伸现象；本文算法对纯纹理和强边缘区域修复的较好，而在纹理和边缘的交界处(如图中手腕处)仍有一点延伸现象。

图 9是对结构相对简单的Peppers彩色破损图像的修复结果。从图 9(c) (e) (f)中可以看出，Criminisi算法^[5]、文献[14]和文献[15]算法在修复平坦区域的过程中产生了错误的填充块；文献[13]算法在修复弱边缘结构时出现了结构断裂的现象，如图 9(d)所示，而本文算法在整体结构修复效果较好。

在对比实验中，也对物体移出进行了仿真实验。图 10是对大物体移除的修复结果图。从图 10(c) (d) (e) (f) 中看出，Criminisi算法^[5]、文献[13]、文献[14]和文献[15]算法在修复大面积的破损区域时容易产生错误的匹配块，导致修复效果不佳；而本文算法克服了上述缺点，使得破损边界与周围区域没有产生块效应。

图 11是对彩色图像window文字去除的修复结果对比。从图 11(d)中看出，文献[13]在部分结构和纹理处产生了错误的延伸；而Criminisi算法^[5]、文献[14]和文献[15]算法在图像的边缘结构处出现了明显的断裂，但对平滑部分都有较理想的修复效果，如图 11(c) (e) (f)所示；而本文算法对图像的平滑区域和边缘结构的修复结果则较为理想。

图 12分别给出了Criminisi算法^[5]、文献[15]算法和本文算法对破损线条的修复过程，其中图 12(c) 图 (d) (e)是3个算法的最终修复结果，图 12(f) (g) (h)为各算法修复过程的中间部分结果(图中修复效果，从左至右按修复步数增加排列)。从图 12(c) (d) (e)中可以看出，Criminisi算法和文献[15]算法的修复效果不佳，出现了明显的断裂现象；而本文算法在修复线性结构时，较好地保持了图像边缘的连续性，其修复结果也符合人眼视觉特性；通过对比图 12(f) (g) (h)可以得出，对破损区域块进行修复时，其正确的填充顺序对整个修复结果至关重要。

表 2和表 3分别给出了各算法的峰值信噪比(PSNR)、结构相似度(SSIM)及运行时间。从表中可以看出，相比于Criminisi算法^[5]、文献[13]算法、文献[14]算法和文献[15]算法，本文算法的PSNR值相对提高了1~3 dB，SSIM也得到改善；但从运行时间来看，本文算法高于其他4种算法，这是因为在修复过程中，采用自适应样本块大小来搜索最佳匹配块以及像素相干因子的计算，都会增加一定的运行时间，从而降低了图像的修复效率。

针对传统Criminisi算法的不足，本文充分利用结构张量能够体现图像的局部结构特性这一特点，利用结构张量对优先权函数进行改进，使得修复顺序更加准确；将结构张量作为自适应选择样本块大小的依据，实现不同区域采用不同样本块大小进行修复，以克服结构不连续和纹理错误延伸的缺点；在最佳匹配准则中引进结构张量特征值信息来减少错误匹配率。实验结果表明，相比于传统Criminisi算法，本文算法对复杂的线性结构和纹理区域都有较理想的修复效果，有效地提高了图像的修复质量，其PSNR值相对提高了1~3 dB；并与其他相关算法相比，其PSNR值和SSIM值都有所提高。但是，改进算法仍在一些不足，如在弱纹理、纹理与边缘的交界处修复不理想，以及修复相对耗时，这些将是今后需要进一步研究解决的问题。