Geometry & Recognition

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

지상에서 포물선 운동을 하는 물체의 속도 벡터를  모두 모아 시작을 같게 만들면 속도 벡터의 끝이 그리는 자취(hodograph라 함)는 직선이 된다. 이는 속도의 차이가 가속도이고 지상에서 중력가속도는 크기와 방향이 일정하기 때문이다.

$$\vec {v}(t) =\vec {v}_0 - g\hat {j} t$$

태양계에서 행성은 태양을 한 초점으로 하는 타원 궤도상에서 움직인다. 그러면 행성의 각 지점에서 속도 벡터의 시작점을 한 지점에 모을 때 벡터의 끝점이 그리는 궤적은 무엇일까? 원 궤도를 그리면 속력이 일정하므로 당연히 원이 될 것으로 예상할 수 있지만,  타원 궤도에서는 속력은 에너지 보존을 고려하면 일정할 수 없다. 그런데 타원궤도를 그리는 경우에도 속도의 hodograph는 원궤도에서와 마찬가지로 벡터 공간에서 원으로 표현된다. 왜 그럴까?

행성의 운동 방정식에서 가속도는 

$$ \frac {d\vec {v}}{dt}= -\frac {GM}{r^2 }\hat {r} \quad \rightarrow\quad |\Delta \vec {v}| =\frac {GM}{r^2} \Delta t$$

이고, 각운동량 보존법칙에서 (각 $\varphi$는 태양을 원점으로 잼)

$$ L =m |\vec{r}\times \vec {v}| = m r^2 \frac {d\varphi}{dt}  \quad\rightarrow\quad \Delta\varphi =\frac {L}{mr^2}\Delta t$$

두 식에서 $r^2$을 소거하면,

$$ |\Delta \vec{v}| = \frac {GMm}{L} \Delta \varphi$$ 

각운동량의 크기가 일정하므로 속도의 증가분이 회전각의 증가분과 비례함을 의미하고, 기하학적으로는 속도 벡터를 모아둔 벡터 공간에서 각이 증가하면서 움직이는 벡터의 끝이 반지름 $GMm/L$인 원주 위에 있을 때만 가능한 관계이다. 벡터 공간에서 벡터의 시작점이 각을 재는 원점이 아니다. hodograph의 반지름은 근일점에서 속력 $v_P$와 원일점에서 속력 $v_A$의 절반임으로 표현되는데 $v_R = (v_A + v_P)/2 = \frac {GMm}{L}$임을 근일점, 원일점에서 역학적 에너지 보존과 각운동량 보존식을  이용해서 확인할 수 있다.

설명 동영상:

youtu.be/xdIjYBtnvZU