Geometry & Recognition

영상에서 추출한 객체(object)의 경계선은 객체의 모양에 대한 많은 정보를 담고 있지만, 어떤 경우에는 불필요하게 많은 정보가 될 수도 있다. 이 경우 원래의 경계를 충분히 닮은 다각형으로 간략화하여서 불필요한 정보를 제거할 수 있다. 이 다각형으로의  근사가  물체의 경계를 얼마나 제대로 표현할 수 있는지에 대한 기준이 있어야 한다. $N$개의 점 $\{(x_i, y_i) | i = 1,... ,N\}$으로 표현된 경계가 있다고 하자. 가장 단순한 근사는 가장 멀리 떨어진 두 점($A, B$)으로 만들어진 선분일 것이다(선분의 길이가 대략적인 물체의 크기를 주고, 선분 방향이 물체의 orientation을 알려준다). 이 선분을 기준으로 다시 가장 멀리 떨어진 점($C$)을 찾아서 일정 거리 이상 떨어져 있으면 꼭짓점으로 추가하여 두 개의 선분으로 분할한다. 두 개의 선분에 대해서 동일한 분할 작업을 재귀적으로 반복하므로 좀 더 정밀하게 물체를 표현하는 다각형을 찾을 수 있다. 선분에서 가장 멀리 떨어진 지점이 기준 거리 이내이면 분할을 멈춘다. 닫힌 다각형 근사일 때는 시작점 = 끝점으로 조건을 부여하면 된다.

** open polyline으로 우리나라 경계를 단순화시키는 예: top-right 쪽이 시작과 끝이므로 simplication 과정에서 변화지 않는다.

// distance square from P to the line segment AB ;

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

// P에서 두 점 A, B을 연결하는 선분까지 거리
double pt2SegDist2(CPoint P, CPoint A, CPoint B) {
    double dx = B.x - A.x, dy = B.y - A.y ;
    double lenAB2 = dx * dx + dy * dy ;
    double du = P.x - A.x, dv = P.y - A.y ;
    if (lenAB2 == 0.) return du * du + dv * dv;
    double dot = dx * du + dy * dv ; 
    if (dot <= 0.) return du * du + dv * dv;
    else if (dot >= lenAB2) {
        du = P.x - B.x; dv = P.y - B.y;
        return du * du + dv * dv;
    } else {
        double slash = du * dy - dv * dx ;
        return slash * slash / lenAB2;
    }
};

//  recursive Douglas-Peucker 알고리즘;

void DouglasPeucker(double tol, std::vector<CPoint>& vtx, int start, int end, 
                   std::vector<bool>& mark) {
    if (end <= start + 1) return; //종료 조건;
    // vtx[start] to vtx[end]을 잇는 선분을 기준으로 분할점을 찾음;
    int break_id = start;			// 선분에서 가장 먼 꼭짓점;
    double maxdist2 = 0;
    const double tolsq = tol * tol;		
    for (int i = start + 1; i < end; i++) {
        double dist2 = pt2SegDist2(vtx[i], vtx[start], vtx[end]);
        if (dist2 <=  maxdist2) continue;
        break_id = i;
        maxdist2 = dist2;
    } 
    if (maxdist2 > tolsq) {		// 가장 먼 꼭짓점까지 거리가 임계값을 넘으면 => 분할;
        mark[break] = true;		// vtx[break_id]를 유지;
        // vtx[break_id] 기준으로 좌/우 polyline을 분할시도;
        DouglasPeucker(tol, vtx, start, break_id, mark);
        DouglasPeucker(tol, vtx, break_id, end, mark);
    }
};

// driver routine;

int polySimplify(double tol, std::vector <CPoint>& V, std::vector <CPoint>& decimated) {
    if (V.size() < 2) return 0;    
    std::vector<bool> mark(V.size(), false); // 꼭짓점 마킹용 버퍼;
    // 알고리즘 호출 전에 너무 붙어있는 꼭짓점을 제거하면 보다 빠르게 동작;
    // 폐곡선이 아닌 경우 처음,끝은 marking; 
    // 폐곡선은 가장 멀리 떨어진 두점중 하나만 해도 충분;
    mark.front() = mark.back() = true; 
    DouglasPeucker( tol, V, 0, V.size()-1, mark );
    // save marked vertices;
    decimated.clear();    
    for (int i = 0; i < V.size(); i++)
        if (mark[i]) decimated.push_back(V[i]);
    return decimated.size();
}

nonrecursive version:

// nonrecursive version;
int DouglasPeucker_nonrec(double tol, std::vector<CPoint>& points,
                            std::vector<CPoint>& decimated) {
    if (points.size() <= 2) return;
    double tol2 = tol * tol; //squaring tolerance;
    // 미리 가까이 충분히 가까이 있는 꼭지점은 제거하면 더 빨라짐;
    std::vector<CPoint> vt;
    vt.push_back(points.front()) ;  // start at the beginning;
    for (int i = 1, pv = 0; i < points.size(); i++) {
        int dx = points[i].x - points[pv].x, dy = points[i].y - points[pv].y;
        double dist2 = dx * dx + dy * dy ;
        if (dist2 < tol2) continue;
        vt.push_back(points[i]); 
        pv = i; // 다시 거리를 재는 기준점;
    }
    if (vt.back()!=points.back()) vt.push_back(points.back());  // finish at the end
    points.swap(vt);
    //
    std::vector<bool> mark(points.size(), false); // vertices marking
    std::vector<int> stack(points.size());  // 분할된 라인 처음-끝 저장;
    // 폐곡선이 아닌 경우 처음,끝은 marking; 
    // 폐곡선은 가장 멀리 떨어진 두점 중 하나만 해도 충분
    mark.front() = mark.back() = true;
    int top = -1;
    // 분할 해야할 폴리라인의 양끝점을 스택에 넣음;
    stack[++top] = 0;
    stack[++top] = points.size() - 1;
    while (top >= 0) {
        // 분할을 시도할 구간의 두 끝점을 꺼냄;
        int end = stack[top--];
        int start = stack[top--];
        // 최대로 멀리 떨어진 점;
        double max_distance2 = 0;
        int break_id = -1;
        for (int i = start + 1; i < end; i++) {
            double dsq = pt2SegDist2(points[i], points[start], points[end]);
            if (dsq > max_distance2) {
                max_distance2 = dsq;
                break_id = i;
            }
        }
        // 주어진 임계값 이상 떨어졌으면 꼭지점을 유지하도록 표시;
        if (max_distance2 > tol2) {
            mark[break_id] = true;       
            // 주어진 꼭지점을 기준으로 양쪽을 다시 검사하기 위해 스택에 쌓음;
            stack[++top] = start;
            stack[++top] = break_id;
            stack[++top] = break_id;
            stack[++top] = end;
        }
    }
    // marking 된 꼭지점을 출력;
    decimated.clear();
    for (int i = 0; i < points.size(); i++)
        if (mark[i]) decimated.push_back(points[i]);
    return decimated.size();
}

내 점으로 만들어진 convex인 사변형에 어떤 점이 포함되어 있는지 않은지를 판별하는 문제는 프로그램을 하다 보면 종종 만나게 된다. 이보다 더 근본적인 문제는 주어진 점이 세 점으로 만들어진 삼각형의 내부점인가를 판단하는 것이다. 세 점이 주어지면, 시계방향이던, 아니면 반시계 방향이던, 한 방향으로 oriented 된 삼각형이 구성이 된다.  삼각형의 내부점들은 이 삼각형은 이루는 변(orientation 때문에 이 변들은 방향을 가지는 벡터임)에 대해서 모두 한쪽으로 치우쳐 상태가 된다. 이것을 수학적으로 표현하면, 삼각형의 세 변 벡터와 각각의 변 벡터의 시작점에서 주어진 점까지 선분이 만드는 벡터와의 외적의 부호가 동일한 값을 갖는다(삼각형의 CW나 CCW인 것에 상관없다)

$$\vec {A} =(A_x, A_y), ~\vec {B}=(B_x, B_y) ~\Longrightarrow ~\vec {A}\times\vec {B}=A_x B_y - A_y B_x$$

// 평면에서 두 벡터의 외적은 숫자를 준다; 3차원으로 생각할 때 z-성분이다. 
BOOL insideTriangle(POINT p, POINT A, POINT POINT B, POINT C) {
    ab = vec(A,B)와 vec(A,P)의 외적; 
    bc = vec(B,C)와 vec(B,P)의 외적; 
    ca = vec(C,A)와 vec(C,P)의 외적; 
    // 점이 변 위에 놓인 경우 외적이 0이 되므로 sign함수를 쓰는 경우 주의! 
    // return SIGN(ab) == SIGN(bc) && SIGN(bc) == SIGN(ca); 
    return (ab * bc >= 0) && (bc * ca >=0) && (ca * ab >= 0) 
}; 

(convex) 사변형은 항상 두 개의 삼각형으로 나누어지므로, 각각의 삼각형에 대해서 테스트하면 된다(이 방식이 직접적으로 4 변에 대해서 모두 다 같은 쪽에 치우쳐 있는가를 테스트하는 것보다 효과적이다. 삼각형을 이용하면 4 변에 대해서 다 테스트하지 않아도 중간에 결과를 얻을 수 있다. 물론, 주어진 사변형이 처음에 어떤 방향으로 orient 되었는지의 정보가 주어졌다면 훨씬 더 빨리 판별할 수 있는 방법이 있다).

// 점이 변 위에 있는 경우는 포함한다; 
// 사변형의 orientation에 무관하다(단, convex이어야 한다); 
BOOL insideQuadrilateral(POINT p, POINT A, POINT B, POINT C, POINT D) { 
    if (insideTriangle(p, A, B, D)) return TRUE; 
    else if (insideTriangle(p, B, C, D)) return TRUE; 
    return FALSE; 
}