segmentation.cpp

00001 /*M///////////////////////////////////////////////////////////////////////////////////////
00002 //
00003 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
00004 //
00005 //  By downloading, copying, installing or using the software you agree to this license.
00006 //  If you do not agree to this license, do not download, install,
00007 //  copy or use the software.
00008 //
00009 //
00010 //                           License Agreement
00011 //                For Open Source Computer Vision Library
00012 //
00013 // Copyright (C) 2000, Intel Corporation, all rights reserved.
00014 // Copyright (C) 2013, OpenCV Foundation, all rights reserved.
00015 // Third party copyrights are property of their respective owners.
00016 //
00017 // Redistribution and use in source and binary forms, with or without modification,
00018 // are permitted provided that the following conditions are met:
00019 //
00020 //   * Redistribution's of source code must retain the above copyright notice,
00021 //     this list of conditions and the following disclaimer.
00022 //
00023 //   * Redistribution's in binary form must reproduce the above copyright notice,
00024 //     this list of conditions and the following disclaimer in the documentation
00025 //     and/or other materials provided with the distribution.
00026 //
00027 //   * The name of the copyright holders may not be used to endorse or promote products
00028 //     derived from this software without specific prior written permission.
00029 //
00030 // This software is provided by the copyright holders and contributors "as is" and
00031 // any express or implied warranties, including, but not limited to, the implied
00032 // warranties of merchantability and fitness for a particular purpose are disclaimed.
00033 // In no event shall the Intel Corporation or contributors be liable for any direct,
00034 // indirect, incidental, special, exemplary, or consequential damages
00035 // (including, but not limited to, procurement of substitute goods or services;
00036 // loss of use, data, or profits; or business interruption) however caused
00037 // and on any theory of liability, whether in contract, strict liability,
00038 // or tort (including negligence or otherwise) arising in any way out of
00039 // the use of this software, even if advised of the possibility of such damage.
00040 //
00041 //M*/
00042 
00043 #include "precomp.hpp"
00044 
00045 /****************************************************************************************\
00046 *                                       Watershed                                        *
00047 \****************************************************************************************/
00048 
00049 namespace cv
00050 {
00051 // A node represents a pixel to label
00052 struct WSNode
00053 {
00054     int next;
00055     int mask_ofs;
00056     int img_ofs;
00057 };
00058 
00059 // Queue for WSNodes
00060 struct WSQueue
00061 {
00062     WSQueue() { first = last = 0; }
00063     int first, last;
00064 };
00065 
00066 
00067 static int
00068 allocWSNodes( std::vector<WSNode>& storage )
00069 {
00070     int sz = (int)storage.size();
00071     int newsz = MAX(128, sz*3/2);
00072 
00073     storage.resize(newsz);
00074     if( sz == 0 )
00075     {
00076         storage[0].next = 0;
00077         sz = 1;
00078     }
00079     for( int i = sz; i < newsz-1; i++ )
00080         storage[i].next = i+1;
00081     storage[newsz-1].next = 0;
00082     return sz;
00083 }
00084 
00085 }
00086 
00087 
00088 void cv::watershed( InputArray _src, InputOutputArray _markers )
00089 {
00090     // Labels for pixels
00091     const int IN_QUEUE = -2; // Pixel visited
00092     const int WSHED = -1; // Pixel belongs to watershed
00093 
00094     // possible bit values = 2^8
00095     const int NQ = 256;
00096 
00097     Mat src = _src.getMat(), dst = _markers.getMat();
00098     Size size = src.size();
00099 
00100     // Vector of every created node
00101     std::vector<WSNode> storage;
00102     int free_node = 0, node;
00103     // Priority queue of queues of nodes
00104     // from high priority (0) to low priority (255)
00105     WSQueue q[NQ];
00106     // Non-empty queue with highest priority
00107     int active_queue;
00108     int i, j;
00109     // Color differences
00110     int db, dg, dr;
00111     int subs_tab[513];
00112 
00113     // MAX(a,b) = b + MAX(a-b,0)
00114     #define ws_max(a,b) ((b) + subs_tab[(a)-(b)+NQ])
00115     // MIN(a,b) = a - MAX(a-b,0)
00116     #define ws_min(a,b) ((a) - subs_tab[(a)-(b)+NQ])
00117 
00118     // Create a new node with offsets mofs and iofs in queue idx
00119     #define ws_push(idx,mofs,iofs)          \
00120     {                                       \
00121         if( !free_node )                    \
00122             free_node = allocWSNodes( storage );\
00123         node = free_node;                   \
00124         free_node = storage[free_node].next;\
00125         storage[node].next = 0;             \
00126         storage[node].mask_ofs = mofs;      \
00127         storage[node].img_ofs = iofs;       \
00128         if( q[idx].last )                   \
00129             storage[q[idx].last].next=node; \
00130         else                                \
00131             q[idx].first = node;            \
00132         q[idx].last = node;                 \
00133     }
00134 
00135     // Get next node from queue idx
00136     #define ws_pop(idx,mofs,iofs)           \
00137     {                                       \
00138         node = q[idx].first;                \
00139         q[idx].first = storage[node].next;  \
00140         if( !storage[node].next )           \
00141             q[idx].last = 0;                \
00142         storage[node].next = free_node;     \
00143         free_node = node;                   \
00144         mofs = storage[node].mask_ofs;      \
00145         iofs = storage[node].img_ofs;       \
00146     }
00147 
00148     // Get highest absolute channel difference in diff
00149     #define c_diff(ptr1,ptr2,diff)           \
00150     {                                        \
00151         db = std::abs((ptr1)[0] - (ptr2)[0]);\
00152         dg = std::abs((ptr1)[1] - (ptr2)[1]);\
00153         dr = std::abs((ptr1)[2] - (ptr2)[2]);\
00154         diff = ws_max(db,dg);                \
00155         diff = ws_max(diff,dr);              \
00156         assert( 0 <= diff && diff <= 255 );  \
00157     }
00158 
00159     CV_Assert( src.type() == CV_8UC3 && dst.type() == CV_32SC1 );
00160     CV_Assert( src.size() == dst.size() );
00161 
00162     // Current pixel in input image
00163     const uchar* img = src.ptr();
00164     // Step size to next row in input image
00165     int istep = int(src.step/sizeof(img[0]));
00166 
00167     // Current pixel in mask image
00168     int* mask = dst.ptr<int>();
00169     // Step size to next row in mask image
00170     int mstep = int(dst.step / sizeof(mask[0]));
00171 
00172     for( i = 0; i < 256; i++ )
00173         subs_tab[i] = 0;
00174     for( i = 256; i <= 512; i++ )
00175         subs_tab[i] = i - 256;
00176 
00177     // draw a pixel-wide border of dummy "watershed" (i.e. boundary) pixels
00178     for( j = 0; j < size.width; j++ )
00179         mask[j] = mask[j + mstep*(size.height-1)] = WSHED;
00180 
00181     // initial phase: put all the neighbor pixels of each marker to the ordered queue -
00182     // determine the initial boundaries of the basins
00183     for( i = 1; i < size.height-1; i++ )
00184     {
00185         img += istep; mask += mstep;
00186         mask[0] = mask[size.width-1] = WSHED; // boundary pixels
00187 
00188         for( j = 1; j < size.width-1; j++ )
00189         {
00190             int* m = mask + j;
00191             if( m[0] < 0 ) m[0] = 0;
00192             if( m[0] == 0 && (m[-1] > 0 || m[1] > 0 || m[-mstep] > 0 || m[mstep] > 0) )
00193             {
00194                 // Find smallest difference to adjacent markers
00195                 const uchar* ptr = img + j*3;
00196                 int idx = 256, t;
00197                 if( m[-1] > 0 )
00198                     c_diff( ptr, ptr - 3, idx );
00199                 if( m[1] > 0 )
00200                 {
00201                     c_diff( ptr, ptr + 3, t );
00202                     idx = ws_min( idx, t );
00203                 }
00204                 if( m[-mstep] > 0 )
00205                 {
00206                     c_diff( ptr, ptr - istep, t );
00207                     idx = ws_min( idx, t );
00208                 }
00209                 if( m[mstep] > 0 )
00210                 {
00211                     c_diff( ptr, ptr + istep, t );
00212                     idx = ws_min( idx, t );
00213                 }
00214 
00215                 // Add to according queue
00216                 assert( 0 <= idx && idx <= 255 );
00217                 ws_push( idx, i*mstep + j, i*istep + j*3 );
00218                 m[0] = IN_QUEUE;
00219             }
00220         }
00221     }
00222 
00223     // find the first non-empty queue
00224     for( i = 0; i < NQ; i++ )
00225         if( q[i].first )
00226             break;
00227 
00228     // if there is no markers, exit immediately
00229     if( i == NQ )
00230         return;
00231 
00232     active_queue = i;
00233     img = src.ptr();
00234     mask = dst.ptr<int>();
00235 
00236     // recursively fill the basins
00237     for(;;)
00238     {
00239         int mofs, iofs;
00240         int lab = 0, t;
00241         int* m;
00242         const uchar* ptr;
00243 
00244         // Get non-empty queue with highest priority
00245         // Exit condition: empty priority queue
00246         if( q[active_queue].first == 0 )
00247         {
00248             for( i = active_queue+1; i < NQ; i++ )
00249                 if( q[i].first )
00250                     break;
00251             if( i == NQ )
00252                 break;
00253             active_queue = i;
00254         }
00255 
00256         // Get next node
00257         ws_pop( active_queue, mofs, iofs );
00258 
00259         // Calculate pointer to current pixel in input and marker image
00260         m = mask + mofs;
00261         ptr = img + iofs;
00262 
00263         // Check surrounding pixels for labels
00264         // to determine label for current pixel
00265         t = m[-1]; // Left
00266         if( t > 0 ) lab = t;
00267         t = m[1]; // Right
00268         if( t > 0 )
00269         {
00270             if( lab == 0 ) lab = t;
00271             else if( t != lab ) lab = WSHED;
00272         }
00273         t = m[-mstep]; // Top
00274         if( t > 0 )
00275         {
00276             if( lab == 0 ) lab = t;
00277             else if( t != lab ) lab = WSHED;
00278         }
00279         t = m[mstep]; // Bottom
00280         if( t > 0 )
00281         {
00282             if( lab == 0 ) lab = t;
00283             else if( t != lab ) lab = WSHED;
00284         }
00285 
00286         // Set label to current pixel in marker image
00287         assert( lab != 0 );
00288         m[0] = lab;
00289 
00290         if( lab == WSHED )
00291             continue;
00292 
00293         // Add adjacent, unlabeled pixels to corresponding queue
00294         if( m[-1] == 0 )
00295         {
00296             c_diff( ptr, ptr - 3, t );
00297             ws_push( t, mofs - 1, iofs - 3 );
00298             active_queue = ws_min( active_queue, t );
00299             m[-1] = IN_QUEUE;
00300         }
00301         if( m[1] == 0 )
00302         {
00303             c_diff( ptr, ptr + 3, t );
00304             ws_push( t, mofs + 1, iofs + 3 );
00305             active_queue = ws_min( active_queue, t );
00306             m[1] = IN_QUEUE;
00307         }
00308         if( m[-mstep] == 0 )
00309         {
00310             c_diff( ptr, ptr - istep, t );
00311             ws_push( t, mofs - mstep, iofs - istep );
00312             active_queue = ws_min( active_queue, t );
00313             m[-mstep] = IN_QUEUE;
00314         }
00315         if( m[mstep] == 0 )
00316         {
00317             c_diff( ptr, ptr + istep, t );
00318             ws_push( t, mofs + mstep, iofs + istep );
00319             active_queue = ws_min( active_queue, t );
00320             m[mstep] = IN_QUEUE;
00321         }
00322     }
00323 }
00324 
00325 
00326 /****************************************************************************************\
00327 *                                         Meanshift                                      *
00328 \****************************************************************************************/
00329 
00330 
00331 void cv::pyrMeanShiftFiltering( InputArray _src, OutputArray _dst,
00332                                 double sp0, double sr, int max_level,
00333                                 TermCriteria termcrit )
00334 {
00335     Mat src0 = _src.getMat();
00336 
00337     if( src0.empty() )
00338         return;
00339 
00340     _dst.create( src0.size(), src0.type() );
00341     Mat dst0 = _dst.getMat();
00342 
00343     const int cn = 3;
00344     const int MAX_LEVELS = 8;
00345 
00346     if( (unsigned)max_level > (unsigned)MAX_LEVELS )
00347         CV_Error( CV_StsOutOfRange, "The number of pyramid levels is too large or negative" );
00348 
00349     std::vector<cv::Mat> src_pyramid(max_level+1);
00350     std::vector<cv::Mat> dst_pyramid(max_level+1);
00351     cv::Mat mask0;
00352     int i, j, level;
00353     //uchar* submask = 0;
00354 
00355     #define cdiff(ofs0) (tab[c0-dptr[ofs0]+255] + \
00356         tab[c1-dptr[(ofs0)+1]+255] + tab[c2-dptr[(ofs0)+2]+255] >= isr22)
00357 
00358     double sr2 = sr * sr;
00359     int isr2 = cvRound(sr2), isr22 = MAX(isr2,16);
00360     int tab[768];
00361 
00362 
00363     if( src0.type() != CV_8UC3 )
00364         CV_Error( CV_StsUnsupportedFormat, "Only 8-bit, 3-channel images are supported" );
00365 
00366     if( src0.type() != dst0.type() )
00367         CV_Error( CV_StsUnmatchedFormats, "The input and output images must have the same type" );
00368 
00369     if( src0.size() != dst0.size() )
00370         CV_Error( CV_StsUnmatchedSizes, "The input and output images must have the same size" );
00371 
00372     if( !(termcrit.type & CV_TERMCRIT_ITER) )
00373         termcrit.maxCount = 5;
00374     termcrit.maxCount = MAX(termcrit.maxCount,1);
00375     termcrit.maxCount = MIN(termcrit.maxCount,100);
00376     if( !(termcrit.type & CV_TERMCRIT_EPS) )
00377         termcrit.epsilon = 1.f;
00378     termcrit.epsilon = MAX(termcrit.epsilon, 0.f);
00379 
00380     for( i = 0; i < 768; i++ )
00381         tab[i] = (i - 255)*(i - 255);
00382 
00383     // 1. construct pyramid
00384     src_pyramid[0] = src0;
00385     dst_pyramid[0] = dst0;
00386     for( level = 1; level <= max_level; level++ )
00387     {
00388         src_pyramid[level].create( (src_pyramid[level-1].rows+1)/2,
00389                         (src_pyramid[level-1].cols+1)/2, src_pyramid[level-1].type() );
00390         dst_pyramid[level].create( src_pyramid[level].rows,
00391                         src_pyramid[level].cols, src_pyramid[level].type() );
00392         cv::pyrDown( src_pyramid[level-1], src_pyramid[level], src_pyramid[level].size() );
00393         //CV_CALL( cvResize( src_pyramid[level-1], src_pyramid[level], CV_INTER_AREA ));
00394     }
00395 
00396     mask0.create(src0.rows, src0.cols, CV_8UC1);
00397     //CV_CALL( submask = (uchar*)cvAlloc( (sp+2)*(sp+2) ));
00398 
00399     // 2. apply meanshift, starting from the pyramid top (i.e. the smallest layer)
00400     for( level = max_level; level >= 0; level-- )
00401     {
00402         cv::Mat src = src_pyramid[level];
00403         cv::Size size = src.size();
00404         const uchar* sptr = src.ptr();
00405         int sstep = (int)src.step;
00406         uchar* mask = 0;
00407         int mstep = 0;
00408         uchar* dptr;
00409         int dstep;
00410         float sp = (float)(sp0 / (1 << level));
00411         sp = MAX( sp, 1 );
00412 
00413         if( level < max_level )
00414         {
00415             cv::Size size1 = dst_pyramid[level+1].size();
00416             cv::Mat m( size.height, size.width, CV_8UC1, mask0.ptr() );
00417             dstep = (int)dst_pyramid[level+1].step;
00418             dptr = dst_pyramid[level+1].ptr() + dstep + cn;
00419             mstep = (int)m.step;
00420             mask = m.ptr() + mstep;
00421             //cvResize( dst_pyramid[level+1], dst_pyramid[level], CV_INTER_CUBIC );
00422             cv::pyrUp( dst_pyramid[level+1], dst_pyramid[level], dst_pyramid[level].size() );
00423             m.setTo(cv::Scalar::all(0));
00424 
00425             for( i = 1; i < size1.height-1; i++, dptr += dstep - (size1.width-2)*3, mask += mstep*2 )
00426             {
00427                 for( j = 1; j < size1.width-1; j++, dptr += cn )
00428                 {
00429                     int c0 = dptr[0], c1 = dptr[1], c2 = dptr[2];
00430                     mask[j*2 - 1] = cdiff(-3) || cdiff(3) || cdiff(-dstep-3) || cdiff(-dstep) ||
00431                         cdiff(-dstep+3) || cdiff(dstep-3) || cdiff(dstep) || cdiff(dstep+3);
00432                 }
00433             }
00434 
00435             cv::dilate( m, m, cv::Mat() );
00436             mask = m.ptr();
00437         }
00438 
00439         dptr = dst_pyramid[level].ptr();
00440         dstep = (int)dst_pyramid[level].step;
00441 
00442         for( i = 0; i < size.height; i++, sptr += sstep - size.width*3,
00443                                           dptr += dstep - size.width*3,
00444                                           mask += mstep )
00445         {
00446             for( j = 0; j < size.width; j++, sptr += 3, dptr += 3 )
00447             {
00448                 int x0 = j, y0 = i, x1, y1, iter;
00449                 int c0, c1, c2;
00450 
00451                 if( mask && !mask[j] )
00452                     continue;
00453 
00454                 c0 = sptr[0], c1 = sptr[1], c2 = sptr[2];
00455 
00456                 // iterate meanshift procedure
00457                 for( iter = 0; iter < termcrit.maxCount; iter++ )
00458                 {
00459                     const uchar* ptr;
00460                     int x, y, count = 0;
00461                     int minx, miny, maxx, maxy;
00462                     int s0 = 0, s1 = 0, s2 = 0, sx = 0, sy = 0;
00463                     double icount;
00464                     int stop_flag;
00465 
00466                     //mean shift: process pixels in window (p-sigmaSp)x(p+sigmaSp)
00467                     minx = cvRound(x0 - sp); minx = MAX(minx, 0);
00468                     miny = cvRound(y0 - sp); miny = MAX(miny, 0);
00469                     maxx = cvRound(x0 + sp); maxx = MIN(maxx, size.width-1);
00470                     maxy = cvRound(y0 + sp); maxy = MIN(maxy, size.height-1);
00471                     ptr = sptr + (miny - i)*sstep + (minx - j)*3;
00472 
00473                     for( y = miny; y <= maxy; y++, ptr += sstep - (maxx-minx+1)*3 )
00474                     {
00475                         int row_count = 0;
00476                         x = minx;
00477                         #if CV_ENABLE_UNROLLED
00478                         for( ; x + 3 <= maxx; x += 4, ptr += 12 )
00479                         {
00480                             int t0 = ptr[0], t1 = ptr[1], t2 = ptr[2];
00481                             if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )
00482                             {
00483                                 s0 += t0; s1 += t1; s2 += t2;
00484                                 sx += x; row_count++;
00485                             }
00486                             t0 = ptr[3], t1 = ptr[4], t2 = ptr[5];
00487                             if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )
00488                             {
00489                                 s0 += t0; s1 += t1; s2 += t2;
00490                                 sx += x+1; row_count++;
00491                             }
00492                             t0 = ptr[6], t1 = ptr[7], t2 = ptr[8];
00493                             if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )
00494                             {
00495                                 s0 += t0; s1 += t1; s2 += t2;
00496                                 sx += x+2; row_count++;
00497                             }
00498                             t0 = ptr[9], t1 = ptr[10], t2 = ptr[11];
00499                             if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )
00500                             {
00501                                 s0 += t0; s1 += t1; s2 += t2;
00502                                 sx += x+3; row_count++;
00503                             }
00504                         }
00505                         #endif
00506                         for( ; x <= maxx; x++, ptr += 3 )
00507                         {
00508                             int t0 = ptr[0], t1 = ptr[1], t2 = ptr[2];
00509                             if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )
00510                             {
00511                                 s0 += t0; s1 += t1; s2 += t2;
00512                                 sx += x; row_count++;
00513                             }
00514                         }
00515                         count += row_count;
00516                         sy += y*row_count;
00517                     }
00518 
00519                     if( count == 0 )
00520                         break;
00521 
00522                     icount = 1./count;
00523                     x1 = cvRound(sx*icount);
00524                     y1 = cvRound(sy*icount);
00525                     s0 = cvRound(s0*icount);
00526                     s1 = cvRound(s1*icount);
00527                     s2 = cvRound(s2*icount);
00528 
00529                     stop_flag = (x0 == x1 && y0 == y1) || std::abs(x1-x0) + std::abs(y1-y0) +
00530                         tab[s0 - c0 + 255] + tab[s1 - c1 + 255] +
00531                         tab[s2 - c2 + 255] <= termcrit.epsilon;
00532 
00533                     x0 = x1; y0 = y1;
00534                     c0 = s0; c1 = s1; c2 = s2;
00535 
00536                     if( stop_flag )
00537                         break;
00538                 }
00539 
00540                 dptr[0] = (uchar)c0;
00541                 dptr[1] = (uchar)c1;
00542                 dptr[2] = (uchar)c2;
00543             }
00544         }
00545     }
00546 }
00547 
00548 
00549 ///////////////////////////////////////////////////////////////////////////////////////////////
00550 
00551 CV_IMPL void cvWatershed( const CvArr* _src, CvArr* _markers )
00552 {
00553     cv::Mat src = cv::cvarrToMat(_src), markers = cv::cvarrToMat(_markers);
00554     cv::watershed(src, markers);
00555 }
00556 
00557 
00558 CV_IMPL void
00559 cvPyrMeanShiftFiltering( const CvArr* srcarr, CvArr* dstarr,
00560                         double sp0, double sr, int max_level,
00561                         CvTermCriteria  termcrit )
00562 {
00563     cv::Mat src = cv::cvarrToMat(srcarr);
00564     const cv::Mat dst = cv::cvarrToMat(dstarr);
00565 
00566     cv::pyrMeanShiftFiltering(src, dst, sp0, sr, max_level, termcrit);
00567 }
00568