Geometry & Recognition

화면이나 그림에서 직선을 그리려면 직선의 방정식을 만족하는 모든 픽셀들을 찾으면 된다. 하지만 픽셀의 위치는 연속적인 값이 아니라  이산적이므로 그림에 그려진 직선을 구성하는 픽셀은 그 직선을 표현을 하는 식이 주는 값에 가장 가까운 정수를 찾아서 그린 것이다. 직선을 그리기 위해서 실수 연산을 한 후 정수로 반올림하는 과정이 들어가게 된다. 

Bresenham 알고리즘은 평면 위의 두 점 $(x_0,y_0),(x_1,y_1)$을 연결하는 직선을 정수 연산만을 사용하여 그릴 수 있음을 보여준다. 직선의 방정식은 $y=mx+b$의 형태로 표현이 되고, 두 점이 주어지는 경우 기울기는 $m=\frac{y_1-y_0}{x_1-x_0}=\frac{\mathrm{\Delta }y}{\mathrm{\Delta }x}$, $y-$절편은 $b=y_0-m{x}_0$로 얻을 수 있다. 직선의 기울기가 1보다 크면 $x$와 $x+1$에서 변할 때 종속변수 $y$값의 차이가 커지므로 선이 연속적으로 그려지지 않게 보인다.  이 경우 $y$를 독립변수로 사용하면 되고, 기울기가 음수면 시작과 끝을 바꾸면 된다. 따라서 기울기가 $0\le m\le 1$인 경우만 살펴보면 충분하다. 

이제, $(x_k,y_k)$ 위치에서 직선의 픽셀이 정해진 경우,  $x_{k+1}=x_k+1$에서 직선식의 값 $y=m(x_k+1)+b$는 어떤 정수 값으로 대체되어야 하는가? 기울기가 1보다 작으므로 $y$는 구간 $[y_k,y_k+1]$에서 값을 취한다. 즉,  오른쪽 옆이나 오른쪽 위 픽셀에 점을 찍어야 한다. 수식으로는 구간의 아래($d_{lower}$)와 위쪽($d_{upper}$)과의 차이를 비교하여 작은 쪽으로 결정하면 된다:

$$d_{lower}=y-y_k=m(x_k+1)+b-y_k$$

$$d_{upper}=(y_k+1)-y=y_k+1-m(x_k+1)-b$$

이 과정을 실수 연산 없이 판별하는 방법을 알아내자. $d_{lower}$와 $d_{upper}$의 차이는 

$$d_{lower}-d_{upper}=2m(x_k+1)-2y_k+2b-1$$

로 표현되고, 나눗셈을 하지 않기 위해서 $\mathrm{\Delta }x$을 곱한 값을 discriminator로 사용하자:

$$p_k=\mathrm{\Delta }x(d_{lower}-d_{upper})=2\mathrm{\Delta }y(x_k+1)-2\mathrm{\Delta }x{y}_k+\mathrm{\Delta }x(2b-1)$$

의 꼴로 표현되므로 직선의 정보($\mathrm{\Delta }x,\mathrm{\Delta }y$)와 현재 점을 찍은 픽셀($x_k,y_k$)을 알면 다음번에 점을 찍을 픽셀을 정수 연산만 가지고 판별할 수 있는 장점이 있다. $k+1$일 때 판별식이

$$p_{k+1}=2\mathrm{\Delta }y(x_k+1)-2\mathrm{\Delta }x{y}_{k+1}+\mathrm{\Delta }x(2b-1)$$

이므로 다음과 같은 recursion relation을 만족한다:

$$\begin{gathered} p_{k+1}=p_k+2\mathrm{\Delta }y-2\mathrm{\Delta }x(y_{k+1}-y_k) \\ {p}_0=2\mathrm{\Delta }y-\mathrm{\Delta }x \end{gathered}$$

이를 이용하면 시작 픽셀 $(x_0,y_0)$에서 $p_0=2\mathrm{\Delta }y-\mathrm{\Delta }x$을 계산하여 0보다 작지 않으면 $(x_k+1,y_k+1)$ 픽셀에 점을 찍고,  그다음 점찍을 픽셀은  $p_{k+1}=p_k+2\mathrm{\Delta }y-2\mathrm{\Delta }x$을 이용해서 판별한다. $p_0<0$ 이면 ($x_k+1,y_k$) 픽셀에 점을 찍고,  그다음 점을 찍을 픽셀은 $p_{k+1}=p_k+2\mathrm{\Delta }y$을 계산하여 판별한다.

void BresenhamLine(int x1, int y1, int x2, int y2, DWORD color) {
    int dx = x2 - x1 ;
    int incx = 1;
    if (dx < 0) {
        dx = -dx;
        incx = -1;
    }
    int dy = y2 - y1;
    int incy = 1;
    if (dy < 0) {
        dy = -dy;
        incy = -1;
    }
    setPixel(x1, y1,color);
    if (dx > dy) {
        int p = (dy << 1) - dx;
        while (x1 != x2) {
            if (p >= 0) {
                y1 += incy;
                p -= (dx << 1);
            }
            x1 += incx;
            p += (dy << 1);
            setPixel(x1, y1, color);
        }
    } else {
        int p = (dx << 1) - dy;
        while (y1 != y2) {
            if (p >= 0) {
                x1 += incx;
                p -= (dy << 1);
            }
            y1 += incy;
            p += (dx << 1);
            setPixel(x1, y1, color);
        }
    }
};

볼록 다각형(convex polygon)에서 한 방향으로(예를 들면 y-축방향) 가장 멀리 떨어진 두 꼭짓점(antipodal points) 통과하는 평행한 두 직선을 생각해보자. 평행 직선의 방향은 이 경우에는 x-축(y-축과 수직)과 나란해야 한다. 이제 두 평행 직선의 방향을 각각의 꼭짓점을 축으로 해서 조금씩 회전시키면, 둘 중 하나는 볼록 다각형에 먼저 접하게 된다. 이 접하는 변을 기준으로 다시 antipodal point을  찾아 (당연히 접하지 않는 평행직선이 통과하는 꼭짓점이 된다) 새로운 평행 직선을 구성한 후 동일한 과정을 반복한다. 이 과정을 N회 반복하면 볼록 다각형의 모든 변에 대해서 변을 공유하는 직선과 그 직선의 antipodal point를 찾을 수 있다. 

이 rotating calipers를 이용하면 convex hull의 bounding box, convex hull의 diameter, convex polygon 사이 거리 등의 계산에 응용할 수 있다.

double COS(CPoint A, CPoint B) {             // directional cosine;
    double norm2A = A.x * A.x + A.y * A.y;
    double norm2B = B.x * B.x + B.y * B.y;
    return (A.x * B.x + A.y * B.y) / sqrt(norm2A * norm2B);
}
int ccw(CPoint A, CPoint B, CPoint C) {     // cross(AB, BC); 반시계방향 > 0;
    return (B.x - A.x) * (C.y - B.y) - (B.y - A.y) * (C.x - B.x) ;
}
// find antipodal points along y-direction;
bool initCalipers(std::vector<CPoint> &Q, 
                 int &minid, double &minCos,
                 int &maxid, double &maxCos) {
    int N = Q.size();
    // ccw check (중복 // 일직선을 고려);
    int orient = 0;
    for (int i = 0; i < N; i++) {
        orient = ccw(Q[(i - 1 + N) % N], Q[i], Q[(i + 1) % N]);
        if (orient != 0) break;
    }
    if (orient == 0) return false;

    minid = maxid = 0;
    for (int i = 1; i < N; i++) {
        if (Q[i].y < Q[minid].y) minid = i;
        if (Q[i].y > Q[maxid].y) maxid = i;
    }

    CPoint n0 = CPoint(1, 0);
    if (orient < 0) n0 = -n0;  // convex hull이 CW인 경우;
    CPoint dir0 = Q[(minid + 1) % N] - Q[minid];
    CPoint dir1 = Q[(maxid + 1) % N] - Q[maxid];
    minCos = COS(dir0, n0);
    maxCos = COS(dir1, -n0);
    return true;
}
int rotatingCalipers(std::vector<CPoint> &Q,
                    int &i0, double &dirCos0,
                    int &i1, double &dirCos1) {
    int N = Q.size();
    int i0next = (i0 + 1) % N ;
    int i1next = (i1 + 1) % N ;
    CPoint dir0 = Q[i0next] - Q[i0];
    CPoint dir1 = Q[i1next] - Q[i1];
    if (dirCos0 > dirCos1) {          // i0을 포함하는 평행선이 다음 edge와 이루는 각이 작음;
        CPoint ndir0 = Q[(i0next + 1) % N] - Q[i0next]; // 각을 재는 기준선은 현재 다다음 edge;
        dirCos0 = COS(ndir0, dir0);    
        dirCos1 = COS(dir1, -dir0);    
        i0 = i0next;                  // i0는 다음 꼭지점으로 이동; i1 = antipodal point;
        return 0;
    } else {                          // i1을 포함하는 평행선이 다음 edge와 이루는 각이 작음
        CPoint ndir1 = Q[(i1next + 1) % N] - Q[i1next];
        dirCos1 = COS(ndir1, dir1);   
        dirCos0 = COS(dir0, -dir1);   
        i1 = i1next;                  // i1은 다음 꼭지점으로 이동; i0 = antipodal point;
        return 1;
    }
}
void runCalipers(std::vector<CPoint> &Q) {
    int i0, i1;
    double dirCos0, dirCos1;
    bool done = false;
    if (!initCalipers(Q, i0, dirCos0, i1, dirCos1)) return ;  // fails
    while (!done) {
    	int i0old = i0, i1old = i1;
        if (rotatingCalipers(Q, i0, dirCos0, i1, dirCos1) == 0) {
            // Q[i1] = antipodal point;
            // do something useful...
        } else {
            // Q[i0] = antipodal point;
            // do something useful...
        }
    }
};

2차원 점집합의 median

• www.cs.princeton.edu/~chazelle/pubs/ConvexLayers.pdf
• arxiv.org/ftp/arxiv/papers/1803/1803.04964.pdf