1 /* Redistribution and use in source and binary forms, with or
2 * without modification, are permitted provided that the following
3 * conditions are met:
4 * Redistributions of source code must retain the above
5 * copyright notice, this list of conditions and the following
6 * disclaimer.
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.
14 *
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
27 * OF SUCH DAMAGE.
28 * Copyright© 2009, Liu Liu All rights reserved.
29 *
30 * OpenCV functions for MSER extraction
31 *
32 * 1. there are two different implementation of MSER, one for gray image, one for color image
33 * 2. the gray 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 gray 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.
40 */
41
42 #include "precomp.hpp"
43 #include "opencv2/imgproc/imgproc_c.h"
44 #include <limits>
45
46 namespace cv
47 {
48
49 using std::vector;
50
51 class MSER_Impl : public MSER
52 {
53 public:
54 struct Params
55 {
Paramscv::MSER_Impl::Params56 Params( int _delta=5, int _min_area=60, int _max_area=14400,
57 double _max_variation=0.25, double _min_diversity=.2,
58 int _max_evolution=200, double _area_threshold=1.01,
59 double _min_margin=0.003, int _edge_blur_size=5 )
60 {
61 delta = _delta;
62 minArea = _min_area;
63 maxArea = _max_area;
64 maxVariation = _max_variation;
65 minDiversity = _min_diversity;
66 maxEvolution = _max_evolution;
67 areaThreshold = _area_threshold;
68 minMargin = _min_margin;
69 edgeBlurSize = _edge_blur_size;
70 pass2Only = false;
71 }
72
73 int delta;
74 int minArea;
75 int maxArea;
76 double maxVariation;
77 double minDiversity;
78 bool pass2Only;
79
80 int maxEvolution;
81 double areaThreshold;
82 double minMargin;
83 int edgeBlurSize;
84 };
85
MSER_Impl(const Params & _params)86 explicit MSER_Impl(const Params& _params) : params(_params) {}
87
~MSER_Impl()88 virtual ~MSER_Impl() {}
89
setDelta(int delta)90 void setDelta(int delta) { params.delta = delta; }
getDelta() const91 int getDelta() const { return params.delta; }
92
setMinArea(int minArea)93 void setMinArea(int minArea) { params.minArea = minArea; }
getMinArea() const94 int getMinArea() const { return params.minArea; }
95
setMaxArea(int maxArea)96 void setMaxArea(int maxArea) { params.maxArea = maxArea; }
getMaxArea() const97 int getMaxArea() const { return params.maxArea; }
98
setPass2Only(bool f)99 void setPass2Only(bool f) { params.pass2Only = f; }
getPass2Only() const100 bool getPass2Only() const { return params.pass2Only; }
101
102 enum { DIR_SHIFT = 29, NEXT_MASK = ((1<<DIR_SHIFT)-1) };
103
104 struct Pixel
105 {
Pixelcv::MSER_Impl::Pixel106 Pixel() : val(0) {}
Pixelcv::MSER_Impl::Pixel107 Pixel(int _val) : val(_val) {}
108
getGraycv::MSER_Impl::Pixel109 int getGray(const Pixel* ptr0, const uchar* imgptr0, int mask) const
110 {
111 return imgptr0[this - ptr0] ^ mask;
112 }
getNextcv::MSER_Impl::Pixel113 int getNext() const { return (val & NEXT_MASK); }
setNextcv::MSER_Impl::Pixel114 void setNext(int next) { val = (val & ~NEXT_MASK) | next; }
115
getDircv::MSER_Impl::Pixel116 int getDir() const { return (int)((unsigned)val >> DIR_SHIFT); }
setDircv::MSER_Impl::Pixel117 void setDir(int dir) { val = (val & NEXT_MASK) | (dir << DIR_SHIFT); }
isVisitedcv::MSER_Impl::Pixel118 bool isVisited() const { return (val & ~NEXT_MASK) != 0; }
119
120 int val;
121 };
122 typedef int PPixel;
123
124 struct WParams
125 {
126 Params p;
127 vector<vector<Point> >* msers;
128 vector<Rect>* bboxvec;
129 Pixel* pix0;
130 int step;
131 };
132
133 // the history of region grown
134 struct CompHistory
135 {
CompHistorycv::MSER_Impl::CompHistory136 CompHistory()
137 {
138 parent_ = child_ = next_ = 0;
139 val = size = 0;
140 var = -1.f;
141 head = 0;
142 checked = false;
143 }
updateTreecv::MSER_Impl::CompHistory144 void updateTree( WParams& wp, CompHistory** _h0, CompHistory** _h1, bool final )
145 {
146 if( var >= 0.f )
147 return;
148 int delta = wp.p.delta;
149
150 CompHistory* h0_ = 0, *h1_ = 0;
151 CompHistory* c = child_;
152 if( size >= wp.p.minArea )
153 {
154 for( ; c != 0; c = c->next_ )
155 {
156 if( c->var < 0.f )
157 c->updateTree(wp, c == child_ ? &h0_ : 0, c == child_ ? &h1_ : 0, final);
158 if( c->var < 0.f )
159 return;
160 }
161 }
162
163 // find h0 and h1 such that:
164 // h0->val >= h->val - delta and (h0->parent == 0 or h0->parent->val < h->val - delta)
165 // h1->val <= h->val + delta and (h1->child == 0 or h1->child->val < h->val + delta)
166 // then we will adjust h0 and h1 as h moves towards latest
167 CompHistory* h0 = this, *h1 = h1_ && h1_->size > size ? h1_ : this;
168 if( h0_ )
169 {
170 for( h0 = h0_; h0 != this && h0->val < val - delta; h0 = h0->parent_ )
171 ;
172 }
173 else
174 {
175 for( ; h0->child_ && h0->child_->val >= val - delta; h0 = h0->child_ )
176 ;
177 }
178
179 for( ; h1->parent_ && h1->parent_->val <= val + delta; h1 = h1->parent_ )
180 ;
181
182 if( _h0 ) *_h0 = h0;
183 if( _h1 ) *_h1 = h1;
184
185 // when we do not well-defined ER(h->val + delta), we stop
186 // the process of computing variances unless we are at the final step
187 if( !final && !h1->parent_ && h1->val < val + delta )
188 return;
189
190 var = (float)(h1->size - h0->size)/size;
191 c = child_;
192 for( ; c != 0; c = c->next_ )
193 c->checkAndCapture(wp);
194 if( final && !parent_ )
195 checkAndCapture(wp);
196 }
197
checkAndCapturecv::MSER_Impl::CompHistory198 void checkAndCapture( WParams& wp )
199 {
200 if( checked )
201 return;
202 checked = true;
203 if( size < wp.p.minArea || size > wp.p.maxArea || var < 0.f || var > wp.p.maxVariation )
204 return;
205 if( child_ )
206 {
207 CompHistory* c = child_;
208 for( ; c != 0; c = c->next_ )
209 {
210 if( c->var >= 0.f && var > c->var )
211 return;
212 }
213 }
214 if( parent_ && parent_->var >= 0.f && var >= parent_->var )
215 return;
216 int xmin = INT_MAX, ymin = INT_MAX, xmax = INT_MIN, ymax = INT_MIN, j = 0;
217 wp.msers->push_back(vector<Point>());
218 vector<Point>& region = wp.msers->back();
219 region.resize(size);
220 const Pixel* pix0 = wp.pix0;
221 int step = wp.step;
222
223 for( PPixel pix = head; j < size; j++, pix = pix0[pix].getNext() )
224 {
225 int y = pix/step;
226 int x = pix - y*step;
227
228 xmin = std::min(xmin, x);
229 xmax = std::max(xmax, x);
230 ymin = std::min(ymin, y);
231 ymax = std::max(ymax, y);
232
233 region[j] = Point(x, y);
234 }
235
236 wp.bboxvec->push_back(Rect(xmin, ymin, xmax - xmin + 1, ymax - ymin + 1));
237 }
238
239 CompHistory* child_;
240 CompHistory* parent_;
241 CompHistory* next_;
242 int val;
243 int size;
244 float var;
245 PPixel head;
246 bool checked;
247 };
248
249 struct ConnectedComp
250 {
ConnectedCompcv::MSER_Impl::ConnectedComp251 ConnectedComp()
252 {
253 init(0);
254 }
255
initcv::MSER_Impl::ConnectedComp256 void init(int gray)
257 {
258 head = tail = 0;
259 history = 0;
260 size = 0;
261 gray_level = gray;
262 }
263
264 // add history chunk to a connected component
growHistorycv::MSER_Impl::ConnectedComp265 void growHistory( CompHistory*& hptr, WParams& wp, int new_gray_level, bool final, bool force=false )
266 {
267 bool update = final;
268 if( new_gray_level < 0 )
269 new_gray_level = gray_level;
270 if( !history || (history->size != size && size > 0 &&
271 (gray_level != history->val || force)))
272 {
273 CompHistory* h = hptr++;
274 h->parent_ = 0;
275 h->child_ = history;
276 h->next_ = 0;
277 if( history )
278 history->parent_ = h;
279 h->val = gray_level;
280 h->size = size;
281 h->head = head;
282
283 history = h;
284 h->var = FLT_MAX;
285 h->checked = true;
286 if( h->size >= wp.p.minArea )
287 {
288 h->var = -1.f;
289 h->checked = false;
290 update = true;
291 }
292 }
293 gray_level = new_gray_level;
294 if( update && history )
295 history->updateTree(wp, 0, 0, final);
296 }
297
298 // merging two connected components
mergecv::MSER_Impl::ConnectedComp299 void merge( ConnectedComp* comp1, ConnectedComp* comp2,
300 CompHistory*& hptr, WParams& wp )
301 {
302 comp1->growHistory( hptr, wp, -1, false );
303 comp2->growHistory( hptr, wp, -1, false );
304
305 if( comp1->size < comp2->size )
306 std::swap(comp1, comp2);
307
308 if( comp2->size == 0 )
309 {
310 gray_level = comp1->gray_level;
311 head = comp1->head;
312 tail = comp1->tail;
313 size = comp1->size;
314 history = comp1->history;
315 return;
316 }
317
318 CompHistory* h1 = comp1->history;
319 CompHistory* h2 = comp2->history;
320
321 gray_level = std::max(comp1->gray_level, comp2->gray_level);
322 history = comp1->history;
323 wp.pix0[comp1->tail].setNext(comp2->head);
324
325 head = comp1->head;
326 tail = comp2->tail;
327 size = comp1->size + comp2->size;
328 bool keep_2nd = h2->size > wp.p.minArea;
329 growHistory( hptr, wp, -1, false, keep_2nd );
330 if( keep_2nd )
331 {
332 h1->next_ = h2;
333 h2->parent_ = history;
334 }
335 }
336
337 PPixel head;
338 PPixel tail;
339 CompHistory* history;
340 int gray_level;
341 int size;
342 };
343
344 void detectRegions( InputArray image,
345 std::vector<std::vector<Point> >& msers,
346 std::vector<Rect>& bboxes );
347 void detect( InputArray _src, vector<KeyPoint>& keypoints, InputArray _mask );
348
preprocess1(const Mat & img,int * level_size)349 void preprocess1( const Mat& img, int* level_size )
350 {
351 memset(level_size, 0, 256*sizeof(level_size[0]));
352
353 int i, j, cols = img.cols, rows = img.rows;
354 int step = cols;
355 pixbuf.resize(step*rows);
356 heapbuf.resize(cols*rows + 256);
357 histbuf.resize(cols*rows);
358 Pixel borderpix;
359 borderpix.setDir(5);
360
361 for( j = 0; j < step; j++ )
362 {
363 pixbuf[j] = pixbuf[j + (rows-1)*step] = borderpix;
364 }
365
366 for( i = 1; i < rows-1; i++ )
367 {
368 const uchar* imgptr = img.ptr(i);
369 Pixel* pptr = &pixbuf[i*step];
370 pptr[0] = pptr[cols-1] = borderpix;
371 for( j = 1; j < cols-1; j++ )
372 {
373 int val = imgptr[j];
374 level_size[val]++;
375 pptr[j].val = 0;
376 }
377 }
378 }
379
preprocess2(const Mat & img,int * level_size)380 void preprocess2( const Mat& img, int* level_size )
381 {
382 int i;
383
384 for( i = 0; i < 128; i++ )
385 std::swap(level_size[i], level_size[255-i]);
386
387 if( !params.pass2Only )
388 {
389 int j, cols = img.cols, rows = img.rows;
390 int step = cols;
391 for( i = 1; i < rows-1; i++ )
392 {
393 Pixel* pptr = &pixbuf[i*step + 1];
394 for( j = 1; j < cols-1; j++ )
395 {
396 pptr[j].val = 0;
397 }
398 }
399 }
400 }
401
pass(const Mat & img,vector<vector<Point>> & msers,vector<Rect> & bboxvec,Size size,const int * level_size,int mask)402 void pass( const Mat& img, vector<vector<Point> >& msers, vector<Rect>& bboxvec,
403 Size size, const int* level_size, int mask )
404 {
405 CompHistory* histptr = &histbuf[0];
406 int step = size.width;
407 Pixel *ptr0 = &pixbuf[0], *ptr = &ptr0[step+1];
408 const uchar* imgptr0 = img.ptr();
409 Pixel** heap[256];
410 ConnectedComp comp[257];
411 ConnectedComp* comptr = &comp[0];
412 WParams wp;
413 wp.p = params;
414 wp.msers = &msers;
415 wp.bboxvec = &bboxvec;
416 wp.pix0 = ptr0;
417 wp.step = step;
418
419 heap[0] = &heapbuf[0];
420 heap[0][0] = 0;
421
422 for( int i = 1; i < 256; i++ )
423 {
424 heap[i] = heap[i-1] + level_size[i-1] + 1;
425 heap[i][0] = 0;
426 }
427
428 comptr->gray_level = 256;
429 comptr++;
430 comptr->gray_level = ptr->getGray(ptr0, imgptr0, mask);
431 ptr->setDir(1);
432 int dir[] = { 0, 1, step, -1, -step };
433 for( ;; )
434 {
435 int curr_gray = ptr->getGray(ptr0, imgptr0, mask);
436 int nbr_idx = ptr->getDir();
437 // take tour of all the 4 directions
438 for( ; nbr_idx <= 4; nbr_idx++ )
439 {
440 // get the neighbor
441 Pixel* ptr_nbr = ptr + dir[nbr_idx];
442 if( !ptr_nbr->isVisited() )
443 {
444 // set dir=1, next=0
445 ptr_nbr->val = 1 << DIR_SHIFT;
446 int nbr_gray = ptr_nbr->getGray(ptr0, imgptr0, mask);
447 if( nbr_gray < curr_gray )
448 {
449 // when the value of neighbor smaller than current
450 // push current to boundary heap and make the neighbor to be the current one
451 // create an empty comp
452 *(++heap[curr_gray]) = ptr;
453 ptr->val = (nbr_idx+1) << DIR_SHIFT;
454 ptr = ptr_nbr;
455 comptr++;
456 comptr->init(nbr_gray);
457 curr_gray = nbr_gray;
458 nbr_idx = 0;
459 continue;
460 }
461 // otherwise, push the neighbor to boundary heap
462 *(++heap[nbr_gray]) = ptr_nbr;
463 }
464 }
465
466 // set dir = nbr_idx, next = 0
467 ptr->val = nbr_idx << DIR_SHIFT;
468 int ptrofs = (int)(ptr - ptr0);
469 CV_Assert(ptrofs != 0);
470
471 // add a pixel to the pixel list
472 if( comptr->tail )
473 ptr0[comptr->tail].setNext(ptrofs);
474 else
475 comptr->head = ptrofs;
476 comptr->tail = ptrofs;
477 comptr->size++;
478 // get the next pixel from boundary heap
479 if( *heap[curr_gray] )
480 {
481 ptr = *heap[curr_gray];
482 heap[curr_gray]--;
483 }
484 else
485 {
486 for( curr_gray++; curr_gray < 256; curr_gray++ )
487 {
488 if( *heap[curr_gray] )
489 break;
490 }
491 if( curr_gray >= 256 )
492 break;
493
494 ptr = *heap[curr_gray];
495 heap[curr_gray]--;
496
497 if( curr_gray < comptr[-1].gray_level )
498 comptr->growHistory(histptr, wp, curr_gray, false);
499 else
500 {
501 // keep merging top two comp in stack until the gray level >= pixel_val
502 for(;;)
503 {
504 comptr--;
505 comptr->merge(comptr, comptr+1, histptr, wp);
506 if( curr_gray <= comptr[0].gray_level )
507 break;
508 if( curr_gray < comptr[-1].gray_level )
509 {
510 comptr->growHistory(histptr, wp, curr_gray, false);
511 break;
512 }
513 }
514 }
515 }
516 }
517
518 for( ; comptr->gray_level != 256; comptr-- )
519 {
520 comptr->growHistory(histptr, wp, 256, true, true);
521 }
522 }
523
524 Mat tempsrc;
525 vector<Pixel> pixbuf;
526 vector<Pixel*> heapbuf;
527 vector<CompHistory> histbuf;
528
529 Params params;
530 };
531
532 /*
533
534 TODO:
535 the color MSER has not been completely refactored yet. We leave it mostly as-is,
536 with just enough changes to convert C structures to C++ ones and
537 add support for color images into MSER_Impl::detectAndLabel.
538 */
539
540 const int TABLE_SIZE = 400;
541
542 static const float chitab3[]=
543 {
544 0.f, 0.0150057f, 0.0239478f, 0.0315227f,
545 0.0383427f, 0.0446605f, 0.0506115f, 0.0562786f,
546 0.0617174f, 0.0669672f, 0.0720573f, 0.0770099f,
547 0.081843f, 0.0865705f, 0.0912043f, 0.0957541f,
548 0.100228f, 0.104633f, 0.108976f, 0.113261f,
549 0.117493f, 0.121676f, 0.125814f, 0.12991f,
550 0.133967f, 0.137987f, 0.141974f, 0.145929f,
551 0.149853f, 0.15375f, 0.15762f, 0.161466f,
552 0.165287f, 0.169087f, 0.172866f, 0.176625f,
553 0.180365f, 0.184088f, 0.187794f, 0.191483f,
554 0.195158f, 0.198819f, 0.202466f, 0.2061f,
555 0.209722f, 0.213332f, 0.216932f, 0.220521f,
556 0.2241f, 0.22767f, 0.231231f, 0.234783f,
557 0.238328f, 0.241865f, 0.245395f, 0.248918f,
558 0.252435f, 0.255947f, 0.259452f, 0.262952f,
559 0.266448f, 0.269939f, 0.273425f, 0.276908f,
560 0.280386f, 0.283862f, 0.287334f, 0.290803f,
561 0.29427f, 0.297734f, 0.301197f, 0.304657f,
562 0.308115f, 0.311573f, 0.315028f, 0.318483f,
563 0.321937f, 0.32539f, 0.328843f, 0.332296f,
564 0.335749f, 0.339201f, 0.342654f, 0.346108f,
565 0.349562f, 0.353017f, 0.356473f, 0.35993f,
566 0.363389f, 0.366849f, 0.37031f, 0.373774f,
567 0.377239f, 0.380706f, 0.384176f, 0.387648f,
568 0.391123f, 0.3946f, 0.39808f, 0.401563f,
569 0.405049f, 0.408539f, 0.412032f, 0.415528f,
570 0.419028f, 0.422531f, 0.426039f, 0.429551f,
571 0.433066f, 0.436586f, 0.440111f, 0.44364f,
572 0.447173f, 0.450712f, 0.454255f, 0.457803f,
573 0.461356f, 0.464915f, 0.468479f, 0.472049f,
574 0.475624f, 0.479205f, 0.482792f, 0.486384f,
575 0.489983f, 0.493588f, 0.4972f, 0.500818f,
576 0.504442f, 0.508073f, 0.511711f, 0.515356f,
577 0.519008f, 0.522667f, 0.526334f, 0.530008f,
578 0.533689f, 0.537378f, 0.541075f, 0.54478f,
579 0.548492f, 0.552213f, 0.555942f, 0.55968f,
580 0.563425f, 0.56718f, 0.570943f, 0.574715f,
581 0.578497f, 0.582287f, 0.586086f, 0.589895f,
582 0.593713f, 0.597541f, 0.601379f, 0.605227f,
583 0.609084f, 0.612952f, 0.61683f, 0.620718f,
584 0.624617f, 0.628526f, 0.632447f, 0.636378f,
585 0.64032f, 0.644274f, 0.648239f, 0.652215f,
586 0.656203f, 0.660203f, 0.664215f, 0.668238f,
587 0.672274f, 0.676323f, 0.680384f, 0.684457f,
588 0.688543f, 0.692643f, 0.696755f, 0.700881f,
589 0.70502f, 0.709172f, 0.713339f, 0.717519f,
590 0.721714f, 0.725922f, 0.730145f, 0.734383f,
591 0.738636f, 0.742903f, 0.747185f, 0.751483f,
592 0.755796f, 0.760125f, 0.76447f, 0.768831f,
593 0.773208f, 0.777601f, 0.782011f, 0.786438f,
594 0.790882f, 0.795343f, 0.799821f, 0.804318f,
595 0.808831f, 0.813363f, 0.817913f, 0.822482f,
596 0.827069f, 0.831676f, 0.836301f, 0.840946f,
597 0.84561f, 0.850295f, 0.854999f, 0.859724f,
598 0.864469f, 0.869235f, 0.874022f, 0.878831f,
599 0.883661f, 0.888513f, 0.893387f, 0.898284f,
600 0.903204f, 0.908146f, 0.913112f, 0.918101f,
601 0.923114f, 0.928152f, 0.933214f, 0.938301f,
602 0.943413f, 0.94855f, 0.953713f, 0.958903f,
603 0.964119f, 0.969361f, 0.974631f, 0.979929f,
604 0.985254f, 0.990608f, 0.99599f, 1.0014f,
605 1.00684f, 1.01231f, 1.01781f, 1.02335f,
606 1.02891f, 1.0345f, 1.04013f, 1.04579f,
607 1.05148f, 1.05721f, 1.06296f, 1.06876f,
608 1.07459f, 1.08045f, 1.08635f, 1.09228f,
609 1.09826f, 1.10427f, 1.11032f, 1.1164f,
610 1.12253f, 1.1287f, 1.1349f, 1.14115f,
611 1.14744f, 1.15377f, 1.16015f, 1.16656f,
612 1.17303f, 1.17954f, 1.18609f, 1.19269f,
613 1.19934f, 1.20603f, 1.21278f, 1.21958f,
614 1.22642f, 1.23332f, 1.24027f, 1.24727f,
615 1.25433f, 1.26144f, 1.26861f, 1.27584f,
616 1.28312f, 1.29047f, 1.29787f, 1.30534f,
617 1.31287f, 1.32046f, 1.32812f, 1.33585f,
618 1.34364f, 1.3515f, 1.35943f, 1.36744f,
619 1.37551f, 1.38367f, 1.39189f, 1.4002f,
620 1.40859f, 1.41705f, 1.42561f, 1.43424f,
621 1.44296f, 1.45177f, 1.46068f, 1.46967f,
622 1.47876f, 1.48795f, 1.49723f, 1.50662f,
623 1.51611f, 1.52571f, 1.53541f, 1.54523f,
624 1.55517f, 1.56522f, 1.57539f, 1.58568f,
625 1.59611f, 1.60666f, 1.61735f, 1.62817f,
626 1.63914f, 1.65025f, 1.66152f, 1.67293f,
627 1.68451f, 1.69625f, 1.70815f, 1.72023f,
628 1.73249f, 1.74494f, 1.75757f, 1.77041f,
629 1.78344f, 1.79669f, 1.81016f, 1.82385f,
630 1.83777f, 1.85194f, 1.86635f, 1.88103f,
631 1.89598f, 1.91121f, 1.92674f, 1.94257f,
632 1.95871f, 1.97519f, 1.99201f, 2.0092f,
633 2.02676f, 2.04471f, 2.06309f, 2.08189f,
634 2.10115f, 2.12089f, 2.14114f, 2.16192f,
635 2.18326f, 2.2052f, 2.22777f, 2.25101f,
636 2.27496f, 2.29966f, 2.32518f, 2.35156f,
637 2.37886f, 2.40717f, 2.43655f, 2.46709f,
638 2.49889f, 2.53206f, 2.56673f, 2.60305f,
639 2.64117f, 2.6813f, 2.72367f, 2.76854f,
640 2.81623f, 2.86714f, 2.92173f, 2.98059f,
641 3.04446f, 3.1143f, 3.19135f, 3.27731f,
642 3.37455f, 3.48653f, 3.61862f, 3.77982f,
643 3.98692f, 4.2776f, 4.77167f, 133.333f
644 };
645
646 struct MSCRNode;
647
648 struct TempMSCR
649 {
650 MSCRNode* head;
651 MSCRNode* tail;
652 double m; // the margin used to prune area later
653 int size;
654 };
655
656 struct MSCRNode
657 {
658 MSCRNode* shortcut;
659 // to make the finding of root less painful
660 MSCRNode* prev;
661 MSCRNode* next;
662 // a point double-linked list
663 TempMSCR* tmsr;
664 // the temporary msr (set to NULL at every re-initialise)
665 TempMSCR* gmsr;
666 // the global msr (once set, never to NULL)
667 int index;
668 // the index of the node, at this point, it should be x at the first 16-bits, and y at the last 16-bits.
669 int rank;
670 int reinit;
671 int size, sizei;
672 double dt, di;
673 double s;
674 };
675
676 struct MSCREdge
677 {
678 double chi;
679 MSCRNode* left;
680 MSCRNode* right;
681 };
682
ChiSquaredDistance(const uchar * x,const uchar * y)683 static double ChiSquaredDistance( const uchar* x, const uchar* y )
684 {
685 return (double)((x[0]-y[0])*(x[0]-y[0]))/(double)(x[0]+y[0]+1e-10)+
686 (double)((x[1]-y[1])*(x[1]-y[1]))/(double)(x[1]+y[1]+1e-10)+
687 (double)((x[2]-y[2])*(x[2]-y[2]))/(double)(x[2]+y[2]+1e-10);
688 }
689
initMSCRNode(MSCRNode * node)690 static void initMSCRNode( MSCRNode* node )
691 {
692 node->gmsr = node->tmsr = NULL;
693 node->reinit = 0xffff;
694 node->rank = 0;
695 node->sizei = node->size = 1;
696 node->prev = node->next = node->shortcut = node;
697 }
698
699 // the preprocess to get the edge list with proper gaussian blur
preprocessMSER_8uC3(MSCRNode * node,MSCREdge * edge,double * total,const Mat & src,Mat & dx,Mat & dy,int Ne,int edgeBlurSize)700 static int preprocessMSER_8uC3( MSCRNode* node,
701 MSCREdge* edge,
702 double* total,
703 const Mat& src,
704 Mat& dx,
705 Mat& dy,
706 int Ne,
707 int edgeBlurSize )
708 {
709 int srccpt = (int)(src.step-src.cols*3);
710 const uchar* srcptr = src.ptr();
711 const uchar* lastptr = srcptr+3;
712 double* dxptr = dx.ptr<double>();
713 for ( int i = 0; i < src.rows; i++ )
714 {
715 for ( int j = 0; j < src.cols-1; j++ )
716 {
717 *dxptr = ChiSquaredDistance( srcptr, lastptr );
718 dxptr++;
719 srcptr += 3;
720 lastptr += 3;
721 }
722 srcptr += srccpt+3;
723 lastptr += srccpt+3;
724 }
725 srcptr = src.ptr();
726 lastptr = srcptr+src.step;
727 double* dyptr = dy.ptr<double>();
728 for ( int i = 0; i < src.rows-1; i++ )
729 {
730 for ( int j = 0; j < src.cols; j++ )
731 {
732 *dyptr = ChiSquaredDistance( srcptr, lastptr );
733 dyptr++;
734 srcptr += 3;
735 lastptr += 3;
736 }
737 srcptr += srccpt;
738 lastptr += srccpt;
739 }
740 // get dx and dy and blur it
741 if ( edgeBlurSize >= 1 )
742 {
743 GaussianBlur( dx, dx, Size(edgeBlurSize, edgeBlurSize), 0 );
744 GaussianBlur( dy, dy, Size(edgeBlurSize, edgeBlurSize), 0 );
745 }
746 dxptr = dx.ptr<double>();
747 dyptr = dy.ptr<double>();
748 // assian dx, dy to proper edge list and initialize mscr node
749 // the nasty code here intended to avoid extra loops
750 MSCRNode* nodeptr = node;
751 initMSCRNode( nodeptr );
752 nodeptr->index = 0;
753 *total += edge->chi = *dxptr;
754 dxptr++;
755 edge->left = nodeptr;
756 edge->right = nodeptr+1;
757 edge++;
758 nodeptr++;
759 for ( int i = 1; i < src.cols-1; i++ )
760 {
761 initMSCRNode( nodeptr );
762 nodeptr->index = i;
763 *total += edge->chi = *dxptr;
764 dxptr++;
765 edge->left = nodeptr;
766 edge->right = nodeptr+1;
767 edge++;
768 nodeptr++;
769 }
770 initMSCRNode( nodeptr );
771 nodeptr->index = src.cols-1;
772 nodeptr++;
773 for ( int i = 1; i < src.rows-1; i++ )
774 {
775 initMSCRNode( nodeptr );
776 nodeptr->index = i<<16;
777 *total += edge->chi = *dyptr;
778 dyptr++;
779 edge->left = nodeptr-src.cols;
780 edge->right = nodeptr;
781 edge++;
782 *total += edge->chi = *dxptr;
783 dxptr++;
784 edge->left = nodeptr;
785 edge->right = nodeptr+1;
786 edge++;
787 nodeptr++;
788 for ( int j = 1; j < src.cols-1; j++ )
789 {
790 initMSCRNode( nodeptr );
791 nodeptr->index = (i<<16)|j;
792 *total += edge->chi = *dyptr;
793 dyptr++;
794 edge->left = nodeptr-src.cols;
795 edge->right = nodeptr;
796 edge++;
797 *total += edge->chi = *dxptr;
798 dxptr++;
799 edge->left = nodeptr;
800 edge->right = nodeptr+1;
801 edge++;
802 nodeptr++;
803 }
804 initMSCRNode( nodeptr );
805 nodeptr->index = (i<<16)|(src.cols-1);
806 *total += edge->chi = *dyptr;
807 dyptr++;
808 edge->left = nodeptr-src.cols;
809 edge->right = nodeptr;
810 edge++;
811 nodeptr++;
812 }
813 initMSCRNode( nodeptr );
814 nodeptr->index = (src.rows-1)<<16;
815 *total += edge->chi = *dxptr;
816 dxptr++;
817 edge->left = nodeptr;
818 edge->right = nodeptr+1;
819 edge++;
820 *total += edge->chi = *dyptr;
821 dyptr++;
822 edge->left = nodeptr-src.cols;
823 edge->right = nodeptr;
824 edge++;
825 nodeptr++;
826 for ( int i = 1; i < src.cols-1; i++ )
827 {
828 initMSCRNode( nodeptr );
829 nodeptr->index = ((src.rows-1)<<16)|i;
830 *total += edge->chi = *dxptr;
831 dxptr++;
832 edge->left = nodeptr;
833 edge->right = nodeptr+1;
834 edge++;
835 *total += edge->chi = *dyptr;
836 dyptr++;
837 edge->left = nodeptr-src.cols;
838 edge->right = nodeptr;
839 edge++;
840 nodeptr++;
841 }
842 initMSCRNode( nodeptr );
843 nodeptr->index = ((src.rows-1)<<16)|(src.cols-1);
844 *total += edge->chi = *dyptr;
845 edge->left = nodeptr-src.cols;
846 edge->right = nodeptr;
847
848 return Ne;
849 }
850
851 class LessThanEdge
852 {
853 public:
operator ()(const MSCREdge & a,const MSCREdge & b) const854 bool operator()(const MSCREdge& a, const MSCREdge& b) const { return a.chi < b.chi; }
855 };
856
857 // to find the root of one region
findMSCR(MSCRNode * x)858 static MSCRNode* findMSCR( MSCRNode* x )
859 {
860 MSCRNode* prev = x;
861 MSCRNode* next;
862 for ( ; ; )
863 {
864 next = x->shortcut;
865 x->shortcut = prev;
866 if ( next == x ) break;
867 prev= x;
868 x = next;
869 }
870 MSCRNode* root = x;
871 for ( ; ; )
872 {
873 prev = x->shortcut;
874 x->shortcut = root;
875 if ( prev == x ) break;
876 x = prev;
877 }
878 return root;
879 }
880
881 // the stable mscr should be:
882 // bigger than minArea and smaller than maxArea
883 // differ from its ancestor more than minDiversity
MSCRStableCheck(MSCRNode * x,const MSER_Impl::Params & params)884 static bool MSCRStableCheck( MSCRNode* x, const MSER_Impl::Params& params )
885 {
886 if ( x->size <= params.minArea || x->size >= params.maxArea )
887 return false;
888 if ( x->gmsr == NULL )
889 return true;
890 double div = (double)(x->size-x->gmsr->size)/(double)x->size;
891 return div > params.minDiversity;
892 }
893
894 static void
extractMSER_8uC3(const Mat & src,vector<vector<Point>> & msers,vector<Rect> & bboxvec,const MSER_Impl::Params & params)895 extractMSER_8uC3( const Mat& src,
896 vector<vector<Point> >& msers,
897 vector<Rect>& bboxvec,
898 const MSER_Impl::Params& params )
899 {
900 bboxvec.clear();
901 MSCRNode* map = (MSCRNode*)cvAlloc( src.cols*src.rows*sizeof(map[0]) );
902 int Ne = src.cols*src.rows*2-src.cols-src.rows;
903 MSCREdge* edge = (MSCREdge*)cvAlloc( Ne*sizeof(edge[0]) );
904 TempMSCR* mscr = (TempMSCR*)cvAlloc( src.cols*src.rows*sizeof(mscr[0]) );
905 double emean = 0;
906 Mat dx( src.rows, src.cols-1, CV_64FC1 );
907 Mat dy( src.rows-1, src.cols, CV_64FC1 );
908 Ne = preprocessMSER_8uC3( map, edge, &emean, src, dx, dy, Ne, params.edgeBlurSize );
909 emean = emean / (double)Ne;
910 std::sort(edge, edge + Ne, LessThanEdge());
911 MSCREdge* edge_ub = edge+Ne;
912 MSCREdge* edgeptr = edge;
913 TempMSCR* mscrptr = mscr;
914 // the evolution process
915 for ( int i = 0; i < params.maxEvolution; i++ )
916 {
917 double k = (double)i/(double)params.maxEvolution*(TABLE_SIZE-1);
918 int ti = cvFloor(k);
919 double reminder = k-ti;
920 double thres = emean*(chitab3[ti]*(1-reminder)+chitab3[ti+1]*reminder);
921 // to process all the edges in the list that chi < thres
922 while ( edgeptr < edge_ub && edgeptr->chi < thres )
923 {
924 MSCRNode* lr = findMSCR( edgeptr->left );
925 MSCRNode* rr = findMSCR( edgeptr->right );
926 // get the region root (who is responsible)
927 if ( lr != rr )
928 {
929 // rank idea take from: N-tree Disjoint-Set Forests for Maximally Stable Extremal Regions
930 if ( rr->rank > lr->rank )
931 {
932 MSCRNode* tmp;
933 CV_SWAP( lr, rr, tmp );
934 } else if ( lr->rank == rr->rank ) {
935 // at the same rank, we will compare the size
936 if ( lr->size > rr->size )
937 {
938 MSCRNode* tmp;
939 CV_SWAP( lr, rr, tmp );
940 }
941 lr->rank++;
942 }
943 rr->shortcut = lr;
944 lr->size += rr->size;
945 // join rr to the end of list lr (lr is a endless double-linked list)
946 lr->prev->next = rr;
947 lr->prev = rr->prev;
948 rr->prev->next = lr;
949 rr->prev = lr;
950 // area threshold force to reinitialize
951 if ( lr->size > (lr->size-rr->size)*params.areaThreshold )
952 {
953 lr->sizei = lr->size;
954 lr->reinit = i;
955 if ( lr->tmsr != NULL )
956 {
957 lr->tmsr->m = lr->dt-lr->di;
958 lr->tmsr = NULL;
959 }
960 lr->di = edgeptr->chi;
961 lr->s = 1e10;
962 }
963 lr->dt = edgeptr->chi;
964 if ( i > lr->reinit )
965 {
966 double s = (double)(lr->size-lr->sizei)/(lr->dt-lr->di);
967 if ( s < lr->s )
968 {
969 // skip the first one and check stablity
970 if ( i > lr->reinit+1 && MSCRStableCheck( lr, params ) )
971 {
972 if ( lr->tmsr == NULL )
973 {
974 lr->gmsr = lr->tmsr = mscrptr;
975 mscrptr++;
976 }
977 lr->tmsr->size = lr->size;
978 lr->tmsr->head = lr;
979 lr->tmsr->tail = lr->prev;
980 lr->tmsr->m = 0;
981 }
982 lr->s = s;
983 }
984 }
985 }
986 edgeptr++;
987 }
988 if ( edgeptr >= edge_ub )
989 break;
990 }
991 for ( TempMSCR* ptr = mscr; ptr < mscrptr; ptr++ )
992 // to prune area with margin less than minMargin
993 if ( ptr->m > params.minMargin )
994 {
995 MSCRNode* lpt = ptr->head;
996 int xmin = INT_MAX, ymin = INT_MAX, xmax = INT_MIN, ymax = INT_MIN;
997 msers.push_back(vector<Point>());
998 vector<Point>& mser = msers.back();
999
1000 for ( int i = 0; i < ptr->size; i++ )
1001 {
1002 Point pt;
1003 pt.x = (lpt->index)&0xffff;
1004 pt.y = (lpt->index)>>16;
1005 xmin = std::min(xmin, pt.x);
1006 xmax = std::max(xmax, pt.x);
1007 ymin = std::min(ymin, pt.y);
1008 ymax = std::max(ymax, pt.y);
1009
1010 lpt = lpt->next;
1011 mser.push_back(pt);
1012 }
1013 bboxvec.push_back(Rect(xmin, ymin, xmax - xmin + 1, ymax - ymin + 1));
1014 }
1015 cvFree( &mscr );
1016 cvFree( &edge );
1017 cvFree( &map );
1018 }
1019
detectRegions(InputArray _src,vector<vector<Point>> & msers,vector<Rect> & bboxes)1020 void MSER_Impl::detectRegions( InputArray _src, vector<vector<Point> >& msers, vector<Rect>& bboxes )
1021 {
1022 Mat src = _src.getMat();
1023 size_t npix = src.total();
1024
1025 msers.clear();
1026 bboxes.clear();
1027
1028 if( npix == 0 )
1029 return;
1030
1031 Size size = src.size();
1032
1033 if( src.type() == CV_8U )
1034 {
1035 int level_size[256];
1036 if( !src.isContinuous() )
1037 {
1038 src.copyTo(tempsrc);
1039 src = tempsrc;
1040 }
1041
1042 // darker to brighter (MSER+)
1043 preprocess1( src, level_size );
1044 if( !params.pass2Only )
1045 pass( src, msers, bboxes, size, level_size, 0 );
1046 // brighter to darker (MSER-)
1047 preprocess2( src, level_size );
1048 pass( src, msers, bboxes, size, level_size, 255 );
1049 }
1050 else
1051 {
1052 CV_Assert( src.type() == CV_8UC3 || src.type() == CV_8UC4 );
1053 extractMSER_8uC3( src, msers, bboxes, params );
1054 }
1055 }
1056
detect(InputArray _image,vector<KeyPoint> & keypoints,InputArray _mask)1057 void MSER_Impl::detect( InputArray _image, vector<KeyPoint>& keypoints, InputArray _mask )
1058 {
1059 vector<Rect> bboxes;
1060 vector<vector<Point> > msers;
1061 Mat mask = _mask.getMat();
1062
1063 detectRegions(_image, msers, bboxes);
1064 int i, ncomps = (int)msers.size();
1065
1066 keypoints.clear();
1067 for( i = 0; i < ncomps; i++ )
1068 {
1069 Rect r = bboxes[i];
1070 // TODO check transformation from MSER region to KeyPoint
1071 RotatedRect rect = fitEllipse(Mat(msers[i]));
1072 float diam = std::sqrt(rect.size.height*rect.size.width);
1073
1074 if( diam > std::numeric_limits<float>::epsilon() && r.contains(rect.center) &&
1075 (mask.empty() || mask.at<uchar>(cvRound(rect.center.y), cvRound(rect.center.x)) != 0) )
1076 keypoints.push_back( KeyPoint(rect.center, diam) );
1077 }
1078 }
1079
create(int _delta,int _min_area,int _max_area,double _max_variation,double _min_diversity,int _max_evolution,double _area_threshold,double _min_margin,int _edge_blur_size)1080 Ptr<MSER> MSER::create( int _delta, int _min_area, int _max_area,
1081 double _max_variation, double _min_diversity,
1082 int _max_evolution, double _area_threshold,
1083 double _min_margin, int _edge_blur_size )
1084 {
1085 return makePtr<MSER_Impl>(
1086 MSER_Impl::Params(_delta, _min_area, _max_area,
1087 _max_variation, _min_diversity,
1088 _max_evolution, _area_threshold,
1089 _min_margin, _edge_blur_size));
1090 }
1091
1092 }
1093