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

Ref: Efficient Graph-Based Image Segmentation, Pedro F. Felzenszwalb and Daniel P. Huttenlocher, International Journal of Computer Vision, 59(2) September 2004.

struct edge {
    float w;	//weight=color_distance
    int a, b;   //vertices=pixel positions;
    edge() {}
    edge(int _a, int _b, float _w): a(_a), b(_b), w(_w) {}
};
bool operator<(const edge &a, const edge &b) {
  return a.w < b.w;
}
#define SQR(x) ((x)*(x)) 
#define DIFF(p,q) \
    (sqrt(SQR(red[(p)]-red[(q)])+SQR(green[(p)]-green[(q)])+SQR(blue[(p)]-blue[(q)])))
int segment_image(CRaster& raster, float c, int min_size, CRaster& out) {
    if (raster.GetBPP() != 24) return 0;
    const CSize sz = raster.GetSize();
    const int width = sz.cx, height = sz.cy;
    std::vector<float> red(width * height);
    std::vector<float> green(width * height);
    std::vector<float> blue(width * height);
    for (int y = 0, curr = 0; y <height; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        for (int x = 0; x < width; x++, curr++) 
            blue[curr] = *p++, green[curr] = *p++, red[curr] = *p++;
    }
    // gaussian smoothing;
    const float sigma = 0.5F;
    smooth(blue,  width, height, sigma);
    smooth(green, width, height, sigma);
    smooth(red,   width, height, sigma);

    std::vector<edge> edges;
    edges.reserve(width * height * 4);
    for (int y = 0, curr = 0; y < height; y++) {
        for (int x = 0; x < width; x++, curr++) {
            if (x < width-1)
                edges.push_back(edge(curr, curr+1, DIFF(curr, curr+1)));
            if (y < height-1) 
                edges.push_back(edge(curr, curr+width, DIFF(curr, curr+width)));
            if ((x < width-1) && (y < height-1))
                edges.push_back(edge(curr, curr+width+1, DIFF(curr, curr+width+1)));
            if ((x < width-1) && (y > 0)) 
                edges.push_back(edge(curr, curr-width+1, DIFF(curr, curr-width+1)));
        }
    }

    if (c < 0) c = 1500; 
    if (min_size < 0) min_size = 10;
    universe *u = segment_graph(width*height, edges, c);
    // join small size regions;
    for (int i = edges.size(); i--> 0;) {
        edge &e = edges[i];
        int a = u->find(e.a), b = u->find(e.b);
        if ((a != b) && ((u->size(a) < min_size) || (u->size(b) < min_size)))
            u->join(a, b);
    }
    int num_rgns = u->count();

    // paint segmented rgns with root pixel color;
    out = raster;
    for (int y = height-1; y-->0;) 
        for (int x = width-1; x-->0;) {
            int a = u->find(y * width + x);
            out.SetPixel0(x, y, RGB(blue[a],green[a],red[a]));
        } 
    delete u;
    return num_rgns;
}
// Segment a graph
// c: constant for treshold function.
universe *segment_graph(int num_vertices, std::vector<edge>& edges, float c) { 
    // sort by weight
    std::sort(edges.begin(), edges.end());
    // disjoint-set
    universe *u = new universe(num_vertices);
    // init thresholds = {c};
    std::vector<float> threshold(num_vertices, c);

    for (int i = 0; i < edges.size(); i++) {
        edge &e = edges[i];
        int a = u->find(e.a), b = u->find(e.b);
        if (a != b)
            if ((e.w <= threshold[a]) && (e.w <= threshold[b])) {
                u->join(a, b);
                a = u->find(a);
                threshold[a] = e.w + c / u->size(a);
            }
    }
    return u;
};

https://kipl.tistory.com/460

주어진 데이터를 잘 피팅하는 직선을 찾기 위해서는 데이터를 이루는 각 점의 $y$ 값과 같은 위치에서 구하려는 직선과의 차이 (residual)의 제곱을 최소화시키는 조건을 사용한다. 그러나 직선의 기울기가 매우 커지는 경우에는  데이터에 들어있는 outlier에 대해서는 그 차이가 매우 커지는 경우가 발생할 수도 있어 올바른 피팅을 방해할 수 있다. 이를 수정하는 간단한 방법 중 하나는 $y$값의 차이를 이용하지 않고 데이터와 직선 간의 최단거리를 이용하는 것이다. 

수평에서 기울어진 각도가 $\theta$이고 원점에서 거리가 $s$인 직선의 방정식은 

$$x\sin \theta -y\cos \theta +s=0$$

이고, 한 점 $(x_i,y_i)$에서 이 직선까지 거리는

$$d_i=|\sin \theta y_i-\cos \theta x_i+s|$$

이다. 따라서 주어진 데이터에서 떨어진 거리 제곱의 합이 최소인 직선을 구하기 위해서는 다음을 최소화시키는 $\theta ,s$을 구해야 한다. $$L=\sum _i(\sin \theta x_i-\cos \theta y_i+s{)}^2$$$$(\theta ,s)=\text{argmin}(L)$$

$\theta$와 $s$에 대한 극값 조건에서 

$$\frac{1}{2}\frac{\partial L}{\partial \theta }=\frac{1}{2}\sin 2\theta \sum _i(x_i^2-y_i^2)-\cos 2\theta \sum x_i{y}_i+s\sin \theta \sum _i{x}_i+s\cos \theta \sum _i{y}_i=0$$

$$\frac{1}{2}\frac{\partial L}{\partial s}=\cos \theta \sum y_i-\sin \theta \sum _i{x}_i-Ns=0$$

주어진 데이터의 질량중심계에서 계산을 수행하면 $\sum _i{x}_i=\sum _i{y}_i=0$ 이므로 데이터의 2차 모멘트를 $$A=\sum _i(x_i^2-y_i^2),B=\sum _i{x}_i{y}_i$$로 놓으면 직선의 파라미터를 결정하는 식은

$$\begin{gathered} \frac{1}{2}A\sin 2\theta -B\cos 2\theta =0\to \tan 2\theta =\frac{2B}{A} \\ s=0 \end{gathered}$$

두 번째 식은 직선이 질량중심(질량중심계에서 원점)을 통과함을 의미한다. 첫번째 식을 풀면

$$\tan \theta =\frac{-A\pm \sqrt{A^2+(2B{)}^2}}{2B}$$

두 해 중에서 극소값 조건을 만족시키는 해가 직선을 결정한다. 그런데

$$\frac{1}{2}\frac{{\partial }^{2}L}{\partial {\theta }^2}=A\cos 2\theta +2B\sin 2\theta =\pm \sqrt{A^2+(2B{)}^2}>0$$

이므로 위쪽 부호로 직선($x\sin \theta =y\cos \theta$)이 정해진다. 질량중심계에서는 원점을 지나지만 원좌표계로 돌아오면 데이터의 질량중심을 통과하도록 평행이동시키면 된다.

$$(-A+\sqrt{A^2+(2B{)}^2})(x-\overline{x})=2B(y-\overline{y})$$

여기서 주어진 데이터의 질량중심은 원좌표계에서

$$\overline{x}=\frac{1}{N}\sum _i{x}_i,\overline{y}=\frac{1}{N}\sum _i{y}_i$$

이다. 또한 원좌표계에서 $A$와 $B$의 계산은 

$$A=\sum _i[(x_i-\overline{x}{)}^2-(y_i-\overline{y}{)}^2],B=\sum (x_i-\overline{x})(y_i-\overline{y})$$

이 결과는 데이터 분포에 PCA를 적용해서 얻은 결과와 동일하다. PCA에서는 공분산 행렬의 고유값이 큰 쪽에 해당하는 고유벡터가 직선의 방향을 결정했다. https://kipl.tistory.com/211.  또한 통계에서는 Deming regression이라고 불린다.