modules/features2d/src/mser.cpp

   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 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.
  40  */
  41
  42 #include "precomp.hpp"
  43
  44 namespace cv
  45 {
  46
  47 const int TABLE_SIZE = 400;
  48
  49 static double chitab3[]={0,  0.0150057,  0.0239478,  0.0315227,
  50                   0.0383427,  0.0446605,  0.0506115,  0.0562786,
  51                   0.0617174,  0.0669672,  0.0720573,  0.0770099,
  52                   0.081843,  0.0865705,  0.0912043,  0.0957541,
  53                   0.100228,  0.104633,  0.108976,  0.113261,
  54                   0.117493,  0.121676,  0.125814,  0.12991,
  55                   0.133967,  0.137987,  0.141974,  0.145929,
  56                   0.149853,  0.15375,  0.15762,  0.161466,
  57                   0.165287,  0.169087,  0.172866,  0.176625,
  58                   0.180365,  0.184088,  0.187794,  0.191483,
  59                   0.195158,  0.198819,  0.202466,  0.2061,
  60                   0.209722,  0.213332,  0.216932,  0.220521,
  61                   0.2241,  0.22767,  0.231231,  0.234783,
  62                   0.238328,  0.241865,  0.245395,  0.248918,
  63                   0.252435,  0.255947,  0.259452,  0.262952,
  64                   0.266448,  0.269939,  0.273425,  0.276908,
  65                   0.280386,  0.283862,  0.287334,  0.290803,
  66                   0.29427,  0.297734,  0.301197,  0.304657,
  67                   0.308115,  0.311573,  0.315028,  0.318483,
  68                   0.321937,  0.32539,  0.328843,  0.332296,
  69                   0.335749,  0.339201,  0.342654,  0.346108,
  70                   0.349562,  0.353017,  0.356473,  0.35993,
  71                   0.363389,  0.366849,  0.37031,  0.373774,
  72                   0.377239,  0.380706,  0.384176,  0.387648,
  73                   0.391123,  0.3946,  0.39808,  0.401563,
  74                   0.405049,  0.408539,  0.412032,  0.415528,
  75                   0.419028,  0.422531,  0.426039,  0.429551,
  76                   0.433066,  0.436586,  0.440111,  0.44364,
  77                   0.447173,  0.450712,  0.454255,  0.457803,
  78                   0.461356,  0.464915,  0.468479,  0.472049,
  79                   0.475624,  0.479205,  0.482792,  0.486384,
  80                   0.489983,  0.493588,  0.4972,  0.500818,
  81                   0.504442,  0.508073,  0.511711,  0.515356,
  82                   0.519008,  0.522667,  0.526334,  0.530008,
  83                   0.533689,  0.537378,  0.541075,  0.54478,
  84                   0.548492,  0.552213,  0.555942,  0.55968,
  85                   0.563425,  0.56718,  0.570943,  0.574715,
  86                   0.578497,  0.582287,  0.586086,  0.589895,
  87                   0.593713,  0.597541,  0.601379,  0.605227,
  88                   0.609084,  0.612952,  0.61683,  0.620718,
  89                   0.624617,  0.628526,  0.632447,  0.636378,
  90                   0.64032,  0.644274,  0.648239,  0.652215,
  91                   0.656203,  0.660203,  0.664215,  0.668238,
  92                   0.672274,  0.676323,  0.680384,  0.684457,
  93                   0.688543,  0.692643,  0.696755,  0.700881,
  94                   0.70502,  0.709172,  0.713339,  0.717519,
  95                   0.721714,  0.725922,  0.730145,  0.734383,
  96                   0.738636,  0.742903,  0.747185,  0.751483,
  97                   0.755796,  0.760125,  0.76447,  0.768831,
  98                   0.773208,  0.777601,  0.782011,  0.786438,
  99                   0.790882,  0.795343,  0.799821,  0.804318,
 100                   0.808831,  0.813363,  0.817913,  0.822482,
 101                   0.827069,  0.831676,  0.836301,  0.840946,
 102                   0.84561,  0.850295,  0.854999,  0.859724,
 103                   0.864469,  0.869235,  0.874022,  0.878831,
 104                   0.883661,  0.888513,  0.893387,  0.898284,
 105                   0.903204,  0.908146,  0.913112,  0.918101,
 106                   0.923114,  0.928152,  0.933214,  0.938301,
 107                   0.943413,  0.94855,  0.953713,  0.958903,
 108                   0.964119,  0.969361,  0.974631,  0.979929,
 109                   0.985254,  0.990608,  0.99599,  1.0014,
 110                   1.00684,  1.01231,  1.01781,  1.02335,
 111                   1.02891,  1.0345,  1.04013,  1.04579,
 112                   1.05148,  1.05721,  1.06296,  1.06876,
 113                   1.07459,  1.08045,  1.08635,  1.09228,
 114                   1.09826,  1.10427,  1.11032,  1.1164,
 115                   1.12253,  1.1287,  1.1349,  1.14115,
 116                   1.14744,  1.15377,  1.16015,  1.16656,
 117                   1.17303,  1.17954,  1.18609,  1.19269,
 118                   1.19934,  1.20603,  1.21278,  1.21958,
 119                   1.22642,  1.23332,  1.24027,  1.24727,
 120                   1.25433,  1.26144,  1.26861,  1.27584,
 121                   1.28312,  1.29047,  1.29787,  1.30534,
 122                   1.31287,  1.32046,  1.32812,  1.33585,
 123                   1.34364,  1.3515,  1.35943,  1.36744,
 124                   1.37551,  1.38367,  1.39189,  1.4002,
 125                   1.40859,  1.41705,  1.42561,  1.43424,
 126                   1.44296,  1.45177,  1.46068,  1.46967,
 127                   1.47876,  1.48795,  1.49723,  1.50662,
 128                   1.51611,  1.52571,  1.53541,  1.54523,
 129                   1.55517,  1.56522,  1.57539,  1.58568,
 130                   1.59611,  1.60666,  1.61735,  1.62817,
 131                   1.63914,  1.65025,  1.66152,  1.67293,
 132                   1.68451,  1.69625,  1.70815,  1.72023,
 133                   1.73249,  1.74494,  1.75757,  1.77041,
 134                   1.78344,  1.79669,  1.81016,  1.82385,
 135                   1.83777,  1.85194,  1.86635,  1.88103,
 136                   1.89598,  1.91121,  1.92674,  1.94257,
 137                   1.95871,  1.97519,  1.99201,  2.0092,
 138                   2.02676,  2.04471,  2.06309,  2.08189,
 139                   2.10115,  2.12089,  2.14114,  2.16192,
 140                   2.18326,  2.2052,  2.22777,  2.25101,
 141                   2.27496,  2.29966,  2.32518,  2.35156,
 142                   2.37886,  2.40717,  2.43655,  2.46709,
 143                   2.49889,  2.53206,  2.56673,  2.60305,
 144                   2.64117,  2.6813,  2.72367,  2.76854,
 145                   2.81623,  2.86714,  2.92173,  2.98059,
 146                   3.04446,  3.1143,  3.19135,  3.27731,
 147                   3.37455,  3.48653,  3.61862,  3.77982,
 148                   3.98692,  4.2776,  4.77167,  133.333 };
 149
 150 typedef struct LinkedPoint
 151 {
 152     struct LinkedPoint* prev;
 153     struct LinkedPoint* next;
 154     Point pt;
 155 }
 156 LinkedPoint;
 157
 158 // the history of region grown
 159 typedef struct MSERGrowHistory
 160 {
 161     struct MSERGrowHistory* shortcut;
 162     struct MSERGrowHistory* child;
 163     int stable; // when it ever stabled before, record the size
 164     int val;
 165     int size;
 166 }
 167 MSERGrowHistory;
 168
 169 typedef struct MSERConnectedComp
 170 {
 171     LinkedPoint* head;
 172     LinkedPoint* tail;
 173     MSERGrowHistory* history;
 174     unsigned long grey_level;
 175     int size;
 176     int dvar; // the derivative of last var
 177     float var; // the current variation (most time is the variation of one-step back)
 178 }
 179 MSERConnectedComp;
 180
 181 // Linear Time MSER claims by using bsf can get performance gain, here is the implementation
 182 // however it seems that will not do any good in real world test
 183 inline void _bitset(unsigned long * a, unsigned long b)
 184 {
 185     *a |= 1<<b;
 186 }
 187 inline void _bitreset(unsigned long * a, unsigned long b)
 188 {
 189     *a &= ~(1<<b);
 190 }
 191
 192 struct MSERParams
 193 {
 194     MSERParams( int _delta, int _minArea, int _maxArea, double _maxVariation,
 195                 double _minDiversity, int _maxEvolution, double _areaThreshold,
 196                 double _minMargin, int _edgeBlurSize )
 197         : delta(_delta), minArea(_minArea), maxArea(_maxArea), maxVariation(_maxVariation),
 198         minDiversity(_minDiversity), maxEvolution(_maxEvolution), areaThreshold(_areaThreshold),
 199         minMargin(_minMargin), edgeBlurSize(_edgeBlurSize)
 200     {}
 201     int delta;
 202     int minArea;
 203     int maxArea;
 204     double maxVariation;
 205     double minDiversity;
 206     int maxEvolution;
 207     double areaThreshold;
 208     double minMargin;
 209     int edgeBlurSize;
 210 };
 211
 212 // clear the connected component in stack
 213 static void
 214 initMSERComp( MSERConnectedComp* comp )
 215 {
 216     comp->size = 0;
 217     comp->var = 0;
 218     comp->dvar = 1;
 219     comp->history = NULL;
 220 }
 221
 222 // add history of size to a connected component
 223 static void
 224 MSERNewHistory( MSERConnectedComp* comp, MSERGrowHistory* history )
 225 {
 226     history->child = history;
 227     if ( NULL == comp->history )
 228     {
 229         history->shortcut = history;
 230         history->stable = 0;
 231     } else {
 232         comp->history->child = history;
 233         history->shortcut = comp->history->shortcut;
 234         history->stable = comp->history->stable;
 235     }
 236     history->val = comp->grey_level;
 237     history->size = comp->size;
 238     comp->history = history;
 239 }
 240
 241 // merging two connected component
 242 static void
 243 MSERMergeComp( MSERConnectedComp* comp1,
 244           MSERConnectedComp* comp2,
 245           MSERConnectedComp* comp,
 246           MSERGrowHistory* history )
 247 {
 248     LinkedPoint* head;
 249     LinkedPoint* tail;
 250     comp->grey_level = comp2->grey_level;
 251     history->child = history;
 252     // select the winner by size
 253     if ( comp1->size >= comp2->size )
 254     {
 255         if ( NULL == comp1->history )
 256         {
 257             history->shortcut = history;
 258             history->stable = 0;
 259         } else {
 260             comp1->history->child = history;
 261             history->shortcut = comp1->history->shortcut;
 262             history->stable = comp1->history->stable;
 263         }
 264         if ( NULL != comp2->history && comp2->history->stable > history->stable )
 265             history->stable = comp2->history->stable;
 266         history->val = comp1->grey_level;
 267         history->size = comp1->size;
 268         // put comp1 to history
 269         comp->var = comp1->var;
 270         comp->dvar = comp1->dvar;
 271         if ( comp1->size > 0 && comp2->size > 0 )
 272         {
 273             comp1->tail->next = comp2->head;
 274             comp2->head->prev = comp1->tail;
 275         }
 276         head = ( comp1->size > 0 ) ? comp1->head : comp2->head;
 277         tail = ( comp2->size > 0 ) ? comp2->tail : comp1->tail;
 278         // always made the newly added in the last of the pixel list (comp1 ... comp2)
 279     } else {
 280         if ( NULL == comp2->history )
 281         {
 282             history->shortcut = history;
 283             history->stable = 0;
 284         } else {
 285             comp2->history->child = history;
 286             history->shortcut = comp2->history->shortcut;
 287             history->stable = comp2->history->stable;
 288         }
 289         if ( NULL != comp1->history && comp1->history->stable > history->stable )
 290             history->stable = comp1->history->stable;
 291         history->val = comp2->grey_level;
 292         history->size = comp2->size;
 293         // put comp2 to history
 294         comp->var = comp2->var;
 295         comp->dvar = comp2->dvar;
 296         if ( comp1->size > 0 && comp2->size > 0 )
 297         {
 298             comp2->tail->next = comp1->head;
 299             comp1->head->prev = comp2->tail;
 300         }
 301         head = ( comp2->size > 0 ) ? comp2->head : comp1->head;
 302         tail = ( comp1->size > 0 ) ? comp1->tail : comp2->tail;
 303         // always made the newly added in the last of the pixel list (comp2 ... comp1)
 304     }
 305     comp->head = head;
 306     comp->tail = tail;
 307     comp->history = history;
 308     comp->size = comp1->size + comp2->size;
 309 }
 310
 311 static float
 312 MSERVariationCalc( MSERConnectedComp* comp, int delta )
 313 {
 314     MSERGrowHistory* history = comp->history;
 315     int val = comp->grey_level;
 316     if ( NULL != history )
 317     {
 318         MSERGrowHistory* shortcut = history->shortcut;
 319         while ( shortcut != shortcut->shortcut && shortcut->val + delta > val )
 320             shortcut = shortcut->shortcut;
 321         MSERGrowHistory* child = shortcut->child;
 322         while ( child != child->child && child->val + delta <= val )
 323         {
 324             shortcut = child;
 325             child = child->child;
 326         }
 327         // get the position of history where the shortcut->val <= delta+val and shortcut->child->val >= delta+val
 328         history->shortcut = shortcut;
 329         return (float)(comp->size-shortcut->size)/(float)shortcut->size;
 330         // here is a small modification of MSER where cal ||R_{i}-R_{i-delta}||/||R_{i-delta}||
 331         // in standard MSER, cal ||R_{i+delta}-R_{i-delta}||/||R_{i}||
 332         // my calculation is simpler and much easier to implement
 333     }
 334     return 1.;
 335 }
 336
 337 static bool MSERStableCheck( MSERConnectedComp* comp, MSERParams params )
 338 {
 339     // tricky part: it actually check the stablity of one-step back
 340     if ( comp->history == NULL || comp->history->size <= params.minArea || comp->history->size >= params.maxArea )
 341         return 0;
 342     float div = (float)(comp->history->size-comp->history->stable)/(float)comp->history->size;
 343     float var = MSERVariationCalc( comp, params.delta );
 344     int dvar = ( comp->var < var || (unsigned long)(comp->history->val + 1) < comp->grey_level );
 345     int stable = ( dvar && !comp->dvar && comp->var < params.maxVariation && div > params.minDiversity );
 346     comp->var = var;
 347     comp->dvar = dvar;
 348     if ( stable )
 349         comp->history->stable = comp->history->size;
 350     return stable != 0;
 351 }
 352
 353 // add a pixel to the pixel list
 354 static void accumulateMSERComp( MSERConnectedComp* comp, LinkedPoint* point )
 355 {
 356     if ( comp->size > 0 )
 357     {
 358         point->prev = comp->tail;
 359         comp->tail->next = point;
 360         point->next = NULL;
 361     } else {
 362         point->prev = NULL;
 363         point->next = NULL;
 364         comp->head = point;
 365     }
 366     comp->tail = point;
 367     comp->size++;
 368 }
 369
 370 // convert the point set to CvSeq
 371 static CvContour* MSERToContour( MSERConnectedComp* comp, CvMemStorage* storage )
 372 {
 373     CvSeq* _contour = cvCreateSeq( CV_SEQ_KIND_GENERIC|CV_32SC2, sizeof(CvContour), sizeof(CvPoint), storage );
 374     CvContour* contour = (CvContour*)_contour;
 375     cvSeqPushMulti( _contour, 0, comp->history->size );
 376     LinkedPoint* lpt = comp->head;
 377     for ( int i = 0; i < comp->history->size; i++ )
 378     {
 379         CvPoint* pt = CV_GET_SEQ_ELEM( CvPoint, _contour, i );
 380         pt->x = lpt->pt.x;
 381         pt->y = lpt->pt.y;
 382         lpt = lpt->next;
 383     }
 384     cvBoundingRect( contour );
 385     return contour;
 386 }
 387
 388 // to preprocess src image to following format
 389 // 32-bit image
 390 // > 0 is available, < 0 is visited
 391 // 17~19 bits is the direction
 392 // 8~11 bits is the bucket it falls to (for BitScanForward)
 393 // 0~8 bits is the color
 394 static int* preprocessMSER_8UC1( CvMat* img,
 395             int*** heap_cur,
 396             CvMat* src,
 397             CvMat* mask )
 398 {
 399     int srccpt = src->step-src->cols;
 400     int cpt_1 = img->cols-src->cols-1;
 401     int* imgptr = img->data.i;
 402     int* startptr;
 403
 404     int level_size[256];
 405     for ( int i = 0; i < 256; i++ )
 406         level_size[i] = 0;
 407
 408     for ( int i = 0; i < src->cols+2; i++ )
 409     {
 410         *imgptr = -1;
 411         imgptr++;
 412     }
 413     imgptr += cpt_1-1;
 414     uchar* srcptr = src->data.ptr;
 415     if ( mask )
 416     {
 417         startptr = 0;
 418         uchar* maskptr = mask->data.ptr;
 419         for ( int i = 0; i < src->rows; i++ )
 420         {
 421             *imgptr = -1;
 422             imgptr++;
 423             for ( int j = 0; j < src->cols; j++ )
 424             {
 425                 if ( *maskptr )
 426                 {
 427                     if ( !startptr )
 428                         startptr = imgptr;
 429                     *srcptr = 0xff-*srcptr;
 430                     level_size[*srcptr]++;
 431                     *imgptr = ((*srcptr>>5)<<8)|(*srcptr);
 432                 } else {
 433                     *imgptr = -1;
 434                 }
 435                 imgptr++;
 436                 srcptr++;
 437                 maskptr++;
 438             }
 439             *imgptr = -1;
 440             imgptr += cpt_1;
 441             srcptr += srccpt;
 442             maskptr += srccpt;
 443         }
 444     } else {
 445         startptr = imgptr+img->cols+1;
 446         for ( int i = 0; i < src->rows; i++ )
 447         {
 448             *imgptr = -1;
 449             imgptr++;
 450             for ( int j = 0; j < src->cols; j++ )
 451             {
 452                 *srcptr = 0xff-*srcptr;
 453                 level_size[*srcptr]++;
 454                 *imgptr = ((*srcptr>>5)<<8)|(*srcptr);
 455                 imgptr++;
 456                 srcptr++;
 457             }
 458             *imgptr = -1;
 459             imgptr += cpt_1;
 460             srcptr += srccpt;
 461         }
 462     }
 463     for ( int i = 0; i < src->cols+2; i++ )
 464     {
 465         *imgptr = -1;
 466         imgptr++;
 467     }
 468
 469     heap_cur[0][0] = 0;
 470     for ( int i = 1; i < 256; i++ )
 471     {
 472         heap_cur[i] = heap_cur[i-1]+level_size[i-1]+1;
 473         heap_cur[i][0] = 0;
 474     }
 475     return startptr;
 476 }
 477
 478 static void extractMSER_8UC1_Pass( int* ioptr,
 479               int* imgptr,
 480               int*** heap_cur,
 481               LinkedPoint* ptsptr,
 482               MSERGrowHistory* histptr,
 483               MSERConnectedComp* comptr,
 484               int step,
 485               int stepmask,
 486               int stepgap,
 487               MSERParams params,
 488               int color,
 489               CvSeq* contours,
 490               CvMemStorage* storage )
 491 {
 492     comptr->grey_level = 256;
 493     comptr++;
 494     comptr->grey_level = (*imgptr)&0xff;
 495     initMSERComp( comptr );
 496     *imgptr |= 0x80000000;
 497     heap_cur += (*imgptr)&0xff;
 498     int dir[] = { 1, step, -1, -step };
 499 #ifdef __INTRIN_ENABLED__
 500     unsigned long heapbit[] = { 0, 0, 0, 0, 0, 0, 0, 0 };
 501     unsigned long* bit_cur = heapbit+(((*imgptr)&0x700)>>8);
 502 #endif
 503     for ( ; ; )
 504     {
 505         // take tour of all the 4 directions
 506         while ( ((*imgptr)&0x70000) < 0x40000 )
 507         {
 508             // get the neighbor
 509             int* imgptr_nbr = imgptr+dir[((*imgptr)&0x70000)>>16];
 510             if ( *imgptr_nbr >= 0 ) // if the neighbor is not visited yet
 511             {
 512                 *imgptr_nbr |= 0x80000000; // mark it as visited
 513                 if ( ((*imgptr_nbr)&0xff) < ((*imgptr)&0xff) )
 514                 {
 515                     // when the value of neighbor smaller than current
 516                     // push current to boundary heap and make the neighbor to be the current one
 517                     // create an empty comp
 518                     (*heap_cur)++;
 519                     **heap_cur = imgptr;
 520                     *imgptr += 0x10000;
 521                     heap_cur += ((*imgptr_nbr)&0xff)-((*imgptr)&0xff);
 522 #ifdef __INTRIN_ENABLED__
 523                     _bitset( bit_cur, (*imgptr)&0x1f );
 524                     bit_cur += (((*imgptr_nbr)&0x700)-((*imgptr)&0x700))>>8;
 525 #endif
 526                     imgptr = imgptr_nbr;
 527                     comptr++;
 528                     initMSERComp( comptr );
 529                     comptr->grey_level = (*imgptr)&0xff;
 530                     continue;
 531                 } else {
 532                     // otherwise, push the neighbor to boundary heap
 533                     heap_cur[((*imgptr_nbr)&0xff)-((*imgptr)&0xff)]++;
 534                     *heap_cur[((*imgptr_nbr)&0xff)-((*imgptr)&0xff)] = imgptr_nbr;
 535 #ifdef __INTRIN_ENABLED__
 536                     _bitset( bit_cur+((((*imgptr_nbr)&0x700)-((*imgptr)&0x700))>>8), (*imgptr_nbr)&0x1f );
 537 #endif
 538                 }
 539             }
 540             *imgptr += 0x10000;
 541         }
 542         int imsk = (int)(imgptr-ioptr);
 543         ptsptr->pt = cvPoint( imsk&stepmask, imsk>>stepgap );
 544         // get the current location
 545         accumulateMSERComp( comptr, ptsptr );
 546         ptsptr++;
 547         // get the next pixel from boundary heap
 548         if ( **heap_cur )
 549         {
 550             imgptr = **heap_cur;
 551             (*heap_cur)--;
 552 #ifdef __INTRIN_ENABLED__
 553             if ( !**heap_cur )
 554                 _bitreset( bit_cur, (*imgptr)&0x1f );
 555 #endif
 556         } else {
 557 #ifdef __INTRIN_ENABLED__
 558             bool found_pixel = 0;
 559             unsigned long pixel_val;
 560             for ( int i = ((*imgptr)&0x700)>>8; i < 8; i++ )
 561             {
 562                 if ( _BitScanForward( &pixel_val, *bit_cur ) )
 563                 {
 564                     found_pixel = 1;
 565                     pixel_val += i<<5;
 566                     heap_cur += pixel_val-((*imgptr)&0xff);
 567                     break;
 568                 }
 569                 bit_cur++;
 570             }
 571             if ( found_pixel )
 572 #else
 573             heap_cur++;
 574             unsigned long pixel_val = 0;
 575             for ( unsigned long i = ((*imgptr)&0xff)+1; i < 256; i++ )
 576             {
 577                 if ( **heap_cur )
 578                 {
 579                     pixel_val = i;
 580                     break;
 581                 }
 582                 heap_cur++;
 583             }
 584             if ( pixel_val )
 585 #endif
 586             {
 587                 imgptr = **heap_cur;
 588                 (*heap_cur)--;
 589 #ifdef __INTRIN_ENABLED__
 590                 if ( !**heap_cur )
 591                     _bitreset( bit_cur, pixel_val&0x1f );
 592 #endif
 593                 if ( pixel_val < comptr[-1].grey_level )
 594                 {
 595                     // check the stablity and push a new history, increase the grey level
 596                     if ( MSERStableCheck( comptr, params ) )
 597                     {
 598                         CvContour* contour = MSERToContour( comptr, storage );
 599                         contour->color = color;
 600                         cvSeqPush( contours, &contour );
 601                     }
 602                     MSERNewHistory( comptr, histptr );
 603                     comptr[0].grey_level = pixel_val;
 604                     histptr++;
 605                 } else {
 606                     // keep merging top two comp in stack until the grey level >= pixel_val
 607                     for ( ; ; )
 608                     {
 609                         comptr--;
 610                         MSERMergeComp( comptr+1, comptr, comptr, histptr );
 611                         histptr++;
 612                         if ( pixel_val <= comptr[0].grey_level )
 613                             break;
 614                         if ( pixel_val < comptr[-1].grey_level )
 615                         {
 616                             // check the stablity here otherwise it wouldn't be an ER
 617                             if ( MSERStableCheck( comptr, params ) )
 618                             {
 619                                 CvContour* contour = MSERToContour( comptr, storage );
 620                                 contour->color = color;
 621                                 cvSeqPush( contours, &contour );
 622                             }
 623                             MSERNewHistory( comptr, histptr );
 624                             comptr[0].grey_level = pixel_val;
 625                             histptr++;
 626                             break;
 627                         }
 628                     }
 629                 }
 630             } else
 631                 break;
 632         }
 633     }
 634 }
 635
 636 static void extractMSER_8UC1( CvMat* src,
 637              CvMat* mask,
 638              CvSeq* contours,
 639              CvMemStorage* storage,
 640              MSERParams params )
 641 {
 642     int step = 8;
 643     int stepgap = 3;
 644     while ( step < src->step+2 )
 645     {
 646         step <<= 1;
 647         stepgap++;
 648     }
 649     int stepmask = step-1;
 650
 651     // to speedup the process, make the width to be 2^N
 652     CvMat* img = cvCreateMat( src->rows+2, step, CV_32SC1 );
 653     int* ioptr = img->data.i+step+1;
 654     int* imgptr;
 655
 656     // pre-allocate boundary heap
 657     int** heap = (int**)cvAlloc( (src->rows*src->cols+256)*sizeof(heap[0]) );
 658     int** heap_start[256];
 659     heap_start[0] = heap;
 660
 661     // pre-allocate linked point and grow history
 662     LinkedPoint* pts = (LinkedPoint*)cvAlloc( src->rows*src->cols*sizeof(pts[0]) );
 663     MSERGrowHistory* history = (MSERGrowHistory*)cvAlloc( src->rows*src->cols*sizeof(history[0]) );
 664     MSERConnectedComp comp[257];
 665
 666     // darker to brighter (MSER-)
 667     imgptr = preprocessMSER_8UC1( img, heap_start, src, mask );
 668     extractMSER_8UC1_Pass( ioptr, imgptr, heap_start, pts, history, comp, step, stepmask, stepgap, params, -1, contours, storage );
 669     // brighter to darker (MSER+)
 670     imgptr = preprocessMSER_8UC1( img, heap_start, src, mask );
 671     extractMSER_8UC1_Pass( ioptr, imgptr, heap_start, pts, history, comp, step, stepmask, stepgap, params, 1, contours, storage );
 672
 673     // clean up
 674     cvFree( &history );
 675     cvFree( &heap );
 676     cvFree( &pts );
 677     cvReleaseMat( &img );
 678 }
 679
 680 struct MSCRNode;
 681
 682 struct TempMSCR
 683 {
 684     MSCRNode* head;
 685     MSCRNode* tail;
 686     double m; // the margin used to prune area later
 687     int size;
 688 };
 689
 690 struct MSCRNode
 691 {
 692     MSCRNode* shortcut;
 693     // to make the finding of root less painful
 694     MSCRNode* prev;
 695     MSCRNode* next;
 696     // a point double-linked list
 697     TempMSCR* tmsr;
 698     // the temporary msr (set to NULL at every re-initialise)
 699     TempMSCR* gmsr;
 700     // the global msr (once set, never to NULL)
 701     int index;
 702     // the index of the node, at this point, it should be x at the first 16-bits, and y at the last 16-bits.
 703     int rank;
 704     int reinit;
 705     int size, sizei;
 706     double dt, di;
 707     double s;
 708 };
 709
 710 struct MSCREdge
 711 {
 712     double chi;
 713     MSCRNode* left;
 714     MSCRNode* right;
 715 };
 716
 717 static double ChiSquaredDistance( uchar* x, uchar* y )
 718 {
 719     return (double)((x[0]-y[0])*(x[0]-y[0]))/(double)(x[0]+y[0]+1e-10)+
 720            (double)((x[1]-y[1])*(x[1]-y[1]))/(double)(x[1]+y[1]+1e-10)+
 721            (double)((x[2]-y[2])*(x[2]-y[2]))/(double)(x[2]+y[2]+1e-10);
 722 }
 723
 724 static void initMSCRNode( MSCRNode* node )
 725 {
 726     node->gmsr = node->tmsr = NULL;
 727     node->reinit = 0xffff;
 728     node->rank = 0;
 729     node->sizei = node->size = 1;
 730     node->prev = node->next = node->shortcut = node;
 731 }
 732
 733 // the preprocess to get the edge list with proper gaussian blur
 734 static int preprocessMSER_8UC3( MSCRNode* node,
 735             MSCREdge* edge,
 736             double* total,
 737             CvMat* src,
 738             CvMat* mask,
 739             CvMat* dx,
 740             CvMat* dy,
 741             int Ne,
 742             int edgeBlurSize )
 743 {
 744     int srccpt = src->step-src->cols*3;
 745     uchar* srcptr = src->data.ptr;
 746     uchar* lastptr = src->data.ptr+3;
 747     double* dxptr = dx->data.db;
 748     for ( int i = 0; i < src->rows; i++ )
 749     {
 750         for ( int j = 0; j < src->cols-1; j++ )
 751         {
 752             *dxptr = ChiSquaredDistance( srcptr, lastptr );
 753             dxptr++;
 754             srcptr += 3;
 755             lastptr += 3;
 756         }
 757         srcptr += srccpt+3;
 758         lastptr += srccpt+3;
 759     }
 760     srcptr = src->data.ptr;
 761     lastptr = src->data.ptr+src->step;
 762     double* dyptr = dy->data.db;
 763     for ( int i = 0; i < src->rows-1; i++ )
 764     {
 765         for ( int j = 0; j < src->cols; j++ )
 766         {
 767             *dyptr = ChiSquaredDistance( srcptr, lastptr );
 768             dyptr++;
 769             srcptr += 3;
 770             lastptr += 3;
 771         }
 772         srcptr += srccpt;
 773         lastptr += srccpt;
 774     }
 775     // get dx and dy and blur it
 776     if ( edgeBlurSize >= 1 )
 777     {
 778         cvSmooth( dx, dx, CV_GAUSSIAN, edgeBlurSize, edgeBlurSize );
 779         cvSmooth( dy, dy, CV_GAUSSIAN, edgeBlurSize, edgeBlurSize );
 780     }
 781     dxptr = dx->data.db;
 782     dyptr = dy->data.db;
 783     // assian dx, dy to proper edge list and initialize mscr node
 784     // the nasty code here intended to avoid extra loops
 785     if ( mask )
 786     {
 787         Ne = 0;
 788         int maskcpt = mask->step-mask->cols+1;
 789         uchar* maskptr = mask->data.ptr;
 790         MSCRNode* nodeptr = node;
 791         initMSCRNode( nodeptr );
 792         nodeptr->index = 0;
 793         *total += edge->chi = *dxptr;
 794         if ( maskptr[0] && maskptr[1] )
 795         {
 796             edge->left = nodeptr;
 797             edge->right = nodeptr+1;
 798             edge++;
 799             Ne++;
 800         }
 801         dxptr++;
 802         nodeptr++;
 803         maskptr++;
 804         for ( int i = 1; i < src->cols-1; i++ )
 805         {
 806             initMSCRNode( nodeptr );
 807             nodeptr->index = i;
 808             if ( maskptr[0] && maskptr[1] )
 809             {
 810                 *total += edge->chi = *dxptr;
 811                 edge->left = nodeptr;
 812                 edge->right = nodeptr+1;
 813                 edge++;
 814                 Ne++;
 815             }
 816             dxptr++;
 817             nodeptr++;
 818             maskptr++;
 819         }
 820         initMSCRNode( nodeptr );
 821         nodeptr->index = src->cols-1;
 822         nodeptr++;
 823         maskptr += maskcpt;
 824         for ( int i = 1; i < src->rows-1; i++ )
 825         {
 826             initMSCRNode( nodeptr );
 827             nodeptr->index = i<<16;
 828             if ( maskptr[0] )
 829             {
 830                 if ( maskptr[-mask->step] )
 831                 {
 832                     *total += edge->chi = *dyptr;
 833                     edge->left = nodeptr-src->cols;
 834                     edge->right = nodeptr;
 835                     edge++;
 836                     Ne++;
 837                 }
 838                 if ( maskptr[1] )
 839                 {
 840                     *total += edge->chi = *dxptr;
 841                     edge->left = nodeptr;
 842                     edge->right = nodeptr+1;
 843                     edge++;
 844                     Ne++;
 845                 }
 846             }
 847             dyptr++;
 848             dxptr++;
 849             nodeptr++;
 850             maskptr++;
 851             for ( int j = 1; j < src->cols-1; j++ )
 852             {
 853                 initMSCRNode( nodeptr );
 854                 nodeptr->index = (i<<16)|j;
 855                 if ( maskptr[0] )
 856                 {
 857                     if ( maskptr[-mask->step] )
 858                     {
 859                         *total += edge->chi = *dyptr;
 860                         edge->left = nodeptr-src->cols;
 861                         edge->right = nodeptr;
 862                         edge++;
 863                         Ne++;
 864                     }
 865                     if ( maskptr[1] )
 866                     {
 867                         *total += edge->chi = *dxptr;
 868                         edge->left = nodeptr;
 869                         edge->right = nodeptr+1;
 870                         edge++;
 871                         Ne++;
 872                     }
 873                 }
 874                 dyptr++;
 875                 dxptr++;
 876                 nodeptr++;
 877                 maskptr++;
 878             }
 879             initMSCRNode( nodeptr );
 880             nodeptr->index = (i<<16)|(src->cols-1);
 881             if ( maskptr[0] && maskptr[-mask->step] )
 882             {
 883                 *total += edge->chi = *dyptr;
 884                 edge->left = nodeptr-src->cols;
 885                 edge->right = nodeptr;
 886                 edge++;
 887                 Ne++;
 888             }
 889             dyptr++;
 890             nodeptr++;
 891             maskptr += maskcpt;
 892         }
 893         initMSCRNode( nodeptr );
 894         nodeptr->index = (src->rows-1)<<16;
 895         if ( maskptr[0] )
 896         {
 897             if ( maskptr[1] )
 898             {
 899                 *total += edge->chi = *dxptr;
 900                 edge->left = nodeptr;
 901                 edge->right = nodeptr+1;
 902                 edge++;
 903                 Ne++;
 904             }
 905             if ( maskptr[-mask->step] )
 906             {
 907                 *total += edge->chi = *dyptr;
 908                 edge->left = nodeptr-src->cols;
 909                 edge->right = nodeptr;
 910                 edge++;
 911                 Ne++;
 912             }
 913         }
 914         dxptr++;
 915         dyptr++;
 916         nodeptr++;
 917         maskptr++;
 918         for ( int i = 1; i < src->cols-1; i++ )
 919         {
 920             initMSCRNode( nodeptr );
 921             nodeptr->index = ((src->rows-1)<<16)|i;
 922             if ( maskptr[0] )
 923             {
 924                 if ( maskptr[1] )
 925                 {
 926                     *total += edge->chi = *dxptr;
 927                     edge->left = nodeptr;
 928                     edge->right = nodeptr+1;
 929                     edge++;
 930                     Ne++;
 931                 }
 932                 if ( maskptr[-mask->step] )
 933                 {
 934                     *total += edge->chi = *dyptr;
 935                     edge->left = nodeptr-src->cols;
 936                     edge->right = nodeptr;
 937                     edge++;
 938                     Ne++;
 939                 }
 940             }
 941             dxptr++;
 942             dyptr++;
 943             nodeptr++;
 944             maskptr++;
 945         }
 946         initMSCRNode( nodeptr );
 947         nodeptr->index = ((src->rows-1)<<16)|(src->cols-1);
 948         if ( maskptr[0] && maskptr[-mask->step] )
 949         {
 950             *total += edge->chi = *dyptr;
 951             edge->left = nodeptr-src->cols;
 952             edge->right = nodeptr;
 953             Ne++;
 954         }
 955     } else {
 956         MSCRNode* nodeptr = node;
 957         initMSCRNode( nodeptr );
 958         nodeptr->index = 0;
 959         *total += edge->chi = *dxptr;
 960         dxptr++;
 961         edge->left = nodeptr;
 962         edge->right = nodeptr+1;
 963         edge++;
 964         nodeptr++;
 965         for ( int i = 1; i < src->cols-1; i++ )
 966         {
 967             initMSCRNode( nodeptr );
 968             nodeptr->index = i;
 969             *total += edge->chi = *dxptr;
 970             dxptr++;
 971             edge->left = nodeptr;
 972             edge->right = nodeptr+1;
 973             edge++;
 974             nodeptr++;
 975         }
 976         initMSCRNode( nodeptr );
 977         nodeptr->index = src->cols-1;
 978         nodeptr++;
 979         for ( int i = 1; i < src->rows-1; i++ )
 980         {
 981             initMSCRNode( nodeptr );
 982             nodeptr->index = i<<16;
 983             *total += edge->chi = *dyptr;
 984             dyptr++;
 985             edge->left = nodeptr-src->cols;
 986             edge->right = nodeptr;
 987             edge++;
 988             *total += edge->chi = *dxptr;
 989             dxptr++;
 990             edge->left = nodeptr;
 991             edge->right = nodeptr+1;
 992             edge++;
 993             nodeptr++;
 994             for ( int j = 1; j < src->cols-1; j++ )
 995             {
 996                 initMSCRNode( nodeptr );
 997                 nodeptr->index = (i<<16)|j;
 998                 *total += edge->chi = *dyptr;
 999                 dyptr++;
1000                 edge->left = nodeptr-src->cols;
1001                 edge->right = nodeptr;
1002                 edge++;
1003                 *total += edge->chi = *dxptr;
1004                 dxptr++;
1005                 edge->left = nodeptr;
1006                 edge->right = nodeptr+1;
1007                 edge++;
1008                 nodeptr++;
1009             }
1010             initMSCRNode( nodeptr );
1011             nodeptr->index = (i<<16)|(src->cols-1);
1012             *total += edge->chi = *dyptr;
1013             dyptr++;
1014             edge->left = nodeptr-src->cols;
1015             edge->right = nodeptr;
1016             edge++;
1017             nodeptr++;
1018         }
1019         initMSCRNode( nodeptr );
1020         nodeptr->index = (src->rows-1)<<16;
1021         *total += edge->chi = *dxptr;
1022         dxptr++;
1023         edge->left = nodeptr;
1024         edge->right = nodeptr+1;
1025         edge++;
1026         *total += edge->chi = *dyptr;
1027         dyptr++;
1028         edge->left = nodeptr-src->cols;
1029         edge->right = nodeptr;
1030         edge++;
1031         nodeptr++;
1032         for ( int i = 1; i < src->cols-1; i++ )
1033         {
1034             initMSCRNode( nodeptr );
1035             nodeptr->index = ((src->rows-1)<<16)|i;
1036             *total += edge->chi = *dxptr;
1037             dxptr++;
1038             edge->left = nodeptr;
1039             edge->right = nodeptr+1;
1040             edge++;
1041             *total += edge->chi = *dyptr;
1042             dyptr++;
1043             edge->left = nodeptr-src->cols;
1044             edge->right = nodeptr;
1045             edge++;
1046             nodeptr++;
1047         }
1048         initMSCRNode( nodeptr );
1049         nodeptr->index = ((src->rows-1)<<16)|(src->cols-1);
1050         *total += edge->chi = *dyptr;
1051         edge->left = nodeptr-src->cols;
1052         edge->right = nodeptr;
1053     }
1054     return Ne;
1055 }
1056
1057 #define cmp_mscr_edge(edge1, edge2) \
1058     ((edge1).chi < (edge2).chi)
1059
1060 static CV_IMPLEMENT_QSORT( QuickSortMSCREdge, MSCREdge, cmp_mscr_edge )
1061
1062 // to find the root of one region
1063 static MSCRNode* findMSCR( MSCRNode* x )
1064 {
1065     MSCRNode* prev = x;
1066     MSCRNode* next;
1067     for ( ; ; )
1068     {
1069         next = x->shortcut;
1070         x->shortcut = prev;
1071         if ( next == x ) break;
1072         prev= x;
1073         x = next;
1074     }
1075     MSCRNode* root = x;
1076     for ( ; ; )
1077     {
1078         prev = x->shortcut;
1079         x->shortcut = root;
1080         if ( prev == x ) break;
1081         x = prev;
1082     }
1083     return root;
1084 }
1085
1086 // the stable mscr should be:
1087 // bigger than minArea and smaller than maxArea
1088 // differ from its ancestor more than minDiversity
1089 static bool MSCRStableCheck( MSCRNode* x, MSERParams params )
1090 {
1091     if ( x->size <= params.minArea || x->size >= params.maxArea )
1092         return 0;
1093     if ( x->gmsr == NULL )
1094         return 1;
1095     double div = (double)(x->size-x->gmsr->size)/(double)x->size;
1096     return div > params.minDiversity;
1097 }
1098
1099 static void
1100 extractMSER_8UC3( CvMat* src,
1101              CvMat* mask,
1102              CvSeq* contours,
1103              CvMemStorage* storage,
1104              MSERParams params )
1105 {
1106     MSCRNode* map = (MSCRNode*)cvAlloc( src->cols*src->rows*sizeof(map[0]) );
1107     int Ne = src->cols*src->rows*2-src->cols-src->rows;
1108     MSCREdge* edge = (MSCREdge*)cvAlloc( Ne*sizeof(edge[0]) );
1109     TempMSCR* mscr = (TempMSCR*)cvAlloc( src->cols*src->rows*sizeof(mscr[0]) );
1110     double emean = 0;
1111     CvMat* dx = cvCreateMat( src->rows, src->cols-1, CV_64FC1 );
1112     CvMat* dy = cvCreateMat( src->rows-1, src->cols, CV_64FC1 );
1113     Ne = preprocessMSER_8UC3( map, edge, &emean, src, mask, dx, dy, Ne, params.edgeBlurSize );
1114     emean = emean / (double)Ne;
1115     QuickSortMSCREdge( edge, Ne, 0 );
1116     MSCREdge* edge_ub = edge+Ne;
1117     MSCREdge* edgeptr = edge;
1118     TempMSCR* mscrptr = mscr;
1119     // the evolution process
1120     for ( int i = 0; i < params.maxEvolution; i++ )
1121     {
1122         double k = (double)i/(double)params.maxEvolution*(TABLE_SIZE-1);
1123         int ti = cvFloor(k);
1124         double reminder = k-ti;
1125         double thres = emean*(chitab3[ti]*(1-reminder)+chitab3[ti+1]*reminder);
1126         // to process all the edges in the list that chi < thres
1127         while ( edgeptr < edge_ub && edgeptr->chi < thres )
1128         {
1129             MSCRNode* lr = findMSCR( edgeptr->left );
1130             MSCRNode* rr = findMSCR( edgeptr->right );
1131             // get the region root (who is responsible)
1132             if ( lr != rr )
1133             {
1134                 // rank idea take from: N-tree Disjoint-Set Forests for Maximally Stable Extremal Regions
1135                 if ( rr->rank > lr->rank )
1136                 {
1137                     MSCRNode* tmp;
1138                     CV_SWAP( lr, rr, tmp );
1139                 } else if ( lr->rank == rr->rank ) {
1140                     // at the same rank, we will compare the size
1141                     if ( lr->size > rr->size )
1142                     {
1143                         MSCRNode* tmp;
1144                         CV_SWAP( lr, rr, tmp );
1145                     }
1146                     lr->rank++;
1147                 }
1148                 rr->shortcut = lr;
1149                 lr->size += rr->size;
1150                 // join rr to the end of list lr (lr is a endless double-linked list)
1151                 lr->prev->next = rr;
1152                 lr->prev = rr->prev;
1153                 rr->prev->next = lr;
1154                 rr->prev = lr;
1155                 // area threshold force to reinitialize
1156                 if ( lr->size > (lr->size-rr->size)*params.areaThreshold )
1157                 {
1158                     lr->sizei = lr->size;
1159                     lr->reinit = i;
1160                     if ( lr->tmsr != NULL )
1161                     {
1162                         lr->tmsr->m = lr->dt-lr->di;
1163                         lr->tmsr = NULL;
1164                     }
1165                     lr->di = edgeptr->chi;
1166                     lr->s = 1e10;
1167                 }
1168                 lr->dt = edgeptr->chi;
1169                 if ( i > lr->reinit )
1170                 {
1171                     double s = (double)(lr->size-lr->sizei)/(lr->dt-lr->di);
1172                     if ( s < lr->s )
1173                     {
1174                         // skip the first one and check stablity
1175                         if ( i > lr->reinit+1 && MSCRStableCheck( lr, params ) )
1176                         {
1177                             if ( lr->tmsr == NULL )
1178                             {
1179                                 lr->gmsr = lr->tmsr = mscrptr;
1180                                 mscrptr++;
1181                             }
1182                             lr->tmsr->size = lr->size;
1183                             lr->tmsr->head = lr;
1184                             lr->tmsr->tail = lr->prev;
1185                             lr->tmsr->m = 0;
1186                         }
1187                         lr->s = s;
1188                     }
1189                 }
1190             }
1191             edgeptr++;
1192         }
1193         if ( edgeptr >= edge_ub )
1194             break;
1195     }
1196     for ( TempMSCR* ptr = mscr; ptr < mscrptr; ptr++ )
1197         // to prune area with margin less than minMargin
1198         if ( ptr->m > params.minMargin )
1199         {
1200             CvSeq* _contour = cvCreateSeq( CV_SEQ_KIND_GENERIC|CV_32SC2, sizeof(CvContour), sizeof(CvPoint), storage );
1201             cvSeqPushMulti( _contour, 0, ptr->size );
1202             MSCRNode* lpt = ptr->head;
1203             for ( int i = 0; i < ptr->size; i++ )
1204             {
1205                 CvPoint* pt = CV_GET_SEQ_ELEM( CvPoint, _contour, i );
1206                 pt->x = (lpt->index)&0xffff;
1207                 pt->y = (lpt->index)>>16;
1208                 lpt = lpt->next;
1209             }
1210             CvContour* contour = (CvContour*)_contour;
1211             cvBoundingRect( contour );
1212             contour->color = 0;
1213             cvSeqPush( contours, &contour );
1214         }
1215     cvReleaseMat( &dx );
1216     cvReleaseMat( &dy );
1217     cvFree( &mscr );
1218     cvFree( &edge );
1219     cvFree( &map );
1220 }
1221
1222 static void
1223 extractMSER( CvArr* _img,
1224            CvArr* _mask,
1225            CvSeq** _contours,
1226            CvMemStorage* storage,
1227            MSERParams params )
1228 {
1229     CvMat srchdr, *src = cvGetMat( _img, &srchdr );
1230     CvMat maskhdr, *mask = _mask ? cvGetMat( _mask, &maskhdr ) : 0;
1231     CvSeq* contours = 0;
1232
1233     CV_Assert(src != 0);
1234     CV_Assert(CV_MAT_TYPE(src->type) == CV_8UC1 || CV_MAT_TYPE(src->type) == CV_8UC3);
1235     CV_Assert(mask == 0 || (CV_ARE_SIZES_EQ(src, mask) && CV_MAT_TYPE(mask->type) == CV_8UC1));
1236     CV_Assert(storage != 0);
1237
1238     contours = *_contours = cvCreateSeq( 0, sizeof(CvSeq), sizeof(CvSeq*), storage );
1239
1240     // choose different method for different image type
1241     // for grey image, it is: Linear Time Maximally Stable Extremal Regions
1242     // for color image, it is: Maximally Stable Colour Regions for Recognition and Matching
1243     switch ( CV_MAT_TYPE(src->type) )
1244     {
1245         case CV_8UC1:
1246             extractMSER_8UC1( src, mask, contours, storage, params );
1247             break;
1248         case CV_8UC3:
1249             extractMSER_8UC3( src, mask, contours, storage, params );
1250             break;
1251     }
1252 }
1253
1254
1255 MSER::MSER( int _delta, int _min_area, int _max_area,
1256       double _max_variation, double _min_diversity,
1257       int _max_evolution, double _area_threshold,
1258       double _min_margin, int _edge_blur_size )
1259     : delta(_delta), minArea(_min_area), maxArea(_max_area),
1260     maxVariation(_max_variation), minDiversity(_min_diversity),
1261     maxEvolution(_max_evolution), areaThreshold(_area_threshold),
1262     minMargin(_min_margin), edgeBlurSize(_edge_blur_size)
1263 {
1264 }
1265
1266 void MSER::operator()( const Mat& image, vector<vector<Point> >& dstcontours, const Mat& mask ) const
1267 {
1268     CvMat _image = image, _mask, *pmask = 0;
1269     if( mask.data )
1270         pmask = &(_mask = mask);
1271     MemStorage storage(cvCreateMemStorage(0));
1272     Seq<CvSeq*> contours;
1273     extractMSER( &_image, pmask, &contours.seq, storage,
1274                  MSERParams(delta, minArea, maxArea, maxVariation, minDiversity,
1275                             maxEvolution, areaThreshold, minMargin, edgeBlurSize));
1276     SeqIterator<CvSeq*> it = contours.begin();
1277     size_t i, ncontours = contours.size();
1278     dstcontours.resize(ncontours);
1279     for( i = 0; i < ncontours; i++, ++it )
1280         Seq<Point>(*it).copyTo(dstcontours[i]);
1281 }
1282
1283
1284 void MserFeatureDetector::detectImpl( const Mat& image, vector<KeyPoint>& keypoints, const Mat& mask ) const
1285 {
1286     vector<vector<Point> > msers;
1287
1288     (*this)(image, msers, mask);
1289
1290     vector<vector<Point> >::const_iterator contour_it = msers.begin();
1291     Rect r(0, 0, image.cols, image.rows);
1292     for( ; contour_it != msers.end(); ++contour_it )
1293     {
1294         // TODO check transformation from MSER region to KeyPoint
1295         RotatedRect rect = fitEllipse(Mat(*contour_it));
1296         float diam = sqrt(rect.size.height*rect.size.width);
1297
1298         if( diam > std::numeric_limits<float>::epsilon() && r.contains(rect.center) )
1299             keypoints.push_back( KeyPoint(rect.center, diam) );
1300     }
1301
1302 }
1303
1304 }