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