MSER算法介绍

Posted on 2013-09-23 11:50 Frisch' Blog 阅读(20652) 评论(5) 收藏举报

MSER代码编译:

matlabroot

%如果是VS2010则解压VS2010MEX支持文件到MATLAB根目录

unzip('E:\Software\develop Tools\VS2010MEXSupport.zip',matlabroot)

mex -setup

%设置代码文件夹编译路径

cd('E:\Koder Quelle\Image process Package\mser-0.5')

mser_compile

%生成PDB文件供调试用

mex -g mser.mex.c -output 'mser.mex.pdb'

%编译MEX文件

mex -g mser.mex.c -output 'mser'

mex -g erfill.mex.c -output 'erfill'

MSER的原理介绍:

图像I为一个映射 $I:D \subset {\mathbb{Z}^2} \to S$ . 极值区域在图像上的定义为:
1. S是全序的，例如存在满足自反性，反对称性和传递性的二元关系“小于等于”。论文里考虑到的范围仅是S={0,1,2,...,255},但是极值区域可以被定义在实值图像上(S=R)。

2.邻域关系的定义为 $A \subset D \times D$ （备注，这里集合的乘积为笛卡尔积），这里只用到了4邻域，例如 $p,q \in D$ 是相邻的（用 $pAq$ 表示），那么当且仅当满足 $\sum\nolimits_{i = 1}^d {|{p_i} - {q_i}| \le 1}$ 的时候它们才是邻域关系。

区域Q是D的连续子集，例如对每个 $p,q \in Q$ 这里存在一个序列 $p,{a_1},{a_2},...,{a_n},q$ 并且 $pA{a_1},...,{a_i}A{a_{i + 1}},{a_n}Aq$

这个性质表示区域Q为单连通的。

区域的边界(外部的) $\partial Q = \{ q \in D\backslash Q:\exists p \in Q:pAq\}$ ，例如Q的边界 $\partial Q$ 是这样一个像素集合，它至少有一个像素是属于Q的邻域像素但是不属于区域Q。

极值区域 $q \subset D$ 是这样的一个区域:对于所有 $p \in Q,q \in \partial Q:I(p) > I(q)$ [最大亮度区域]或者 $I(p) < I(q)$ [最小亮度区域]。

最大稳定极值区域（MSER）。令 ${Q_1},...,{Q_{i - 1}},{Q_i},...$ 为嵌套的极值区域，例如 ${Q_{i - 1}} \subset {Q_i}$ ，当且仅当 ${q_i} = |{Q_{i + \Delta }}\backslash {Q_{i - \Delta }}|/|{Q_i}|$ 在 $i*$ 时存在局部最小值，极值区域 ${Q_{i*}}$ 才是最大稳定的（|.|表示集合的势,对于离散图像来说就是像素面积）， $\Delta \in S$ 属于方法参数

Component Tree原理

当前元素和所有周围邻近当前区域的元素并且小于当前元素的区域都属于该区域

从最小值Xi像素开始建立一个region，如果Xi的邻域亮度值等于该像素的亮度值，则将Xi合并进邻域像素所在的region（且父节点是邻域像素父节点），如果Xi邻域的像素值小于Xi的像素值，则将邻域像素归入当前region，并将当前region作为邻域像素region的根节点（且邻域像素父节点指向当前像素的父节点），直到遍历完所有像素为止。

Component Tree 代码解读:

4    4    4    4    4    4    4
4    2    2    4    2    1    4
4    2    2    4    1    2    4
4    3    3    4    3    3    4
4    2    2    4    4    3    4
4    1    1    2    4    3    4
4    4    4    4    4    4    4

Extremal regions: 8

ner: 0. value:1,area: 2, parent: 3
ner: 1. value:1,area: 2, parent: 4
ner: 2. value:2,area: 4, parent: 5
ner: 3. value:2,area: 5, parent: 5
ner: 4. value:2,area: 4, parent: 6
ner: 5. value:3,area: 11, parent: 7
ner: 6. value:3,area: 8, parent: 7
ner: 7. value:4,area: 49, parent: 7

从像素值是1的像素开始:

节点0，1属于节点3

节点2，3属于节点5

节点4属于节点6

节点5，6属于节点7

反之则有：

1 1 1 1 1 1 1
1 3 3 2 3 4 1
1 3 3 2 3 4 1
1 1 1 1 1 3 1
1 3 4 2 1 1 1
1 4 3 2 2 2 1
1 1 1 1 1 1 1

ner: 0. value:1,area: 30, parent: 1
ner: 1. value:2,area: 36, parent: 2
ner: 2. value:3,area: 45, parent: 3
ner: 3. value:4,area: 49, parent: 3

从像素值是1的像素开始，刚开始只有节点0一个区域，将像素2的区域并入0后形成节点1，将像素3并入节点1后形成节点2，最后全部合并形成整个区域作为父节点.

Component Tree的步骤:

1.在计算开始之前首先对所有像素按照亮度大小进行排序。

for (i = 0; i < pixel_numbers; i++)
{
     1.弹出当前像素的索引和值
     2.将邻结点的索引指向当前索引
     3.将当前像素结点插入到树中
     4.while(依次处理以该像素为中心的8邻域像素)
     {
          a.邻域像素在图像边界以内
          b.邻域像素不是自身
          c.邻域像素已经在树中了，也就是说这个像素已经被处理过了
          if(满足a,b,c)
          {
               1.将当前此节点，以及当前此节点的上一层父节点,..., 直到根节点的所有链路节点的shortcut都置为根节点
               2.将此邻域节点，以及此邻域节点的上一层父节点,..., 直到根节点的所有链路节点的shortcut都置为根节点

               i. 如果当前节点的根节点等于邻域像素的根节点，这种情况下两个树已经加入了同一个树当中了
               ii. 如果当前节点的根节点的亮度等于邻域像素的根节点的亮度，此种情况索引将被扩展到具有相同亮度的一个极值区域，因为邻域像素的根节点不会是整个图像的一个极值区域；当前索引的根节点能够安全的被添加到邻域像素的根节点作为邻域像素的孩子结点。
               iii.如果当前节点的根节点的亮度大于邻域像素的根节点的亮度并且当前节点是一个极值区域，这时仅仅增加它的层数。这种情况下邻域像素根节点将会是最终整个图像的一个极值区域并且只有唯一的可能性作为添加到当前像素所属的树中当前像素的孩子

          }
     }
}

MSER区域的椭圆拟合(注意！这里是二值图像):

X,X correlation[X,X]

X,Y correlation[X,Y]

Y,Y correlation[Y,Y]

$E(X{Y^T}) = \frac{1}{{|R|}}\sum\limits_{x \in R} {X{X^T}}$

$E(X) = \frac{1}{{|R|}}\sum\limits_{x,y \in R} X$

$E(x) = \frac{1}{{|R|}}\sum\limits_{x \in R} x$

$E(y) = \frac{1}{{|R|}}\sum\limits_{y \in R} y$

$E({x^2}) = \frac{1}{{|R|}}\sum\limits_{x \in R} {{x^2}}$

$E({y^2}) = \frac{1}{{|R|}}\sum\limits_{x \in R} {{y^2}}$

$E(xy) = \frac{1}{{|R|}}\sum\limits_{x,y \in R} {xy}$

其中E(x),E(y)分别为拟合椭圆的中心点坐标，D(x),D(y),COV(x,y)分别为该区域内所有点的横坐标方差，纵坐标方差和横纵坐标协方差

椭圆拟合的代码:

 1   /* -----------------------------------------------------------------
 2    *                                                      Fit ellipses
 3    * -------------------------------------------------------------- */
 4   ell_pt = 0 ;
 5   if (nout >= 1) {
 6     int midx = 1 ;
 7     int d, index, j ;
 8     
 9     verbose && mexPrintf("Fitting ellipses...\n") ;
10 
11     /* enumerate maxstable regions */
12     for(i = 0 ; i < ner ; ++i) {      
13       if(! regions_pt [i].maxstable) continue ;
14       regions_pt [i].maxstable = midx++ ;
15     }
16     
17     /* allocate space */
18     acc_pt = mxMalloc(sizeof(acc_t) * nel) ;
19     ell_pt = mxMalloc(sizeof(acc_t) * gdl * nmer) ; 
20     
21     /* clear accumulators */
22     memset(ell_pt, 0, sizeof(int) * gdl * nmer) ;
23     for(index = 0 ; index < nel ; ++ index){
24           mexPrintf("%d\t",joins_pt[index]) ;
25           if( 0 == (index+1)%dims[0]) mexPrintf("\n");
26       }
27       
28      for(index = 0 ; index < nel ; ++ index){
29           mexPrintf("%d\t",forest_pt [ index ].parent) ;
30           if( 0 == (index+1)%dims[0]) mexPrintf("\n");
31       }
32       
33     /* for each gdl */
34     for(d = 0 ; d < gdl ; ++d) {
35       /* initalize parameter */
36       memset(subs_pt, 0, sizeof(int) * ndims) ;
37       
38       //注意这里的0和1下标不是0阶矩和1阶矩，而是表示X和Y之间的均值，相关系数
39       //例如(0,0)表示corr(X,X), (0,1)表示corr(X,Y),(1,0)表示corr(Y,X),(1,1)表示corr(Y,Y)不要搞混！！！ by frisch
40       if(d < ndims) {
41         //图像每个像素的行[0], 列坐标[1]坐标, x方向和y方向的均值，by frisch
42         verbose && mexPrintf(" mean %d\n",d) ;
43         for(index = 0 ; index < nel ; ++ index) {
44           acc_pt[index] = subs_pt[d] ;
45           adv(dims, ndims, subs_pt) ;      
46           
47         }
48 
49       } else {
50 
51       //图像每个像素的行[0], 列坐标[1]的相关度计算, by frisch
52       //corr[0,0](行与行), corr[1,0](行与列), corr[0,1](列与行), corr[1,1](列与列), by frisch
53         /* decode d-ndims into a (i,j) pair */
54         i = d-ndims ; 
55         j = 0 ;
56         while(i > j) {
57           i -= j + 1 ;
58           j ++ ;
59         }
60 
61         verbose && mexPrintf(" corr (%d,%d)\n",i,j) ;
62 
63         /* add x_i * x_j */
64         //计算相关值，将其放入二维数组
65         for(index = 0 ; index < nel ; ++ index){
66           acc_pt[index] = subs_pt[i]*subs_pt[j] ;
67           mexPrintf("(%d\t%d\t)",subs_pt[i],subs_pt[j]) ;
68           if( 0 == (index+1)%dims[0]) mexPrintf("\n");
69           adv(dims, ndims, subs_pt) ;      
70         }
71       }
72 
73       for(index = 0 ; index < nel ; ++ index){
74            mexPrintf("%d\t",acc_pt[index]) ;
75           if( 0 == (index+1)%dims[0]) mexPrintf("\n");
76       }
77 
78       /* integrate parameter */
79       //对每个区域进行均值或者相关值积分,by frisch
80       for(i = 0 ; i < njoins ; ++i) {      
81         idx_t index  = joins_pt[i] ;
82         idx_t parent = forest_pt [ index ].parent ;
83         acc_pt[parent] += acc_pt[index] ;
84       }
85     
86       /* save back to ellpises */
87       //保存椭圆参数, 五个参数(包含椭圆的中心点，长短轴幅值，方向向量) by frisch
88       for(i = 0 ; i < ner ; ++i) {
89         idx_t region = regions_pt [i].maxstable ;
90 
91         /* skip if not extremal region */
92         if(region-- == 0) continue ;
93         ell_pt [d + gdl*region] = acc_pt [ regions_pt[i].index ] ;
94       }
95 
96       /* next gdl */
97     }
98     mxFree(acc_pt) ;
99   }

reference:

椭圆拟合部分参考《多通道图像MSER局部不变特征》国防科技大学，柳涛

刷新页面返回顶部

frischzenger

公告

MSER算法介绍