Geometry & Recognition

세 점 $A$, $B$, $C$가 만드는 삼각형의 외접원의 반지름을 구해보자. 세 변을 나타내는 벡터는

$$\overrightarrow{a}=\overrightarrow{OC}-\overrightarrow{OB},\overrightarrow{b}=\overrightarrow{OA}-\overrightarrow{OC},\overrightarrow{c}=\overrightarrow{OB}-\overrightarrow{OA}.$$

세 점이 일직선 위에 놓이지 않기 위해서는 

$$\overrightarrow{a}\times \overrightarrow{b}\ne 0$$

를 만족해야 한다. 

그림의 삼각형에서 $\mathrm{\angle }(BOC)$가 $\mathrm{\angle }(BAC)$의 두 배이므로,

$$\sin \theta =\frac{a/2}{R}.$$

이제 삼각형 면적은 세 변의 길이 ($a,b,c$)와 외접원의 반지름($R$)으로 표현된다:

$$\text{area}=\frac{1}{2}bc\sin \theta ⟶\frac{abc}{4R}.$$

또한, 삼각형의 면적을 세 점이 일직선 위에 있는가를 테스트하는 식으로 표현하면 

$$\text{area}=\frac{1}{2}|\overrightarrow{a}\times \overrightarrow{b}|$$

처럼 쓸 수 있다. 따라서 삼각형의 외접원의 반지름은 아래처럼 주어진다.

$$R=\frac{abc}{4\cdot \text{area}}=\frac{|\overrightarrow{a}||\overrightarrow{b}||\overrightarrow{c}|}{2|\overrightarrow{a}\times \overrightarrow{b}|}.$$

double circumradius(CPoint A, CPoint B, CPoint C) {
    double ax = C.x - B.x, ay = C.y - B.y;
    double bx = A.x - C.x, by = A.y - C.y;
    double crossab = ax * by - ay * bx;
    if (crossab != 0) { 
        double a = hypot(ax, ay);
        double b = hypot(bx, by); 
        double cx = B.x - A.x, cy = B.y - A.y;       
        double c = hypot(cx, cy);
        return (0.5 * a * b * c / fabs(crossab));
    } else 
        return 0;   //collinear
};

static 
int ccw(const CPoint& A, const CPoint& B, const CPoint& C) {
    // cross (C-A, B-A) >= 0 이면 AB의 오른쪽에 C가 놓임;
    int cax = C.x - A.x;
    int cay = C.y - A.y;
    int bax = B.x - A.x;
    int bay = B.y - A.y;
    return (cax * bay - cay * bax) >= 0; 
}

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

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

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

std::vector<CPoint> quickHull(const std::vector<CPoint>& pts) {
    if (pts.size() < 3) std::vector<CPoint> ();//null_vector;
    //Find left and right most points, say A & B, and add A & B to convex hull ;
    findMaxMin(pts);
    std::vector<CPoint> hull;
    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);
    return 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를 이용하거나 복사복은 이용하여야 한다.