爱乐眼底图像分析

?找回密码
?立即注册

QQ登录

只需一步,快速开始

搜索
查看: 9069|回复: 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算法基于的思想是,先处理待修复区域边缘上的像素点,然后层层向内推进,直到修复完所有的像素点。
下面以灰度图为例,我们只需要计算出像素新的灰度值即可。对于彩色图像,分别用同样的方法处理各个通道即可。
一、先说一下如何修复一个像素点的。
参考上图,Ω区域是待修复的区域;δΩ指Ω的边界);要修复Ω中的像素,就需要计算出新的像素值来代替原值。
现在假设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 的贡献越大。

本帖子中包含更多资源

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

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. }
复制代码
5#
?楼主| 发表于 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-2 19:08 , Processed in 0.346465 second(s), 27 queries .

Powered by Discuz! X3.1

? 2001-2013 Comsenz Inc.

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