opencv之分水岭算法watershed源码注释

转载自：http://blog.csdn.net/u011375993/article/details/46793655

为了研究分水岭算法，阅读了OpenCV 2.4.9 中watershed函数的源码实现部分，代码位于 opencv\sources\modules\imgproc\src\segmentation.cpp 文件中。先贴出加了注解的代码，以后补充对分水岭算法的解释。

#include "precomp.hpp"/*******************************************************                                    Watershed                                    **************************************************************************************/// 结点，用于存储原始图img中像素的偏移量和输出图mask中像素的偏移量typedef struct CvWSNode{    struct CvWSNode* next;    int mask_ofs;    int img_ofs;}CvWSNode;// 队列，用于存储结点 CvWSNodetypedef struct CvWSQueue{    CvWSNode* first;    CvWSNode* last;}CvWSQueue;// 分配空间static CvWSNode*icvAllocWSNodes( CvMemStorage* storage ){    CvWSNode* n = 0;    int i, count = (storage->block_size - sizeof(CvMemBlock))/sizeof(*n) - 1;    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;    return n;}CV_IMPL voidcvWatershed( const CvArr* srcarr, CvArr* dstarr ){    const int IN_QUEUE = -2;        // 加入到队列q中的点定义为 -2    const int WSHED = -1;           // “分水岭”在mask中定义为 -1     const int NQ = 256;             // 队列的数量 256，其实是对应灰度的数量    cv::Ptr<CvMemStorage> storage;    CvMat sstub, *src;    CvMat dstub, *dst;    CvSize size;    CvWSNode* free_node = 0, *node;    CvWSQueue q[NQ];                // 长度为256的CvWSQueue数组，注意数组中每个元素都是一个队列，队列中每个元素是一个节点    int active_queue;               // 指明当前处理的队列，q[active_queue]    int i, j;    int db, dg, dr;    int* mask;                      // 指向标记图像的指针    uchar* img;                     // 指向原始图像的指针    int mstep, istep;               // mstep是mask对应的一行像素数（不是字节数），istep是img对应的一行像素数    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 )            \            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;       \    }    // 求出 ptr1 和 ptr2 指向的像素 r,g,b 差值的最大值    #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 ); \    }    src = cvGetMat( srcarr, &sstub );    dst = cvGetMat( dstarr, &dstub );    // 对参数做检查，要求图像src的类型是8UC3,dst的类型是32SC1，src和dst size相同    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);       // 获取图像的size    storage = cvCreateMemStorage();    // 步长 = 一行字节数 / sizeof(像素数据类型)    istep = src->step;            // img是uchar型， sizeof(uchar) = 1，所以忽略除数      img = src->data.ptr;          // 获取 uchar类型指针    mstep = dst->step / sizeof(mask[0]);    // mask是int（32SC1）型，sizeof(mask[0]) = 4    mask = dst->data.i;           // 获取 int类型指针    memset( q, 0, NQ*sizeof(q[0]) );    // 初始化队列q    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    // 把图像四个边的像素画成分水岭    // mask的首行和末行画成分水岭    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    // 初始阶段：把每个标记的所有邻居像素放到有序队列中去，以确定聚水盆的初始边界    // 即每个标记（种子，全为正值，1,2,3...）都是一个初始聚水盆，标记的周围一圈的邻居像素就是聚水盆的初始边界    // 这里用的是一种逆向思维，不是找标记点，而是判断每一个点是否为标记点的邻居,若是，则该点也被扩充为与标记点同类型的标记点    // 若是多个标记点的邻居，选择梯度最小的标记点的类型，作为该点的标记点类型    for( i = 1; i < size.height-1; i++ )    {        img += istep; mask += mstep;            // 逐行扫描        mask[0] = mask[size.width-1] = WSHED;   // 每一行的首列和末列画成分水岭，加上前面的首行和末行，mask被分水岭方框围起来        for( j = 1; j < size.width-1; j++ )     // 逐列        {            int* m = mask + j;                  // mask的每个像素            if( m[0] < 0 ) m[0] = 0;            // 该点若为负值，先置为零（初始状态下除了四边是分水岭(-1)其余点不应该存在负值？）            if( m[0] == 0 && (m[-1] > 0 || m[1] > 0 || m[-mstep] > 0 || m[mstep] > 0) ) // 若该点为非标记点（0），且四邻域存在标记点（>0）            {                // 求出原图中该点到有标记点的四邻域中，梯度值最小（idx）方向的点，将该点和对应的最小梯度值放入q[idex]队列中                // 两个像素的r,g,b 三个通道中相差最大的值作为像素间的梯度值                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 );     // 将该点在img和mask中的坐标（一维表示）存储在q[idx]队列中                m[0] = IN_QUEUE;         // 在mask中标记该点已入队            }        }    }    // find the first non-empty queue    // 定位到第一个非空的队列    for( i = 0; i < NQ; i++ )        if( q[i].first )            break;    // if there is no markers, exit immediately    // 若i=256，说明数组q中所有队列为空    if( i == NQ )        return;    active_queue = i;    img = src->data.ptr;    mask = dst->data.i;    // recursively fill the basins    // 递归地填满聚水盆    for(;;)    {        int mofs, iofs;         // 将二维图像线性化后图像像素的坐标 mask_offset 和 img_offset 的缩写        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 );     // 从q[active_queue]队列中取出一个结点数据        // 找到这个结点记录的img和mask中的像素点，比较该点在mask中的邻居点        // 邻居点中如果有标记点：该点与邻居点的标记类型不同，则该点为分水岭；该点与邻居点标记类型相同，则该点不变        // 如果有非标记点：将非标记点扩充为标记点        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;            // 如果该像素点的标记类型和邻居像素标记类型都 > 0 且不同，则为分水岭        }        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;        }         // 因为标记点要么是初始种子点，要么是初始阶段延伸的种子点的邻接点         // 该点一定存在一个邻接点是标记点，所以lab一定会赋值一次，不为 0        assert( lab != 0 );          // 若lab > 0 ，则该点被周围的标记点扩充；若lab = -1（WSHED）,则该点定义为分水岭，继续下一个循环              m[0] = lab;                         if( lab == WSHED )            continue;        // lab > 0 的情况，确定该点为标记点，且邻居点中存在未标记点的情况，将未标记点扩充为标记点        if( m[-1] == 0 )        {            c_diff( ptr, ptr - 3, t );                  // 计算梯度t            ws_push( t, mofs - 1, iofs - 3 );           // 将m[-1]这一未标记的点扩充为标记点，进队            active_queue = ws_min( active_queue, t );   // 判断，若t < 当前处理的队列active_queue值，则下一次循环中处理q[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 + istep, t );            ws_push( t, mofs + mstep, iofs + istep );            active_queue = ws_min( active_queue, t );            m[mstep] = IN_QUEUE;        }    }}void cv::watershed( InputArray _src, InputOutputArray markers ){    Mat src = _src.getMat();    CvMat c_src = _src.getMat(), c_markers = markers.getMat();    cvWatershed( &c_src, &c_markers );}