Watershed Segmentation

Watershed Segmentation은 영상을 지형(topographic surface)으로 해석하여 영역을 분할하는 대표적인 영상 분할 알고리즘이다. 일반적으로 영상의 밝기 값이나 에지의 크기를 지형의 높이로 간주하며, 밝기가 낮은 영역은 계곡(basin), 밝기가 높은 영역은 능선(ridge)에 해당한다고 생각한다.

지형의 가장 낮은 지점인 local minimum에서부터 물이 차오른다고 가정하면, 물은 수위가 높아짐에 따라 각 계곡을 따라 점차 퍼져 나간다. 이 과정에서 서로 다른 계곡에서 흘러온 물이 만나려고 하면 두 영역이 합쳐지는 것을 막기 위해 그 경계에 둑을 쌓는다. 물이 영상 전체를 덮을 때까지 이러한 과정을 반복하면, 최종적으로 형성된 둑이 서로 다른 영역을 구분하는 경계, 즉 watershed가 된다.

Watershed Segmentation은 닫힌 경계를 자동으로 생성하므로 서로 맞닿아 있는 물체를 분리하는 데 매우 효과적이다. 그러나 잡음이나 작은 밝기 변화에도 많은 국소 최소점이 생성되어 수많은 작은 영역으로 분할되는 over-segmentation 문제가 발생하기 쉽다. 이를 해결하기 위해 실제 응용에서는 영상을 먼저 평활화(smoothing)하여 불필요한 국소 최소점을 줄이거나, 사용자가 지정한 marker만을 시작점으로 사용하는 marker-controlled watershed 기법을 이용한다. 이 방법은 과분할을 크게 줄이면서도 원하는 물체를 안정적으로 분할할 수 있다. 

아래의 예제는 이진 영상의 distance map을 계산한 후 이를 입력으로 watershed 알고리즘을 적용하여, 서로 맞닿아 있는 객체를 분리한 결과를 보여준다.

segmentation result

1. gradient 대신 negate시킨 distance map을 사용한다.
2. 경계영역 처리를 보완해야 한다.
3. distance-map을 충분히 smooth하게 만들어야 한다: mean-filter 반복 적용.
4. 참고: https://kipl.tistory.com/66의 구현보다 단순하지만 성능은 떨어진다.

#define insideROI(x, y) ((x) >= 0 && (x) < w && (y) >= 0 && (y) < h)
#define BACKGND 0xFF
int watershed(float *distMap, int w, int h, int *rgnMap) {
    const int xend = w - 1, yend = h - 1;
    const int nx[] =  {-1, 0, 1, -1, 0, 1, -1, 0, 1};
    const int ny[] =  {-1, -1, -1, 0, 0, 0, 1, 1, 1};
    // processing image;
    std::vector<float>   wbuffer(w * h);
    std::vector<float* > procImg(h); 
    for (int y = 0; y < h; y++)
        procImg[y] = &wbuffer[y * w];
        
    // presmooth; 10 iterations;
    mean_filter(distMap, w, h, 10, procImg); 

    // 입력 영상에서 경계는 global maximum으로 설정;
    for (int y = 0; y < h; y++ ) procImg[y][0] = procImg[y][xend] = BACKGND;
    for (int x = 0; x < w; x++ ) procImg[0][x] = procImg[yend][x] = BACKGND;
    // find local minima for each pixel;
    int minCnt = 1;
    for (int y = 0, pos = 0; y < h; y++) {
        for (int x = 0; x < w; x++, pos++) {
            float minVal = procImg[y][x];
            if (minVal == BACKGND) {
                rgnMap[pos] = 1;	// background(white);
                continue;
            }
            int minIdx = 4; //current position;
            for (int k = 0; k < 9; k++) {
                int xx = x + nx[k], yy = y + ny[k];
                if (insideROI(xx, yy) && procImg[yy][xx] <= minVal) {
                    minVal = procImg[yy][xx];
                    minIdx = k;
                }
            }
            if (minIdx == 4) rgnMap[pos] = ++minCnt; // center(x, y) is a new local minimum;
            else             rgnMap[pos] = -minIdx;  // direction of local-min =  "-(dir)"
        }
    }
    // follow gradient downhill for each pixel;
    // reset rgnMap to have the id's of local minimum connected with current pixel;
    for (int y = 1, pos = y * w; y < yend; y++ ) {
        pos++;
        for (int x = 1; x < xend; x++, pos++ ) {
            int minpos = pos;
            while ( rgnMap[minpos] <= 0 ) { //it is not a local minimum.
                switch ( rgnMap[minpos] ) {
                case  0: minpos += -w - 1; break; // top-lef = downhill direction;
                case -1: minpos += -w;	   break; // top
                case -2: minpos += -w + 1; break; // top-right;
                case -3: minpos--;         break; // left;
                case -5: minpos++;         break; // right;
                case -6: minpos += w - 1;  break; // bot-left;
                case -7: minpos += w;      break; // bot;
                case -8: minpos += w + 1;  break; // bot-right;
                }
            }
            // speed-up: use a stack to record downhill path.
            //assign the id of a local min connected with current pixel;
            rgnMap[pos] = rgnMap[minpos]; 
        }
        pos++;
    }
    // rgnMap는 각 픽셀과 연결되는 국소점의 id을 알려주는 lookup table이다;
    return ( minCnt );
}