Geometry & Recognition

/* A->B->C가 반시계방향으로 정렬되거나(< 0), 일직선상에 있으면(=) 참이다*/
static int ccw(POINT A, POINT B, POINT C) {
    int a = A.x - B.x, b = A.y - B.y, c = C.x - B.x, d = C.y - B.y;
    return (a * d - b * c) <= 0; /* true if points A,B,C counterclockwise */
}

static int cmpl(const void *a, const void *b) {
    POINT *A = (POINT* )a, *B = (POINT* )b;
    int v = A->x - B->x ;       //lower-hull은 x 증가: y감소 순으로 정렬;
    if (v > 0) return 1; 
    if (v < 0) return -1;
    v = B->y - A->y ;
    if (v > 0) return 1; 
    if (v < 0) return -1;
    return 0;
}
static int cmph(const void *a, const void *b) {
    return cmpl(b, a);      // upper-hull은 x감소, y-증가 순으로 정렬;
}

static int makeChain(POINT *V, int n, int (*cmp)(const void*, const void*)) {
    qsort(V, n, sizeof(POINT), cmp);
    // 꼭짓점(i(>j))가 현재까지 만들어진 hull(s까지)의 에지보다 더 오른쪽에 
    // 있으면(ccw(i,j,j-1)), 왼쪽에 위치할 때까지 기존 hull을 줄인다. 
    // 이 점을 hull에 추가.
    int s = 1;
    for (int i = 2; i < n; i++) {
        int j = s;
        while (j >= 1 && ccw(V[i], V[j], V[j - 1])) j--; // ccw의 =때문에 일직선상의 
                                                         // 마지막 점만 들어감;
        s = j + 1; // j = last index of hull;
        // (j + 1)에 새 hull-point을 추가;
        POINT t = V[i]; V[i] = V[s]; V[s] = t;         // swap: V[i] <-> V[j + 1]
    }
    return s;
}
int chainHull(std::vector<POINT>& pts, std::vector<POINT>& chull) {
    int n = pts.size();
    if (n < 3) return 0;
    // copy to working array;
    chull.resize(n + 1);
    for (int i = 0; i < n; i++) chull[i] = pts[i];
    /* make lower hull */
    int u = makeChain(&chull[0], n, cmpl);      
    /* make upper hull */   
    chull[n] = chull[0];
    u += makeChain(&chull[u], n - u + 1, cmph); 
    chull.resize(u);
    return u;
};

Input = a set S of n points 
    Assume that there are at least 2 points in the input set S of points
QuickHull (S) { 
    // Find convex hull from the set S of n points
    Convex Hull := {} 
    가장 왼쪽(A)과 오른쪽(B)에 있는 두 점을 찾아 convec hull에 추가;
    선분 AB을 기준으로 나머지 (n-2)개 점들을 두 그룹 S1과 S2로 나눔;
        S1 = 선분 AB의 오른쪽에 놓인 점; 
        S2 = 선분 BA의 오른쪽에 놓인 점;
    FindHull (S1, A, B) 
    FindHull (S2, B, A) 
};
FindHull (Sk, P, Q) { 
   If Sk has no point, 
        then  return. 
   선분 PQ에서 가장 멀리 떨어진 점(C)을 찾아 convex hull에서 P와 Q의 사이에 추가;
   세점 P, C, Q는 나머지 Sk의 나머지 점들을 세 그룹S0, S1, and S2 
        S0 = 삼각형 PCQ에 포함되는 점;
        S1 = 선분 PC의 오른쪽에 놓인 점;
        S2 = 선분 CQ의 오른쪽에 놓인 점; 
    FindHull(S1, P, C) 
    FindHull(S2, C, Q) 
}
Output = convex hull

참고: http://www.cse.yorku.ca/~aaw/Hang/quick_hull/Algorithm.html

// P가 line(A, B)의 완전히 왼쪽에 있는가?

static bool leftSide(CPoint P, CPoint A, CPoint B) ;

static bool leftSide(CPoint P, CPoint A, CPoint B) {
   int ax = P.x - A.x,  ay = P.y - A.y;
   int bx = B.x - A.x,  by = B.y - A.y;
   return (ax * by - ay * bx) >= 0;
};

double dist2ToLine(CPoint P, CPoint A, CPoint B);

double dist2ToLine(CPoint P, CPoint A, CPoint B) {
    double dx = (B.x - A.x), dy = (B.y - A.y);
    double den = (dx * dx + dy * dy);
    if (den == 0.) return 0; //degenerate;
    double du = (P.x - A.x), dv = (P.y - A.y);
    double dist = (du * dy - dv * dx);
    return dist * dist / den; //distance-square!
};

#define SWAP_POINT(A, B)

static void findHull(CPoint *S, int n, CPoint A, CPoint B, std::vector <CPoint>& hull) ;

static void findHull(CPoint *S, int n, CPoint A, CPoint B, std::vector<CPoint>& hull) {
    if (n < 1) return ;
    //Find furtherst point 
    double maxdist = 0;
    for (int i = 0; i < n; i++) {
        double dist = dist2ToLine(S[i], A, B);
        if (dist > maxdist) {
            SWAP_POINT(S[0], S[i]);
            maxdist = dist;
        }
    }
    if (maxdist == 0.) return ; //collinear case;
    hull.push_back(S[0]);
    int j = 1;
    for (int i = 1; i < n; i++) {
        if (leftSide(S[i], A, S[0])) {
            SWAP_POINT(S[i], S[j]);
            j++;
        };
    };
    int k = j;
    for (int i = j; i < n; i++) {
        if (leftSide(S[i], S[0], B)) {
            SWAP_POINT(S[i], S[k]);
            k++;
        }
    };
    //1,...,j-1;
    findHull(&S[1], j-1, A, S[0], hull);
    //j,...,k-1;
    findHull(&S[j], k-j, S[0], B, hull);
};

void findMaxMin(std::vector<CPoint>& pts) ;

void findMaxMin(std::vector<CPoint>& pts) {
    int minx = pts[0].x, maxx = minx;
    int minid = 0, maxid = 0;
    for (int i = pts.size(); i-->1;) {
        if (pts[i].x > maxx) maxid = i, maxx = pts[i].x;
        if (pts[i].x < minx) minid = i, minx = pts[i].x;
    }
    SWAP_POINT(pts[0], pts[minid]); if (maxid == 0) maxid = minid;
    SWAP_POINT(pts[1], pts[maxid]);
};

void hullToPolyline(std::vector<CPoint>& hull);

void quickHull(std::vector<CPoint>& pts, std::vector<CPoint>& hull) {
    if (pts.size() < 3) return;
    hull.clear();
    //Find left and right most points, say A & B, and add A & B to convex hull ;
    findMaxMin(pts);
    hull.push_back(pts[0]);
    hull.push_back(pts[1]);
    int j = 2;
    for (int i = 2; i < pts.size(); i++) {
        if (leftSide(pts[i], pts[0], pts[1])) {
            SWAP_POINT(pts[i], pts[j]);
            j++;
        }
    }
    //2,3,...,j-1;
    findHull(&pts[2], j-2, pts[0], pts[1], hull);
    //j,j+1,...,pts.size()-1;
    if (j < pts.size()) // in order to avoid dereferencing &pts[pts.size()];
        findHull(&pts[j], pts.size()-j, pts[1], pts[0], hull);
    hullToPolyline(hull);
};

출력 hull은 단순히 convex hull을 구성하는 정렬이 안된 점들만 주므로, hull의 에지를 구하고 싶으면 추가적인 수정이 필요함.

static int cmph(const void *a, const void *b) {
    CPoint *A = (CPoint *)a , *B = (CPoint *)b ;
    int v = (A->x - B->x) ;
    if (v > 0 ) return 1;
    if (v < 0 ) return -1;
    v = B->y - A->y ;
    if (v > 0 ) return 1;
    if (v < 0 ) return -1;
    return 0;
};
static int cmpl(const void * a, const void *b) {
    return cmph(b, a);
};
void hullToPolyline(std::vector<CPoint>& hull) {
    CPoint A = hull[0];
    CPoint B = hull[1];
    // 가장 멀리 떨어진 두 점(hull[0], hull[1])을 연결하는 직선을 기준으로 프로젝션을 구하여서 순차적으로 
    int k = 2;
    for (int i = 2; i < hull.size(); i++) {
        if (leftSide(hull[i], A, B)) {
            SWAP_POINT(hull[i], hull[k]); k++;
        };
    };
    // k-1; last index of hull left side of line(A,B);
    // upper part reordering: 
    qsort(&hull[0], k, sizeof(CPoint), cmph);
    //lower part reordering;
    if (k < hull.size())
        qsort(&hull[k], hull.size()-k, sizeof(CPoint), cmpl);
    }
};

또한 입력점들의 순서를 그대로 유지하고 싶으면, double pointer를 이용하거나 복사복은 이용하여야 한다.

3차의 Hermite spline을 쓰는 interpolation에서 입력 데이터 점에서의 기울기 값을 수정하면 입력 데이터의 monotonic 한 성질을 그대로 유지하는 보간을 할 수 있다. 주어진 데이터가 $\{(x_1, y_1), (x_2, y_2), \dots, (x_n, y_n)\}$이고, $x/y$-좌표가 증가하는 순으로 정렬되었다고 하자. 구간 $[x_k, x_{k+1} ]$의 시작과 끝에서 값 ($p_0 = y_k $, $p_1=y_{k+1}$)과 기울기 ($m_0, m_1$) 가 주어지면 Cubic Hermite spline 보간 함수는 다음과 같이 매개변수 $t \in [0,1]$의 함수로 표현된다: 

$$f(t) = p_0 h_{00}(t) + m_1 h_{10}(t) + p_1 h_{01}(t) + m_1 h_{11}(t),            \quad \quad t= \frac {x-x_k}{x_{k+1}-x_k} $$

$$f(0) = p_0,\quad f'(0) = m_0,\quad f(1) = p_1,\quad f'(1)=m_1$$

Hermite 기저 함수의 정의는 아래와 같다:

$$ h_{00}(t) = 2t^3 - 3t^2+1, \quad h_{10}(t) = t^3 - 2t^2+1,\quad h_{01}(t)=-2t^3+3t^2,\quad h_{11}(t)=t^3-t^2$$

보간된 spline이 단조성(monotonicity)을 유지하도록 각 입력점에서 기울기 값($m_{k}$)을 조절해주어야 한다. 그럼 기울기 $m_{k}$는 어떻게 설정해야 인접 구간의 보간 함수와 smooth 하게 연결될 수 있을까? 입력점의 각 위치에서 기울기는 각각 왼쪽과 오른쪽 인접 입력점과의 평균 변화율($\delta_{k-1}, \delta_{k}$)을 구한 후 다시 이들의 평균값을 기울기로 사용하면 된다. (양 끝점에서는 오른쪽이나 왼쪽 구간의 평균 변화율로 대체)

$$ m_k = \frac {1}{2}(\delta_{k-1}+\delta_{k} )= \frac {1}{2}\left(\frac { y_k-y_{k-1}}{x_k - x_{k-1}}+\frac {y_{k+1}-y_k}{x_{k+1} + x_k}\right)$$

이것만으로는 단조성을 보증을 할 수 없으므로 다음 조건을 추가로 부여해야 됨이 관련 논문에서 증명이 되었다. 

$$1. \text{For}~~\delta_{k} = 0 \quad\Longrightarrow\quad m_k =m_{k+1} = 0$$

$$2. \text{For}~~\epsilon=\sqrt {\left(\frac {m_k }{\delta_k }\right)^2 + \left(\frac {m_{k+1}}{\delta_{k} }\right)^2 } > 3\\ \quad\Longrightarrow\quad m_{k} \leftarrow3\frac {m_{k}}{\epsilon}, ~m_{k+1} \leftarrow 3\frac {m_{k+1}}{\epsilon}$$

아래의 코드는 위키피디아에 제시된 알고리즘을 구현한 것이다

(입력 데이터가 한 방향으로 monotonic 할 때만 동작함)

참고: http://en.wikipedia.org/wiki/Monotone_cubic_interpolation

void estimateTangents(std::vector <double>& x, std::vector <double>& y, std::vector<double>& m);

더보기

static void estimateTangents(std::vector<double>& x, std::vector<double>& y, 
                             std::vector<double>& m) {
    int n = x.size() ;
    std::vector<double> delta (n) ;
    // slopes;
    for (int k = 0; k < n - 1; k++)
        delta[k] = (y[k + 1] - y[k]) / (x[k + 1] - x[k]) ;
        
    // average tangent at data points;
    m[0]   = delta[0];
    m[n-1] = delta[n-2];
    for (int k = 1; k < n - 1; k++)
        m[k] = (delta[k-1] + delta[k]) / 2.;  //error corrected.
        
    for (int k = 0; k < n - 1; k++) {
        if (delta[k] == 0.) m[k] = m[k+1] = 0;
        else {
            double alpha = m[k] / delta[k];
            double beta  = m[k + 1] / delta[k];
            double epsilon = sqrt(alpha * alpha + beta * beta);
            if (epsilon > 3) {
                double tau = 3 * delta[k] / epsilon ;
                m[k]     = tau * alpha;
                m[k + 1] = tau * beta;
            }
        }
    }
};

void monotoneCubicSplineInterp(std::vector<CPoint>& points, std::vector<CPoint>& curve) {
    curve.clear();
    if (points.size() < 2) return ;
    //if (!isMonotonicPolyline(points)) return ;
    std::vector<double> xx(points.size()), yy(points.size());
    for (int i = points.size(); i-->0;) {
        xx[i] = points[i].x ;
        yy[i] = points[i].y;
    };
    //tangent table;
    std::vector<double> m(points.size(), 0);
    estimateTangents(xx, yy, m);
    //hermite spline;
#define STEPS 12
    curve.push_back(points[0]); //add first point;
    for (int i = 0; i < points.size() - 1; i++) {
        double h = xx[i + 1] - xx[i];
        for (int k = 1; k <= STEPS; k++) {
            double t = double(k) / STEPS;         // 0 < t <= 1;
            double t2 = t * t ;
            double t3 = t2 * t ;
            double h00 =  2. * t3 - 3. * t2     + 1.;
            double h10 =       t3 - 2. * t2 + t;
            double h01 = -2. * t3 + 3. * t2;
            double h11 =       t3      - t2;
            double yf = h00 * yy[i] + h10 * h * m[i] + h01 * yy[i + 1] + h11 * h * m[i + 1];
            curve.push_back(CPoint(int(xx[i] + t * h + 0.5), int(yf + 0.5)));
        }
    }
};

Natural cubic spline과의 비교:

Convex hull은 집합으로 주어진 점이나 영역을 포함하는 가장 작은 볼록 집합이다. 알려진 알고리즘으로는 

1. Jarvis March (Gift Wrap)

2. Graham scan

3. Quick hull (Chain hull)

4. Monotone chain

5.......

여기서는 소개하는 convex hull 알고리즘은 가장 기초적인 알고리즘이다. 두 점에 의해서 만들어지는 선분이 convex hull의 일부이면 나머지 점들은 항상 왼편(설정에 따라 오른편으로 해도 된다)에 있게 된다는 사실을 이용한다. N개의 점이 있을 떄 가능한 선분의 수는 N(N-1) (방향까지 고려함)이고, 나머지 점에 대해서 테스트를 해야 하므로 총 비교 횟수는 N(N-1)(N-2) ~ O(N^3) 정도이다. 알고리즘의 구현이 간단하고, 단순 비교연산이므로 매우 robust해서 점집합의 크기가 작을 때 효과적이다.

static int ccw(CPoint A, CPoint B, CPoint P) ; // cross(AB, AP) > 0 if <ABP> is ccw;

static int ccw(CPoint A, CPoint B, CPoint P) { //cross(AB, AP) > 0 if <ABP> is ccw;
    return ((B.x - A.x) * (P.y - A.y) - (B.y - A.y) * (P.x - A.x));
}

int BruteForceConvexHull(std::vector<CPoint>& pts, std::vector<CPoint>& hull) {
    // O(n^3) ;
    if (pts.size() < 3) return 0;
    hull.clear() ;    
    std::vector<int> head; //hull을 구성하는 에지의 head index;
    std::vector<int> tail; //hull을 구성하는 에지의 tail index;
    for (int i = 0; i < pts.size(); i++) {
        for (int j = 0; j < pts.size(); j++) {
            if (pts[j] == pts[i]) continue;     //에지는 중복점이 아닌 경우만(점을 비교해야 함)
            bool onHull = true;
            for (int k = 0; k < pts.size(); k++) {
                if (pts[k] == pts[i] || pts[k] == pts[j]) continue;
                if (ccw(pts[i], pts[j], pts[k]) < 0) {  // cw;
                    onHull = false;                     // (i,j) is not an edge;
                    break ;
                }
            }
            if (onHull) { //(i->j) edge is on the hull;
                tail.push_back(i); head.push_back(j);
            }
        }
    }
    if (head.size() < 1) return 0;
    // hull을 구성하는 에지정보에서 단순폴리곤 만들기(bug 수정);
    int currIdx = 0;
    hull.push_back(pts[tail[currIdx]]);
    hull.push_back(pts[head[currIdx]]);
    while (head[currIdx] != tail[0]) {
        for (int j = 0; j < head.size(); j++)
            if (tail[j] == head[currIdx]) {
                currIdx = j;
                hull.push_back(pts[head[currIdx]]);
                break;
            }
    }
    return hull.size();
};

convex hull의 edge가 중간에 (collinear) 점을 포함하지 않도록 하기 위해서는, ccw == 0일 때 ABC가 접혀있는지를 체크하면 된다.

주변에 보이는 많은 선 중에는 직선도 있고 휘어진 곡선도 있다. 그럼 수학적으로 곡선이 휘어진 정도를 어떻게 정의할까? 곡선의 각각의 부분에서 휘어진 정도가 다  다르므로, 휘어진 정도는 곡선의 위치에 따라 달라지는 값이 될 것이다. 평면에 놓인 곡선의 기술할 때 보통은 $x, y$ 좌표를 많이 사용하는데, $x$ 방향으로 움직이면서 $y$의 값이 변하는 정도를 재는 것으로 곡선의 휘어짐을 정의할 수 있지 않을까? 직선을 살펴보면 $x$ 축으로 움직이면 $y$의 값이 변하는데, 통상 직선은 휘어져 있다고 이야기하지 않으므로, 이러한 방법은 휘어짐을 정의하기에는 적당하지 않다. 곡선 위의 임의 지점에서 곡선의 휘어짐은 그 지점에서 접선 방향으로 움직였을 때, 접점에서 많이 벗어날 수로 접선과의 차이가 크면 곡선은 많이 휘어졌다는 것은 직관적으로 알 수 있다. 이것을 좀 더 수학적으로 이야기하면 곡선을 따라가면 접선의 방향이 얼마나 변했는가를 재면 곡선이 휘어진 정도를 알 수 있다는 것을 의미한다. 그런데 곡선을 따라갈 때 빨리 갈 수도 또는 느리게 갈 수도 있으므로, 곡선을 따라가는 비율을 고정해야 한다. 가장 간단한 방법은 곡선을 따라 움직이는 거리를 기준으로 접선의 변화가 얼마나 생기는 가를 재면 된다. 곡선의 길이($s$)를 매개변수로 사용하면 곡선 위의 임의 지점에서 접선 벡터는 크기는 항상 1로 주어진다:

$$ds^2=dx^2 + dy^2 \quad \rightarrow\quad \sqrt{ \Big(\frac{dx}{ds}\Big)^2 + \Big(\frac{dy}{ds}\Big)^2 } = 1.$$

곡선의 휘어짐은 곡률이라는 용어를 사용하는데, 엄밀하게 정의하면 곡선의 한 지점에서 곡률은 그 지점에서 접선 벡터

$$\vec{T} = \left(\frac{dx}{ds},\frac{dy}{ds} \right)$$ 

의 미분계수의 크기로 주어진다:

$$\kappa = \Big| \frac{d\vec{T}}{ds}\Big| = \sqrt{\Big(\frac{d^2x}{ds^2} \Big)^2 + \Big(\frac{d^2 y }{ds^2} \Big)^2 }.$$

접선 벡터 $\vec{T}$는 길이가 1이므로 미분한 값은 $\vec{T}$에 수직하고, 곡선이 안쪽으로 휘어지는 쪽을 가리킨다.  $\vec{T}$에 수직인 방향의 단위 벡터를 $\vec{N}$이라 하면(보통 $\vec{T}$에서 반시계 방향으로 회전된 방향, $\vec{T}\cdot \vec{N}=0$) $d\vec{T}/ds$는 $\vec{N}$에 비례하고 비례항의 크기가 그 지점에서 곡률이다. 

$$\frac{d\vec{T}}{ds} =k \vec{N},\quad \kappa = |k|. $$

즉, 곡률이 클수록 같은 길이를 옮겨갈 때 $T$의 방향 변화가 심하게 일어난다. 주어진 $T$방향에 대해서 곡선이 휘어지는 방향이 오른쪽이나 왼쪽이 될 수 있으므로 $k$는 부호를 가지게 된다. 반지름이 $R$인 원이 있다고 하자(중심=원점, $(R,0)$에서 거리를 재면)

,$$x(s)=R \cos(\theta)=R\cos (s/R), \quad y(s) = R \sin (\theta)=R \sin (s/R).$$

$$ \therefore \kappa =\frac{1}{R}$$

이므로 곡률은 반지름의 역수로 주어진다. 반지름이 작을수록 같은 거리를 움직일 때 접선의 방향 변화가 심하므로 곡률이 더 크게 나타날 것이라는 것은 쉽게 예상할 수 있다. 직선의 경우에는 접선 벡터가 일정하므로 곡률은 당연히 0이다.

일반적인 매개변수($t$)인 경우 

$$\vec{T} = \left(\frac{dx}{ds},\frac{dy}{ds} \right) = \frac{1}{\dot{s}}   ( \dot{x}, \dot{y} ) = \frac{1}{\sqrt{\dot{x}^2 + \dot{y}^2 }} (\dot{x}, \dot{y})$$ 이므로 (overdot은 매개변수 $t$-에 대한 미분)

$$\kappa = \frac{|\dot{x}\ddot{y} - \dot{y} \ddot{x}|}{( \dot{x}^2 + \dot{y}^2 )^{3/2}}.$$

하나의 예로 catenary의 경우를 계산해 보자. 매개변수로 표현된 catenary 곡선은 

$$ x(t) = a \sinh^{-1}(t),\quad y(t) = a\sqrt{1+t^2}\quad  (y=a\cosh(x/a))$$

이고, $a$는 장력의 수평성분을 줄의 선밀도와 중력가속도의 곱으로 나눈 양이다. $a$가 증가할수록 줄이 평평해지려는 경향이 있으므로 곡률이 전반적으로 줄어들 것으로 예상할 수 있다. 실제로 곡률을 계산하면 

$$ \kappa = \frac{1}{a(1+t^2)}=\frac{a}{y^2}$$

임을 알 수 있다.