Watershed Algorithm 구현

Watershed 알고리즘은 이미지를 여러 개의 영역으로 분할(segmentation)하는 대표적인 방법이다. 이름 그대로 지형의 유역(watershed)을 찾는 과정을 이미지에 적용한 알고리즘이다.

이미지의 밝기값을 3차원 지형으로 생각하면, 밝은 픽셀은 높은 산, 어두운 픽셀은 낮은 골짜기에 해당한다. 이 지형에 물을 채운다고 가정하면, 물은 각 골짜기에서부터 차오르기 시작한다. 서로 다른 골짜기에서 시작된 물이 만나기 직전에 둑을 쌓으면, 이 둑이 각 영역을 구분하는 경계가 된다. 이것이 Watershed 알고리즘의 기본 원리이다.

실제 영상에서는 원본 영상보다 에지의 크기(gradient magnitude)를 지형으로 사용하는 경우가 많다. 이 경우 물체 내부는 하나의 분지(basin)를 이루고, 에지 부분은 높은 능선(ridge)이 되어 물체의 경계를 효과적으로 찾을 수 있다.

그러나 잡음이 많은 영상에서는 작은 골짜기가 많이 생겨 하나의 물체가 여러 영역으로 나누어지는 과분할(over-segmentation) 문제가 발생한다. 따라서 일반적으로는 영상을 평활화하거나 marker를 이용하는 Marker-controlled Watershed 기법을 함께 사용한다.

Watershed 알고리즘은 서로 붙어 있는 물체를 분리하는 데 매우 효과적이며, 세포 분할, 의료 영상, 산업용 비전 등 다양한 영상 처리 분야에서 널리 활용되고 있다. 특히 Distance Transform과 함께 사용하면 서로 맞닿아 있는 원형 물체들을 효과적으로 분리할 수 있어 실제 응용에서 매우 자주 사용되는 영상 분할 기법이다.

적용 예는 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;
}