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