vcs.maemo.org Git - opencv/blob - src/cv/cvsegmentation.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
  42 #include "_cv.h"
  43
  44 /****************************************************************************************\
  45 *                                       Watershed                                        *
  46 \****************************************************************************************/
  47
  48 typedef struct CvWSNode
  49 {
  50     struct CvWSNode* next;
  51     int mask_ofs;
  52     int img_ofs;
  53 }
  54 CvWSNode;
  55
  56 typedef struct CvWSQueue
  57 {
  58     CvWSNode* first;
  59     CvWSNode* last;
  60 }
  61 CvWSQueue;
  62
  63 static CvWSNode*
  64 icvAllocWSNodes( CvMemStorage* storage )
  65 {
  66     CvWSNode* n = 0;
  67
  68     CV_FUNCNAME( "icvAllocWSNodes" );
  69
  70     __BEGIN__;
  71
  72     int i, count = (storage->block_size - sizeof(CvMemBlock))/sizeof(*n) - 1;
  73
  74     CV_CALL( n = (CvWSNode*)cvMemStorageAlloc( storage, count*sizeof(*n) ));
  75     for( i = 0; i < count-1; i++ )
  76         n[i].next = n + i + 1;
  77     n[count-1].next = 0;
  78
  79     __END__;
  80
  81     return n;
  82 }
  83
  84
  85 CV_IMPL void
  86 cvWatershed( const CvArr* srcarr, CvArr* dstarr )
  87 {
  88     const int IN_QUEUE = -2;
  89     const int WSHED = -1;
  90     const int NQ = 256;
  91     CvMemStorage* storage = 0;
  92
  93     CV_FUNCNAME( "cvWatershed" );
  94
  95     __BEGIN__;
  96
  97     CvMat sstub, *src;
  98     CvMat dstub, *dst;
  99     CvSize size;
 100     CvWSNode* free_node = 0, *node;
 101     CvWSQueue q[NQ];
 102     int active_queue;
 103     int i, j;
 104     int db, dg, dr;
 105     int* mask;
 106     uchar* img;
 107     int mstep, istep;
 108     int subs_tab[513];
 109
 110     // MAX(a,b) = b + MAX(a-b,0)
 111     #define ws_max(a,b) ((b) + subs_tab[(a)-(b)+NQ])
 112     // MIN(a,b) = a - MAX(a-b,0)
 113     #define ws_min(a,b) ((a) - subs_tab[(a)-(b)+NQ])
 114
 115     #define ws_push(idx,mofs,iofs)  \
 116     {                               \
 117         if( !free_node )            \
 118             CV_CALL( free_node = icvAllocWSNodes( storage ));\
 119         node = free_node;           \
 120         free_node = free_node->next;\
 121         node->next = 0;             \
 122         node->mask_ofs = mofs;      \
 123         node->img_ofs = iofs;       \
 124         if( q[idx].last )           \
 125             q[idx].last->next=node; \
 126         else                        \
 127             q[idx].first = node;    \
 128         q[idx].last = node;         \
 129     }
 130
 131     #define ws_pop(idx,mofs,iofs)   \
 132     {                               \
 133         node = q[idx].first;        \
 134         q[idx].first = node->next;  \
 135         if( !node->next )           \
 136             q[idx].last = 0;        \
 137         node->next = free_node;     \
 138         free_node = node;           \
 139         mofs = node->mask_ofs;      \
 140         iofs = node->img_ofs;       \
 141     }
 142
 143     #define c_diff(ptr1,ptr2,diff)      \
 144     {                                   \
 145         db = abs((ptr1)[0] - (ptr2)[0]);\
 146         dg = abs((ptr1)[1] - (ptr2)[1]);\
 147         dr = abs((ptr1)[2] - (ptr2)[2]);\
 148         diff = ws_max(db,dg);           \
 149         diff = ws_max(diff,dr);         \
 150         assert( 0 <= diff && diff <= 255 ); \
 151     }
 152
 153     CV_CALL( src = cvGetMat( srcarr, &sstub ));
 154     CV_CALL( dst = cvGetMat( dstarr, &dstub ));
 155
 156     if( CV_MAT_TYPE(src->type) != CV_8UC3 )
 157         CV_ERROR( CV_StsUnsupportedFormat, "Only 8-bit, 3-channel input images are supported" );
 158
 159     if( CV_MAT_TYPE(dst->type) != CV_32SC1 )
 160         CV_ERROR( CV_StsUnsupportedFormat,
 161             "Only 32-bit, 1-channel output images are supported" );
 162
 163     if( !CV_ARE_SIZES_EQ( src, dst ))
 164         CV_ERROR( CV_StsUnmatchedSizes, "The input and output images must have the same size" );
 165
 166     size = cvGetMatSize(src);
 167
 168     CV_CALL( storage = cvCreateMemStorage() );
 169
 170     istep = src->step;
 171     img = src->data.ptr;
 172     mstep = dst->step / sizeof(mask[0]);
 173     mask = dst->data.i;
 174
 175     memset( q, 0, NQ*sizeof(q[0]) );
 176
 177     for( i = 0; i < 256; i++ )
 178         subs_tab[i] = 0;
 179     for( i = 256; i <= 512; i++ )
 180         subs_tab[i] = i - 256;
 181
 182     // draw a pixel-wide border of dummy "watershed" (i.e. boundary) pixels
 183     for( j = 0; j < size.width; j++ )
 184         mask[j] = mask[j + mstep*(size.height-1)] = WSHED;
 185
 186     // initial phase: put all the neighbor pixels of each marker to the ordered queue -
 187     // determine the initial boundaries of the basins
 188     for( i = 1; i < size.height-1; i++ )
 189     {
 190         img += istep; mask += mstep;
 191         mask[0] = mask[size.width-1] = WSHED;
 192
 193         for( j = 1; j < size.width-1; j++ )
 194         {
 195             int* m = mask + j;
 196             if( m[0] < 0 ) m[0] = 0;
 197             if( m[0] == 0 && (m[-1] > 0 || m[1] > 0 || m[-mstep] > 0 || m[mstep] > 0) )
 198             {
 199                 uchar* ptr = img + j*3;
 200                 int idx = 256, t;
 201                 if( m[-1] > 0 )
 202                     c_diff( ptr, ptr - 3, idx );
 203                 if( m[1] > 0 )
 204                 {
 205                     c_diff( ptr, ptr + 3, t );
 206                     idx = ws_min( idx, t );
 207                 }
 208                 if( m[-mstep] > 0 )
 209                 {
 210                     c_diff( ptr, ptr - istep, t );
 211                     idx = ws_min( idx, t );
 212                 }
 213                 if( m[mstep] > 0 )
 214                 {
 215                     c_diff( ptr, ptr + istep, t );
 216                     idx = ws_min( idx, t );
 217                 }
 218                 assert( 0 <= idx && idx <= 255 );
 219                 ws_push( idx, i*mstep + j, i*istep + j*3 );
 220                 m[0] = IN_QUEUE;
 221             }
 222         }
 223     }
 224
 225     // find the first non-empty queue
 226     for( i = 0; i < NQ; i++ )
 227         if( q[i].first )
 228             break;
 229
 230     // if there is no markers, exit immediately
 231     if( i == NQ )
 232         EXIT;
 233
 234     active_queue = i;
 235     img = src->data.ptr;
 236     mask = dst->data.i;
 237
 238     // recursively fill the basins
 239     for(;;)
 240     {
 241         int mofs, iofs;
 242         int lab = 0, t;
 243         int* m;
 244         uchar* ptr;
 245
 246         if( q[active_queue].first == 0 )
 247         {
 248             for( i = active_queue+1; i < NQ; i++ )
 249                 if( q[i].first )
 250                     break;
 251             if( i == NQ )
 252                 break;
 253             active_queue = i;
 254         }
 255
 256         ws_pop( active_queue, mofs, iofs );
 257
 258         m = mask + mofs;
 259         ptr = img + iofs;
 260         t = m[-1];
 261         if( t > 0 ) lab = t;
 262         t = m[1];
 263         if( t > 0 )
 264         {
 265             if( lab == 0 ) lab = t;
 266             else if( t != lab ) lab = WSHED;
 267         }
 268         t = m[-mstep];
 269         if( t > 0 )
 270         {
 271             if( lab == 0 ) lab = t;
 272             else if( t != lab ) lab = WSHED;
 273         }
 274         t = m[mstep];
 275         if( t > 0 )
 276         {
 277             if( lab == 0 ) lab = t;
 278             else if( t != lab ) lab = WSHED;
 279         }
 280         assert( lab != 0 );
 281         m[0] = lab;
 282         if( lab == WSHED )
 283             continue;
 284
 285         if( m[-1] == 0 )
 286         {
 287             c_diff( ptr, ptr - 3, t );
 288             ws_push( t, mofs - 1, iofs - 3 );
 289             active_queue = ws_min( active_queue, t );
 290             m[-1] = IN_QUEUE;
 291         }
 292         if( m[1] == 0 )
 293         {
 294             c_diff( ptr, ptr + 3, t );
 295             ws_push( t, mofs + 1, iofs + 3 );
 296             active_queue = ws_min( active_queue, t );
 297             m[1] = IN_QUEUE;
 298         }
 299         if( m[-mstep] == 0 )
 300         {
 301             c_diff( ptr, ptr - istep, t );
 302             ws_push( t, mofs - mstep, iofs - istep );
 303             active_queue = ws_min( active_queue, t );
 304             m[-mstep] = IN_QUEUE;
 305         }
 306         if( m[mstep] == 0 )
 307         {
 308             c_diff( ptr, ptr + 3, t );
 309             ws_push( t, mofs + mstep, iofs + istep );
 310             active_queue = ws_min( active_queue, t );
 311             m[mstep] = IN_QUEUE;
 312         }
 313     }
 314
 315     __END__;
 316
 317     cvReleaseMemStorage( &storage );
 318 }
 319
 320
 321 void cv::watershed( const Mat& src, Mat& markers )
 322 {
 323     CvMat _src = src, _markers = markers;
 324     cvWatershed( &_src, &_markers );
 325 }
 326
 327
 328 /****************************************************************************************\
 329 *                                         Meanshift                                      *
 330 \****************************************************************************************/
 331
 332 CV_IMPL void
 333 cvPyrMeanShiftFiltering( const CvArr* srcarr, CvArr* dstarr,
 334                          double sp0, double sr, int max_level,
 335                          CvTermCriteria termcrit )
 336 {
 337     const int cn = 3;
 338     const int MAX_LEVELS = 8;
 339     CvMat* src_pyramid[MAX_LEVELS+1];
 340     CvMat* dst_pyramid[MAX_LEVELS+1];
 341     CvMat* mask0 = 0;
 342     int i, j, level;
 343     //uchar* submask = 0;
 344
 345     #define cdiff(ofs0) (tab[c0-dptr[ofs0]+255] + \
 346         tab[c1-dptr[(ofs0)+1]+255] + tab[c2-dptr[(ofs0)+2]+255] >= isr22)
 347
 348     memset( src_pyramid, 0, sizeof(src_pyramid) );
 349     memset( dst_pyramid, 0, sizeof(dst_pyramid) );
 350
 351     CV_FUNCNAME( "cvPyrMeanShiftFiltering" );
 352
 353     __BEGIN__;
 354
 355     double sr2 = sr * sr;
 356     int isr2 = cvRound(sr2), isr22 = MAX(isr2,16);
 357     int tab[768];
 358     CvMat sstub0, *src0;
 359     CvMat dstub0, *dst0;
 360
 361     CV_CALL( src0 = cvGetMat( srcarr, &sstub0 ));
 362     CV_CALL( dst0 = cvGetMat( dstarr, &dstub0 ));
 363
 364     if( CV_MAT_TYPE(src0->type) != CV_8UC3 )
 365         CV_ERROR( CV_StsUnsupportedFormat, "Only 8-bit, 3-channel images are supported" );
 366
 367     if( !CV_ARE_TYPES_EQ( src0, dst0 ))
 368         CV_ERROR( CV_StsUnmatchedFormats, "The input and output images must have the same type" );
 369
 370     if( !CV_ARE_SIZES_EQ( src0, dst0 ))
 371         CV_ERROR( CV_StsUnmatchedSizes, "The input and output images must have the same size" );
 372
 373     if( (unsigned)max_level > (unsigned)MAX_LEVELS )
 374         CV_ERROR( CV_StsOutOfRange, "The number of pyramid levels is too large or negative" );
 375
 376     if( !(termcrit.type & CV_TERMCRIT_ITER) )
 377         termcrit.max_iter = 5;
 378     termcrit.max_iter = MAX(termcrit.max_iter,1);
 379     termcrit.max_iter = MIN(termcrit.max_iter,100);
 380     if( !(termcrit.type & CV_TERMCRIT_EPS) )
 381         termcrit.epsilon = 1.f;
 382     termcrit.epsilon = MAX(termcrit.epsilon, 0.f);
 383
 384     for( i = 0; i < 768; i++ )
 385         tab[i] = (i - 255)*(i - 255);
 386
 387     // 1. construct pyramid
 388     src_pyramid[0] = src0;
 389     dst_pyramid[0] = dst0;
 390     for( level = 1; level <= max_level; level++ )
 391     {
 392         CV_CALL( src_pyramid[level] = cvCreateMat( (src_pyramid[level-1]->rows+1)/2,
 393                         (src_pyramid[level-1]->cols+1)/2, src_pyramid[level-1]->type ));
 394         CV_CALL( dst_pyramid[level] = cvCreateMat( src_pyramid[level]->rows,
 395                         src_pyramid[level]->cols, src_pyramid[level]->type ));
 396         CV_CALL( cvPyrDown( src_pyramid[level-1], src_pyramid[level] ));
 397         //CV_CALL( cvResize( src_pyramid[level-1], src_pyramid[level], CV_INTER_AREA ));
 398     }
 399
 400     CV_CALL( mask0 = cvCreateMat( src0->rows, src0->cols, CV_8UC1 ));
 401     //CV_CALL( submask = (uchar*)cvAlloc( (sp+2)*(sp+2) ));
 402
 403     // 2. apply meanshift, starting from the pyramid top (i.e. the smallest layer)
 404     for( level = max_level; level >= 0; level-- )
 405     {
 406         CvMat* src = src_pyramid[level];
 407         CvSize size = cvGetMatSize(src);
 408         uchar* sptr = src->data.ptr;
 409         int sstep = src->step;
 410         uchar* mask = 0;
 411         int mstep = 0;
 412         uchar* dptr;
 413         int dstep;
 414         float sp = (float)(sp0 / (1 << level));
 415         sp = MAX( sp, 1 );
 416
 417         if( level < max_level )
 418         {
 419             CvSize size1 = cvGetMatSize(dst_pyramid[level+1]);
 420             CvMat m = cvMat( size.height, size.width, CV_8UC1, mask0->data.ptr );
 421             dstep = dst_pyramid[level+1]->step;
 422             dptr = dst_pyramid[level+1]->data.ptr + dstep + cn;
 423             mstep = m.step;
 424             mask = m.data.ptr + mstep;
 425             //cvResize( dst_pyramid[level+1], dst_pyramid[level], CV_INTER_CUBIC );
 426             cvPyrUp( dst_pyramid[level+1], dst_pyramid[level] );
 427             cvZero( &m );
 428
 429             for( i = 1; i < size1.height-1; i++, dptr += dstep - (size1.width-2)*3, mask += mstep*2 )
 430             {
 431                 for( j = 1; j < size1.width-1; j++, dptr += cn )
 432                 {
 433                     int c0 = dptr[0], c1 = dptr[1], c2 = dptr[2];
 434                     mask[j*2 - 1] = cdiff(-3) || cdiff(3) || cdiff(-dstep-3) || cdiff(-dstep) ||
 435                         cdiff(-dstep+3) || cdiff(dstep-3) || cdiff(dstep) || cdiff(dstep+3);
 436                 }
 437             }
 438
 439             cvDilate( &m, &m, 0, 1 );
 440             mask = m.data.ptr;
 441         }
 442
 443         dptr = dst_pyramid[level]->data.ptr;
 444         dstep = dst_pyramid[level]->step;
 445
 446         for( i = 0; i < size.height; i++, sptr += sstep - size.width*3,
 447                                           dptr += dstep - size.width*3,
 448                                           mask += mstep )
 449         {
 450             for( j = 0; j < size.width; j++, sptr += 3, dptr += 3 )
 451             {
 452                 int x0 = j, y0 = i, x1, y1, iter;
 453                 int c0, c1, c2;
 454
 455                 if( mask && !mask[j] )
 456                     continue;
 457
 458                 c0 = sptr[0], c1 = sptr[1], c2 = sptr[2];
 459
 460                 // iterate meanshift procedure
 461                 for( iter = 0; iter < termcrit.max_iter; iter++ )
 462                 {
 463                     uchar* ptr;
 464                     int x, y, count = 0;
 465                     int minx, miny, maxx, maxy;
 466                     int s0 = 0, s1 = 0, s2 = 0, sx = 0, sy = 0;
 467                     double icount;
 468                     int stop_flag;
 469
 470                     //mean shift: process pixels in window (p-sigmaSp)x(p+sigmaSp)
 471                     minx = cvRound(x0 - sp); minx = MAX(minx, 0);
 472                     miny = cvRound(y0 - sp); miny = MAX(miny, 0);
 473                     maxx = cvRound(x0 + sp); maxx = MIN(maxx, size.width-1);
 474                     maxy = cvRound(y0 + sp); maxy = MIN(maxy, size.height-1);
 475                     ptr = sptr + (miny - i)*sstep + (minx - j)*3;
 476
 477                     for( y = miny; y <= maxy; y++, ptr += sstep - (maxx-minx+1)*3 )
 478                     {
 479                         int row_count = 0;
 480                         x = minx;
 481                         for( ; x + 3 <= maxx; x += 4, ptr += 12 )
 482                         {
 483                             int t0 = ptr[0], t1 = ptr[1], t2 = ptr[2];
 484                             if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )
 485                             {
 486                                 s0 += t0; s1 += t1; s2 += t2;
 487                                 sx += x; row_count++;
 488                             }
 489                             t0 = ptr[3], t1 = ptr[4], t2 = ptr[5];
 490                             if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )
 491                             {
 492                                 s0 += t0; s1 += t1; s2 += t2;
 493                                 sx += x+1; row_count++;
 494                             }
 495                             t0 = ptr[6], t1 = ptr[7], t2 = ptr[8];
 496                             if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )
 497                             {
 498                                 s0 += t0; s1 += t1; s2 += t2;
 499                                 sx += x+2; row_count++;
 500                             }
 501                             t0 = ptr[9], t1 = ptr[10], t2 = ptr[11];
 502                             if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )
 503                             {
 504                                 s0 += t0; s1 += t1; s2 += t2;
 505                                 sx += x+3; row_count++;
 506                             }
 507                         }
 508
 509                         for( ; x <= maxx; x++, ptr += 3 )
 510                         {
 511                             int t0 = ptr[0], t1 = ptr[1], t2 = ptr[2];
 512                             if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )
 513                             {
 514                                 s0 += t0; s1 += t1; s2 += t2;
 515                                 sx += x; row_count++;
 516                             }
 517                         }
 518                         count += row_count;
 519                         sy += y*row_count;
 520                     }
 521
 522                     if( count == 0 )
 523                         break;
 524
 525                     icount = 1./count;
 526                     x1 = cvRound(sx*icount);
 527                     y1 = cvRound(sy*icount);
 528                     s0 = cvRound(s0*icount);
 529                     s1 = cvRound(s1*icount);
 530                     s2 = cvRound(s2*icount);
 531
 532                     stop_flag = (x0 == x1 && y0 == y1) || abs(x1-x0) + abs(y1-y0) +
 533                         tab[s0 - c0 + 255] + tab[s1 - c1 + 255] +
 534                         tab[s2 - c2 + 255] <= termcrit.epsilon;
 535
 536                     x0 = x1; y0 = y1;
 537                     c0 = s0; c1 = s1; c2 = s2;
 538
 539                     if( stop_flag )
 540                         break;
 541                 }
 542
 543                 dptr[0] = (uchar)c0;
 544                 dptr[1] = (uchar)c1;
 545                 dptr[2] = (uchar)c2;
 546             }
 547         }
 548     }
 549
 550     __END__;
 551
 552     for( i = 1; i <= MAX_LEVELS; i++ )
 553     {
 554         cvReleaseMat( &src_pyramid[i] );
 555         cvReleaseMat( &dst_pyramid[i] );
 556     }
 557     cvReleaseMat( &mask0 );
 558 }
 559