mser函数包 matlab,Maximally Stable Extremal Regions (matlab code)

说明： 不太熟悉简书的模式，图表也没时间整理了。请看原网址 github.com/acoder-fin/…

谢谢！

[1] Robust wide baseline stereo from maximally stable extremal regions

[2] Linear time maximally stable extremal regions

MSER 算法有J. Matas等人于2002年发表在BMVC上。对我个人而言，这个算法很容易就可以理解，但是却不知道如何去实现。

当然现在已经可以实现了。

这篇文章提出一种叫extremal regions的新的元素集合，其实质也就是一种image blob。

该特征的优点可以从原文中看到，不再赘述。

1 经典回顾 [1]

Reliable extraction of a manageable number of potentially corresponding image elements is necessary but certainly not a sufficient preerquisite for successful wide-baseline matching.

In watershed computation, focus is on the thresholds where regions merge (and two watersheds touch). Such threshold are of little interest here, since they are highly unstable - after merge, the region area jumps.

2 算法实现 MSER pure matlab code

在获得图像最小值和最大值后，就可以不断地把阈值从最小值升到最大值。每次阈值增加1。

不同的阈值下，可以得到不同的联通的区域，这些区域就是extremal regions。根据作者的说法，这些区域不一是最稳定的，即maximally stable。

然后作者就提出了一种度量方法，用以度量一个区域的稳定度。那怎么度量呢？作者是根据一个区域的面积的变化情况来度量。

在MSER算法中，区域只有两种可能性，即增长或者与其他区域合并。要知道一个区域面积的变化情况，就需要记录区域增长或与其他区域合并的情况。

这样，就可以得到区域面积的变化情况。

Minimum intensity region

For threshold t:

bw = ones(size(im)); # which is fast than *bw = 0*im+1;*

bw(im>=t) = 0; # be careful that it is '>='

复制代码

以lena图像为例：

[图片上传失败...(image-14d028-1564134107562)]

设阈值为200，则二值化以后的图像为

[图片上传失败...(image-73f032-1564134107562)]

2.2 得到区域

在Matlab中，有一个‘bwconncomp’的函数，通过这个函数，我们就可以得到在某个阈值下，图像的区域个数及其位置。

cc = bwconncomp(bw,4); # here, 4 means the neighbors we need to consider, usually it is 4 or 8 for 2D images

复制代码

2.3 记录区域关系

要想记录区域关系，就得知道最开始的区域，然后才可以考虑该区域的增长和与其他区域合并情况。

因此，这里对一个突然出现的区域和 在突然出现的区域的基础上增长而来(只和一个区域有交集)或 由多个突然出现的区域合并而来(和多个区域有交集) 的区域进行区分。

把区域面积的变化情况，全部记录在那些突然出现的区域里。同时记录不同面积对应的区域编号。

seed_list = []; # preserve the seed for each 'pop out' region

# if a 'pop out' region occurs, we need to append a seed to seed_list

# Just choose the first one as the seed, it does not matter

# if a new region has a intersection with the seed_list, it means that

# this region is not a 'pop out' region and it must from one 'pop out' region or from merging by

# more than one 'pop out' regions.

idx = cc.PixelIdxList;

if isempty(intersect(idx,seed_list))

复制代码

2.4 判断区域是否Maximally stable

penelope 对判断区域是否稳定的条件有详细地说明。

主要有6个阈值，一个delta用来计算区域的变化，然后得到面积变化为局部最小的区域，可能为maximally stable。

max_variation也是用来限定面积变化，即面积变化值不能超过该值。

区域的面积也不可以过大或过小。

min_diversity 用来去除太过相似的区域。

2.5 结果

用上边方法得到的结果为(coins.png)

[图片上传失败...(image-a0cc1b-1564134107563)]

Matlab函数detectMSERFeatures检测结果

[图片上传失败...(image-870bc4-1564134107563)]

可以用其他图像进行测试，也可以去除那段去重复的程序(效果一般)，则可以得到更多的区域。

2.6 上边方法的优缺点

对图像灰度值范围没有限制(0-255,0-65535,etc.)

添加了一段去重复的程序，在一定程度上去除了一些非常临近的区域

耗时长！

2.7 完整代码

谢谢阅读！

2.8 补充1-突然出现的区域 / pop out regions

假设有一幅4*4的图像，如下图所示：

51

51

51

51

51

50

50

51

51

50

50

51

51

51

51

51

当图像阈值为51时，其二值化图像为：

0

0

0

0

0

1

1

0

0

1

1

0

0

0

0

0

其中，四个为1的像素点构成一个区域，该区域是突然出现的，或pop out region。将第一个像素的索引，即6，记入seed_list中。当阈值为52时，图像二值化，得到一个4*4的区域。

该区域与seed_list求交集，非空。说明该区域不是突然出现的，因为是由其他区域增长/合并而来的。面积记录为[4,16]。合并的情况类似。

2.8 补充2 去重问题

去重思路：

(1) 根据得到区域的椭圆拟合的中心，得到不同中心间的距离；

(2) 根据不同椭圆的短轴，判断另一个区域是否在该区域内；

(3) 可直接去除，选择面积最大的，也可以结合图像像素值大小进行去除。

2.8 补充3 耗时问题

要不断地改变阈值，若图像的灰度阶不小，则非常耗时。该情况是否可避免呢？Linear time maximally stable extremal regions

很好地回答了这个问题。需要说明的是，耗时长是一个主观感受，所用电脑的不同，编程软件的不同等都会造成同一个算法的运算时间不一样。因此，需要考虑计算复杂度来分析算法的效率。

该文发表于2008年ECCV。相对于MSER要求阈值不断上升，该文提出了一种更符合物理现实的实现方法。即水平面不是同时上涨的，而是在一个区域填满后，水溢出到其他区域。

LTMSER不需要不断地将图像二值化，也不需要union find 来获取二值图像里面的区域信息。

该算法先找到一个区域的最小值，然后不断地抬高水平面，直到该区域的最大值。这里的最小值和最大值指的都是灰度值。

2.9.1 三种数据结构在Matlab中如何实现 [2]

A binary mask of accessible pixels: accessible = zeros(size(im))

A priority quere of boundary pixels, where priority is minus the grey-level: boundary = cell(priority,1)

A stack C of component information: stack = struct();

2.9.2 具体实现 step by step [2]

clear the accessible pixel mask, the heap of boundary pixels and the component stack. Push a dummy-component on the stack,

with grey-level higher than any allowed in the image.

step 1

accessible = zeros(size(im)); # faster than 0*im

boundary = cell(priority,1);

regions = struct(); # this one is used to record the regions

复制代码Make the source pixel (with its first edge) the current pixel, mark it as accessible and store

the gray level of it in the variable current_level

step 2

current_index = 2; # start from the second pixel

current_level = im(current_index);

accessible(current_index) = 1;

复制代码Push an empty component with current_level onto the component stack

step 3

stack = initialize_regions(current_level,[]);

复制代码Explore the remaining edges to the neighbors of the current pixel, in order,

as follows: For each neighbor, check if the neighbor is already accessible. If it

is not, mark it as accessible and retrieve its grey-level. If the grey-level is not

lower than the current one, push it onto the heap of boundary pixels. If on

the other hand the grey-level is lower than the current one, enter the current

pixel back into the queue of boundary pixels for later processing (with the

next edge number), consider the new pixel and its grey-level and go to 3.

step 4

在Matlab中，没有goto的功能，那么如何实现go to step 3呢？在请教了朋友之后，才发现可以用循环来实现。

这里设置一个flag，如何找到一个比current_level小的像素，则改变flag的值，然后再回到step 3。

neighbor_index = current_index + neighbor_template;

neighbor_index(neighbor_index<1|neighbor_index>rows*cols) = [];

neighbor_length = length(neighbor_index);

for i = 1 : neighbor_length

flag = 0;

if accessible(neighbor_index(i))~=1

accessible(neighbor_index(i)) = 1;

gray_level = im(neighbor_index(i));

if gray_level >= current_level

boundary{gray_level} = [boundary{gray_level},neighbor_index(i)];

if gray_level < priority

priority = gray_level;

end

else

boundary{current_level} = [boundary{current_level},current_index];

if current_level < priority

priority = current_level;

end

current_index = neighbor_index(i);

current_level = gray_level;

flag = 1;

break;

end

end

end

if flag == 1

break;

end

复制代码Accumulate the current pixel to the component at the top of the stack (water saturates the current pixel).

step 5

add_pixel_to_region();

复制代码Pop the heap of the boundary pixels. If the heap is empty, we are done. If the returned pixel is at the same

grey-level as the previous, go to 4.

The returned pixel is at a higher grey-level, so we must now process all

components on the component stack until we reach the higher grey-level.

This is done with the ProcessStack sub-routine, see below. Then go to 4.

step 6

寻找可能属于当前区域的其他像素点。

if priority == max(im(:)) + 1

stable_regions = is_stable(delta,min_area,max_area,max_variation,min_diversity);

return;

end

current_index = boundary{priority}(end);

boundary{priority}(end) = [];

while isempty(boundary{priority})&&priority

priority = priority + 1;

end

next_level = im(current_index);

% if current region level is not equal to the boundary, which means

% we should make the region level higher (region merge not accumulate or add pixel to region)

if next_level ~= current_level

current_level = next_level;

% step 7 合并区域

process_stack();

end

复制代码

2.9.3 完整代码

code

谢谢阅读！