1 /* Redistribution and use in source and binary forms, with or
2 * without modification, are permitted provided that the following
4 * Redistributions of source code must retain the above
5 * copyright notice, this list of conditions and the following
7 * Redistributions in binary form must reproduce the above
8 * copyright notice, this list of conditions and the following
9 * disclaimer in the documentation and/or other materials
10 * provided with the distribution.
11 * The name of Contributor may not be used to endorse or
12 * promote products derived from this software without
13 * specific prior written permission.
15 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
16 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
17 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
18 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19 * DISCLAIMED. IN NO EVENT SHALL THE CONTRIBUTORS BE LIABLE FOR ANY
20 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
21 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
22 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
23 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
24 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
25 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
26 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
28 * Copyright© 2009, Liu Liu All rights reserved.
30 * OpenCV functions for MSER extraction
32 * 1. there are two different implementation of MSER, one for grey image, one for color image
33 * 2. the grey image algorithm is taken from: Linear Time Maximally Stable Extremal Regions;
34 * the paper claims to be faster than union-find method;
35 * it actually get 1.5~2m/s on my centrino L7200 1.2GHz laptop.
36 * 3. the color image algorithm is taken from: Maximally Stable Colour Regions for Recognition and Match;
37 * it should be much slower than grey image method ( 3~4 times );
38 * the chi_table.h file is taken directly from paper's source code which is distributed under GPL.
39 * 4. though the name is *contours*, the result actually is a list of point set.
42 #include "precomp.hpp"
43 #include "opencv2/imgproc/imgproc_c.h"
49 const int TABLE_SIZE = 400;
51 static double chitab3[]={0, 0.0150057, 0.0239478, 0.0315227,
52 0.0383427, 0.0446605, 0.0506115, 0.0562786,
53 0.0617174, 0.0669672, 0.0720573, 0.0770099,
54 0.081843, 0.0865705, 0.0912043, 0.0957541,
55 0.100228, 0.104633, 0.108976, 0.113261,
56 0.117493, 0.121676, 0.125814, 0.12991,
57 0.133967, 0.137987, 0.141974, 0.145929,
58 0.149853, 0.15375, 0.15762, 0.161466,
59 0.165287, 0.169087, 0.172866, 0.176625,
60 0.180365, 0.184088, 0.187794, 0.191483,
61 0.195158, 0.198819, 0.202466, 0.2061,
62 0.209722, 0.213332, 0.216932, 0.220521,
63 0.2241, 0.22767, 0.231231, 0.234783,
64 0.238328, 0.241865, 0.245395, 0.248918,
65 0.252435, 0.255947, 0.259452, 0.262952,
66 0.266448, 0.269939, 0.273425, 0.276908,
67 0.280386, 0.283862, 0.287334, 0.290803,
68 0.29427, 0.297734, 0.301197, 0.304657,
69 0.308115, 0.311573, 0.315028, 0.318483,
70 0.321937, 0.32539, 0.328843, 0.332296,
71 0.335749, 0.339201, 0.342654, 0.346108,
72 0.349562, 0.353017, 0.356473, 0.35993,
73 0.363389, 0.366849, 0.37031, 0.373774,
74 0.377239, 0.380706, 0.384176, 0.387648,
75 0.391123, 0.3946, 0.39808, 0.401563,
76 0.405049, 0.408539, 0.412032, 0.415528,
77 0.419028, 0.422531, 0.426039, 0.429551,
78 0.433066, 0.436586, 0.440111, 0.44364,
79 0.447173, 0.450712, 0.454255, 0.457803,
80 0.461356, 0.464915, 0.468479, 0.472049,
81 0.475624, 0.479205, 0.482792, 0.486384,
82 0.489983, 0.493588, 0.4972, 0.500818,
83 0.504442, 0.508073, 0.511711, 0.515356,
84 0.519008, 0.522667, 0.526334, 0.530008,
85 0.533689, 0.537378, 0.541075, 0.54478,
86 0.548492, 0.552213, 0.555942, 0.55968,
87 0.563425, 0.56718, 0.570943, 0.574715,
88 0.578497, 0.582287, 0.586086, 0.589895,
89 0.593713, 0.597541, 0.601379, 0.605227,
90 0.609084, 0.612952, 0.61683, 0.620718,
91 0.624617, 0.628526, 0.632447, 0.636378,
92 0.64032, 0.644274, 0.648239, 0.652215,
93 0.656203, 0.660203, 0.664215, 0.668238,
94 0.672274, 0.676323, 0.680384, 0.684457,
95 0.688543, 0.692643, 0.696755, 0.700881,
96 0.70502, 0.709172, 0.713339, 0.717519,
97 0.721714, 0.725922, 0.730145, 0.734383,
98 0.738636, 0.742903, 0.747185, 0.751483,
99 0.755796, 0.760125, 0.76447, 0.768831,
100 0.773208, 0.777601, 0.782011, 0.786438,
101 0.790882, 0.795343, 0.799821, 0.804318,
102 0.808831, 0.813363, 0.817913, 0.822482,
103 0.827069, 0.831676, 0.836301, 0.840946,
104 0.84561, 0.850295, 0.854999, 0.859724,
105 0.864469, 0.869235, 0.874022, 0.878831,
106 0.883661, 0.888513, 0.893387, 0.898284,
107 0.903204, 0.908146, 0.913112, 0.918101,
108 0.923114, 0.928152, 0.933214, 0.938301,
109 0.943413, 0.94855, 0.953713, 0.958903,
110 0.964119, 0.969361, 0.974631, 0.979929,
111 0.985254, 0.990608, 0.99599, 1.0014,
112 1.00684, 1.01231, 1.01781, 1.02335,
113 1.02891, 1.0345, 1.04013, 1.04579,
114 1.05148, 1.05721, 1.06296, 1.06876,
115 1.07459, 1.08045, 1.08635, 1.09228,
116 1.09826, 1.10427, 1.11032, 1.1164,
117 1.12253, 1.1287, 1.1349, 1.14115,
118 1.14744, 1.15377, 1.16015, 1.16656,
119 1.17303, 1.17954, 1.18609, 1.19269,
120 1.19934, 1.20603, 1.21278, 1.21958,
121 1.22642, 1.23332, 1.24027, 1.24727,
122 1.25433, 1.26144, 1.26861, 1.27584,
123 1.28312, 1.29047, 1.29787, 1.30534,
124 1.31287, 1.32046, 1.32812, 1.33585,
125 1.34364, 1.3515, 1.35943, 1.36744,
126 1.37551, 1.38367, 1.39189, 1.4002,
127 1.40859, 1.41705, 1.42561, 1.43424,
128 1.44296, 1.45177, 1.46068, 1.46967,
129 1.47876, 1.48795, 1.49723, 1.50662,
130 1.51611, 1.52571, 1.53541, 1.54523,
131 1.55517, 1.56522, 1.57539, 1.58568,
132 1.59611, 1.60666, 1.61735, 1.62817,
133 1.63914, 1.65025, 1.66152, 1.67293,
134 1.68451, 1.69625, 1.70815, 1.72023,
135 1.73249, 1.74494, 1.75757, 1.77041,
136 1.78344, 1.79669, 1.81016, 1.82385,
137 1.83777, 1.85194, 1.86635, 1.88103,
138 1.89598, 1.91121, 1.92674, 1.94257,
139 1.95871, 1.97519, 1.99201, 2.0092,
140 2.02676, 2.04471, 2.06309, 2.08189,
141 2.10115, 2.12089, 2.14114, 2.16192,
142 2.18326, 2.2052, 2.22777, 2.25101,
143 2.27496, 2.29966, 2.32518, 2.35156,
144 2.37886, 2.40717, 2.43655, 2.46709,
145 2.49889, 2.53206, 2.56673, 2.60305,
146 2.64117, 2.6813, 2.72367, 2.76854,
147 2.81623, 2.86714, 2.92173, 2.98059,
148 3.04446, 3.1143, 3.19135, 3.27731,
149 3.37455, 3.48653, 3.61862, 3.77982,
150 3.98692, 4.2776, 4.77167, 133.333 };
152 typedef struct LinkedPoint
154 struct LinkedPoint* prev;
155 struct LinkedPoint* next;
160 // the history of region grown
161 typedef struct MSERGrowHistory
163 struct MSERGrowHistory* shortcut;
164 struct MSERGrowHistory* child;
165 int stable; // when it ever stabled before, record the size
171 typedef struct MSERConnectedComp
175 MSERGrowHistory* history;
176 unsigned long grey_level;
178 int dvar; // the derivative of last var
179 float var; // the current variation (most time is the variation of one-step back)
183 // Linear Time MSER claims by using bsf can get performance gain, here is the implementation
184 // however it seems that will not do any good in real world test
185 inline void _bitset(unsigned long * a, unsigned long b)
189 inline void _bitreset(unsigned long * a, unsigned long b)
196 MSERParams( int _delta, int _minArea, int _maxArea, double _maxVariation,
197 double _minDiversity, int _maxEvolution, double _areaThreshold,
198 double _minMargin, int _edgeBlurSize )
199 : delta(_delta), minArea(_minArea), maxArea(_maxArea), maxVariation(_maxVariation),
200 minDiversity(_minDiversity), maxEvolution(_maxEvolution), areaThreshold(_areaThreshold),
201 minMargin(_minMargin), edgeBlurSize(_edgeBlurSize)
209 double areaThreshold;
214 // clear the connected component in stack
216 initMSERComp( MSERConnectedComp* comp )
221 comp->history = NULL;
224 // add history of size to a connected component
226 MSERNewHistory( MSERConnectedComp* comp, MSERGrowHistory* history )
228 history->child = history;
229 if ( NULL == comp->history )
231 history->shortcut = history;
234 comp->history->child = history;
235 history->shortcut = comp->history->shortcut;
236 history->stable = comp->history->stable;
238 history->val = comp->grey_level;
239 history->size = comp->size;
240 comp->history = history;
243 // merging two connected component
245 MSERMergeComp( MSERConnectedComp* comp1,
246 MSERConnectedComp* comp2,
247 MSERConnectedComp* comp,
248 MSERGrowHistory* history )
252 comp->grey_level = comp2->grey_level;
253 history->child = history;
254 // select the winner by size
255 if ( comp1->size >= comp2->size )
257 if ( NULL == comp1->history )
259 history->shortcut = history;
262 comp1->history->child = history;
263 history->shortcut = comp1->history->shortcut;
264 history->stable = comp1->history->stable;
266 if ( NULL != comp2->history && comp2->history->stable > history->stable )
267 history->stable = comp2->history->stable;
268 history->val = comp1->grey_level;
269 history->size = comp1->size;
270 // put comp1 to history
271 comp->var = comp1->var;
272 comp->dvar = comp1->dvar;
273 if ( comp1->size > 0 && comp2->size > 0 )
275 comp1->tail->next = comp2->head;
276 comp2->head->prev = comp1->tail;
278 head = ( comp1->size > 0 ) ? comp1->head : comp2->head;
279 tail = ( comp2->size > 0 ) ? comp2->tail : comp1->tail;
280 // always made the newly added in the last of the pixel list (comp1 ... comp2)
282 if ( NULL == comp2->history )
284 history->shortcut = history;
287 comp2->history->child = history;
288 history->shortcut = comp2->history->shortcut;
289 history->stable = comp2->history->stable;
291 if ( NULL != comp1->history && comp1->history->stable > history->stable )
292 history->stable = comp1->history->stable;
293 history->val = comp2->grey_level;
294 history->size = comp2->size;
295 // put comp2 to history
296 comp->var = comp2->var;
297 comp->dvar = comp2->dvar;
298 if ( comp1->size > 0 && comp2->size > 0 )
300 comp2->tail->next = comp1->head;
301 comp1->head->prev = comp2->tail;
303 head = ( comp2->size > 0 ) ? comp2->head : comp1->head;
304 tail = ( comp1->size > 0 ) ? comp1->tail : comp2->tail;
305 // always made the newly added in the last of the pixel list (comp2 ... comp1)
309 comp->history = history;
310 comp->size = comp1->size + comp2->size;
314 MSERVariationCalc( MSERConnectedComp* comp, int delta )
316 MSERGrowHistory* history = comp->history;
317 int val = comp->grey_level;
318 if ( NULL != history )
320 MSERGrowHistory* shortcut = history->shortcut;
321 while ( shortcut != shortcut->shortcut && shortcut->val + delta > val )
322 shortcut = shortcut->shortcut;
323 MSERGrowHistory* child = shortcut->child;
324 while ( child != child->child && child->val + delta <= val )
327 child = child->child;
329 // get the position of history where the shortcut->val <= delta+val and shortcut->child->val >= delta+val
330 history->shortcut = shortcut;
331 return (float)(comp->size-shortcut->size)/(float)shortcut->size;
332 // here is a small modification of MSER where cal ||R_{i}-R_{i-delta}||/||R_{i-delta}||
333 // in standard MSER, cal ||R_{i+delta}-R_{i-delta}||/||R_{i}||
334 // my calculation is simpler and much easier to implement
339 static bool MSERStableCheck( MSERConnectedComp* comp, MSERParams params )
341 // tricky part: it actually check the stablity of one-step back
342 if ( comp->history == NULL || comp->history->size <= params.minArea || comp->history->size >= params.maxArea )
344 float div = (float)(comp->history->size-comp->history->stable)/(float)comp->history->size;
345 float var = MSERVariationCalc( comp, params.delta );
346 int dvar = ( comp->var < var || (unsigned long)(comp->history->val + 1) < comp->grey_level );
347 int stable = ( dvar && !comp->dvar && comp->var < params.maxVariation && div > params.minDiversity );
351 comp->history->stable = comp->history->size;
355 // add a pixel to the pixel list
356 static void accumulateMSERComp( MSERConnectedComp* comp, LinkedPoint* point )
358 if ( comp->size > 0 )
360 point->prev = comp->tail;
361 comp->tail->next = point;
372 // convert the point set to CvSeq
373 static CvContour* MSERToContour( MSERConnectedComp* comp, CvMemStorage* storage )
375 CvSeq* _contour = cvCreateSeq( CV_SEQ_KIND_GENERIC|CV_32SC2, sizeof(CvContour), sizeof(CvPoint), storage );
376 CvContour* contour = (CvContour*)_contour;
377 cvSeqPushMulti( _contour, 0, comp->history->size );
378 LinkedPoint* lpt = comp->head;
379 for ( int i = 0; i < comp->history->size; i++ )
381 CvPoint* pt = CV_GET_SEQ_ELEM( CvPoint, _contour, i );
386 cvBoundingRect( contour );
390 // to preprocess src image to following format
392 // > 0 is available, < 0 is visited
393 // 17~19 bits is the direction
394 // 8~11 bits is the bucket it falls to (for BitScanForward)
395 // 0~8 bits is the color
396 static int* preprocessMSER_8UC1( CvMat* img,
401 int srccpt = src->step-src->cols;
402 int cpt_1 = img->cols-src->cols-1;
403 int* imgptr = img->data.i;
407 for ( int i = 0; i < 256; i++ )
410 for ( int i = 0; i < src->cols+2; i++ )
416 uchar* srcptr = src->data.ptr;
420 uchar* maskptr = mask->data.ptr;
421 for ( int i = 0; i < src->rows; i++ )
425 for ( int j = 0; j < src->cols; j++ )
431 *srcptr = 0xff-*srcptr;
432 level_size[*srcptr]++;
433 *imgptr = ((*srcptr>>5)<<8)|(*srcptr);
447 startptr = imgptr+img->cols+1;
448 for ( int i = 0; i < src->rows; i++ )
452 for ( int j = 0; j < src->cols; j++ )
454 *srcptr = 0xff-*srcptr;
455 level_size[*srcptr]++;
456 *imgptr = ((*srcptr>>5)<<8)|(*srcptr);
465 for ( int i = 0; i < src->cols+2; i++ )
472 for ( int i = 1; i < 256; i++ )
474 heap_cur[i] = heap_cur[i-1]+level_size[i-1]+1;
480 static void extractMSER_8UC1_Pass( int* ioptr,
484 MSERGrowHistory* histptr,
485 MSERConnectedComp* comptr,
492 CvMemStorage* storage )
494 comptr->grey_level = 256;
496 comptr->grey_level = (*imgptr)&0xff;
497 initMSERComp( comptr );
498 *imgptr |= 0x80000000;
499 heap_cur += (*imgptr)&0xff;
500 int dir[] = { 1, step, -1, -step };
501 #ifdef __INTRIN_ENABLED__
502 unsigned long heapbit[] = { 0, 0, 0, 0, 0, 0, 0, 0 };
503 unsigned long* bit_cur = heapbit+(((*imgptr)&0x700)>>8);
507 // take tour of all the 4 directions
508 while ( ((*imgptr)&0x70000) < 0x40000 )
511 int* imgptr_nbr = imgptr+dir[((*imgptr)&0x70000)>>16];
512 if ( *imgptr_nbr >= 0 ) // if the neighbor is not visited yet
514 *imgptr_nbr |= 0x80000000; // mark it as visited
515 if ( ((*imgptr_nbr)&0xff) < ((*imgptr)&0xff) )
517 // when the value of neighbor smaller than current
518 // push current to boundary heap and make the neighbor to be the current one
519 // create an empty comp
523 heap_cur += ((*imgptr_nbr)&0xff)-((*imgptr)&0xff);
524 #ifdef __INTRIN_ENABLED__
525 _bitset( bit_cur, (*imgptr)&0x1f );
526 bit_cur += (((*imgptr_nbr)&0x700)-((*imgptr)&0x700))>>8;
530 initMSERComp( comptr );
531 comptr->grey_level = (*imgptr)&0xff;
534 // otherwise, push the neighbor to boundary heap
535 heap_cur[((*imgptr_nbr)&0xff)-((*imgptr)&0xff)]++;
536 *heap_cur[((*imgptr_nbr)&0xff)-((*imgptr)&0xff)] = imgptr_nbr;
537 #ifdef __INTRIN_ENABLED__
538 _bitset( bit_cur+((((*imgptr_nbr)&0x700)-((*imgptr)&0x700))>>8), (*imgptr_nbr)&0x1f );
544 int imsk = (int)(imgptr-ioptr);
545 ptsptr->pt = cvPoint( imsk&stepmask, imsk>>stepgap );
546 // get the current location
547 accumulateMSERComp( comptr, ptsptr );
549 // get the next pixel from boundary heap
554 #ifdef __INTRIN_ENABLED__
556 _bitreset( bit_cur, (*imgptr)&0x1f );
559 #ifdef __INTRIN_ENABLED__
560 bool found_pixel = 0;
561 unsigned long pixel_val;
562 for ( int i = ((*imgptr)&0x700)>>8; i < 8; i++ )
564 if ( _BitScanForward( &pixel_val, *bit_cur ) )
568 heap_cur += pixel_val-((*imgptr)&0xff);
576 unsigned long pixel_val = 0;
577 for ( unsigned long i = ((*imgptr)&0xff)+1; i < 256; i++ )
591 #ifdef __INTRIN_ENABLED__
593 _bitreset( bit_cur, pixel_val&0x1f );
595 if ( pixel_val < comptr[-1].grey_level )
597 // check the stablity and push a new history, increase the grey level
598 if ( MSERStableCheck( comptr, params ) )
600 CvContour* contour = MSERToContour( comptr, storage );
601 contour->color = color;
602 cvSeqPush( contours, &contour );
604 MSERNewHistory( comptr, histptr );
605 comptr[0].grey_level = pixel_val;
608 // keep merging top two comp in stack until the grey level >= pixel_val
612 MSERMergeComp( comptr+1, comptr, comptr, histptr );
614 if ( pixel_val <= comptr[0].grey_level )
616 if ( pixel_val < comptr[-1].grey_level )
618 // check the stablity here otherwise it wouldn't be an ER
619 if ( MSERStableCheck( comptr, params ) )
621 CvContour* contour = MSERToContour( comptr, storage );
622 contour->color = color;
623 cvSeqPush( contours, &contour );
625 MSERNewHistory( comptr, histptr );
626 comptr[0].grey_level = pixel_val;
638 static void extractMSER_8UC1( CvMat* src,
641 CvMemStorage* storage,
646 while ( step < src->step+2 )
651 int stepmask = step-1;
653 // to speedup the process, make the width to be 2^N
654 CvMat* img = cvCreateMat( src->rows+2, step, CV_32SC1 );
655 int* ioptr = img->data.i+step+1;
658 // pre-allocate boundary heap
659 int** heap = (int**)cvAlloc( (src->rows*src->cols+256)*sizeof(heap[0]) );
660 int** heap_start[256];
661 heap_start[0] = heap;
663 // pre-allocate linked point and grow history
664 LinkedPoint* pts = (LinkedPoint*)cvAlloc( src->rows*src->cols*sizeof(pts[0]) );
665 MSERGrowHistory* history = (MSERGrowHistory*)cvAlloc( src->rows*src->cols*sizeof(history[0]) );
666 MSERConnectedComp comp[257];
668 // darker to brighter (MSER-)
669 imgptr = preprocessMSER_8UC1( img, heap_start, src, mask );
670 extractMSER_8UC1_Pass( ioptr, imgptr, heap_start, pts, history, comp, step, stepmask, stepgap, params, -1, contours, storage );
671 // brighter to darker (MSER+)
672 imgptr = preprocessMSER_8UC1( img, heap_start, src, mask );
673 extractMSER_8UC1_Pass( ioptr, imgptr, heap_start, pts, history, comp, step, stepmask, stepgap, params, 1, contours, storage );
679 cvReleaseMat( &img );
688 double m; // the margin used to prune area later
695 // to make the finding of root less painful
698 // a point double-linked list
700 // the temporary msr (set to NULL at every re-initialise)
702 // the global msr (once set, never to NULL)
704 // the index of the node, at this point, it should be x at the first 16-bits, and y at the last 16-bits.
719 static double ChiSquaredDistance( uchar* x, uchar* y )
721 return (double)((x[0]-y[0])*(x[0]-y[0]))/(double)(x[0]+y[0]+1e-10)+
722 (double)((x[1]-y[1])*(x[1]-y[1]))/(double)(x[1]+y[1]+1e-10)+
723 (double)((x[2]-y[2])*(x[2]-y[2]))/(double)(x[2]+y[2]+1e-10);
726 static void initMSCRNode( MSCRNode* node )
728 node->gmsr = node->tmsr = NULL;
729 node->reinit = 0xffff;
731 node->sizei = node->size = 1;
732 node->prev = node->next = node->shortcut = node;
735 // the preprocess to get the edge list with proper gaussian blur
736 static int preprocessMSER_8UC3( MSCRNode* node,
746 int srccpt = src->step-src->cols*3;
747 uchar* srcptr = src->data.ptr;
748 uchar* lastptr = src->data.ptr+3;
749 double* dxptr = dx->data.db;
750 for ( int i = 0; i < src->rows; i++ )
752 for ( int j = 0; j < src->cols-1; j++ )
754 *dxptr = ChiSquaredDistance( srcptr, lastptr );
762 srcptr = src->data.ptr;
763 lastptr = src->data.ptr+src->step;
764 double* dyptr = dy->data.db;
765 for ( int i = 0; i < src->rows-1; i++ )
767 for ( int j = 0; j < src->cols; j++ )
769 *dyptr = ChiSquaredDistance( srcptr, lastptr );
777 // get dx and dy and blur it
778 if ( edgeBlurSize >= 1 )
780 cvSmooth( dx, dx, CV_GAUSSIAN, edgeBlurSize, edgeBlurSize );
781 cvSmooth( dy, dy, CV_GAUSSIAN, edgeBlurSize, edgeBlurSize );
785 // assian dx, dy to proper edge list and initialize mscr node
786 // the nasty code here intended to avoid extra loops
790 int maskcpt = mask->step-mask->cols+1;
791 uchar* maskptr = mask->data.ptr;
792 MSCRNode* nodeptr = node;
793 initMSCRNode( nodeptr );
795 *total += edge->chi = *dxptr;
796 if ( maskptr[0] && maskptr[1] )
798 edge->left = nodeptr;
799 edge->right = nodeptr+1;
806 for ( int i = 1; i < src->cols-1; i++ )
808 initMSCRNode( nodeptr );
810 if ( maskptr[0] && maskptr[1] )
812 *total += edge->chi = *dxptr;
813 edge->left = nodeptr;
814 edge->right = nodeptr+1;
822 initMSCRNode( nodeptr );
823 nodeptr->index = src->cols-1;
826 for ( int i = 1; i < src->rows-1; i++ )
828 initMSCRNode( nodeptr );
829 nodeptr->index = i<<16;
832 if ( maskptr[-mask->step] )
834 *total += edge->chi = *dyptr;
835 edge->left = nodeptr-src->cols;
836 edge->right = nodeptr;
842 *total += edge->chi = *dxptr;
843 edge->left = nodeptr;
844 edge->right = nodeptr+1;
853 for ( int j = 1; j < src->cols-1; j++ )
855 initMSCRNode( nodeptr );
856 nodeptr->index = (i<<16)|j;
859 if ( maskptr[-mask->step] )
861 *total += edge->chi = *dyptr;
862 edge->left = nodeptr-src->cols;
863 edge->right = nodeptr;
869 *total += edge->chi = *dxptr;
870 edge->left = nodeptr;
871 edge->right = nodeptr+1;
881 initMSCRNode( nodeptr );
882 nodeptr->index = (i<<16)|(src->cols-1);
883 if ( maskptr[0] && maskptr[-mask->step] )
885 *total += edge->chi = *dyptr;
886 edge->left = nodeptr-src->cols;
887 edge->right = nodeptr;
895 initMSCRNode( nodeptr );
896 nodeptr->index = (src->rows-1)<<16;
901 *total += edge->chi = *dxptr;
902 edge->left = nodeptr;
903 edge->right = nodeptr+1;
907 if ( maskptr[-mask->step] )
909 *total += edge->chi = *dyptr;
910 edge->left = nodeptr-src->cols;
911 edge->right = nodeptr;
920 for ( int i = 1; i < src->cols-1; i++ )
922 initMSCRNode( nodeptr );
923 nodeptr->index = ((src->rows-1)<<16)|i;
928 *total += edge->chi = *dxptr;
929 edge->left = nodeptr;
930 edge->right = nodeptr+1;
934 if ( maskptr[-mask->step] )
936 *total += edge->chi = *dyptr;
937 edge->left = nodeptr-src->cols;
938 edge->right = nodeptr;
948 initMSCRNode( nodeptr );
949 nodeptr->index = ((src->rows-1)<<16)|(src->cols-1);
950 if ( maskptr[0] && maskptr[-mask->step] )
952 *total += edge->chi = *dyptr;
953 edge->left = nodeptr-src->cols;
954 edge->right = nodeptr;
958 MSCRNode* nodeptr = node;
959 initMSCRNode( nodeptr );
961 *total += edge->chi = *dxptr;
963 edge->left = nodeptr;
964 edge->right = nodeptr+1;
967 for ( int i = 1; i < src->cols-1; i++ )
969 initMSCRNode( nodeptr );
971 *total += edge->chi = *dxptr;
973 edge->left = nodeptr;
974 edge->right = nodeptr+1;
978 initMSCRNode( nodeptr );
979 nodeptr->index = src->cols-1;
981 for ( int i = 1; i < src->rows-1; i++ )
983 initMSCRNode( nodeptr );
984 nodeptr->index = i<<16;
985 *total += edge->chi = *dyptr;
987 edge->left = nodeptr-src->cols;
988 edge->right = nodeptr;
990 *total += edge->chi = *dxptr;
992 edge->left = nodeptr;
993 edge->right = nodeptr+1;
996 for ( int j = 1; j < src->cols-1; j++ )
998 initMSCRNode( nodeptr );
999 nodeptr->index = (i<<16)|j;
1000 *total += edge->chi = *dyptr;
1002 edge->left = nodeptr-src->cols;
1003 edge->right = nodeptr;
1005 *total += edge->chi = *dxptr;
1007 edge->left = nodeptr;
1008 edge->right = nodeptr+1;
1012 initMSCRNode( nodeptr );
1013 nodeptr->index = (i<<16)|(src->cols-1);
1014 *total += edge->chi = *dyptr;
1016 edge->left = nodeptr-src->cols;
1017 edge->right = nodeptr;
1021 initMSCRNode( nodeptr );
1022 nodeptr->index = (src->rows-1)<<16;
1023 *total += edge->chi = *dxptr;
1025 edge->left = nodeptr;
1026 edge->right = nodeptr+1;
1028 *total += edge->chi = *dyptr;
1030 edge->left = nodeptr-src->cols;
1031 edge->right = nodeptr;
1034 for ( int i = 1; i < src->cols-1; i++ )
1036 initMSCRNode( nodeptr );
1037 nodeptr->index = ((src->rows-1)<<16)|i;
1038 *total += edge->chi = *dxptr;
1040 edge->left = nodeptr;
1041 edge->right = nodeptr+1;
1043 *total += edge->chi = *dyptr;
1045 edge->left = nodeptr-src->cols;
1046 edge->right = nodeptr;
1050 initMSCRNode( nodeptr );
1051 nodeptr->index = ((src->rows-1)<<16)|(src->cols-1);
1052 *total += edge->chi = *dyptr;
1053 edge->left = nodeptr-src->cols;
1054 edge->right = nodeptr;
1062 bool operator()(const MSCREdge& a, const MSCREdge& b) const { return a.chi < b.chi; }
1065 // to find the root of one region
1066 static MSCRNode* findMSCR( MSCRNode* x )
1074 if ( next == x ) break;
1083 if ( prev == x ) break;
1089 // the stable mscr should be:
1090 // bigger than minArea and smaller than maxArea
1091 // differ from its ancestor more than minDiversity
1092 static bool MSCRStableCheck( MSCRNode* x, MSERParams params )
1094 if ( x->size <= params.minArea || x->size >= params.maxArea )
1096 if ( x->gmsr == NULL )
1098 double div = (double)(x->size-x->gmsr->size)/(double)x->size;
1099 return div > params.minDiversity;
1103 extractMSER_8UC3( CvMat* src,
1106 CvMemStorage* storage,
1109 MSCRNode* map = (MSCRNode*)cvAlloc( src->cols*src->rows*sizeof(map[0]) );
1110 int Ne = src->cols*src->rows*2-src->cols-src->rows;
1111 MSCREdge* edge = (MSCREdge*)cvAlloc( Ne*sizeof(edge[0]) );
1112 TempMSCR* mscr = (TempMSCR*)cvAlloc( src->cols*src->rows*sizeof(mscr[0]) );
1114 CvMat* dx = cvCreateMat( src->rows, src->cols-1, CV_64FC1 );
1115 CvMat* dy = cvCreateMat( src->rows-1, src->cols, CV_64FC1 );
1116 Ne = preprocessMSER_8UC3( map, edge, &emean, src, mask, dx, dy, Ne, params.edgeBlurSize );
1117 emean = emean / (double)Ne;
1118 std::sort(edge, edge + Ne, LessThanEdge());
1119 MSCREdge* edge_ub = edge+Ne;
1120 MSCREdge* edgeptr = edge;
1121 TempMSCR* mscrptr = mscr;
1122 // the evolution process
1123 for ( int i = 0; i < params.maxEvolution; i++ )
1125 double k = (double)i/(double)params.maxEvolution*(TABLE_SIZE-1);
1126 int ti = cvFloor(k);
1127 double reminder = k-ti;
1128 double thres = emean*(chitab3[ti]*(1-reminder)+chitab3[ti+1]*reminder);
1129 // to process all the edges in the list that chi < thres
1130 while ( edgeptr < edge_ub && edgeptr->chi < thres )
1132 MSCRNode* lr = findMSCR( edgeptr->left );
1133 MSCRNode* rr = findMSCR( edgeptr->right );
1134 // get the region root (who is responsible)
1137 // rank idea take from: N-tree Disjoint-Set Forests for Maximally Stable Extremal Regions
1138 if ( rr->rank > lr->rank )
1141 CV_SWAP( lr, rr, tmp );
1142 } else if ( lr->rank == rr->rank ) {
1143 // at the same rank, we will compare the size
1144 if ( lr->size > rr->size )
1147 CV_SWAP( lr, rr, tmp );
1152 lr->size += rr->size;
1153 // join rr to the end of list lr (lr is a endless double-linked list)
1154 lr->prev->next = rr;
1155 lr->prev = rr->prev;
1156 rr->prev->next = lr;
1158 // area threshold force to reinitialize
1159 if ( lr->size > (lr->size-rr->size)*params.areaThreshold )
1161 lr->sizei = lr->size;
1163 if ( lr->tmsr != NULL )
1165 lr->tmsr->m = lr->dt-lr->di;
1168 lr->di = edgeptr->chi;
1171 lr->dt = edgeptr->chi;
1172 if ( i > lr->reinit )
1174 double s = (double)(lr->size-lr->sizei)/(lr->dt-lr->di);
1177 // skip the first one and check stablity
1178 if ( i > lr->reinit+1 && MSCRStableCheck( lr, params ) )
1180 if ( lr->tmsr == NULL )
1182 lr->gmsr = lr->tmsr = mscrptr;
1185 lr->tmsr->size = lr->size;
1186 lr->tmsr->head = lr;
1187 lr->tmsr->tail = lr->prev;
1196 if ( edgeptr >= edge_ub )
1199 for ( TempMSCR* ptr = mscr; ptr < mscrptr; ptr++ )
1200 // to prune area with margin less than minMargin
1201 if ( ptr->m > params.minMargin )
1203 CvSeq* _contour = cvCreateSeq( CV_SEQ_KIND_GENERIC|CV_32SC2, sizeof(CvContour), sizeof(CvPoint), storage );
1204 cvSeqPushMulti( _contour, 0, ptr->size );
1205 MSCRNode* lpt = ptr->head;
1206 for ( int i = 0; i < ptr->size; i++ )
1208 CvPoint* pt = CV_GET_SEQ_ELEM( CvPoint, _contour, i );
1209 pt->x = (lpt->index)&0xffff;
1210 pt->y = (lpt->index)>>16;
1213 CvContour* contour = (CvContour*)_contour;
1214 cvBoundingRect( contour );
1216 cvSeqPush( contours, &contour );
1218 cvReleaseMat( &dx );
1219 cvReleaseMat( &dy );
1226 extractMSER( CvArr* _img,
1229 CvMemStorage* storage,
1232 CvMat srchdr, *src = cvGetMat( _img, &srchdr );
1233 CvMat maskhdr, *mask = _mask ? cvGetMat( _mask, &maskhdr ) : 0;
1234 CvSeq* contours = 0;
1236 CV_Assert(src != 0);
1237 CV_Assert(CV_MAT_TYPE(src->type) == CV_8UC1 || CV_MAT_TYPE(src->type) == CV_8UC3);
1238 CV_Assert(mask == 0 || (CV_ARE_SIZES_EQ(src, mask) && CV_MAT_TYPE(mask->type) == CV_8UC1));
1239 CV_Assert(storage != 0);
1241 contours = *_contours = cvCreateSeq( 0, sizeof(CvSeq), sizeof(CvSeq*), storage );
1243 // choose different method for different image type
1244 // for grey image, it is: Linear Time Maximally Stable Extremal Regions
1245 // for color image, it is: Maximally Stable Colour Regions for Recognition and Matching
1246 switch ( CV_MAT_TYPE(src->type) )
1249 extractMSER_8UC1( src, mask, contours, storage, params );
1252 extractMSER_8UC3( src, mask, contours, storage, params );
1258 MSER::MSER( int _delta, int _min_area, int _max_area,
1259 double _max_variation, double _min_diversity,
1260 int _max_evolution, double _area_threshold,
1261 double _min_margin, int _edge_blur_size )
1262 : delta(_delta), minArea(_min_area), maxArea(_max_area),
1263 maxVariation(_max_variation), minDiversity(_min_diversity),
1264 maxEvolution(_max_evolution), areaThreshold(_area_threshold),
1265 minMargin(_min_margin), edgeBlurSize(_edge_blur_size)
1269 void MSER::operator()( InputArray image, std::vector<std::vector<Point> >& dstcontours, InputArray mask ) const
1271 CvMat _image = image.getMat(), _mask, *pmask = 0;
1273 pmask = &(_mask = mask.getMat());
1274 MemStorage storage(cvCreateMemStorage(0));
1275 Seq<CvSeq*> contours;
1276 extractMSER( &_image, pmask, &contours.seq, storage,
1277 MSERParams(delta, minArea, maxArea, maxVariation, minDiversity,
1278 maxEvolution, areaThreshold, minMargin, edgeBlurSize));
1279 SeqIterator<CvSeq*> it = contours.begin();
1280 size_t i, ncontours = contours.size();
1281 dstcontours.resize(ncontours);
1282 for( i = 0; i < ncontours; i++, ++it )
1283 Seq<Point>(*it).copyTo(dstcontours[i]);
1287 void MserFeatureDetector::detectImpl( InputArray _image, std::vector<KeyPoint>& keypoints, InputArray _mask ) const
1289 Mat image = _image.getMat(), mask = _mask.getMat();
1290 std::vector<std::vector<Point> > msers;
1292 (*this)(image, msers, mask);
1294 std::vector<std::vector<Point> >::const_iterator contour_it = msers.begin();
1295 Rect r(0, 0, image.cols, image.rows);
1296 for( ; contour_it != msers.end(); ++contour_it )
1298 // TODO check transformation from MSER region to KeyPoint
1299 RotatedRect rect = fitEllipse(Mat(*contour_it));
1300 float diam = std::sqrt(rect.size.height*rect.size.width);
1302 if( diam > std::numeric_limits<float>::epsilon() && r.contains(rect.center) )
1303 keypoints.push_back( KeyPoint(rect.center, diam) );