Geometry & Recognition

각 전경 픽셀에 대해서 모든 배경 픽셀까지 거리를 구한 다음 최솟값을 할당하면 된다. 픽셀 수가 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;
}

Distance Transform은 이미지 상의 각 픽셀에서 가장 가까이 놓인 물체(전경)의 경계까지 거리를 구하는 연산이다. 이진 영상의 경우 각 배경 픽셀에서 전경까지 거리를 구하는 과정이 된다. 이렇게 얻어진 거리 정보를 이용해서 영상 속 전경 물체의 skeleton을 구하거나 (OCR), 로봇의 팔 등을 자동 제어하는 데 있어서 움직이는 경로에 놓인 방해물의 위치를 예측할 수 있으며, 또한 물체까지의 경로 찾기(path finding) 등에 이용이 된다. distance transform은 사용되는 거리 측도(measure)에 따라 그 특성이 달라진다. 우선 가장 단순한 Euclidean 측도를 이용하도록 하자.

distance transform은 전경의 경계선을 모두 찾아서 배경의 각 픽셀에서 거리의 최솟값을 할당하면 된다. 그러나 경계가 많은 픽셀로 구성이 되는 이미지에서는 연산 비용이 매우 커진다. 여기서는 이미지의 크기에 선형으로 비례하는 알고리즘에 대해서 살펴보도록 하겠다. 픽셀 수에 선형이기 위해서는 알고리즘이 raster 스캔 방식을 따라야 한다.

Euclidean 측도를 사용하는 distance transform은 한 점 $(x, y)$에서 전경 픽셀 $(i, j)$까지 거리의 제곱이 $z = (x-i)^2 + (y-j)^2$로 표현된다. 각 전경 픽셀 꼭짓점으로 하는 paraboloid를 생각하면 distance transform은 주어진 배경 픽셀에서 볼 때 가장 밑에 위치한 paraboloid를 찾는 문제로 환원된다. 그러나 이차원 영역에서 이를 찾기는 쉽지 않은 작업이므로 이를 수평과 수직 방향으로 분리해서 적용할 수 알고리즘을 먼저 만들어야 한다. 우선 $x$-방향으로 스캔하면서 $x$-축으로 사영된 거리를 저장하면서 다시 $y$-축 방향으로 스캔하면서 저장된 $x$-방향의 거리의 제곱의 정보를 이용하면 된다. 이것은  $x$-방향으로 스캔할 때 거리의 제곱 값을 이미지의 그레이 값으로 저장하면 보다 쉽게 해결이 된다. 그러면 $y$-방향으로 스캔할 때는 꼭짓점이 이전에 $x$-방향 스캔 시 저장된 거리의 제곱만큼 offset 된 이차곡선들 중 맨 아래에 놓인 것을 찾으면 된다.  2차원 distance transform은 1차원의 distance transform을 행방향 먼저 수행하고, 그 결과를 다시 열방향으로 시행하면 얻을 수 있다.

1차원의 distance transform은 일직선에 놓인 각 전경 픽셀을 꼭짓점으로 하는 이차 곡선들 중에서 가장 아래에 놓인 곡선을 선택하여 그 값을 취하면 된다. 배경 픽셀은 고려대상이 아니므로 꼭지점이 배경 픽셀이면 이차곡선에는 매우 큰 (수학적으로는 무한대) offset을 준다:

$$\text {1-dim distance transform}(x)=\underset {p\in\text {pixels}}{\text {min}}  \left( (x-p)^2+ F(p) \right)$$

초기 offset $F(p)$ 값은 전경 픽셀은 0, 배경 픽셀은 +INFINITY로 설정한다. 전경의 경계에서 형성이 된 2차 곡선 가장 아래에 놓이므로 그것에 의해서 거리가 정해진다.

일반적인 F(q)의 값이 주어지는 경우 각 픽셀에서 전경까지 거리를 구하는 것은 그 픽셀에서 볼 때 가장 아래에 있는 이차곡선이 어느 것이고, 어는 영역까지 유효한지를 찾는 문제로 귀결이 된다. 픽셀 위치 $p$, $q$에 꼭짓점이 있는 두 이차곡선의 교점이 

$$s=\frac { (F(p) + p^2)-(F(q) + q^2) }{ 2(p-q)}$$

로 표현되므로  주어진 2차 곡선이 가장 낮게 되는 영역은 쉽게 찾을 수 있다. 아래 그림에서 $p<x \le s$ 에서는 꼭짓점이 $p$인 이차곡선이, $s <x <q$ 영역에서는 꼭짓점이 $q$인 이차곡선이 거리를 준다.

그리고 다음 그림은 일차원 이미지에 대한 distance transform의 결과이다: 1차원 이진영상에 대해 적용한 것으로 붉은 선은 이미지(0:전경, 255:배경)이고, 녹색 선은 거리 값을 나타낸다.

이 일차원 distance transform을 각각의 행에 대해서 적용하여 거리를 구한 다음 그 결과를 영상의 픽셀 값으로 저장한다. 다시 각각의 열에 대해서 적용하는데, 이때 주어진 픽셀 $p$에서 영상 값이 $F(p)$의 역할을 하게 된다. 아래의 그림은 글씨에 대해서 skeleton을 구하기 위해서 distance transform을 적용하였다(글씨의 skeleton을 얻기 위해서 negate 영상에 대해서 적용했다: 흰색=배경, 검은색=전경)

** 참고: en.wikipedia.org/wiki/Distance_transform

** cs.brown.edu/people/pfelzens/dt/에서 관련 논문과 원 코드를 찾을 수 있다;

#define SQR(a)	((a) * (a))
/* distanceTransform of 1d function using the Euclidean distance */
struct para_envelop {
    int x;           // parabola의 꼭지점 좌표;
    double start;    // 현재 parabola가 가장 아래 놓인 구간의 시작점;
};
double para_intersection(int q, int p, double *data) {
    return ((data[q] + SQR(q)) - (data[p] + SQR(p))) / (q - p) / 2;
}
// dblimg 이미지의 초기화: 배경 = +INF, 전경 = 0 으로 설정한다.
void DistanceTransform_1D(double *data, int n, double* dist) {
    std::vector<para_envelop> env(n + 1);
    // 가장 낮게 형성이 되는 parabola의 꼭지점과 시작점(유효한 영역 시작점)을 env에 등록한다.
    env[0].x = 0; env[0].start = -FLT_MAX; //첫번째 parabola를 모든 영역에서 가장 낮다고 설정.
    env[1].start = +FLT_MAX;
    int k = 0;
    for (int x = 1; x < n; ++x) {
        // 현 parabola와 이전 오른쪽 가장 낮은 env와 교차점을 찾음;
        double s = para_intersection(x, env[k].x, data);
        // 교차점이 이전 env 꼭지점보다 왼쪽에 있으면 현재 parabola가 더 낮으므로 더 이전의 
        // env와 교차점 검사해서 꼭지점보다 오른쪽에 있을 때까지 계속 인덱스를 내린다.
        while (s <= env[k].start) s = para_intersection(x, env[--k].x, data);
        // env[k]까지는 가장 낮은 parabola꼭지점과 그 영역의 시작점이 기록되어 있음.
        // x > s 이후는 현 parabola가 가장 낮으므로 새로이 env등록을 함;
        env[++k].x = x ;    
        env[k].start = s ;
        env[k + 1].start = +FLT_MAX;
    }
    // 등록된 env 정보를 이용해서 거리정보를 업데이트 함;
    // {env[k]}는 0<=x<n사이에서 가장 낮은 parabola의 꼭지점과 그것이 성립하는 구간의 시작점
    // 정보를 담고 있음. 긱 픽셀에서 가장 낮은 env.start가 그 픽셀보다 왼쪽에 있는 env를 찾아
    // env의 꼭지점까지 거리를 구하면 된다.
    k = 0;
    for (int x = 0; x < n; ++x) {
        // x가 env.x 오른쪽, 다음 env.start의 왼쪽에 놓여야 한다;
        while (env[k + 1].start < x) k++;
        // x에서 가장 가까이 있는 전경 픽셀 위치가 env[k].x 이므로 
        // Euclidean 거리 정보를 업데이트; 
        dist[x] = SQR(x - env[k].x) + data[env[k].x];
    }
};

void DistanceTransform_2D(double *img, int width, int height) {
    int sz = max(width, height);
    std::vector<double> dblimg(sz), dist(sz);
    // transform along rows
    for (int y = 0; y < height; y++) {
        double *ptr = &img[y * width];
        for (int x = 0; x < width; x++)
            dblimg[x] = *ptr++ ? 0: INF;    //전경(white)=0, 배경(black)=INF;
        distanceTransform(&dblimg[0], width, &dist[0]);
        memcpy(&img[y * width], &dist[0], width * sizeof(double));
    };
    // transform along columns
    for (int x = 0; x < width; x++) {
        double *ptr = &img[x];  // copy column;
        for (int y = 0; y < height; y++, ptr += width)
            dblimg[y] = *ptr;
        distanceTransform(&dblimg[0], height, &dist[0]);
        ptr = &img[x];
        for (int y = 0; y < height; y++, ptr += width)
            *ptr = dist[y];
    }
};