opencv图像修复算法cvInpaint(Telea的FMM算法）

cqh

在opencv中实现修复有两种算法，这里只介绍Telea的算法，即基于快速行进（FMM）的修复算法。

首先看c++接口中，函数的定义。

1. cv::inpaint( const Mat& src, const Mat& mask, Mat& dst,  double inpaintRange, int flags )
   
2. //src:要修复的图像；
   
3. //mask：修复模板，必须是单通道图像；
   
4. //dst：目标图像；
   
5. //inpaintRange：选取邻域半径;
   
6. //flags:要使用的方法，可以是CV INPAINT NS或CV INPAINT TELEA（本文介绍的方法）。

复制代码

其实c++接口实现的inpaint方法，只是调用了一下c接口中的cvinpaint。

1. cvInpaint( const vArr*_input_img,const CvArr* _inpaint_mask,CvArr* _output_img, double inpaintRange, int flags )

复制代码

cqh

最近项目中要用到修复技术，看opencv里有两种算法，我先尝试了一下Telea在2004年提出的基于快速行进的修复算法（后面简称FMM算法）。找到作者的原文看了一下，对算法有了一定了解,记录一下。

论文题目：An Image Inpainting Technique Based on the Fast Marching Method （2004）

作者主页：http://www.cs.rug.nl/~alext/

论文下载： http://www.cs.rug.nl/~alext/PAPERS/index.html  （编号36的那篇）
在opencv中实现修复有两种算法，这里只介绍Telea的算法，即基于快速行进（FMM）的修复算法。

cqh

首先，FMM算法基于的思想是，先处理待修复区域边缘上的像素点，然后层层向内推进，直到修复完所有的像素点。

下面以灰度图为例，我们只需要计算出像素新的灰度值即可。对于彩色图像，分别用同样的方法处理各个通道即可。

一、先说一下如何修复一个像素点的。

[/url]

参考上图，Ω区域是待修复的区域；δΩ指Ω的边界）；要修复Ω中的像素，就需要计算出新的像素值来代替原值。

现在假设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）：

[url=http://upload.gbs-cqh.net/2012/04/2.png]

这里的w(p, q)就是权值函数，是用来限定邻域中各像素的贡献大小的。

w(p, q) = dir(p, q) · dst(p, q) · lev(p, q)

[url=http://upload.gbs-cqh.net/2012/04/3.png][/url]

其中，d0和 T0分别为距离参数和水平集参数，一般都取为 1。方向因子 dir(p,q)保证了越靠近法线方向 N =  ∇T的像素点对 p 点的贡献最大；几何距离因子 dst(p,q)保证了离 p 点越近的像素点对p 点贡献越大；水平集距离因子lev(p,q)保证了离经过点 p 的待修复区域的轮廓线越近的已知像素点对点 p 的贡献越大。

cqh

二、下面就是考虑是什么样的顺序来处理待修复区域中的所有像素。作者采用的是快速行进方法Fast Marching Method。

行进算法伪代码：

1. δΩi = boundary of region to inpaint//修复区域的边缘
   
2. δΩ = δΩi
   
3. while (δΩ not empty)
   
4. {
   
5.     p = pixel of δΩ closest to δΩi//修复距离边缘最近的像素
   
6.     inpaint p using Eqn.2//利用公式2修复p点
   
7.     advance δΩ into Ω//把边缘向里行进
   
8. }

复制代码

看到这里，就会有疑惑了，怎么确定像素与边缘的距离呢。

算法中为待修复区域边缘构建了一个窄边（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. while (NarrowBand not empty)
   
2. {
   
3.     extract P(i,j) = head(NarrowBand); /* STEP 1 *//*提取T值最小像素*/
   
4.     for (k,l) in (i-1,j),(i,j-1),(i+1,j),(i,j+1)/*依次处理该像素的四邻域像素*/
   
5.     {
   
6.         if (f(k,l)!=KNOWN)
   
7.         {
   
8.             if (f(k,l)==INSIDE)
   
9.             {
   
10.                 f(k,l)=BAND; /* STEP 2修改(k,l)像素点的lag*/
    
11.                 inpaint(k,l); /* STEP 3修复(k,l)像素点*/
    
12.             }
    
13.             T (k,l) = min(solve(k-1,l,k,l-1), /* STEP 4 更新(k,l)像素点的值*/
    
14.             solve(k+1,l,k,l-1),
    
15.             solve(k-1,l,k,l+1),
    
16.             solve(k+1,l,k,l+1));
    
17.             insert(k,l) in NarrowBand; /* STEP 5 将(k,l)像素点加入NarrowBand*/
    
18.         }
    
19.     }
    
20. }

复制代码

下面是具体计算T的代码个地方，原文中的代码有错误，下面是参考opencv代码修改好的）

1. 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));
   
2. float solve(int i1,int j1,int i2,int j2)
   
3. {
   
4.     float sol = 1.0e6;
   
5.     if (f(i1,j1)==KNOWN)
   
6.         if (f(i2,j2)==KNOWN)
   
7.         {
   
8.             float r = sqrt(2-(T(i1,j1)-T(i2,j2))*(T(i1,j1)-T(i2,j2)));
   
9.             float s = (T(i1,j1)+T(i2,j2)-r)/2;
   
10.             if (s>=T(i1,j1) && s>=T(i2,j2)) sol = s;
    
11.             else
    
12.             { s += r; if (s>=T(i1,j1) && s>=T(i2,j2)) sol = s; }
    
13.         }
    
14.     else sol = 1+T(i1,j1));
    
15.     else if (f(i2,j2)==KNOWN) sol = 1+T(i2,j2));
    
16.     return sol;
    
17. }

复制代码

cqh

贴一个原文中的结果图：

[url=http://upload.gbs-cqh.net/2012/04/result.jpg][/url]