X-Git-Url: http://git.maemo.org/git/?a=blobdiff_plain;ds=sidebyside;f=src%2Fcv%2Fcvsegmentation.cpp;fp=src%2Fcv%2Fcvsegmentation.cpp;h=e25dce2914874b68e9c7334e9d6e584da920a4df;hb=e4c14cdbdf2fe805e79cd96ded236f57e7b89060;hp=0000000000000000000000000000000000000000;hpb=454138ff8a20f6edb9b65a910101403d8b520643;p=opencv

diff --git a/src/cv/cvsegmentation.cpp b/src/cv/cvsegmentation.cpp
new file mode 100644
index 0000000..e25dce2
--- /dev/null
+++ b/src/cv/cvsegmentation.cpp
@@ -0,0 +1,559 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "_cv.h"
+
+/****************************************************************************************\
+*                                       Watershed                                        *
+\****************************************************************************************/
+
+typedef struct CvWSNode
+{
+    struct CvWSNode* next;
+    int mask_ofs;
+    int img_ofs;
+}
+CvWSNode;
+
+typedef struct CvWSQueue
+{
+    CvWSNode* first;
+    CvWSNode* last;
+}
+CvWSQueue;
+
+static CvWSNode*
+icvAllocWSNodes( CvMemStorage* storage )
+{
+    CvWSNode* n = 0;
+    
+    CV_FUNCNAME( "icvAllocWSNodes" );
+
+    __BEGIN__;
+
+    int i, count = (storage->block_size - sizeof(CvMemBlock))/sizeof(*n) - 1;
+
+    CV_CALL( n = (CvWSNode*)cvMemStorageAlloc( storage, count*sizeof(*n) ));
+    for( i = 0; i < count-1; i++ )
+        n[i].next = n + i + 1;
+    n[count-1].next = 0;
+
+    __END__;
+
+    return n;
+}
+
+
+CV_IMPL void
+cvWatershed( const CvArr* srcarr, CvArr* dstarr )
+{
+    const int IN_QUEUE = -2;
+    const int WSHED = -1;
+    const int NQ = 256;
+    CvMemStorage* storage = 0;
+    
+    CV_FUNCNAME( "cvWatershed" );
+
+    __BEGIN__;
+
+    CvMat sstub, *src;
+    CvMat dstub, *dst;
+    CvSize size;
+    CvWSNode* free_node = 0, *node;
+    CvWSQueue q[NQ];
+    int active_queue;
+    int i, j;
+    int db, dg, dr;
+    int* mask;
+    uchar* img;
+    int mstep, istep;
+    int subs_tab[513];
+
+    // MAX(a,b) = b + MAX(a-b,0)
+    #define ws_max(a,b) ((b) + subs_tab[(a)-(b)+NQ])
+    // MIN(a,b) = a - MAX(a-b,0)
+    #define ws_min(a,b) ((a) - subs_tab[(a)-(b)+NQ])
+
+    #define ws_push(idx,mofs,iofs)  \
+    {                               \
+        if( !free_node )            \
+            CV_CALL( free_node = icvAllocWSNodes( storage ));\
+        node = free_node;           \
+        free_node = free_node->next;\
+        node->next = 0;             \
+        node->mask_ofs = mofs;      \
+        node->img_ofs = iofs;       \
+        if( q[idx].last )           \
+            q[idx].last->next=node; \
+        else                        \
+            q[idx].first = node;    \
+        q[idx].last = node;         \
+    }
+
+    #define ws_pop(idx,mofs,iofs)   \
+    {                               \
+        node = q[idx].first;        \
+        q[idx].first = node->next;  \
+        if( !node->next )           \
+            q[idx].last = 0;        \
+        node->next = free_node;     \
+        free_node = node;           \
+        mofs = node->mask_ofs;      \
+        iofs = node->img_ofs;       \
+    }
+
+    #define c_diff(ptr1,ptr2,diff)      \
+    {                                   \
+        db = abs((ptr1)[0] - (ptr2)[0]);\
+        dg = abs((ptr1)[1] - (ptr2)[1]);\
+        dr = abs((ptr1)[2] - (ptr2)[2]);\
+        diff = ws_max(db,dg);           \
+        diff = ws_max(diff,dr);         \
+        assert( 0 <= diff && diff <= 255 ); \
+    }
+
+    CV_CALL( src = cvGetMat( srcarr, &sstub ));
+    CV_CALL( dst = cvGetMat( dstarr, &dstub ));
+
+    if( CV_MAT_TYPE(src->type) != CV_8UC3 )
+        CV_ERROR( CV_StsUnsupportedFormat, "Only 8-bit, 3-channel input images are supported" );
+
+    if( CV_MAT_TYPE(dst->type) != CV_32SC1 )
+        CV_ERROR( CV_StsUnsupportedFormat,
+            "Only 32-bit, 1-channel output images are supported" );
+    
+    if( !CV_ARE_SIZES_EQ( src, dst ))
+        CV_ERROR( CV_StsUnmatchedSizes, "The input and output images must have the same size" );
+
+    size = cvGetMatSize(src);
+
+    CV_CALL( storage = cvCreateMemStorage() );
+
+    istep = src->step;
+    img = src->data.ptr;
+    mstep = dst->step / sizeof(mask[0]);
+    mask = dst->data.i;
+
+    memset( q, 0, NQ*sizeof(q[0]) );
+
+    for( i = 0; i < 256; i++ )
+        subs_tab[i] = 0;
+    for( i = 256; i <= 512; i++ )
+        subs_tab[i] = i - 256;
+
+    // draw a pixel-wide border of dummy "watershed" (i.e. boundary) pixels
+    for( j = 0; j < size.width; j++ )
+        mask[j] = mask[j + mstep*(size.height-1)] = WSHED;
+
+    // initial phase: put all the neighbor pixels of each marker to the ordered queue -
+    // determine the initial boundaries of the basins
+    for( i = 1; i < size.height-1; i++ )
+    {
+        img += istep; mask += mstep;
+        mask[0] = mask[size.width-1] = WSHED;
+
+        for( j = 1; j < size.width-1; j++ )
+        {
+            int* m = mask + j;
+            if( m[0] < 0 ) m[0] = 0;
+            if( m[0] == 0 && (m[-1] > 0 || m[1] > 0 || m[-mstep] > 0 || m[mstep] > 0) )
+            {
+                uchar* ptr = img + j*3;
+                int idx = 256, t;
+                if( m[-1] > 0 )
+                    c_diff( ptr, ptr - 3, idx );
+                if( m[1] > 0 )
+                {
+                    c_diff( ptr, ptr + 3, t );
+                    idx = ws_min( idx, t );
+                }
+                if( m[-mstep] > 0 )
+                {
+                    c_diff( ptr, ptr - istep, t );
+                    idx = ws_min( idx, t );
+                }
+                if( m[mstep] > 0 )
+                {
+                    c_diff( ptr, ptr + istep, t );
+                    idx = ws_min( idx, t );
+                }
+                assert( 0 <= idx && idx <= 255 );
+                ws_push( idx, i*mstep + j, i*istep + j*3 );
+                m[0] = IN_QUEUE;
+            }
+        }
+    }
+
+    // find the first non-empty queue
+    for( i = 0; i < NQ; i++ )
+        if( q[i].first )
+            break;
+
+    // if there is no markers, exit immediately
+    if( i == NQ )
+        EXIT;
+
+    active_queue = i;
+    img = src->data.ptr;
+    mask = dst->data.i;
+
+    // recursively fill the basins
+    for(;;)
+    {
+        int mofs, iofs;
+        int lab = 0, t;
+        int* m;
+        uchar* ptr;
+        
+        if( q[active_queue].first == 0 )
+        {
+            for( i = active_queue+1; i < NQ; i++ )
+                if( q[i].first )
+                    break;
+            if( i == NQ )
+                break;
+            active_queue = i;
+        }
+
+        ws_pop( active_queue, mofs, iofs );
+
+        m = mask + mofs;
+        ptr = img + iofs;
+        t = m[-1];
+        if( t > 0 ) lab = t;
+        t = m[1];
+        if( t > 0 )
+        {
+            if( lab == 0 ) lab = t;
+            else if( t != lab ) lab = WSHED;
+        }
+        t = m[-mstep];
+        if( t > 0 )
+        {
+            if( lab == 0 ) lab = t;
+            else if( t != lab ) lab = WSHED;
+        }
+        t = m[mstep];
+        if( t > 0 )
+        {
+            if( lab == 0 ) lab = t;
+            else if( t != lab ) lab = WSHED;
+        }
+        assert( lab != 0 );
+        m[0] = lab;
+        if( lab == WSHED )
+            continue;
+
+        if( m[-1] == 0 )
+        {
+            c_diff( ptr, ptr - 3, t );
+            ws_push( t, mofs - 1, iofs - 3 );
+            active_queue = ws_min( active_queue, t );
+            m[-1] = IN_QUEUE;
+        }
+        if( m[1] == 0 )
+        {
+            c_diff( ptr, ptr + 3, t );
+            ws_push( t, mofs + 1, iofs + 3 );
+            active_queue = ws_min( active_queue, t );
+            m[1] = IN_QUEUE;
+        }
+        if( m[-mstep] == 0 )
+        {
+            c_diff( ptr, ptr - istep, t );
+            ws_push( t, mofs - mstep, iofs - istep );
+            active_queue = ws_min( active_queue, t );
+            m[-mstep] = IN_QUEUE;
+        }
+        if( m[mstep] == 0 )
+        {
+            c_diff( ptr, ptr + 3, t );
+            ws_push( t, mofs + mstep, iofs + istep );
+            active_queue = ws_min( active_queue, t );
+            m[mstep] = IN_QUEUE;
+        }
+    }
+
+    __END__;
+
+    cvReleaseMemStorage( &storage );
+}
+
+
+void cv::watershed( const Mat& src, Mat& markers )
+{
+    CvMat _src = src, _markers = markers;
+    cvWatershed( &_src, &_markers );
+}
+
+
+/****************************************************************************************\
+*                                         Meanshift                                      *
+\****************************************************************************************/
+
+CV_IMPL void
+cvPyrMeanShiftFiltering( const CvArr* srcarr, CvArr* dstarr, 
+                         double sp0, double sr, int max_level,
+                         CvTermCriteria termcrit )
+{
+    const int cn = 3;
+    const int MAX_LEVELS = 8;
+    CvMat* src_pyramid[MAX_LEVELS+1];
+    CvMat* dst_pyramid[MAX_LEVELS+1];
+    CvMat* mask0 = 0;
+    int i, j, level;
+    //uchar* submask = 0;
+
+    #define cdiff(ofs0) (tab[c0-dptr[ofs0]+255] + \
+        tab[c1-dptr[(ofs0)+1]+255] + tab[c2-dptr[(ofs0)+2]+255] >= isr22)
+
+    memset( src_pyramid, 0, sizeof(src_pyramid) );
+    memset( dst_pyramid, 0, sizeof(dst_pyramid) );
+    
+    CV_FUNCNAME( "cvPyrMeanShiftFiltering" );
+
+    __BEGIN__;
+
+    double sr2 = sr * sr;
+    int isr2 = cvRound(sr2), isr22 = MAX(isr2,16);
+    int tab[768];
+    CvMat sstub0, *src0;
+    CvMat dstub0, *dst0;
+
+    CV_CALL( src0 = cvGetMat( srcarr, &sstub0 ));
+    CV_CALL( dst0 = cvGetMat( dstarr, &dstub0 ));
+
+    if( CV_MAT_TYPE(src0->type) != CV_8UC3 )
+        CV_ERROR( CV_StsUnsupportedFormat, "Only 8-bit, 3-channel images are supported" );
+    
+    if( !CV_ARE_TYPES_EQ( src0, dst0 ))
+        CV_ERROR( CV_StsUnmatchedFormats, "The input and output images must have the same type" );
+
+    if( !CV_ARE_SIZES_EQ( src0, dst0 ))
+        CV_ERROR( CV_StsUnmatchedSizes, "The input and output images must have the same size" );
+
+    if( (unsigned)max_level > (unsigned)MAX_LEVELS )
+        CV_ERROR( CV_StsOutOfRange, "The number of pyramid levels is too large or negative" );
+
+    if( !(termcrit.type & CV_TERMCRIT_ITER) )
+        termcrit.max_iter = 5;
+    termcrit.max_iter = MAX(termcrit.max_iter,1);
+    termcrit.max_iter = MIN(termcrit.max_iter,100);
+    if( !(termcrit.type & CV_TERMCRIT_EPS) )
+        termcrit.epsilon = 1.f;
+    termcrit.epsilon = MAX(termcrit.epsilon, 0.f);
+
+    for( i = 0; i < 768; i++ )
+        tab[i] = (i - 255)*(i - 255);
+
+    // 1. construct pyramid
+    src_pyramid[0] = src0;
+    dst_pyramid[0] = dst0;
+    for( level = 1; level <= max_level; level++ )
+    {
+        CV_CALL( src_pyramid[level] = cvCreateMat( (src_pyramid[level-1]->rows+1)/2,
+                        (src_pyramid[level-1]->cols+1)/2, src_pyramid[level-1]->type ));
+        CV_CALL( dst_pyramid[level] = cvCreateMat( src_pyramid[level]->rows,
+                        src_pyramid[level]->cols, src_pyramid[level]->type ));
+        CV_CALL( cvPyrDown( src_pyramid[level-1], src_pyramid[level] ));
+        //CV_CALL( cvResize( src_pyramid[level-1], src_pyramid[level], CV_INTER_AREA ));
+    }
+
+    CV_CALL( mask0 = cvCreateMat( src0->rows, src0->cols, CV_8UC1 ));
+    //CV_CALL( submask = (uchar*)cvAlloc( (sp+2)*(sp+2) ));
+
+    // 2. apply meanshift, starting from the pyramid top (i.e. the smallest layer)
+    for( level = max_level; level >= 0; level-- )
+    {
+        CvMat* src = src_pyramid[level];
+        CvSize size = cvGetMatSize(src);
+        uchar* sptr = src->data.ptr;
+        int sstep = src->step;
+        uchar* mask = 0;
+        int mstep = 0;
+        uchar* dptr;
+        int dstep;
+        float sp = (float)(sp0 / (1 << level));
+        sp = MAX( sp, 1 );
+
+        if( level < max_level )
+        {
+            CvSize size1 = cvGetMatSize(dst_pyramid[level+1]);
+            CvMat m = cvMat( size.height, size.width, CV_8UC1, mask0->data.ptr );
+            dstep = dst_pyramid[level+1]->step;
+            dptr = dst_pyramid[level+1]->data.ptr + dstep + cn;
+            mstep = m.step;
+            mask = m.data.ptr + mstep;
+            //cvResize( dst_pyramid[level+1], dst_pyramid[level], CV_INTER_CUBIC );
+            cvPyrUp( dst_pyramid[level+1], dst_pyramid[level] );
+            cvZero( &m );
+
+            for( i = 1; i < size1.height-1; i++, dptr += dstep - (size1.width-2)*3, mask += mstep*2 )
+            {
+                for( j = 1; j < size1.width-1; j++, dptr += cn )
+                {
+                    int c0 = dptr[0], c1 = dptr[1], c2 = dptr[2];
+                    mask[j*2 - 1] = cdiff(-3) || cdiff(3) || cdiff(-dstep-3) || cdiff(-dstep) ||
+                        cdiff(-dstep+3) || cdiff(dstep-3) || cdiff(dstep) || cdiff(dstep+3);
+                }
+            }
+
+            cvDilate( &m, &m, 0, 1 );
+            mask = m.data.ptr;
+        }
+
+        dptr = dst_pyramid[level]->data.ptr;
+        dstep = dst_pyramid[level]->step;
+
+        for( i = 0; i < size.height; i++, sptr += sstep - size.width*3,
+                                          dptr += dstep - size.width*3,
+                                          mask += mstep )
+        {   
+            for( j = 0; j < size.width; j++, sptr += 3, dptr += 3 )
+            {               
+                int x0 = j, y0 = i, x1, y1, iter;
+                int c0, c1, c2;
+
+                if( mask && !mask[j] )
+                    continue;
+
+                c0 = sptr[0], c1 = sptr[1], c2 = sptr[2];
+
+                // iterate meanshift procedure
+                for( iter = 0; iter < termcrit.max_iter; iter++ )
+                {
+                    uchar* ptr;
+                    int x, y, count = 0;
+                    int minx, miny, maxx, maxy;
+                    int s0 = 0, s1 = 0, s2 = 0, sx = 0, sy = 0;
+                    double icount;
+                    int stop_flag;
+
+                    //mean shift: process pixels in window (p-sigmaSp)x(p+sigmaSp)
+                    minx = cvRound(x0 - sp); minx = MAX(minx, 0);
+                    miny = cvRound(y0 - sp); miny = MAX(miny, 0);
+                    maxx = cvRound(x0 + sp); maxx = MIN(maxx, size.width-1);
+                    maxy = cvRound(y0 + sp); maxy = MIN(maxy, size.height-1);
+                    ptr = sptr + (miny - i)*sstep + (minx - j)*3; 
+
+                    for( y = miny; y <= maxy; y++, ptr += sstep - (maxx-minx+1)*3 )
+                    {
+                        int row_count = 0;
+                        x = minx;
+                        for( ; x + 3 <= maxx; x += 4, ptr += 12 )
+                        {
+                            int t0 = ptr[0], t1 = ptr[1], t2 = ptr[2];
+                            if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )
+                            {                        
+                                s0 += t0; s1 += t1; s2 += t2;
+                                sx += x; row_count++;
+                            }
+                            t0 = ptr[3], t1 = ptr[4], t2 = ptr[5];
+                            if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )
+                            {                        
+                                s0 += t0; s1 += t1; s2 += t2;
+                                sx += x+1; row_count++;
+                            }
+                            t0 = ptr[6], t1 = ptr[7], t2 = ptr[8];
+                            if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )
+                            {                        
+                                s0 += t0; s1 += t1; s2 += t2;
+                                sx += x+2; row_count++;
+                            }
+                            t0 = ptr[9], t1 = ptr[10], t2 = ptr[11];
+                            if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )
+                            {                        
+                                s0 += t0; s1 += t1; s2 += t2;
+                                sx += x+3; row_count++;
+                            }
+                        }
+                        
+                        for( ; x <= maxx; x++, ptr += 3 )
+                        {      
+                            int t0 = ptr[0], t1 = ptr[1], t2 = ptr[2];
+                            if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )
+                            {                        
+                                s0 += t0; s1 += t1; s2 += t2;
+                                sx += x; row_count++;
+                            }
+                        }
+                        count += row_count;
+                        sy += y*row_count;
+                    }
+
+                    if( count == 0 )
+                        break;
+
+                    icount = 1./count;
+                    x1 = cvRound(sx*icount);
+                    y1 = cvRound(sy*icount);
+                    s0 = cvRound(s0*icount);
+                    s1 = cvRound(s1*icount);
+                    s2 = cvRound(s2*icount);
+
+                    stop_flag = (x0 == x1 && y0 == y1) || abs(x1-x0) + abs(y1-y0) +
+                        tab[s0 - c0 + 255] + tab[s1 - c1 + 255] +
+                        tab[s2 - c2 + 255] <= termcrit.epsilon;
+                
+                    x0 = x1; y0 = y1;
+                    c0 = s0; c1 = s1; c2 = s2;
+
+                    if( stop_flag )
+                        break;
+                }
+
+                dptr[0] = (uchar)c0;
+                dptr[1] = (uchar)c1;
+                dptr[2] = (uchar)c2;
+            }
+        }
+    }   
+
+    __END__;
+
+    for( i = 1; i <= MAX_LEVELS; i++ )
+    {
+        cvReleaseMat( &src_pyramid[i] );
+        cvReleaseMat( &dst_pyramid[i] );
+    }
+    cvReleaseMat( &mask0 );
+}
+