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