Watershed Algorithm 구현

적용 예는 http://kipl.tistory.com/1         

std::vector<int> sort_pixels(std::vector<BYTE> &img,   std::vector<int> &sorted) ;

// img의 픽셀 위치를 gray값이 커지는 순서로 정렬하고, cumulative histogram chist를 반함함.

// chist은 sorted img에서 특정한 gray값을 갖는 픽셀이 포함된 위치 구간을 알 수 있게 한다.

// 8-way connectivity case only;
// last update: 2024-7-22;
#define MASK -2
#define WSHED 0
#define INIT -1
#define RANGE(i,j,h,w) (((i)>=0) && ((j)>=0) && ((i)<h) && ((j)<w))
// main procedure;
int watershed8(std::vector<BYTE> &img, int width, int height, 
               std::vector<int*> &lab)
{
    int cur_label = 0 ;
    std::queue<int> q;
    std::vector<int> dmap(img.size(), 0);  
    std::vector<int> sorted(img.size());
    std::vector<int> chist = sort_pixels(img, sorted); 
    int hmin = -1, hmax = 255;
    while (!chist[++hmin]);
    while (chist[--hmax]==chist[255]); hmax++; 
    // prepare label image;
    for (int y = 0; y < height; y++)
        for (int x = 0; x < width; x++) lab[y][x] = INIT;
    //
    int *plab = lab[0]; //1-dim representative of label image
    for (int h = hmin; h <= hmax; h++) {
        int *bin_left = &sorted[0] + (h==0 ? 0: chist[h-1]);
        int *bin_right = &sorted[0] + chist[h];
        for (int *curr = bin_left; curr < bin_right; curr++) {
            int pos = *curr;
            plab[pos]=MASK;      
            int x=pos % width;
            int y=pos / width;      
            for (int k=-1;k<=1;k++)
                for (int m=-1;m<=1;m++) 
                    if (RANGE(y+k, x+m, height, width)) 
                        if ((k!=0) || (m!=0))
                            if ((lab[y+k][x+m]>0) || (lab[y+k][x+m]==WSHED)) {
                                q.push(pos);
                                dmap[pos]=1;
                                k=1; break;
                            }
        }      

        int cur_dist=1;
        q.push(-1);
        while (1) {
            int pos = q.front(); q.pop();
            if (pos==-1) {
                if (q.empty()) break;
                else {
                    q.push(-1);
                    cur_dist++;
                    pos = q.front(); q.pop();
                }
            }
            int x=pos % width;
            int y=pos / width;
            for (int k=-1;k<=1;k++)
                for (int m=-1;m<=1;m++) {
                    if (RANGE(y+k, x+m, height, width)) 
                        if ((k!=0) || (m!=0)) {
                            int nn = (y+k)*width + (x+m);// nearest neightbor;
                            if ((dmap[nn]<cur_dist) && ((plab[nn]>0) || (plab[nn]==WSHED))) { 
                                if (plab[nn]>0) {
                                    if ((plab[pos]==MASK) || (plab[pos]==WSHED))
                                        plab[pos]=plab[nn];
                                    else if(plab[pos]!=plab[nn])
                                        plab[pos]=WSHED;
                                }
                                else if(plab[pos]==MASK) plab[pos]=WSHED;
                            } else if((plab[nn]==MASK) && (dmap[nn]==0)) {
                                dmap[nn]=cur_dist+1;
                                q.push(nn);
                            }
                        }
                }
        }  // end while(1)	    	  
        for (int *curr = bin_left; curr < bin_right; curr++) {
            int pos = *curr;
            dmap[pos]=0;	
            if (plab[pos]==MASK) {
                cur_label++;
                q.push(pos);
                plab[pos] = cur_label;
                while (!q.empty()) {
                    int pos1=q.front(); q.pop();
                    int x=pos1 % width;
                    int y=pos1 / width;
                    for (int k=-1;k<=1;k++)
                        for (int m=-1;m<=1;m++)
                            if (RANGE(y+k, x+m, height, width)) 
                                if ((k!=0) || (m!=0)) {
                                    int nn = (y+k)*width + (x+m);
                                    if (plab[nn]==MASK) {
                                        q.push(nn);
                                        plab[nn]=cur_label;
                                    }
                                }
                }
            }
        }
    }  //endfor h    	    	
    
    // convert the result(labeled rgns + boundary) to a suitable form;
    mark_boundary8(lab, width, height);
    return cur_label;
}