最近项目中要用到修复技术，看opencv里有两种算法，我先尝试了一下Telea在2004年提出的基于快速行进的修复算法（后面简称FMM算法）。找到作者的原文看了一下，对算法有了一定了解,记录一下。
论文题目：An Image Inpainting Technique Based on the Fast Marching Method （2004）
论文下载： http://www.cs./~alext/PAPERS/index.html （编号36的那篇）
|
1
2
3
4
5
6 |
void cv::inpaint( const Mat& src, const Mat& mask, Mat& dst, double inpaintRange, int flags ) //src:要修复的图像； //mask：修复模板，必须是单通道图像； //dst：目标图像； //inpaintRange：选取邻域半径; //flags:要使用的方法，可以是CV INPAINT NS或CV INPAINT TELEA（本文介绍的方法）。 |
|
1 |
cvInpaint( const vArr*_input_img,const CvArr* _inpaint_mask,CvArr* _output_img, double inpaintRange, int flags ) |
首先，FMM算法基于的思想是，先处理待修复区域边缘上的像素点，然后层层向内推进，直到修复完所有的像素点。
下面以灰度图为例，我们只需要计算出像素新的灰度值即可。对于彩色图像，分别用同样的方法处理各个通道即可。
一、先说一下如何修复一个像素点的。
参考上图，Ω区域是待修复的区域；δΩ指Ω的边界）；要修复Ω中的像素，就需要计算出新的像素值来代替原值。
现在假设p点是我们要修复的像素。以p为中心选取一个小邻域B(ε），该邻域中的点像素值都是已知的（只要已知的）。（这个ε就是opencv函数中参数 inpaintRadius）
q为 Bε（p）中的一点，由q点计算P的灰度值公式如下：
I q (p) = I(q) + ? I(q)(p ? q)（? I(q)是q点的亮度梯度值）
显然，我们需要的是用邻域Bε（p）中的所有点计算p点的新灰度值。显然，各个像素点所起的作用应该是不同的，也就引入了权值函数来决定哪些像素的值对新像素值影响更大，哪些比较小。采用下面的公式（公式2）：
这里的w(p, q)就是权值函数，是用来限定邻域中各像素的贡献大小的。
w(p, q) = dir(p, q) · dst(p, q) · lev(p, q)
其中，d0和 T0分别为距离参数和水平集参数，一般都取为 1。方向因子 dir(p,q)保证了越靠近法线方向 N = ?T的像素点对 p 点的贡献最大；几何距离因子 dst(p,q)保证了离 p 点越近的像素点对p 点贡献越大；水平集距离因子lev(p,q)保证了离经过点 p 的待修复区域的轮廓线越近的已知像素点对点 p 的贡献越大。
二、下面就是考虑是什么样的顺序来处理待修复区域中的所有像素。作者采用的是快速行进方法Fast Marching Method。
行进算法伪代码：
|
1
2
3
4
5
6
7
8 |
δΩi = boundary of region to inpaint//修复区域的边缘 δΩ = δΩi while (δΩ not empty) { p = pixel of δΩ closest to δΩi//修复距离边缘最近的像素 inpaint p using Eqn.2//利用公式2修复p点 advance δΩ into Ω//把边缘向里行进 } |
看到这里，就会有疑惑了，怎么确定像素与边缘的距离呢。
算法中为待修复区域边缘构建了一个窄边（narrowBand），就是上面所说的δΩ。在opencv里是利用先将mask膨胀得到mask2（结构元素是长为2*ε+1的十字形，以中心点为原点），再用mask2减去mask得到band图，则band中非0元素即narrowBand）。从这里可以看出最初的narrowBand（即δΩ1）是不需要修复的。确定窄边的目的就是为了找到下面要修复的像素。
首先将像素分为三类，用flag标识记录：BAND：其实就是δΩ上的像素； KNOWN:就是δΩ外部不需要修复的像素；INSIDE:就是δΩ内部的等待修复的像素。
另外，每个像素还需要存储两个值：T（该像素离到边缘 δΩ的距离）；I（灰度值）。
下面先说一下处理像素是按怎样的行进方式的：
1． 初始化。首先按上面说的方法找到narrowBand，flag记为BAND;窄边内部的待修复区域记为INSIDE，已知像素flag设为KNOWN。BAND和KNOWN类型的像素T值初始化为0（这里看opencv代码里好像把KNOWN也设为106），INSIDE类型像素T值设为无限大（实际中设为106）。
2． 定义一个数据结构NarrowBand（opencv中采用双向链表实现），将窄边中的像素按T值升序排列，依次加入到NarrowBand中，先处理T最小的像素。假设为p点，将p点类型改为KNOWN,然后依次处理p点的四邻域点Pi。如果Pi类型为INSIDE，若是则重新计算I，修复该点，并更新其T值,修改该点类型为BAND,加入NarrowBand（这里仍按顺序，即始终保持NarrowBand是按升序排列的）。依次进行，每次处理的都是NarrowBand中T最小的像素，直到NarrowBand中没有像素。
代码如下：
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20 |
while (NarrowBand not empty) { extract P(i,j) = head(NarrowBand); /* STEP 1 *//*提取T值最小像素*/ for (k,l) in (i-1,j),(i,j-1),(i+1,j),(i,j+1)/*依次处理该像素的四邻域像素*/ { if (f(k,l)!=KNOWN) { if (f(k,l)==INSIDE) { f(k,l)=BAND; /* STEP 2修改(k,l)像素点的lag*/ inpaint(k,l); /* STEP 3修复(k,l)像素点*/ } T (k,l) = min(solve(k-1,l,k,l-1), /* STEP 4 更新(k,l)像素点的值*/ solve(k+1,l,k,l-1), solve(k-1,l,k,l+1), solve(k+1,l,k,l+1)); insert(k,l) in NarrowBand; /* STEP 5 将(k,l)像素点加入NarrowBand*/ } } } |
下面是具体计算T的代码个地方，原文中的代码有错误，下面是参考opencv代码修改好的）
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 |
T(k,l)=min(solve(k-1,l,k,l-1),solve(k+1,l,k,l-1),solve(k-1,l,k,l+1),solve(k+1,l,k,l+1)); 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(i2,j2)); return sol; } |
贴一个原文中的结果图：
