文章目录
  1. 1. 分水岭分割简介
  2. 2. Vincent-Soille分水岭算法
  3. 3. OpenCV分水岭算法

本文详细介绍了Vincent-Soille分水岭算法,并对OpenCV实现的分水岭算法源码进行了分析。

分水岭分割简介

在对处理后的图像进行数据分析之前,图像分割时最重要的步骤之一,其主要目的就是为了有效提取感兴趣的目标物体。这种分割可以是基于图像特征直方图,或者是基于边缘的分割,另外一种是基于区域的分割。分水岭算法就是一种基于区域的分割算法。
分水岭(watershed)和集水盆地(catchment basins)的概念在地形学中是人所共知的。分水线分开了每个积水盆地。我们可以将图像数据理解成地形表面,其中灰度值的大小就表示地貌的高度。因此,分水岭分割的过程可以大致上这么理解,在全图范围内降水,盆地部分慢慢被填满,当不同集水盆将要相互连通是,筑起水坝,水坝也就是分水线。当所有的区域都被水淹过后,分水岭也就都建完成了,整个算法结束。
上面的只是很笼统的介绍这个算法的思路,真正算法的实现还是有很多有意思的地方。

Vincent-Soille分水岭算法

首先先来介绍一下Vincent-Soille的分水岭算法
[1]Vincent L. and Soille P. Watersheds in digital spaces: An efficient algorithm based on immersion simulations, 1991.
论文的前面大概表达了这么个意思:“这篇论文是迄今(1991)为止效率最高,精度最好的分水岭分割算法实现”。然后描述了前人的一些工作,还有哪些不足,本篇论文进行了一一改进云云。接下来作者才甩出大招————一段伪代码,全部的代码注释已经添加在伪代码旁边:




伪代码


如果你按照这个伪代码原封不动的写完,算法已经实现了95%了。但是作者在论文并没有给出完全正确的代码,在关键性的步骤里面耍了几个花枪(十分恶劣)。以下给出代码运行结果




药片图像及其距离变换图像




论文代码直译后运行结果(8连通)




修改代码后运行结果(8连通)

上面的处理步骤先对原图二值化,然后计算距离变换图像,接着需要反转距离图像,最后运用分水岭算法进行分割。
在这里并不打算贴出修改后的源码,但是可以提示读者错误出在“设置水坝”的那一段代码中,有兴趣的可以找一找。

上述的分水岭分割并不能直接对原图处理,一般来说都需要对原图进行二值化,然后计算它的距离变换图像再进行处理,即使这样,也是非常容易分割出一些细小的区域。下面介绍另外一种分水岭填充算法。

OpenCV分水岭算法

继续看下上面的例子,如果直接对原图使用分水岭分割算法会有什么结果呢————过分割现象,这也是分水岭算法的一大弊端:




过分割现象

在OpenCV中运行的watershed算法是参考
[2]F. Meyer.“Color image segmentation”1992的论文,很遗憾没有下载到该论文,不过从OpenCV的源码中可以基本了解到算法的原理。不过OpenCV中只支持彩色图像,下面的代码中修改成支持灰度图像:

void watershedImpl( const CvArr* srcarr, CvArr*                 dstarr )
{
    const int IN_QUEUE = -2;
    const int WSHED = -1;
    const int NQ = 256;
    cv::Ptr<CvMemStorage> storage;  

    CvMat sstub, *src;
    CvMat dstub, *dst;
    CvSize size;
    CvWSNode* free_node = 0, *node;
    CvWSQueue q[NQ];
    int active_queue;
    int i, j, ptr_i;
    int db, dg, dr;
    int* mask;
    uchar* img;
    int mstep, istep, icnl;
    int subs_tab[513];//用于查表计算MAX(a-b,0)    

    // 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 )            \
            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;       \
    }   

    // 计算3个通道的时候像素差
    #define c_diff3c(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 ); \
    }   

    // 计算1个通道的时候像素差
    #define c_diff1c(ptr1,ptr2,diff)    \
    {                                   \
        diff = abs((ptr1)[0] - (ptr2)[0]);  \
        assert( 0 <= diff && diff <= 255 ); \
    }   

    #define c_diff(ptr1,ptr2,diff,cnl)  \
    {                                   \
        if (cnl == 1)                   \
            c_diff1c(ptr1,ptr2,diff)    \
        else                            \
            c_diff3c(ptr1,ptr2,diff);   \
    }   

    src = cvGetMat( srcarr, &sstub );
    dst = cvGetMat( dstarr, &dstub );   

    if( CV_MAT_TYPE(src->type) != CV_8UC1 || CV_MAT_TYPE(src->type) != CV_8UC3)
        CV_Error( CV_StsUnsupportedFormat, "Only 8-bit, 1-channel or 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);
    storage = cvCreateMemStorage(); 

    icnl = src->type == CV_8UC1 ? 1 : 3;
    istep = src->step;
    img = src->data.ptr;
    mstep = dst->step / sizeof(mask[0]);//除以4
    mask = dst->data.i;//CvMat中data结构为union类型   

    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;// 将上下两条边设置成分水线   

    // 初始化过程:将未标记的,但是和已标记的像素相邻的像素点加入队列,队列的个数为256,根据距离排列
    // 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 * icnl;
                int idx = 256, t;
                if( m[-1] > 0 )
                    c_diff( ptr, ptr - icnl, idx, icnl );
                if( m[1] > 0 )
                {
                    c_diff( ptr, ptr + icnl, t, icnl );
                    idx = ws_min( idx, t );
                }
                if( m[-mstep] > 0 )
                {
                    c_diff( ptr, ptr - istep, t, icnl );
                    idx = ws_min( idx, t );
                }
                if( m[mstep] > 0 )
                {
                    c_diff( ptr, ptr + istep, t, icnl );
                    idx = ws_min( idx, t );
                }
                assert( 0 <= idx && idx <= 255 );
                ws_push( idx, i*mstep + j, i*istep + j*icnl );
                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 )
        return; 

    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 - icnl, t , icnl);
            ws_push( t, mofs - 1, iofs - icnl );
            active_queue = ws_min( active_queue, t );
            m[-1] = IN_QUEUE;
        }
        if( m[1] == 0 )
        {
            c_diff( ptr, ptr + icnl, t , icnl);
            ws_push( t, mofs + 1, iofs + icnl );
            active_queue = ws_min( active_queue, t );
            m[1] = IN_QUEUE;
        }
        if( m[-mstep] == 0 )
        {
            c_diff( ptr, ptr - istep, t , icnl);
            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 + istep, t , icnl);
            ws_push( t, mofs + mstep, iofs + istep );
            active_queue = ws_min( active_queue, t );
            m[mstep] = IN_QUEUE;
        }
    }
}

OpenCV的这套代码需要传入带分割图像和”markers”图像,markers是标记图像,里面大于0的值表示设置的标记。也就是说,在分割之前,用户需要事先手动设置大概要把图像分割成几块。这里需要注意的是背景区域也需要设置一个标记。这个事先设置好标记个数的算法能有效避免过分割现象。

当传入图像后,首先代码中先把markers的上下左右边界设置成分水岭,这实际上是为了避免这种情况:由于图像数据一行一行存储,那么在图像边界上的点比如mdata[step]这个像素(mdata是markers的data指针,step是一行的长度),它的一个邻接像素mdata[step+1]在实际图像中是另起一行,如果不把它设置成WSHED,就会误认为是邻接像素。

做好处理工作后,接着是选取已经标记的像素邻接的像素点加入到队列queue中,这里的queue是一个256大小的指针数组,按照邻接像素和标记像素的像素差值来加入到对应的指针中。比如像素差5,就加入到queue[4]中。可以看到这种算法不是严格意义上的分水岭算法,它利用的是像素之间的差来作为了标记的依据(像素差小的邻接像素会优先处理)。

初始化完队列,就进入最关键的填充步骤了。这个过程中,按照距离的大小,从0到255遍历队列数组中的所有像素。如果当前像素的邻接像素中只有一个像素已经被标记,那么就把它标记为相同的标记。如果有两个及以上的标记,就把这个像素标记为WSHED,即分水岭。如果这个像素不是分水岭,需要继续处理这个像素的邻接像素。这里可以看出OpenCV的算法是深度优先遍历,而之前的那个算法是广度优先遍历。直到所有的像素都处理完,算法结束。

使用这种算法的优点是可以避免过分割现象,同时考虑了邻近的像素差值这个信息。但是需要事先手动设置好标记的个数,同时深度优先的算法有时会有意向不到的结果。

关于OpenCV的代码运行结果就不在这里上图片了,在官方的文档、实例代码中都可以找到说明,另外很多的博客中都有详细的案例。

文章目录
  1. 1. 分水岭分割简介
  2. 2. Vincent-Soille分水岭算法
  3. 3. OpenCV分水岭算法