设为首页收藏本站

爱乐眼底图像分析

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 9169|回复: 5

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

[复制链接]
发表于 2012-6-3 19:21:14 | 显示全部楼层 |阅读模式
在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 )
复制代码
回复

使用道具 举报

 楼主| 发表于 2012-6-3 19:49:07 | 显示全部楼层
最近项目中要用到修复技术,看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)的修复算法。
回复 支持 反对

使用道具 举报

 楼主| 发表于 2012-6-3 19:49:37 | 显示全部楼层
首先,FMM算法基于的思想是,先处理待修复区域边缘上的像素点,然后层层向内推进,直到修复完所有的像素点。
下面以灰度图为例,我们只需要计算出像素新的灰度值即可。对于彩色图像,分别用同样的方法处理各个通道即可。
一、先说一下如何修复一个像素点的。
这里的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 的贡献越大。

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

 楼主| 发表于 2012-6-3 19:51:55 | 显示全部楼层
二、下面就是考虑是什么样的顺序来处理待修复区域中的所有像素。作者采用的是快速行进方法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. }
复制代码
回复 支持 反对

使用道具 举报

 楼主| 发表于 2012-6-3 19:52:23 | 显示全部楼层
贴一个原文中的结果图:
[url=http://upload.gbs-cqh.net/2012/04/result.jpg][/url]

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则



QQ|Archiver|手机版|小黑屋|爱乐眼底图像分析 ( 京ICP备1201155号          

GMT+8, 2019-10-22 04:12 , Processed in 0.318731 second(s), 30 queries .

Powered by Discuz! X3.1

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表