曲率驱动扩散模型的影像修复

数学基础

基于TV模型影像修复的局限性

在TV模型1中，扩散强度只取决于对比度或等照度线的强度，在表达式中表现为：
$$\hat{D}=\frac{1}{|\Delta u|}\hspace{1cm}(1)$$
根据上式2，扩散强度并不取决于等照线的集合信息，对于一个曲面，其几何形状通过一个标量，曲率$k$来描述，对于文章3中图2TV模型失效的原因在于，在四个交点，其曲率为$\infty$相反门根据心理学校正在图中的校正结果，所有等照度线应该被处理的尽可能伸展和平整，因此曲率应该尽可能的小才对。

曲率扩散校正模型的数学基础

正是由于基于全局变差模型的校正方法存在以上问题，因此需要对模型进行改进将原始TV模型求$\hat{D}$的方法：
$$\hat{D}=f(|\Delta u|)\,,\hspace{1cm}(2)$$
修改为：
$$\hat{D}=\frac{g(|k|)}{\Delta u}\,,\hspace{1cm}(3)$$
其中$g(s)$是一个去除曲率极大值，稳定曲率极小值的函数，通过此函数，在曲率具有较大值的时候有比较强的扩散作用。因此对于文章图2中所示模型，能够更好的符合人类感知。在文章中函数$g(s)$为：
$$g(s)=s^p\,,\quad s>0\,,p\ge1\hspace{1cm}(4)$$
而在点$x$上的曲率$k$为等照线的量化曲率方向通过如下公式获取：
$$k=\Delta[\frac{\Delta u}{|\Delta u|}]\,.\hspace{1cm}(5)$$
因此CDD模型为：
$$\left\{\begin{array}{lr} \frac{\delta u}{\delta t}(or 0)=\Delta\cdot[\frac{g(|k|)}{|\Delta u|}\Delta u]\,,\quad x\in D \\ u=u^0\,,\quad x\in D^c \end{array}\hspace{1cm}(7) \right.$$
则求解时间进行的方程，初始值满足：$u(x,0)=u^0(x)\,,x\in D^c\,.$则曲率驱动的扩散场流为：
$$\textbf{j}=-\hat{D}\Delta u=-\frac{g(|k|)}{|\Delta u|}\Delta u\,,\hspace{1cm}(8)$$
在大多数情况下原始影像是存在噪声的，因此如果直接使用公式（7）进行修复，则对噪声比较敏感，文章提出了两种方式解决噪声的问题：
1.在使用CDD模型进行修复之前，先进行去噪处理；
2.在实用CDD模型的过程中进行去噪，对于此过程分为两个自然阶段，在修复区域内部使用公式（7）所示的CDD影像修复模型，而在正常的区域采用如下模型：
$$\frac{\delta u}{\delta t}(or 0)=\Delta\cdot[\frac{\Delta u}{\|\Delta u\|}]+\lambda(u-u^0)\quad(8)$$

数值解法

同样是以公式（7）为例，其数值解法为：
$$\frac{\delta u}{\delta t}=-\Delta \cdot\textbf{j}$$
则

TV模型数学理论

$\hspace{1cm}$TV模型与基于偏微分方程的主要用于修复较小的图像缺失区域和噪声，假设$u^0$是一个定义在二维平面$\Omega$上的光滑影像函数， $D$是一个待修复的区域，其区域直径为$d$，令$u^0|_d$为在区域$D$上$u^0$的限制，则影像的修复过程可以描述为获取$u_D$在$u^0|_d$上的拟合。
$\hspace{1cm}$如果区域$D$的直径$d$收缩到0，则可以采用线性插值的方法获取任意光滑的测试数据$u^0$：
$$\lVert u_D-d^0|_D\rVert_{\infty}=O(d^2)\hspace{1cm}(1)$$ 同样，影像修复过程可以使用k阶修复方式：
$$\lVert u_D-d^0|_D\rVert_{\infty}=O(d^{k+1})\hspace{1cm}(2)$$

二阶格林公式1影像修复

数学基础

假设$\Delta$表示拉普拉斯算子： $$\Delta u:=\frac{\delta^2u}{\delta x^2}+\frac{\delta^2u}{\delta y^2}\,.\hspace{1cm}(3)$$ 在区域$D$上的二阶格林公式为： $$\int_D(u\Delta v-v\Delta u)dxdy=\int_{\Gamma}(u\frac{\delta v}{\delta \textbf{n}}-v\frac{\delta u}{\delta\textbf{n}})ds\,,\hspace{1cm}(4)$$ 其中：$u$和$v$为在闭合环$D$上定义的二维函数，n为区域边界$\Gamma$上定义的朝外的法向量，$s$为边界的长度。
$\hspace{1cm}$假设$G(z_0,z)$是一个基于在$D$上基于珀松方程的格林函数，格林函数的意义为：$z_0$处所产生的脉冲对$z$处的场的作用，他是距离的函数，则对于任意一个$z_0=(x_0,y_0)\in D$，作为对应的点$z=(x,y)\in D$则$G(z_0,z)$求解为：
$$-\Delta G=\delta(z-z_0)\,,G|_{\Gamma}=0\,.\hspace{1cm}(5)$$ 则应用格林二阶函数（$u=u^0(z),v=-G(z_0,z)$）则讲$u,v$代入公式（4）中，可以得到：
$$u^0(z_0)=\int_\Gamma u^0(z(s))\frac{\delta(-G(z_0,z))}{\delta\textbf{n}}ds+\int_DG(z_0,z)(-\Delta u^0(z))dz\hspace{1cm}(6)$$ 其中$dz=dxdy\,.$公式(6)中的第一项表示为$u^h(z_0)$表示$f=u^0|_\Gamma$的光滑的扩展。

---

到公式（6）开始我就开始看看不懂了，简单的来说公式（4）就是格林公式，然后将$u,v$代入到公式（4）中得到公式（6），但是具体怎么得到的由于数学基础不够没法进行数学证明

---

$$dw_{z_0}=\frac{\delta(-G(z_0,z))}{\delta\textbf{n}}ds\hspace{1cm}(7)$$ 为边界$\Gamma$与源点$z_0$关联的和谐程度的度量，而非和谐程度可以表示为$u^a=u^0-u^h$且满足泊松方程2: $$\Delta u^a(z)=\Delta u^0(z)\,,z\in D\,,\qquad and\,\,u^a|_\Gamma=0\hspace{1cm}(8)$$
假设$d$为待修复影像区域$D$的直径$G(z_0,z)$是泊松分布的格林函数，则有以下条件成立：
$$\int_DG(z_0,z)dxdy\le\frac{d^2}{4}\,.\hspace{1cm}(9)$$
公式（9）的证明有兴趣的可以参看原文3 ，基于上述理论我们可以很简单的建立准确的基于二阶格林公式的影像修复方法：
（a）Linear inpainting via harmonic extension 假设我们只通过和谐扩展修复$u^0|_D$，则线性修复的表达为：
$$\lVert u^h-u^0|_D\rVert_\infty=O(d^2)\hspace{1cm}(10)$$
当$d\to0$时公式（10）成立，根据公式（6）均匀性修复误差正好就是非均匀性组份$u^a$。由于$u^0$为固定的光滑函数，因此存在常量$M$使得满足:
$$|\Delta u^0(z)|\le M\hspace{1cm}(11)$$
对于所有的$z\in D\,.$上式都成立，且对于$z_0\in D$根据公式（9）则有：
$$|u^a(z_0)|\le M\int_DG(z_0,z)dz\le \frac{Md^2}{4}\,.\hspace{1cm}(12)$$
（b）Cubic inpainting via Green’s formula 为了提高影像修复的精度，除了需要对整体结构进行修复外还需要对细节组份$u^a$进行修复。假设$u^\Delta_D$为线性修复的结果，然后我们修复$u^a|_D$通过$u^a_D$根据如下积分方程：
$$u^a_D(z_0)=\int_DG(z_0,z)(-u^\Delta_D(z))dz\,,\hspace{1cm}(13)$$
或者通过求解已知的泊松方程：
$$-\Delta u^a_D(z)=-u^\Delta_D(z)\,,\quad z\in D;\quad u^a_D|\Gamma=0\hspace{1cm}(14)$$
最后将这些新的细节信息加到均匀性修复过程中，可疑得到更加准确的修复结果$u_D$
$$u_D(z)=u^h(z)+u^a_D(z)\,.\hspace{1cm}(15)$$

局限性

实际上上述模型在大部分情况下是不适用的，因为：
1.影像通常情况下来说并不是平滑的函数，而是包含各种边界信息和非连续性的
2.影像在统计上也不能完全一致，由于存在影像噪声

全局变差模型的影像修复

模型数学基础

如上图所示，其中$D$为待修复区域，待修复区域具有分段光滑的边界$\Gamma$，$E$为待修复边界的外环。假设$u^0|_E$是被白噪声污染的区域，则基于变差修复模型为在扩展域找到变差函数$u$使得如下规则函数取得最小值：
$$R[u]:=\int_{E\bigcup D}r(|\Delta u|)dxdy\,,\hspace{1cm}(16)$$
在以下约束下：
$$\frac{1}{Area(E)\int_E|u-u^0|^2dxdy=\delta^2\hspace{1cm}(17)}$$
对式（16）（17）求解其欧拉拉格朗日方程为：
$$J_\lambda[u]=\int_{E\bigcup D}|\Delta u|dxdy+\frac{\lambda}{2}\int_E|u-u^0|^2dxdy\,,\hspace{1cm}(18)$$
其中$\lambda$为拉格朗日乘数子，则剩下的问题为求解这个欧拉拉格朗日方程。

FMM影像修复方法

数学描述

$\hspace{1cm}$以上为FMM方法的示意图，图中$p$为待修复的点，$\delta\Omega$为待修复区域的边界$\Omega$为待修复的区域，对于一个待修复点$p$周围的的邻域小区域$B_\epsilon(p)$，区域大小为$\epsilon$，如Figure2图（a）所示，则未知点$p$的像素值应该由已知的区域$B_\epsilon(p)$得到，首先考虑灰度影像，彩色影像是灰度影像的扩展。如果区域$\epsilon$足够小的话，我们认为使用一个一阶导的近似可以进行修复，进行修复的公式为：
$$I_{q}(p)=I(q)+\Delta I(q)(p-q)\,.(1)$$
$\hspace{1cm}$上述公式表明在位置$p$的校正结果与$q$点位置和$q$点位置和$p$点位置距离差值的一阶近似，也就是$q$点的1阶近似（在区域比较小的时候比较适用）。接下来我们修复$p$点实用区域$B_\epsilon(p)$内所有像素进行处理，对每一个像素都有一个权值，则$p$的修复后的像素值为：
$$I(p)=\frac{\Sigma_{q\in{B_\epsilon(p)}}\omega(p,q)[I(q)+\Delta I()q(p-q)]}{\Sigma_{q\in{B_\epsilon(p)}}\omega(p,q)}\,.(2)$$
其中$\omega(p,q)$为权函数，权函数的设计既要考虑到$p$点信息的传递也要考虑到在$B_\epsilon(p)$区域影像的细节。

FMM影像修复

$\hspace{1cm}$上一章节描述了影像修复的基本原理和数学基础，则对整个区域的$\Omega$的修复方法为，对处于$\delta\Omega$边界上的点根据公式（2）进行修复，然后进行迭代，直到整个区域都修复完成，在这里有一个修复顺序的问题，在处理过程中使用待修复边界点与原始边界点的距离进行判断，距离原始边界越越近的点首先进行修复，这是因为距离边界越近的点修复的可靠性越大。对于上述操作可以采用FMM方法进行处理，简单的说FMM方法就是求解 $Eikonal$方程1

$$|\Delta T|=1\hspace{1cm}on \hspace{3pt}\Omega,\hspace{1cm}with\,T=0\,on\,\delta\Omega\,.(3)$$
其中方程（3）中对$T$的求解就是$\Omega$中像素距离初始边界$\delta\Omega$的距离，整个求解过程的伪代码为：

$\delta\Omega_i\,$=boundary of region to inpaint
$\delta\Omega\,$= $\,\delta\Omega_i$
while($\delta\Omega\,$not empty)
{
$p\,$=$\,$pixel of $\delta\Omega\,$closest to $\delta\Omega_i\,$
inpaint $\,p$ using Eqn (2)
advance $\delta\Omega\,$into $\Omega$
}

FMM方法的主要优势在于此方法明确平衡了区分已知和未知区域窄带，确定了下一个需要进行修复的像素，这里的像素窄带表示的就是我们的边界$\delta\Omega$，则对于每一个像素除了储存像素的灰度值之外还储存一个标志变量，这个标志变量包含三个值

• BAND：像素值属于窄带，其值经过更新。
• KNOWN：在边界$\delta\Omega$之外的点，也就是在已知区域的点，其T和灰度值I是已知的。
• INSIDE：在边界的像素值点，在边界区域中，T和I的值目前是未知的。

则FMM算法的初始值和传播相位如下：首先设置T在边界$\delta\Omega$上或边界外值为0，在待修复边界内设置极大值，对所有的影像都按上述条件设置初始值

while (NarrowBand not empty)
{
extract P(i,j) = head(NarrowBand); /* STEP 1 */
f(i,j) = KNOWN;
for (k,l) in (i1,j),(i,j1),(i+1,j),(i,j+1)
if (f(k,l)!=KNOWN)
{
if (f(k,l)==INSIDE)
{
f(k,l)=BAND; /* STEP 2 */
inpaint(k,l); /* STEP 3 */
}
T (k,l) = min(solve(k1,l,k,l1), /* STEP 4 */
solve(k+1,l,k,l1),
solve(k1,l,k,l+1),
solve(k+1,l,k,l+1));
insert(k,l) in NarrowBand; /* STEP 5 */
}
float solve(int i1,int j1,int i2,int j2)
{
float sol = 1.0e6;
if (f(i1,j1)==KNOWN)
if (f(i2,j2)==KNOWN)
{
float r = sqrt(2(T(i1,j1)T(i2,j2))*(T(i1,j1)T(i2,j2)));
float s = (T(i1,j1)+T(i2,j2)r)/2;
if (s>=T(i1,j1) && s>=T(i2,j2)) sol = s;
else
{ s += r; if (s>=T(i1,j1) && s>=T(i2,j2)) sol = s; }
}
else sol = 1+T(i1,j1));
else if (f(i2,j2)==KNOWN) sol = 1+T(i1,j2));
return sol;
}

对于单个像素的修复伪代码如下：

void inpaint(int i,int j)
{
for (all (k,l) in Bε(i,j) such that f(k,l)!=OUTSIDE)
{
r = vector from (i,j) to (k,l);
dir = r * gradT(i,j)/length(r);
dst = 1/(length(r)*length(r));
lev = 1/(1+fabs(T(k,l)T(i,j)));
w = dir*dst*lev;
if (f(k+1,l)!=OUTSIDE && f(k1,l)!=OUTSIDE &&
f(k,l+1)!=OUTSIDE && f(k,l1)!=OUTSIDE)
gradI = (I(k+1,l)I(k1,l),I(k,l+1)I(k,l1));
Ia += w * (I(k,l) + gradI * r);
s += w;
}
I(i,j) = Ia/s;
}

FMM算法细节

首先使用FMM方法计算初始待修复边界外的区域，获取待修复边界外部距离。由于我们只需要距离边界点小于$\epsilon$的点，因此只在边界外小于$\epsilon$的区域内计算T，这样我们就提高了计算速度，另外我们对边界内的点运行FMM进行处理，因此，整个影像上的T可以表现为：
$$\begin{equation} T(p)= \left\{ \begin{array}{lr} &T_{in}(p)\, & if \,p\in\Omega \\ -&T_{out(p)}\,&if\,p\not\in\Omega\, \end{array} \right. \,.(5) \end{equation}$$
然后通过$3\times 3$滤波算子进行滤波，然后计算$\Delta T$通过中心差分。

以上为文修复的效果，从上图可以看出，与BSCB方法相比，文章提出的方法在细节上表现并没有比BSCB方法好，但是文章的方法处理效率要远高于BSCB方法2。

局限性

对于以上采用微分校正的方法只能校正比较小的，或者线性的缺失区域，缺失范围较大的区域，校正方法会造成过度平滑，并不能很好进行处理，因此对于比较大的缺失需要进一步处理，但是采用微分校正的方法可疑获取校正趋势3。

结合纹理合成的多分辨率影像修复方法

引言

$\hspace{1cm}$原则上，使用纹理合成的方法能够对数字影像进行修复，特别是在一些纹理结构比较强的区域，通过纹理合成的方法进行修复可疑获取比较好的影像修复效果，然而如果待修复的区域包含额外的颜色或者是强度梯度信息，则采用纹理合成的方法就不能达到理想的效果了。影像修复的方法可以作为纹理合成方法的补充，影像修复方法一般是通过待修复区域周围像素信息传递对待修复区域进行修复的，与基于纹理合成的方法不同，影像修复的方法能够正确处理影像强度信息，但是对于纹理细节并不能够进行很好的恢复。

修复方法

整个修复方法的处理流程为：

fig 1

$\hspace{1cm}$从Fig2的处理流程我们可以看出整个处理过程为：首先对影像进行分解得到高频和低频的分量，对高频和低频的分量分别进行处理。低频分量主要包含影像的强度信息，包含较少的细节信息，因此对低频的修复着重修复影像强度，而基于PDE1的影像修复方法对于纹理细节比较少的区域能够较好的保留其强度信息，而对于高频分量，将其通过高斯金字塔分解为n+1层，假设每一层为$H_i\,(i=0,...,n)\,.$从第$H_n$层开始采用多分辨率纹理合成方法2对每一层的待修复区域进行修复。将修复后的低频成分和高频成分结合得到修正后的影像。

频率分解

$\hspace{1cm}$处理过程中第一步就是要获取影像$I$分解后得到高频和低频分量。首先通过离散DCT3变换将原始影像变换为一系列的子带，然后我们选取前$k$个子带（前$k$个子带包含了主要信息），然后通过前$k$个子带进行DCT逆变换，通过逆变换后就可以得到低频影像$L_k$。而相应的高频影像则通过原始影像减去低频影像得到：$H_k:=I-L_k\,.$在这个过程中低频组份选取的$k$值确定了低频组份的模糊程度，$k$值越大则低频组份中的细节信息越丰富，$k$值越小则说明低频组份越模糊，实际上$k$值的选取对影像修复效果有比较大的影响。在校正过程中我们希望能够尽量大的去除低频成分中的细节信息只保留强度信息。而我们面临的问题在于如何确定$k$值的问题，对于这个问题文章给出了一个解决方案4。首先获取高频成分后计算高频成分$H_k$的自相关矩阵5$A_k$：$A_k:=DCT^{-1}(DCT(H_k)\cdot DCT(H_k))\,.$如果输入影像不是一个矩形，则非矩形部分通过0值补全，求得$A_k$之后剩下的工作为求取$A_k$的标准差，取$k$值使得$A_k$的标准差小于0.001，则最小的$k$值为选取的值，文章认为此时能够取得比较好的校正效果。

高斯影像金字塔

$\hspace{1cm}$通过上述过程分别得到高频和低频分量之后剩下的步骤为对高频分量进行高斯金字塔分解，实际上高斯金字塔分解是一个很简单的过程：
$$\begin{eqnarray*} && \omega(u,v)=\widehat{\omega}(u)\,\widehat{\omega}(v)\\ && \widehat{\omega}(0)=0.4,\widehat{\omega}(\pm1)=0.25,\widehat{\omega}(\pm1)=0.05 \end{eqnarray*}$$
高斯金字塔6必须要求影像大小满足$(p2^N+1)\times(q2^N+1)$像素个数且（$p\,,q\,,N\in\mathbb N$）在这样大小的影像下可以保证影像高斯金字塔具有$N+1$层。在此情况下我们将影像分解为$n$层且$n\ll N+1\,.$整个过程可以杯描述为：在纹理合成的第$i$层像素信息可以从更高层$(i+1,...,n)$的核中获取。在实际处理过程中进行三层金字塔分解就能获取比较好的效果。

纹理合成

纹理合成的难点有三点7：1）如何定义纹理单元和纹理内容；2）如何重构纹理概率分布；3）如何在二维合成过程中进行线性化。单像素合成：假设$I$是一个通过纹理样本影像合成的影像$I_{smp}\subset I_{real}$,其中$I_{real}$是真实的无限的纹理。假设$p\in I$是一个像素，$w(p)\subset I$是一个以$p$为中心的一个模版，定义$d(w_1,w_2)$为两个块之间的距离测度，假设除$p$以外$w(p)$中其他像素值都是已知的，在基于MRF的纹理合成模型为：
$$\Omega(p)= { w^`\in I_real : d(w^`,w(p))=0}$$
通过模板$w(p)$在整个影像上做匹配，找到最佳匹配点，采用距离测度则为找到一个匹配区域使得距离最小，而距离测度描述为：$d=d_{SSD}*G$，其中$G$为二维高斯核。找到最佳匹配点之后可以通过纹理合成算法对点$p$进行纹理修复。

处理结果

$\hspace{1cm}$从上图结果可以看出文章算法具有比较好的校正效果，整个校正过程有两个参数起重要作用：1.DCT变换后子带个数的选取，子带个数的选取决定了低频组份$L_k$；2.对高频成分进行高斯金字塔分解。实际上通过标准差小于0.001的范围选取的$k$值能够获取比较好的处理结果。这样对于子带数目$k$的选取就可以自动化进行。对高频组份的金字塔分解选取合适的分解层数是一个比较困难的问题，通常在处理中使用两层或三层金字塔。