Geometry & Recognition

배경의 각 픽셀에 인접 픽셀까지의 수평/수직 방향으로는 Euclidean distance가 1 만큼 떨어져 있고, 대각선 방향으로는 $\sqrt{2}$만큼 떨어져 있다. 이 거리에 $\sqrt{8}\approx 2.83\approx 3$을 곱하면 거리를 각각 $3$과 $4$로 근사할 수 있어서 정수 연산 만으로 distance transform을 구할 수 있다.

#define BACKGROUND (0)
void distanceTransform_approx(BYTE **image, int w, int h, int **map) {
    // initialization;
    for (int y = h; y-->0;)
        for (int x = w; x-->0;) map[y][x] = 0;
    // forward_pass
    // y=x=0;
    if (image[0][0] == BACKGROUND) map[0][0] = 0xFFFF;//  INFINITY;
    // y=0;
    for (int x=1; x<w; x++)
        if (image[0][x] == BACKGROUND) map[0][x] = 3 + map[0][x-1];
    for (int y=1; y<h; y++) {
        // x=0;
        if (image[y][0] == BACKGROUND) 
            map[y][0] = min(3 + map[y-1][0], 4 + map[y-1][1]);
        for (int x=1; x<w-1; x++) 
            if (image[y][x] == BACKGROUND)
                map[y][x] = min(4 + map[y-1][x-1], min(3 + map[y-1][x], 
                           min(4 + map[y-1][x+1], 3 + map[y][x-1])));
        // x=w-1;
        if (image[y][w-1] == BACKGROUND) 
            map[y][w-1] = min(4 + map[y-1][w-2], min(3 + map[y-1][w-1], 3 + map[y][w-2]));
    }
    // backward_pass
    // y=h-1;
    for (int x = w-1; x-->0;)
        if (image[h-1][x] == BACKGROUND) 
            map[h-1][x] = min(map[h-1][x], 3 + map[h-1][x+1]);

    for (int y = h-1; y-->0;) {
        // x=w-1;
        if (image[y][w-1] == BACKGROUND)
            map[y][w-1] = min(map[y][w-1], min(3 + map[y+1][w-1], 4 + map[y+1][w-2]));
        for (int x=w-1; x-->1;) 
            if (image[y][x] == BACKGROUND)
                map[y][x] = min(map[y][x], min(4 + map[y+1][x+1], 
                           min(3 + map[y+1][x], min(4 + map[y+1][x-1], 3 + map[y][x+1]))));
        // x=0;
        if (image[y][0] == BACKGROUND)
            map[y][0] = min(map[y][0], min(4 + map[y+1][1], 
                            min(3 + map[y+1][0], 3 + map[y][1])));
    }
};

각 전경 픽셀에 대해서 모든 배경 픽셀까지 거리를 구한 다음 최솟값을 할당하면 된다. 픽셀 수가 n = w * h 개면 time complexity는 O(n^2)이다 (이미지 폭이나 높이의 4승이므로 연산량이 상당하다.) linear time Euclidean distance transform은 kipl.tistory.com/268.

ListPlot3D[ Reverse@ImageData[pic], PlotRange -> Full]

int brute_force_edt(double *img, int w, int h) {
    int *list = new int [w * h];
    int fg_cnt = 0;
    int bg_ptr = w * h;               //배경 픽셀 위치는 끝에서 역으로 채움;
    //전경과 배경 픽셀 위치를 분리;
    for (int i = w * h; i-->0;)
        if (img[i]) list[fg_cnt++] = i;  // foreground;
        else        list[--bg_ptr] = i;  // background;

    for (int i = 0; i < fg_cnt; ++i) {     // 전경 픽셀;
        int xi = list[i] % w, yi = list[i] / w;
        int d2min = INT_MAX;
        for (int j = w * h; j--> fg_cnt;) { // 배경 픽셀까지 거리를 계산해서 최소값을 할당;
                                            // 배경이 list에 역순으로 저장되으로 역방향 서치;
            int dx = (list[j] % w) - xi, dy = (list[j] / w) - yi;
            int dst = dx * dx + dy * dy;
            if (dst == 1) {
            	d2min = 1; break;
            }
            if (d2min > dst) d2min = dst;
        }
        img[list[i]] = d2min;
    }
    delete [] list;
    return fg_cnt;
}

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