vcs.maemo.org Git - opencv/blob - src/cv/cvpyrsegmentation.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                        Intel License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
  14 // Third party copyrights are property of their respective owners.
  15 //
  16 // Redistribution and use in source and binary forms, with or without modification,
  17 // are permitted provided that the following conditions are met:
  18 //
  19 //   * Redistribution's of source code must retain the above copyright notice,
  20 //     this list of conditions and the following disclaimer.
  21 //
  22 //   * Redistribution's in binary form must reproduce the above copyright notice,
  23 //     this list of conditions and the following disclaimer in the documentation
  24 //     and/or other materials provided with the distribution.
  25 //
  26 //   * The name of Intel Corporation may not be used to endorse or promote products
  27 //     derived from this software without specific prior written permission.
  28 //
  29 // This software is provided by the copyright holders and contributors "as is" and
  30 // any express or implied warranties, including, but not limited to, the implied
  31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  32 // In no event shall the Intel Corporation or contributors be liable for any direct,
  33 // indirect, incidental, special, exemplary, or consequential damages
  34 // (including, but not limited to, procurement of substitute goods or services;
  35 // loss of use, data, or profits; or business interruption) however caused
  36 // and on any theory of liability, whether in contract, strict liability,
  37 // or tort (including negligence or otherwise) arising in any way out of
  38 // the use of this software, even if advised of the possibility of such damage.
  39 //
  40 //M*/
  41 #include "_cv.h"
  42
  43 typedef struct _CvRGBf
  44 {   float blue;
  45     float green;
  46     float red;
  47 }
  48 _CvRGBf;
  49
  50 typedef struct _CvRect16u
  51 {
  52     ushort x1, y1, x2, y2;
  53 }
  54 _CvRect16u;
  55
  56 typedef struct _CvPyramid
  57 {
  58     float c;
  59     struct _CvPyramid *p;
  60     int a;
  61     _CvRect16u rect;      /*  ROI for the connected component    */
  62 } _CvPyramid;
  63
  64 /* element of base layer */
  65 typedef struct _CvPyramidBase
  66 {
  67     float c;
  68     struct _CvPyramid *p;
  69 }
  70 _CvPyramidBase;
  71
  72 typedef struct _CvPyramidC3
  73 {
  74     _CvRGBf c;
  75     struct _CvPyramidC3 *p;
  76     int a;
  77     _CvRect16u rect;      /*  ROI for the connected component    */
  78 } _CvPyramidC3;
  79
  80 /* element of base layer */
  81 typedef struct _CvPyramidBaseC3
  82 {
  83     _CvRGBf c;
  84     struct _CvPyramidC3 *p;
  85 }
  86 _CvPyramidBaseC3;
  87
  88 typedef struct _CvListNode
  89 {
  90     struct _CvListNode* next;
  91     void* data;
  92 }
  93 _CvListNode;
  94
  95
  96 static CvStatus  icvSegmentClusterC1( CvSeq* cmp_seq, CvSeq* res_seq,
  97                                  double threshold,
  98                                  _CvPyramid* first_level_end,
  99                                  CvSize first_level_size );
 100
 101 static CvStatus  icvSegmentClusterC3( CvSeq* cmp_seq, CvSeq* res_seq,
 102                                  double threshold,
 103                                  _CvPyramidC3* first_level_end,
 104                                  CvSize first_level_size );
 105
 106 typedef void (CV_CDECL * CvWriteNodeFunction)(void* seq,void* node);
 107
 108 static CvStatus icvUpdatePyrLinks_8u_C1
 109     (int layer, void *layer_data, CvSize size, void *parent_layer,
 110      void *_writer, float threshold, int is_last_iter, void *_stub, CvWriteNodeFunction /*func*/);
 111
 112 static CvStatus icvUpdatePyrLinks_8u_C3
 113     (int layer, void *layer_data, CvSize size, void *parent_layer,
 114      void *_writer, float threshold, int is_last_iter, void *_stub, CvWriteNodeFunction /*func*/);
 115
 116 static void icvMaxRoi( _CvRect16u *max_rect, _CvRect16u* cur_rect );
 117 static void icvMaxRoi1( _CvRect16u *max_rect, int x, int y );
 118
 119
 120 #define _CV_CHECK( icvFun )                                             \
 121   {                                                                     \
 122     if( icvFun != CV_OK )                                               \
 123      goto M_END;                                                        \
 124   }
 125
 126
 127 #define _CV_MAX3( a, b, c) ((a)>(b) ? ((a)>(c) ? (a) : (c)) : ((b)>(c) ? (b) : (c)))
 128
 129 /*#define _CV_RGB_DIST(a, b)  _CV_MAX3((float)fabs((a).red - (b).red),      \
 130                                        (float)fabs((a).green - (b).green),  \
 131                                        (float)fabs((a).blue - (b).blue))*/
 132
 133 #define _CV_NEXT_BASE_C1(p,n) (_CvPyramid*)((char*)(p) + (n)*sizeof(_CvPyramidBase))
 134 #define _CV_NEXT_BASE_C3(p,n) (_CvPyramidC3*)((char*)(p) + (n)*sizeof(_CvPyramidBaseC3))
 135
 136
 137 CV_INLINE float icvRGBDist_Max( const _CvRGBf& a, const _CvRGBf& b )
 138 {
 139     float tr = (float)fabs(a.red - b.red);
 140     float tg = (float)fabs(a.green - b.green);
 141     float tb = (float)fabs(a.blue - b.blue);
 142
 143     return _CV_MAX3( tr, tg, tb );
 144 }
 145
 146 CV_INLINE float icvRGBDist_Sum( const _CvRGBf& a, const _CvRGBf& b )
 147 {
 148     float tr = (float)fabs(a.red - b.red);
 149     float tg = (float)fabs(a.green - b.green);
 150     float tb = (float)fabs(a.blue - b.blue);
 151
 152     return (tr + tg + tb);
 153 }
 154
 155 #if 1
 156 #define _CV_RGB_DIST  icvRGBDist_Max
 157 #define _CV_RGB_THRESH_SCALE   1
 158 #else
 159 #define _CV_RGB_DIST  icvRGBDist_Sum
 160 #define _CV_RGB_THRESH_SCALE   3
 161 #endif
 162
 163 #define _CV_INV_TAB_SIZE   32
 164
 165 static const float icvInvTab[ /*_CV_INV_TAB_SIZE*/ ] =
 166 {
 167     1.00000000f, 0.50000000f, 0.33333333f, 0.25000000f, 0.20000000f, 0.16666667f,
 168     0.14285714f, 0.12500000f, 0.11111111f, 0.10000000f, 0.09090909f, 0.08333333f,
 169     0.07692308f, 0.07142857f, 0.06666667f, 0.06250000f, 0.05882353f, 0.05555556f,
 170     0.05263158f, 0.05000000f, 0.04761905f, 0.04545455f, 0.04347826f, 0.04166667f,
 171     0.04000000f, 0.03846154f, 0.03703704f, 0.03571429f, 0.03448276f, 0.03333333f,
 172     0.03225806f, 0.03125000f
 173 };
 174
 175 static void
 176 icvWritePyrNode( void *elem, void *writer )
 177 {
 178     CV_WRITE_SEQ_ELEM( *(_CvListNode *) elem, *(CvSeqWriter *) writer );
 179 }
 180
 181
 182 static CvStatus
 183 icvPyrSegmentation8uC1R( uchar * src_image, int src_step,
 184                          uchar * dst_image, int dst_step,
 185                          CvSize roi, CvFilter filter,
 186                          CvSeq ** dst_comp, CvMemStorage * storage,
 187                          int level, int threshold1, int threshold2 )
 188 {
 189     int i, j, l;
 190     int step;
 191     const int max_iter = 3;     /* maximum number of iterations */
 192     int cur_iter = 0;           /* current iteration */
 193
 194     _CvPyramid *pyram[16];      /* pointers to the pyramid down up to level */
 195
 196     float *pyramida = 0;
 197     _CvPyramid stub;
 198
 199     _CvPyramid *p_cur;
 200     _CvPyramidBase *p_base;
 201     _CvListNode cmp_node;
 202
 203     CvSeq *cmp_seq = 0;
 204     CvSeq *res_seq = 0;
 205     CvMemStorage *temp_storage = 0;
 206     CvSize size;
 207     CvStatus status;
 208     CvSeqWriter writer;
 209
 210     int buffer_size;
 211     char *buffer = 0;
 212
 213     status = CV_OK;
 214
 215     /* clear pointer to resultant sequence */
 216     if( dst_comp )
 217         *dst_comp = 0;
 218
 219     /* check args */
 220     if( !src_image || !dst_image || !storage || !dst_comp )
 221         return CV_NULLPTR_ERR;
 222     if( roi.width <= 0 || roi.height <= 0 || src_step < roi.width || dst_step < roi.width )
 223         return CV_BADSIZE_ERR;
 224     if( filter != CV_GAUSSIAN_5x5 )
 225         return CV_BADRANGE_ERR;
 226     if( threshold1 < 0 || threshold2 < 0 )
 227         return CV_BADRANGE_ERR;
 228     if( level <= 0 )
 229         return CV_BADRANGE_ERR;
 230
 231     if( ((roi.width | roi.height) & ((1 << level) - 1)) != 0 )
 232         return CV_BADCOEF_ERR;
 233
 234     temp_storage = cvCreateChildMemStorage( storage );
 235
 236     /* sequence for temporary components */
 237     cmp_seq = cvCreateSeq( 0, sizeof( CvSeq ), sizeof( _CvListNode ), temp_storage );
 238     assert( cmp_seq != 0 );
 239
 240     res_seq = cvCreateSeq( CV_SEQ_CONNECTED_COMP, sizeof( CvSeq ),
 241                            sizeof( CvConnectedComp ), storage );
 242     assert( res_seq != 0 );
 243
 244     /* calculate buffer size */
 245     buffer_size = roi.width * roi.height * (sizeof( float ) + sizeof( _CvPyramidBase ));
 246
 247     for( l = 1; l <= level; l++ )
 248         buffer_size += ((roi.width >> l) + 1) * ((roi.height >> l) + 1) * sizeof(_CvPyramid);
 249
 250     /* allocate buffer */
 251     buffer = (char *) cvAlloc( buffer_size );
 252     if( !buffer )
 253     {
 254         status = CV_OUTOFMEM_ERR;
 255         goto M_END;
 256     }
 257
 258     pyramida = (float *) buffer;
 259
 260     /* initialization pyramid-linking properties down up to level */
 261     step = roi.width * sizeof( float );
 262
 263     {
 264         CvMat _src;
 265         CvMat _pyramida;
 266         cvInitMatHeader( &_src, roi.height, roi.width, CV_8UC1, src_image, src_step );
 267         cvInitMatHeader( &_pyramida, roi.height, roi.width, CV_32FC1, pyramida, step );
 268         cvConvert( &_src, &_pyramida );
 269         /*_CV_CHECK( icvCvtTo_32f_C1R( src_image, src_step, pyramida, step, roi, CV_8UC1 ));*/
 270     }
 271     p_base = (_CvPyramidBase *) (buffer + step * roi.height);
 272     pyram[0] = (_CvPyramid *) p_base;
 273
 274     /* fill base level of pyramid */
 275     for( i = 0; i < roi.height; i++ )
 276     {
 277         for( j = 0; j < roi.width; j++, p_base++ )
 278         {
 279             p_base->c = pyramida[i * roi.width + j];
 280             p_base->p = &stub;
 281         }
 282     }
 283
 284     p_cur = (_CvPyramid *) p_base;
 285     size = roi;
 286
 287     /* calculate initial pyramid */
 288     for( l = 1; l <= level; l++ )
 289     {
 290         CvSize dst_size = { size.width/2+1, size.height/2+1 };
 291         CvMat prev_level = cvMat( size.height, size.width, CV_32FC1 );
 292         CvMat next_level = cvMat( dst_size.height, dst_size.width, CV_32FC1 );
 293
 294         cvSetData( &prev_level, pyramida, step );
 295         cvSetData( &next_level, pyramida, step );
 296         cvPyrDown( &prev_level, &next_level );
 297
 298         //_CV_CHECK( icvPyrDown_Gauss5x5_32f_C1R( pyramida, step, pyramida, step, size, buff ));
 299         //_CV_CHECK( icvPyrDownBorder_32f_CnR( pyramida, step, size, pyramida, step, dst_size, 1 ));
 300         pyram[l] = p_cur;
 301
 302         size.width = dst_size.width - 1;
 303         size.height = dst_size.height - 1;
 304
 305         /* fill layer #l */
 306         for( i = 0; i <= size.height; i++ )
 307         {
 308             for( j = 0; j <= size.width; j++, p_cur++ )
 309             {
 310                 p_cur->c = pyramida[i * roi.width + j];
 311                 p_cur->p = &stub;
 312                 p_cur->a = 0;
 313                 p_cur->rect.x2 = 0;
 314             }
 315         }
 316     }
 317
 318     cvStartAppendToSeq( cmp_seq, &writer );
 319
 320     /* do several iterations to determine son-father links */
 321     for( cur_iter = 0; cur_iter < max_iter; cur_iter++ )
 322     {
 323         int is_last_iter = cur_iter == max_iter - 1;
 324
 325         size = roi;
 326
 327         /* build son-father links down up to level */
 328         for( l = 0; l < level; l++ )
 329         {
 330             icvUpdatePyrLinks_8u_C1( l, pyram[l], size, pyram[l + 1], &writer,
 331                                       (float) threshold1, is_last_iter, &stub,
 332                                       icvWritePyrNode );
 333
 334             /* clear last border row */
 335             if( l > 0 )
 336             {
 337                 p_cur = pyram[l] + (size.width + 1) * size.height;
 338                 for( j = 0; j <= size.width; j++ )
 339                     p_cur[j].c = 0;
 340             }
 341
 342             size.width >>= 1;
 343             size.height >>= 1;
 344         }
 345
 346 /*  clear the old c value for the last level     */
 347         p_cur = pyram[level];
 348         for( i = 0; i <= size.height; i++, p_cur += size.width + 1 )
 349             for( j = 0; j <= size.width; j++ )
 350                 p_cur[j].c = 0;
 351
 352         size = roi;
 353         step = roi.width;
 354
 355 /* calculate average c value for the 0 < l <=level   */
 356         for( l = 0; l < level; l++, step = (step >> 1) + 1 )
 357         {
 358             _CvPyramid *p_prev, *p_row_prev;
 359
 360             stub.c = 0;
 361
 362             /* calculate average c value for the next level   */
 363             if( l == 0 )
 364             {
 365                 p_base = (_CvPyramidBase *) pyram[0];
 366                 for( i = 0; i < roi.height; i++, p_base += size.width )
 367                 {
 368                     for( j = 0; j < size.width; j += 2 )
 369                     {
 370                         _CvPyramid *p1 = p_base[j].p;
 371                         _CvPyramid *p2 = p_base[j + 1].p;
 372
 373                         p1->c += p_base[j].c;
 374                         p2->c += p_base[j + 1].c;
 375                     }
 376                 }
 377             }
 378             else
 379             {
 380                 p_cur = pyram[l];
 381                 for( i = 0; i < size.height; i++, p_cur += size.width + 1 )
 382                 {
 383                     for( j = 0; j < size.width; j += 2 )
 384                     {
 385                         _CvPyramid *p1 = p_cur[j].p;
 386                         _CvPyramid *p2 = p_cur[j + 1].p;
 387
 388                         float t0 = (float) p_cur[j].a * p_cur[j].c;
 389                         float t1 = (float) p_cur[j + 1].a * p_cur[j + 1].c;
 390
 391                         p1->c += t0;
 392                         p2->c += t1;
 393
 394                         if( !is_last_iter )
 395                             p_cur[j].a = p_cur[j + 1].a = 0;
 396                     }
 397                     if( !is_last_iter )
 398                         p_cur[size.width].a = 0;
 399                 }
 400                 if( !is_last_iter )
 401                 {
 402                     for( j = 0; j <= size.width; j++ )
 403                     {
 404                         p_cur[j].a = 0;
 405                     }
 406                 }
 407             }
 408
 409             /* assign random values of the next level null c   */
 410             p_cur = pyram[l + 1];
 411             p_row_prev = p_prev = pyram[l];
 412
 413             size.width >>= 1;
 414             size.height >>= 1;
 415
 416             for( i = 0; i <= size.height; i++, p_cur += size.width + 1 )
 417             {
 418                 if( i < size.height || !is_last_iter )
 419                 {
 420                     for( j = 0; j < size.width; j++ )
 421                     {
 422                         int a = p_cur[j].a;
 423
 424                         if( a != 0 )
 425                         {
 426                             if( a <= _CV_INV_TAB_SIZE )
 427                             {
 428                                 p_cur[j].c *= icvInvTab[a - 1];
 429                             }
 430                             else
 431                             {
 432                                 p_cur[j].c /= a;
 433                             }
 434                         }
 435                         else
 436                         {
 437                             p_cur[j].c = p_prev->c;
 438                         }
 439
 440                         if( l == 0 )
 441                             p_prev = _CV_NEXT_BASE_C1(p_prev,2);
 442                         else
 443                             p_prev += 2;
 444                     }
 445
 446                     if( p_cur[size.width].a == 0 )
 447                     {
 448                         p_cur[size.width].c = p_prev[(l != 0) - 1].c;
 449                     }
 450                     else
 451                     {
 452                         p_cur[size.width].c /= p_cur[size.width].a;
 453                         if( is_last_iter )
 454                         {
 455                             cmp_node.data = p_cur + size.width;
 456                             CV_WRITE_SEQ_ELEM( cmp_node, writer );
 457                         }
 458                     }
 459                 }
 460                 else
 461                 {
 462                     for( j = 0; j <= size.width; j++ )
 463                     {
 464                         int a = p_cur[j].a;
 465
 466                         if( a != 0 )
 467                         {
 468                             if( a <= _CV_INV_TAB_SIZE )
 469                             {
 470                                 p_cur[j].c *= icvInvTab[a - 1];
 471                             }
 472                             else
 473                             {
 474                                 p_cur[j].c /= a;
 475                             }
 476
 477                             cmp_node.data = p_cur + j;
 478                             CV_WRITE_SEQ_ELEM( cmp_node, writer );
 479                         }
 480                         else
 481                         {
 482                             p_cur[j].c = p_prev->c;
 483                         }
 484
 485                         if( l == 0 )
 486                         {
 487                             p_prev = _CV_NEXT_BASE_C1(p_prev, (j * 2 < step - 2 ? 2 : 1));
 488                         }
 489                         else
 490                         {
 491                             p_prev++;
 492                         }
 493                     }
 494                 }
 495
 496                 if( l + 1 == level && !is_last_iter )
 497                     for( j = 0; j <= size.width; j++ )
 498                         p_cur[j].a = 0;
 499
 500                 if( !(i & 1) )
 501                 {
 502                     p_prev = p_row_prev;
 503                 }
 504                 else
 505                 {
 506                     p_prev = (_CvPyramid*)((char*)p_row_prev + step *
 507                         (l == 0 ? sizeof(_CvPyramidBase) : sizeof(_CvPyramid)));
 508                 }
 509             }
 510         }
 511     }                           /*  end of the iteration process  */
 512
 513     /* construct a connected  components   */
 514     size.width = roi.width >> level;
 515     size.height = roi.height >> level;
 516
 517     p_cur = pyram[level];
 518
 519     for( i = 0; i < size.height; i++, p_cur += size.width + 1 )
 520     {
 521         for( j = 0; j < size.width; j++ )
 522         {
 523             if( p_cur[j].a != 0 )
 524             {
 525                 cmp_node.data = p_cur + j;
 526                 CV_WRITE_SEQ_ELEM( cmp_node, writer );
 527             }
 528         }
 529     }
 530
 531     cvEndWriteSeq( &writer );
 532
 533 /* clusterization segmented components and construction
 534    output connected components                            */
 535     icvSegmentClusterC1( cmp_seq, res_seq, threshold2, pyram[1], roi );
 536
 537 /* convert (inplace) resultant segment values to int (top level) */
 538
 539 /* propagate segment values top down */
 540     for( l = level - 1; l >= 0; l-- )
 541     {
 542         p_cur = pyram[l];
 543
 544         size.width <<= 1;
 545         size.height <<= 1;
 546
 547         if( l == 0 )
 548         {
 549             size.width--;
 550             size.height--;
 551         }
 552
 553         for( i = 0; i <= size.height; i++ )
 554         {
 555             for( j = 0; j <= size.width; j++ )
 556             {
 557                 _CvPyramid *p = p_cur->p;
 558
 559                 assert( p != 0 );
 560                 if( p != &stub )
 561                     p_cur->c = p->c;
 562
 563                 if( l == 0 )
 564                 {
 565                     Cv32suf _c;
 566                     /* copy the segmented values to destination image */
 567                     _c.f = p_cur->c; dst_image[j] = (uchar)_c.i;
 568                     p_cur = _CV_NEXT_BASE_C1(p_cur, 1);
 569                 }
 570                 else
 571                 {
 572                     p_cur++;
 573                 }
 574             }
 575             if( l == 0 )
 576                 dst_image += dst_step;
 577         }
 578     }
 579   M_END:
 580
 581     cvFree( &buffer );
 582     cvReleaseMemStorage( &temp_storage );
 583
 584     if( status == CV_OK )
 585         *dst_comp = res_seq;
 586
 587     return status;
 588 }
 589
 590
 591
 592 /****************************************************************************************\
 593     color!!!  image segmentation by pyramid-linking
 594 \****************************************************************************************/
 595 static CvStatus
 596 icvPyrSegmentation8uC3R( uchar * src_image, int src_step,
 597                          uchar * dst_image, int dst_step,
 598                          CvSize roi, CvFilter filter,
 599                          CvSeq ** dst_comp, CvMemStorage * storage,
 600                          int level, int threshold1, int threshold2 )
 601 {
 602     int i, j, l;
 603
 604     int step;
 605     const int max_iter = 3;     /* maximum number of iterations */
 606     int cur_iter = 0;           /* current iteration */
 607
 608     _CvPyramidC3 *pyram[16];    /* pointers to the pyramid down up to level */
 609
 610     float *pyramida = 0;
 611     _CvPyramidC3 stub;
 612
 613     _CvPyramidC3 *p_cur;
 614     _CvPyramidBaseC3 *p_base;
 615     _CvListNode cmp_node;
 616
 617     CvSeq *cmp_seq = 0;
 618     CvSeq *res_seq = 0;
 619     CvMemStorage *temp_storage = 0;
 620     CvSize size;
 621     CvStatus status;
 622     CvSeqWriter writer;
 623
 624     int buffer_size;
 625     char *buffer = 0;
 626
 627     status = CV_OK;
 628
 629     threshold1 *= _CV_RGB_THRESH_SCALE;
 630     threshold2 *= _CV_RGB_THRESH_SCALE;
 631
 632     /* clear pointer to resultant sequence */
 633     if( dst_comp )
 634         *dst_comp = 0;
 635
 636     /* check args */
 637     if( !src_image || !dst_image || !storage || !dst_comp )
 638         return CV_NULLPTR_ERR;
 639     if( roi.width <= 0 || roi.height <= 0 ||
 640         src_step < roi.width * 3 || dst_step < roi.width * 3 ) return CV_BADSIZE_ERR;
 641     if( filter != CV_GAUSSIAN_5x5 )
 642         return CV_BADRANGE_ERR;
 643     if( threshold1 < 0 || threshold2 < 0 )
 644         return CV_BADRANGE_ERR;
 645     if( level <= 0 )
 646         return CV_BADRANGE_ERR;
 647
 648     if( ((roi.width | roi.height) & ((1 << level) - 1)) != 0 )
 649         return CV_BADCOEF_ERR;
 650
 651     temp_storage = cvCreateChildMemStorage( storage );
 652
 653     /* sequence for temporary components */
 654     cmp_seq = cvCreateSeq( 0, sizeof( CvSeq ), sizeof( _CvListNode ), temp_storage );
 655     assert( cmp_seq != 0 );
 656
 657     res_seq = cvCreateSeq( CV_SEQ_CONNECTED_COMP, sizeof( CvSeq ),
 658                            sizeof( CvConnectedComp ), storage );
 659     assert( res_seq != 0 );
 660
 661     /* calculate buffer size */
 662     buffer_size = roi.width * roi.height * (sizeof( _CvRGBf ) + sizeof( _CvPyramidBaseC3 ));
 663
 664     for( l = 1; l <= level; l++ )
 665         buffer_size += ((roi.width >> l) + 1) * ((roi.height >> l) + 1) * sizeof(_CvPyramidC3);
 666
 667     /* allocate buffer */
 668     buffer = (char *) cvAlloc( buffer_size );
 669     if( !buffer )
 670     {
 671         status = CV_OUTOFMEM_ERR;
 672         goto M_END;
 673     }
 674
 675     pyramida = (float *) buffer;
 676
 677     /* initialization pyramid-linking properties down up to level */
 678     step = roi.width * sizeof( _CvRGBf );
 679
 680     {
 681         CvMat _src;
 682         CvMat _pyramida;
 683         cvInitMatHeader( &_src, roi.height, roi.width, CV_8UC3, src_image, src_step );
 684         cvInitMatHeader( &_pyramida, roi.height, roi.width, CV_32FC3, pyramida, step );
 685         cvConvert( &_src, &_pyramida );
 686         /*_CV_CHECK( icvCvtTo_32f_C1R( src_image, src_step, pyramida, step,
 687                                  cvSize( roi.width * 3, roi.height ), CV_8UC1 ));*/
 688     }
 689
 690     p_base = (_CvPyramidBaseC3 *) (buffer + step * roi.height);
 691     pyram[0] = (_CvPyramidC3 *) p_base;
 692
 693     /* fill base level of pyramid */
 694     for( i = 0; i < roi.height; i++ )
 695     {
 696         for( j = 0; j < roi.width; j++, p_base++ )
 697         {
 698             p_base->c = ((_CvRGBf *) pyramida)[i * roi.width + j];
 699             p_base->p = &stub;
 700         }
 701     }
 702
 703     p_cur = (_CvPyramidC3 *) p_base;
 704     size = roi;
 705
 706     /* calculate initial pyramid */
 707     for( l = 1; l <= level; l++ )
 708     {
 709         CvSize dst_size = { size.width/2 + 1, size.height/2 + 1 };
 710         CvMat prev_level = cvMat( size.height, size.width, CV_32FC3 );
 711         CvMat next_level = cvMat( dst_size.height, dst_size.width, CV_32FC3 );
 712
 713         cvSetData( &prev_level, pyramida, step );
 714         cvSetData( &next_level, pyramida, step );
 715         cvPyrDown( &prev_level, &next_level );
 716
 717         //_CV_CHECK( icvPyrDown_Gauss5x5_32f_C3R( pyramida, step, pyramida, step, size, buff ));
 718         //_CV_CHECK( icvPyrDownBorder_32f_CnR( pyramida, step, size, pyramida, step, dst_size, 3 ));
 719         pyram[l] = p_cur;
 720
 721         size.width = dst_size.width - 1;
 722         size.height = dst_size.height - 1;
 723
 724         /* fill layer #l */
 725         for( i = 0; i <= size.height; i++ )
 726         {
 727             assert( (char*)p_cur - buffer < buffer_size );
 728             for( j = 0; j <= size.width; j++, p_cur++ )
 729             {
 730                 p_cur->c = ((_CvRGBf *) pyramida)[i * roi.width + j];
 731                 p_cur->p = &stub;
 732                 p_cur->a = 0;
 733                 p_cur->rect.x2 = 0;
 734             }
 735         }
 736     }
 737
 738     cvStartAppendToSeq( cmp_seq, &writer );
 739
 740     /* do several iterations to determine son-father links */
 741     for( cur_iter = 0; cur_iter < max_iter; cur_iter++ )
 742     {
 743         int is_last_iter = cur_iter == max_iter - 1;
 744
 745         size = roi;
 746
 747         /* build son-father links down up to level */
 748         for( l = 0; l < level; l++ )
 749         {
 750             icvUpdatePyrLinks_8u_C3( l, pyram[l], size, pyram[l + 1], &writer,
 751                                       (float) threshold1, is_last_iter, &stub,
 752                                       icvWritePyrNode );
 753
 754             /* clear last border row */
 755             if( l > 0 )
 756             {
 757                 p_cur = pyram[l] + (size.width + 1) * size.height;
 758                 for( j = 0; j <= size.width; j++ )
 759                     p_cur[j].c.blue = p_cur[j].c.green = p_cur[j].c.red = 0;
 760             }
 761
 762             size.width >>= 1;
 763             size.height >>= 1;
 764         }
 765
 766 /*  clear the old c value for the last level     */
 767         p_cur = pyram[level];
 768         for( i = 0; i <= size.height; i++, p_cur += size.width + 1 )
 769             for( j = 0; j <= size.width; j++ )
 770                 p_cur[j].c.blue = p_cur[j].c.green = p_cur[j].c.red = 0;
 771
 772         size = roi;
 773         step = roi.width;
 774
 775 /* calculate average c value for the 0 < l <=level   */
 776         for( l = 0; l < level; l++, step = (step >> 1) + 1 )
 777         {
 778             _CvPyramidC3 *p_prev, *p_row_prev;
 779
 780             stub.c.blue = stub.c.green = stub.c.red = 0;
 781
 782             /* calculate average c value for the next level   */
 783             if( l == 0 )
 784             {
 785                 p_base = (_CvPyramidBaseC3 *) pyram[0];
 786                 for( i = 0; i < roi.height; i++, p_base += size.width )
 787                 {
 788                     for( j = 0; j < size.width; j++ )
 789                     {
 790                         _CvPyramidC3 *p = p_base[j].p;
 791
 792                         p->c.blue += p_base[j].c.blue;
 793                         p->c.green += p_base[j].c.green;
 794                         p->c.red += p_base[j].c.red;
 795                     }
 796                 }
 797             }
 798             else
 799             {
 800                 p_cur = pyram[l];
 801                 for( i = 0; i < size.height; i++, p_cur += size.width + 1 )
 802                 {
 803                     for( j = 0; j < size.width; j++ )
 804                     {
 805                         _CvPyramidC3 *p = p_cur[j].p;
 806                         float a = (float) p_cur[j].a;
 807
 808                         p->c.blue += a * p_cur[j].c.blue;
 809                         p->c.green += a * p_cur[j].c.green;
 810                         p->c.red += a * p_cur[j].c.red;
 811
 812                         if( !is_last_iter )
 813                             p_cur[j].a = 0;
 814                     }
 815                     if( !is_last_iter )
 816                         p_cur[size.width].a = 0;
 817                 }
 818                 if( !is_last_iter )
 819                 {
 820                     for( j = 0; j <= size.width; j++ )
 821                     {
 822                         p_cur[j].a = 0;
 823                     }
 824                 }
 825             }
 826
 827             /* assign random values of the next level null c   */
 828             p_cur = pyram[l + 1];
 829             p_row_prev = p_prev = pyram[l];
 830
 831             size.width >>= 1;
 832             size.height >>= 1;
 833
 834             for( i = 0; i <= size.height; i++, p_cur += size.width + 1 )
 835             {
 836                 if( i < size.height || !is_last_iter )
 837                 {
 838                     for( j = 0; j < size.width; j++ )
 839                     {
 840                         int a = p_cur[j].a;
 841
 842                         if( a != 0 )
 843                         {
 844                             float inv_a;
 845
 846                             if( a <= _CV_INV_TAB_SIZE )
 847                             {
 848                                 inv_a = icvInvTab[a - 1];
 849                             }
 850                             else
 851                             {
 852                                 inv_a = 1.f / a;
 853                             }
 854                             p_cur[j].c.blue *= inv_a;
 855                             p_cur[j].c.green *= inv_a;
 856                             p_cur[j].c.red *= inv_a;
 857                         }
 858                         else
 859                         {
 860                             p_cur[j].c = p_prev->c;
 861                         }
 862
 863                         if( l == 0 )
 864                             p_prev = _CV_NEXT_BASE_C3( p_prev, 2 );
 865                         else
 866                             p_prev += 2;
 867                     }
 868
 869                     if( p_cur[size.width].a == 0 )
 870                     {
 871                         p_cur[size.width].c = p_prev[(l != 0) - 1].c;
 872                     }
 873                     else
 874                     {
 875                         p_cur[size.width].c.blue /= p_cur[size.width].a;
 876                         p_cur[size.width].c.green /= p_cur[size.width].a;
 877                         p_cur[size.width].c.red /= p_cur[size.width].a;
 878                         if( is_last_iter )
 879                         {
 880                             cmp_node.data = p_cur + size.width;
 881                             CV_WRITE_SEQ_ELEM( cmp_node, writer );
 882                         }
 883                     }
 884                 }
 885                 else
 886                 {
 887                     for( j = 0; j <= size.width; j++ )
 888                     {
 889                         int a = p_cur[j].a;
 890
 891                         if( a != 0 )
 892                         {
 893                             float inv_a;
 894
 895                             if( a <= _CV_INV_TAB_SIZE )
 896                             {
 897                                 inv_a = icvInvTab[a - 1];
 898                             }
 899                             else
 900                             {
 901                                 inv_a = 1.f / a;
 902                             }
 903                             p_cur[j].c.blue *= inv_a;
 904                             p_cur[j].c.green *= inv_a;
 905                             p_cur[j].c.red *= inv_a;
 906
 907                             cmp_node.data = p_cur + j;
 908                             CV_WRITE_SEQ_ELEM( cmp_node, writer );
 909                         }
 910                         else
 911                         {
 912                             p_cur[j].c = p_prev->c;
 913                         }
 914
 915                         if( l == 0 )
 916                         {
 917                             p_prev = _CV_NEXT_BASE_C3( p_prev, (j * 2 < step - 2 ? 2 : 1));
 918                         }
 919                         else
 920                         {
 921                             p_prev++;
 922                         }
 923                     }
 924                 }
 925
 926                 if( l + 1 == level && !is_last_iter )
 927                     for( j = 0; j <= size.width; j++ )
 928                         p_cur[j].a = 0;
 929
 930                 if( !(i & 1) )
 931                 {
 932                     p_prev = p_row_prev;
 933                 }
 934                 else
 935                 {
 936                     p_prev = (_CvPyramidC3*)((char*)p_row_prev + step *
 937                         (l == 0 ? sizeof( _CvPyramidBaseC3 ) : sizeof( _CvPyramidC3 )));
 938                 }
 939             }
 940         }
 941     }                           /*  end of the iteration process  */
 942
 943     /* construct a connected  components   */
 944     size.width = roi.width >> level;
 945     size.height = roi.height >> level;
 946
 947     p_cur = pyram[level];
 948
 949     for( i = 0; i < size.height; i++, p_cur += size.width + 1 )
 950     {
 951         for( j = 0; j < size.width; j++ )
 952         {
 953             if( p_cur[j].a != 0 )
 954             {
 955                 cmp_node.data = p_cur + j;
 956                 CV_WRITE_SEQ_ELEM( cmp_node, writer );
 957             }
 958         }
 959     }
 960
 961     cvEndWriteSeq( &writer );
 962
 963 /* clusterization segmented components and construction
 964    output connected components                            */
 965     icvSegmentClusterC3( cmp_seq, res_seq, threshold2, pyram[1], roi );
 966
 967 /* convert (inplace) resultant segment values to int (top level) */
 968
 969 /* propagate segment values top down */
 970     for( l = level - 1; l >= 0; l-- )
 971     {
 972         p_cur = pyram[l];
 973
 974         size.width <<= 1;
 975         size.height <<= 1;
 976
 977         if( l == 0 )
 978         {
 979             size.width--;
 980             size.height--;
 981         }
 982
 983         for( i = 0; i <= size.height; i++ )
 984         {
 985             for( j = 0; j <= size.width; j++ )
 986             {
 987                 _CvPyramidC3 *p = p_cur->p;
 988
 989                 assert( p != 0 );
 990                 if( p != &stub )
 991                 {
 992                     p_cur->c = p->c;
 993                 }
 994
 995                 if( l == 0 )
 996                 {
 997                     Cv32suf _c;
 998                     /* copy the segmented values to destination image */
 999                     _c.f = p_cur->c.blue; dst_image[j*3] = (uchar)_c.i;
1000                     _c.f = p_cur->c.green; dst_image[j*3+1] = (uchar)_c.i;
1001                     _c.f = p_cur->c.red; dst_image[j*3+2] = (uchar)_c.i;
1002                     p_cur = _CV_NEXT_BASE_C3(p_cur,1);
1003                 }
1004                 else
1005                 {
1006                     p_cur++;
1007                 }
1008             }
1009             if( l == 0 )
1010                 dst_image += dst_step;
1011         }
1012     }
1013
1014   M_END:
1015
1016     cvFree( &buffer );
1017     cvReleaseMemStorage( &temp_storage );
1018
1019     if( status == CV_OK )
1020         *dst_comp = res_seq;
1021
1022     return status;
1023 }
1024
1025
1026 static CvStatus icvUpdatePyrLinks_8u_C1
1027     (int layer, void *layer_data, CvSize size, void *parent_layer,
1028      void *_writer, float threshold, int is_last_iter, void *_stub, CvWriteNodeFunction /*func*/)
1029 {
1030     int i, j;
1031     _CvListNode cmp_node;
1032
1033     _CvPyramid *stub = (_CvPyramid *) _stub;
1034     _CvPyramid *p_cur = (_CvPyramid *) layer_data;
1035     _CvPyramid *p_next1 = (_CvPyramid *) parent_layer;
1036     _CvPyramid *p_next3 = p_next1 + (size.width >> 1) + 1;
1037
1038     CvSeqWriter & writer = *(CvSeqWriter *) _writer;
1039
1040     for( i = 0; i < size.height; i++ )
1041     {
1042         for( j = 0; j < size.width; j += 2 )
1043         {
1044             float c0, c1, c2, c3, c4;
1045             _CvPyramid *p;
1046
1047 /* son-father threshold linking for the current node establish */
1048             c0 = p_cur->c;
1049
1050 /* find pointer for the first pixel */
1051             c1 = (float) fabs( c0 - p_next1[0].c );
1052             c2 = (float) fabs( c0 - p_next1[1].c );
1053             c3 = (float) fabs( c0 - p_next3[0].c );
1054             c4 = (float) fabs( c0 - p_next3[1].c );
1055
1056             p = p_next1;
1057
1058             if( c1 > c2 )
1059             {
1060                 p = p_next1 + 1;
1061                 c1 = c2;
1062             }
1063             if( c1 > c3 )
1064             {
1065                 p = p_next3;
1066                 c1 = c3;
1067             }
1068             if( c1 > c4 )
1069             {
1070                 p = p_next3 + 1;
1071                 c1 = c4;
1072             }
1073
1074             if( c1 <= threshold )
1075             {
1076                 p_cur->p = p;
1077
1078                 if( layer == 0 )
1079                 {
1080                     p->a++;
1081                     p_cur = (_CvPyramid*)((char*)p_cur + sizeof(_CvPyramidBase));
1082                     if( is_last_iter )
1083                         icvMaxRoi1( &(p->rect), j, i );
1084                 }
1085                 else
1086                 {
1087                     int a = p_cur->a;
1088
1089                     p->a += a;
1090                     p_cur->c = 0;
1091                     p_cur++;
1092                     if( is_last_iter && a != 0 )
1093                         icvMaxRoi( &(p->rect), &(p_cur[-1].rect) );
1094                 }
1095             }
1096             else
1097             {
1098                 p_cur->p = stub;
1099                 if( is_last_iter )
1100                 {
1101                     cmp_node.data = p_cur;
1102                     CV_WRITE_SEQ_ELEM( cmp_node, writer );
1103                 }
1104                 if( layer == 0 )
1105                 {
1106                     p_cur = _CV_NEXT_BASE_C1(p_cur,1);
1107                 }
1108                 else
1109                 {
1110                     p_cur->c = 0;
1111                     p_cur++;
1112                 }
1113             }
1114
1115             /* find pointer for the second pixel */
1116             c0 = p_cur->c;
1117
1118             c1 = (float) fabs( c0 - p_next1[0].c );
1119             c2 = (float) fabs( c0 - p_next1[1].c );
1120             c3 = (float) fabs( c0 - p_next3[0].c );
1121             c4 = (float) fabs( c0 - p_next3[1].c );
1122
1123             p = p_next1;
1124             p_next1++;
1125
1126             if( c1 > c2 )
1127             {
1128                 p = p_next1;
1129                 c1 = c2;
1130             }
1131             if( c1 > c3 )
1132             {
1133                 p = p_next3;
1134                 c1 = c3;
1135             }
1136
1137             p_next3++;
1138             if( c1 > c4 )
1139             {
1140                 p = p_next3;
1141                 c1 = c4;
1142             }
1143
1144             if( c1 <= threshold )
1145             {
1146                 p_cur->p = p;
1147
1148                 if( layer == 0 )
1149                 {
1150                     p->a++;
1151                     p_cur = _CV_NEXT_BASE_C1(p_cur,1);
1152                     if( is_last_iter )
1153                         icvMaxRoi1( &(p->rect), j + 1, i );
1154                 }
1155                 else
1156                 {
1157                     int a = p_cur->a;
1158
1159                     p->a += a;
1160                     p_cur->c = 0;
1161                     p_cur++;
1162                     if( is_last_iter && a != 0 )
1163                         icvMaxRoi( &(p->rect), &(p_cur[-1].rect) );
1164                 }
1165             }
1166             else
1167             {
1168                 p_cur->p = stub;
1169                 if( is_last_iter )
1170                 {
1171                     cmp_node.data = p_cur;
1172                     CV_WRITE_SEQ_ELEM( cmp_node, writer );
1173                 }
1174                 if( layer == 0 )
1175                 {
1176                     p_cur = _CV_NEXT_BASE_C1(p_cur,1);
1177                 }
1178                 else
1179                 {
1180                     p_cur->c = 0;
1181                     p_cur++;
1182                 }
1183             }
1184         }
1185
1186         /* clear c's */
1187         if( layer > 0 )
1188         {
1189             p_cur->c = 0;
1190             p_cur++;
1191         }
1192
1193         if( !(i & 1) )
1194         {
1195             p_next1 -= size.width >> 1;
1196             p_next3 -= size.width >> 1;
1197         }
1198         else
1199         {
1200             p_next1++;
1201             p_next3++;
1202         }
1203     }
1204
1205     return CV_OK;
1206 }
1207
1208
1209 static CvStatus icvUpdatePyrLinks_8u_C3
1210     (int layer, void *layer_data, CvSize size, void *parent_layer,
1211      void *_writer, float threshold, int is_last_iter, void *_stub, CvWriteNodeFunction /*func*/)
1212 {
1213     int i, j;
1214     _CvListNode cmp_node;
1215
1216     _CvPyramidC3 *stub = (_CvPyramidC3 *) _stub;
1217     _CvPyramidC3 *p_cur = (_CvPyramidC3 *) layer_data;
1218     _CvPyramidC3 *p_next1 = (_CvPyramidC3 *) parent_layer;
1219     _CvPyramidC3 *p_next3 = p_next1 + (size.width >> 1) + 1;
1220
1221     CvSeqWriter & writer = *(CvSeqWriter *) _writer;
1222
1223     for( i = 0; i < size.height; i++ )
1224     {
1225         for( j = 0; j < size.width; j += 2 )
1226         {
1227             float c1, c2, c3, c4;
1228             _CvPyramidC3 *p;
1229
1230 /* find pointer for the first pixel */
1231             c1 = _CV_RGB_DIST( p_cur->c, p_next1[0].c );
1232             c2 = _CV_RGB_DIST( p_cur->c, p_next1[1].c );
1233             c3 = _CV_RGB_DIST( p_cur->c, p_next3[0].c );
1234             c4 = _CV_RGB_DIST( p_cur->c, p_next3[1].c );
1235
1236             p = p_next1;
1237
1238             if( c1 > c2 )
1239             {
1240                 p = p_next1 + 1;
1241                 c1 = c2;
1242             }
1243             if( c1 > c3 )
1244             {
1245                 p = p_next3;
1246                 c1 = c3;
1247             }
1248             if( c1 > c4 )
1249             {
1250                 p = p_next3 + 1;
1251                 c1 = c4;
1252             }
1253
1254             if( c1 < threshold )
1255             {
1256                 p_cur->p = p;
1257
1258                 if( layer == 0 )
1259                 {
1260                     p->a++;
1261                     p_cur = _CV_NEXT_BASE_C3(p_cur,1);
1262                     if( is_last_iter )
1263                         icvMaxRoi1( &(p->rect), j, i );
1264                 }
1265                 else
1266                 {
1267                     int a = p_cur->a;
1268
1269                     p->a += a;
1270                     p_cur->c.blue = p_cur->c.green = p_cur->c.red = 0;
1271                     p_cur++;
1272                     if( is_last_iter && a != 0 )
1273                         icvMaxRoi( &(p->rect), &(p_cur[-1].rect) );
1274                 }
1275             }
1276             else
1277             {
1278                 p_cur->p = stub;
1279                 if( is_last_iter /* && ( == 0 || p_cur->a != 0) */  )
1280                 {
1281                     cmp_node.data = p_cur;
1282                     CV_WRITE_SEQ_ELEM( cmp_node, writer );
1283                 }
1284
1285                 if( layer == 0 )
1286                 {
1287                     p_cur = _CV_NEXT_BASE_C3(p_cur,1);
1288                 }
1289                 else
1290                 {
1291                     p_cur->c.blue = p_cur->c.green = p_cur->c.red = 0;
1292                     p_cur++;
1293                 }
1294             }
1295
1296             /* find pointer for the second pixel */
1297             c1 = _CV_RGB_DIST( p_cur->c, p_next1[0].c );
1298             c2 = _CV_RGB_DIST( p_cur->c, p_next1[1].c );
1299             c3 = _CV_RGB_DIST( p_cur->c, p_next3[0].c );
1300             c4 = _CV_RGB_DIST( p_cur->c, p_next3[1].c );
1301
1302             p = p_next1;
1303             p_next1++;
1304
1305             if( c1 > c2 )
1306             {
1307                 p = p_next1;
1308                 c1 = c2;
1309             }
1310             if( c1 > c3 )
1311             {
1312                 p = p_next3;
1313                 c1 = c3;
1314             }
1315
1316             p_next3++;
1317             if( c1 > c4 )
1318             {
1319                 p = p_next3;
1320                 c1 = c4;
1321             }
1322
1323             if( c1 < threshold )
1324             {
1325                 p_cur->p = p;
1326
1327                 if( layer == 0 )
1328                 {
1329                     p->a++;
1330                     p_cur = _CV_NEXT_BASE_C3(p_cur,1);
1331                     if( is_last_iter )
1332                         icvMaxRoi1( &(p->rect), j + 1, i );
1333                 }
1334                 else
1335                 {
1336                     int a = p_cur->a;
1337
1338                     p->a += a;
1339                     p_cur->c.blue = p_cur->c.green = p_cur->c.red = 0;
1340                     p_cur++;
1341                     if( is_last_iter && a != 0 )
1342                         icvMaxRoi( &(p->rect), &(p_cur[-1].rect) );
1343                 }
1344             }
1345             else
1346             {
1347                 p_cur->p = stub;
1348                 if( is_last_iter /* && ( == 0 || p_cur->a != 0) */  )
1349                 {
1350                     cmp_node.data = p_cur;
1351                     CV_WRITE_SEQ_ELEM( cmp_node, writer );
1352                 }
1353                 if( layer == 0 )
1354                 {
1355                     p_cur = _CV_NEXT_BASE_C3(p_cur,1);
1356                 }
1357                 else
1358                 {
1359                     p_cur->c.blue = p_cur->c.green = p_cur->c.red = 0;
1360                     p_cur++;
1361                 }
1362             }
1363         }
1364
1365         /* clear c's */
1366         if( layer > 0 )
1367         {
1368             p_cur->c.blue = p_cur->c.green = p_cur->c.red = 0;
1369             p_cur++;
1370         }
1371
1372         if( !(i & 1) )
1373         {
1374             p_next1 -= size.width >> 1;
1375             p_next3 -= size.width >> 1;
1376         }
1377         else
1378         {
1379             p_next1++;
1380             p_next3++;
1381         }
1382     }
1383
1384     return CV_OK;
1385 }
1386
1387
1388
1389 /****************************************************************************************\
1390
1391     clusterization segmented components
1392
1393 \****************************************************************************************/
1394 static void
1395 icvExpandBaseLevelC1( _CvPyramid * base_p, _CvPyramid * p, _CvPyramidBase * start, int width )
1396 {
1397     int x = (int)((_CvPyramidBase *) base_p - start);
1398     int y = x / width;
1399
1400     x -= y * width;
1401     p->a = 1;
1402     p->rect.x1 = (ushort) x;
1403     p->rect.y1 = (ushort) y;
1404     p->rect.x2 = (ushort) (x + 1);
1405     p->rect.y2 = (ushort) (y + 1);
1406     p->c = base_p->c;
1407 }
1408
1409 CvStatus
1410 icvSegmentClusterC1( CvSeq * cmp_seq, CvSeq * res_seq,
1411                      double threshold, _CvPyramid * first_level_end, CvSize first_level_size )
1412 {
1413     const double eps = 1.;
1414     CvSeqWriter writer;
1415     CvSeqReader reader;
1416     _CvPyramid temp_cmp;
1417     _CvPyramidBase *first_level_start = (_CvPyramidBase *) first_level_end -
1418         first_level_size.width * first_level_size.height;
1419     int c, i, count = cmp_seq->total;
1420
1421     cvStartReadSeq( cmp_seq, &reader, 0 );
1422     cvStartAppendToSeq( res_seq, &writer );
1423
1424     if( threshold < eps )
1425     {
1426         /* if threshold is too small then simply copy all
1427            the components to the output sequence */
1428         for( i = 0; i < count; i++ )
1429         {
1430             CvConnectedComp comp;
1431             _CvPyramid *cmp = (_CvPyramid *) (((_CvListNode *) reader.ptr)->data);
1432             Cv32suf _c;
1433
1434             if( cmp < first_level_end )
1435             {
1436                 icvExpandBaseLevelC1( cmp, &temp_cmp, first_level_start,
1437                                       first_level_size.width );
1438                 cmp = &temp_cmp;
1439             }
1440
1441             _c.i = cvRound( cmp->c );
1442             cmp->c = _c.f;
1443             comp.value = cvRealScalar(_c.i);
1444             comp.area = cmp->a;
1445             comp.rect.x = cmp->rect.x1;
1446             comp.rect.y = cmp->rect.y1;
1447             comp.rect.width = cmp->rect.x2 - cmp->rect.x1;
1448             comp.rect.height = cmp->rect.y2 - cmp->rect.y1;
1449             comp.contour = 0;
1450
1451             CV_WRITE_SEQ_ELEM( comp, writer );
1452             CV_NEXT_SEQ_ELEM( sizeof( _CvListNode ), reader );
1453         }
1454     }
1455     else
1456     {
1457         _CvListNode stub_node;
1458         _CvListNode *prev = &stub_node;
1459
1460         stub_node.next = 0;
1461
1462         for( i = 0; i < count; i++ )
1463         {
1464             _CvListNode *node = (_CvListNode *) reader.ptr;
1465
1466             prev->next = node;
1467             prev = node;
1468             CV_NEXT_SEQ_ELEM( sizeof( _CvListNode ), reader );
1469         }
1470         prev->next = 0;
1471         prev = stub_node.next;
1472
1473         while( prev )
1474         {
1475             _CvListNode *node = prev->next;
1476             _CvListNode *acc = prev;
1477             _CvPyramid *cmp = (_CvPyramid *) (acc->data);
1478             CvConnectedComp comp;
1479             float c0 = cmp->c;
1480
1481             if( cmp < first_level_end )
1482             {
1483                 icvExpandBaseLevelC1( cmp, &temp_cmp, first_level_start,
1484                                       first_level_size.width );
1485             }
1486             else
1487             {
1488                 temp_cmp = *cmp;
1489                 temp_cmp.c *= temp_cmp.a;
1490             }
1491
1492             acc->next = 0;
1493             stub_node.next = 0;
1494             prev = &stub_node;
1495
1496             while( node )
1497             {
1498                 cmp = (_CvPyramid *) (node->data);
1499                 if( fabs( c0 - cmp->c ) < threshold )
1500                 {
1501                     _CvPyramid temp;
1502
1503                     /* exclude from global list and add to list of joint component */
1504                     prev->next = node->next;
1505                     node->next = acc;
1506                     acc = node;
1507
1508                     if( cmp < first_level_end )
1509                     {
1510                         icvExpandBaseLevelC1( cmp, &temp, first_level_start,
1511                                               first_level_size.width );
1512                         cmp = &temp;
1513                     }
1514
1515                     temp_cmp.a += cmp->a;
1516                     temp_cmp.c += cmp->c * cmp->a;
1517                     icvMaxRoi( &(temp_cmp.rect), &(cmp->rect) );
1518                 }
1519                 else
1520                 {
1521                     if( prev == &stub_node )
1522                     {
1523                         stub_node.next = node;
1524                     }
1525                     prev = node;
1526                 }
1527                 node = prev->next;
1528             }
1529
1530             if( temp_cmp.a != 0 )
1531             {
1532                 c = cvRound( temp_cmp.c / temp_cmp.a );
1533             }
1534             else
1535             {
1536                 c = cvRound( c0 );
1537             }
1538             node = acc;
1539
1540             while( node )
1541             {
1542                 Cv32suf _c;
1543                 cmp = (_CvPyramid *) (node->data);
1544                 _c.i = c; cmp->c = _c.f;
1545                 node = node->next;
1546             }
1547
1548             comp.value = cvRealScalar(c);
1549             comp.area = temp_cmp.a;
1550             comp.rect.x = temp_cmp.rect.x1;
1551             comp.rect.y = temp_cmp.rect.y1;
1552             comp.rect.width = temp_cmp.rect.x2 - temp_cmp.rect.x1;
1553             comp.rect.height = temp_cmp.rect.y2 - temp_cmp.rect.y1;
1554             comp.contour = 0;
1555
1556             CV_WRITE_SEQ_ELEM( comp, writer );
1557             prev = stub_node.next;
1558         }
1559     }
1560
1561     cvEndWriteSeq( &writer );
1562     return CV_OK;
1563 }
1564
1565 /****************************************************************************************\
1566
1567     clusterization segmented components
1568
1569 \****************************************************************************************/
1570 static void
1571 icvExpandBaseLevelC3( _CvPyramidC3 * base_p, _CvPyramidC3 * p,
1572                       _CvPyramidBaseC3 * start, int width )
1573 {
1574     int x = (int)((_CvPyramidBaseC3 *) base_p - start);
1575     int y = x / width;
1576
1577     x -= y * width;
1578     p->a = 1;
1579     p->rect.x1 = (ushort) x;
1580     p->rect.y1 = (ushort) y;
1581     p->rect.x2 = (ushort) (x + 1);
1582     p->rect.y2 = (ushort) (y + 1);
1583     p->c = base_p->c;
1584 }
1585
1586 CvStatus
1587 icvSegmentClusterC3( CvSeq * cmp_seq, CvSeq * res_seq,
1588                      double threshold,
1589                      _CvPyramidC3 * first_level_end, CvSize first_level_size )
1590 {
1591     const double eps = 1.;
1592     CvSeqWriter writer;
1593     CvSeqReader reader;
1594     _CvPyramidC3 temp_cmp;
1595     _CvPyramidBaseC3 *first_level_start = (_CvPyramidBaseC3 *) first_level_end -
1596         first_level_size.width * first_level_size.height;
1597     int i, count = cmp_seq->total;
1598     int c_blue, c_green, c_red;
1599
1600     cvStartReadSeq( cmp_seq, &reader, 0 );
1601     cvStartAppendToSeq( res_seq, &writer );
1602
1603     if( threshold < eps )
1604     {
1605         /* if threshold is too small then simply copy all
1606            the components to the output sequence */
1607         for( i = 0; i < count; i++ )
1608         {
1609             CvConnectedComp comp;
1610             _CvPyramidC3 *cmp = (_CvPyramidC3 *) (((_CvListNode *) reader.ptr)->data);
1611             Cv32suf _c;
1612
1613             if( cmp < first_level_end )
1614             {
1615                 icvExpandBaseLevelC3( cmp, &temp_cmp, first_level_start,
1616                                       first_level_size.width );
1617                 cmp = &temp_cmp;
1618             }
1619
1620             c_blue = cvRound( cmp->c.blue );
1621             c_green = cvRound( cmp->c.green );
1622             c_red = cvRound( cmp->c.red );
1623             _c.i = c_blue; cmp->c.blue = _c.f;
1624             _c.i = c_green; cmp->c.green = _c.f;
1625             _c.i = c_red; cmp->c.red = _c.f;
1626             comp.value = cvScalar( c_blue, c_green, c_red );
1627             comp.area = cmp->a;
1628             comp.rect.x = cmp->rect.x1;
1629             comp.rect.y = cmp->rect.y1;
1630             comp.rect.width = cmp->rect.x2 - cmp->rect.x1;
1631             comp.rect.height = cmp->rect.y2 - cmp->rect.y1;
1632             comp.contour = 0;
1633
1634             CV_WRITE_SEQ_ELEM( comp, writer );
1635             CV_NEXT_SEQ_ELEM( sizeof( _CvListNode ), reader );
1636         }
1637     }
1638     else
1639     {
1640         _CvListNode stub_node;
1641         _CvListNode *prev = &stub_node;
1642
1643         stub_node.next = 0;
1644
1645         for( i = 0; i < count; i++ )
1646         {
1647             _CvListNode *node = (_CvListNode *) reader.ptr;
1648
1649             prev->next = node;
1650             prev = node;
1651             CV_NEXT_SEQ_ELEM( sizeof( _CvListNode ), reader );
1652         }
1653         prev->next = 0;
1654         prev = stub_node.next;
1655
1656         while( prev )
1657         {
1658             _CvListNode *node = prev->next;
1659             _CvListNode *acc = prev;
1660             _CvPyramidC3 *cmp = (_CvPyramidC3 *) (acc->data);
1661             CvConnectedComp comp;
1662             _CvRGBf c0 = cmp->c;
1663
1664             if( cmp < first_level_end )
1665             {
1666                 icvExpandBaseLevelC3( cmp, &temp_cmp, first_level_start,
1667                                       first_level_size.width );
1668             }
1669             else
1670             {
1671                 temp_cmp = *cmp;
1672                 temp_cmp.c.blue *= temp_cmp.a;
1673                 temp_cmp.c.green *= temp_cmp.a;
1674                 temp_cmp.c.red *= temp_cmp.a;
1675             }
1676
1677             acc->next = 0;
1678             stub_node.next = 0;
1679             prev = &stub_node;
1680
1681             while( node )
1682             {
1683                 cmp = (_CvPyramidC3 *) (node->data);
1684                 if( _CV_RGB_DIST( c0, cmp->c ) < threshold )
1685                 {
1686                     _CvPyramidC3 temp;
1687
1688                     /* exclude from global list and add to list of joint component */
1689                     prev->next = node->next;
1690                     node->next = acc;
1691                     acc = node;
1692
1693                     if( cmp < first_level_end )
1694                     {
1695                         icvExpandBaseLevelC3( cmp, &temp, first_level_start,
1696                                               first_level_size.width );
1697                         cmp = &temp;
1698                     }
1699
1700                     temp_cmp.a += cmp->a;
1701                     temp_cmp.c.blue += cmp->c.blue * cmp->a;
1702                     temp_cmp.c.green += cmp->c.green * cmp->a;
1703                     temp_cmp.c.red += cmp->c.red * cmp->a;
1704                     icvMaxRoi( &(temp_cmp.rect), &(cmp->rect) );
1705                 }
1706                 else
1707                 {
1708                     if( prev == &stub_node )
1709                     {
1710                         stub_node.next = node;
1711                     }
1712                     prev = node;
1713                 }
1714                 node = prev->next;
1715             }
1716
1717             if( temp_cmp.a != 0 )
1718             {
1719                 c_blue = cvRound( temp_cmp.c.blue / temp_cmp.a );
1720                 c_green = cvRound( temp_cmp.c.green / temp_cmp.a );
1721                 c_red = cvRound( temp_cmp.c.red / temp_cmp.a );
1722             }
1723             else
1724             {
1725                 c_blue = cvRound( c0.blue );
1726                 c_green = cvRound( c0.green );
1727                 c_red = cvRound( c0.red );
1728             }
1729             node = acc;
1730
1731             while( node )
1732             {
1733                 Cv32suf _c;
1734                 cmp = (_CvPyramidC3 *) (node->data);
1735                 _c.i = c_blue; cmp->c.blue = _c.f;
1736                 _c.i = c_green; cmp->c.green = _c.f;
1737                 _c.i = c_red; cmp->c.red = _c.f;
1738                 node = node->next;
1739             }
1740
1741             comp.value = cvScalar( c_blue, c_green, c_red );
1742             comp.area = temp_cmp.a;
1743             comp.rect.x = temp_cmp.rect.x1;
1744             comp.rect.y = temp_cmp.rect.y1;
1745             comp.rect.width = temp_cmp.rect.x2 - temp_cmp.rect.x1;
1746             comp.rect.height = temp_cmp.rect.y2 - temp_cmp.rect.y1;
1747             comp.contour = 0;
1748
1749             CV_WRITE_SEQ_ELEM( comp, writer );
1750             prev = stub_node.next;
1751         }
1752     }
1753
1754     cvEndWriteSeq( &writer );
1755     return CV_OK;
1756 }
1757
1758 /****************************************************************************************\
1759
1760                  definition of the maximum roi size
1761
1762 \****************************************************************************************/
1763 void
1764 icvMaxRoi( _CvRect16u * max_rect, _CvRect16u * cur_rect )
1765 {
1766     if( max_rect->x2 == 0 )
1767         *max_rect = *cur_rect;
1768     else
1769     {
1770         if( max_rect->x1 > cur_rect->x1 )
1771             max_rect->x1 = cur_rect->x1;
1772         if( max_rect->y1 > cur_rect->y1 )
1773             max_rect->y1 = cur_rect->y1;
1774
1775         if( max_rect->x2 < cur_rect->x2 )
1776             max_rect->x2 = cur_rect->x2;
1777         if( max_rect->y2 < cur_rect->y2 )
1778             max_rect->y2 = cur_rect->y2;
1779     }
1780 }
1781
1782 void
1783 icvMaxRoi1( _CvRect16u * max_rect, int x, int y )
1784 {
1785     if( max_rect->x2 == 0 )
1786     {
1787         max_rect->x1 = (ushort) x;
1788         max_rect->y1 = (ushort) y;
1789
1790         ++x;
1791         ++y;
1792
1793         max_rect->x2 = (ushort) x;
1794         max_rect->y2 = (ushort) y;
1795     }
1796     else
1797     {
1798         if( max_rect->x1 > x )
1799             max_rect->x1 = (ushort) x;
1800         if( max_rect->y1 > y )
1801             max_rect->y1 = (ushort) y;
1802
1803         ++x;
1804         ++y;
1805
1806         if( max_rect->x2 < x )
1807             max_rect->x2 = (ushort) x;
1808         if( max_rect->y2 < y )
1809             max_rect->y2 = (ushort) y;
1810     }
1811 }
1812
1813
1814 /*F///////////////////////////////////////////////////////////////////////////////////////
1815 //    Name:    cvPyrSegmentation
1816 //    Purpose:
1817 //      segments an image using pyramid-linking technique
1818 //    Context:
1819 //    Parameters:
1820 //      src - source image
1821 //      dst - destination image
1822 //      comp - pointer to returned connected component sequence
1823 //      storage - where the sequence is stored
1824 //      level - maximal pyramid level
1825 //      threshold1 - first threshold, affecting on detalization level when pyramid
1826 //                   is built.
1827 //      threshold2 - second threshold - affects on final components merging.
1828 //    Returns:
1829 //    Notes:
1830 //      Source and destination image must be equal types and channels
1831 //F*/
1832 CV_IMPL void
1833 cvPyrSegmentation( IplImage * src,
1834                    IplImage * dst,
1835                    CvMemStorage * storage,
1836                    CvSeq ** comp, int level, double threshold1, double threshold2 )
1837 {
1838     CvSize src_size, dst_size;
1839     uchar *src_data = 0;
1840     uchar *dst_data = 0;
1841     int src_step = 0, dst_step = 0;
1842     int thresh1 = cvRound( threshold1 );
1843     int thresh2 = cvRound( threshold2 );
1844
1845     CV_FUNCNAME( "cvPyrSegmentation" );
1846
1847     __BEGIN__;
1848
1849     if( src->depth != IPL_DEPTH_8U )
1850         CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1851
1852     if( src->depth != dst->depth || src->nChannels != dst->nChannels )
1853         CV_ERROR( CV_StsBadArg, "src and dst have different formats" );
1854
1855     cvGetRawData( src, &src_data, &src_step, &src_size );
1856     cvGetRawData( dst, &dst_data, &dst_step, &dst_size );
1857
1858     if( src_size.width != dst_size.width ||
1859         src_size.height != dst_size.height )
1860         CV_ERROR( CV_StsBadArg, "src and dst have different ROIs" );
1861
1862     switch (src->nChannels)
1863     {
1864     case 1:
1865         IPPI_CALL( icvPyrSegmentation8uC1R( src_data, src_step,
1866                                             dst_data, dst_step,
1867                                             src_size,
1868                                             CV_GAUSSIAN_5x5,
1869                                             comp, storage, level, thresh1, thresh2 ));
1870         break;
1871     case 3:
1872         IPPI_CALL( icvPyrSegmentation8uC3R( src_data, src_step,
1873                                             dst_data, dst_step,
1874                                             src_size,
1875                                             CV_GAUSSIAN_5x5,
1876                                             comp, storage, level, thresh1, thresh2 ));
1877         break;
1878     default:
1879         CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1880     }
1881     __END__;
1882 }
1883
1884
1885 /* End of file. */