주어진 점집합에서 세 점을 뽑아서 만든 삼각형의 외접원에 다른 점이 하나도 포함하지 않으면 triangulation의 기본 삼각형 cell이 된다. 주어진 점으로 만들 수 있는 삼각형의 총개수가 ${}_nC_3$ 이므로, 기본 삼각형을 찾기 위해서는 이들 각각의 삼각형에 대서 나머지 점을 가지고 incircle 테스트를 수행해야 한다. 따라서 이 알고리즘은 ${\cal O}(n^4)$의 스텝이 필요하게 된다.
/*brute force attack*/
foreach point p1
foreach point p2
foreach point p3
foreach point p4
if(incircle(p1,p2,p3,p4))
iscell=false;
break ;
endfor;
if(iscell)
add(triangle(p1,p2,p3));
endfor;
endfor;
endfor;
세 점에 의해서 형성이 되는 외접원은 대수적으로 쉽게 구할 수 있다. 여기서는 좀 더 기하학적인 접근을 쓰면, 평면의 점은
$$(x, y)\rightarrow (x, y, z=x^2 + y^2)$$
의 mapping에 의해서 3차원 paraboloid 곡면의 점으로 옮길 수 있다. paraboloid 위의 세 점이 형성하는 3차원에서 평면이 paraboloid를 절단하는 곡선을 $x-y$ 평면으로 정사영하면 원이 된다는 것을 쉽게 알 수 있다.(incircle 포스팅 참조). 따라서 주어진 점이 세 점의 외접원에 포함되는가를 테스트하는 작업은 이 점을 paraboloid로 올렸을 때의 점과 (paraboloid로 올려진) 외접원의 3점이 형성하는 3차에서의 평면과 관계를 조사하는 것으로 바꿀 수 있다.
주어진 점이 외접원에 포함되면 paraboloid로 변환된 점은 평면의 아래에 놓이고, 외접원 밖의 점이면 평면 위에 놓이게 된다. 물론 외접원 위의 점은 평면에 놓인다. 따라서 평면의 법선 벡터를 구하고, 삼각형의 한 꼭짓점을 기준한 주어진 점의 변위 벡터와 내적을 구하면 내적의 값은 평면 위인지, 아래인지 또는 평면에 놓인 점인가에 따라서 부호가 달라진다. 평면의 수직 벡터를 고정하면(예제는 아래 방향: $n_z < 0$), 평면 위에 놓인 점과의 내적은 음수, 평면 아래에 놓인 점과의 내적은 양수가 되고, 평면의 점과의 내적은 0이다.
주어진 세 점이 만드는 외접원 내부(and 경계)에 들어가는 점이 없으면 이 삼각형을 선택한다.
** 참고 : Computational Geometry in C(2nd Edition) by Joseph O'Rourke
std::vector<Triple> dt4(const std::vector<double>& x, const std::vector<double>& y) {
const int n = x.size();
if (n < 3) return std::vector<Triple> (); // null_vec;
std::vector<double> z(n);
for (int i = 0; i < n; i++)
z[i] = x[i] * x[i] + y[i] * y[i] ;
std::vector<Triple> triples;
/* For each triple (i,j,k) */
for (int i = 0; i < n - 2; i++ )
for (int j = i + 1; j < n; j++ )
for (int k = i + 1; k < n; k++ )
if ( j != k ) {
/* Compute normal to triangle (i,j,k):: outter_product(j-i, k-i)*/
double nx = (y[j] - y[i]) * (z[k] - z[i]) - (y[k] - y[i]) * (z[j] - z[i]);
double ny = (x[k] - x[i]) * (z[j] - z[i]) - (x[j] - x[i]) * (z[k] - z[i]);
double nz = (x[j] - x[i]) * (y[k] - y[i]) - (x[k] - x[i]) * (y[j] - y[i]);
/* Only examine faces on bottom of paraboloid: nz < 0. */
int flag = (nz < 0);
if (flag) {
/* For each other point m */
for (int m = 0; m < n; m++) {
/* Check if m above (i,j,k)::i점을 기준으로 m 과
** normal 간의 내적으로 체크(내적이 양수이면
** m이 원의 내부에 완전히 들어 있는 경우가 된다.
** 0인 경우는 원주상에 놓인 경우이므로 배제한다
*/
flag &= ((x[m]-x[i])*nx + (y[m]-y[i])*ny + (z[m]-z[i])*nz <= 0);
}
}
if (flag) {
// (i, j, k)의 외접원이 다른 점을 포함하지 않으므로 이 삼각형은
// 삼각분할의 한 면을 형성하게 된다.
triples.push_back(Triple(i, j, k));
}
}
return triples;
}
단순 다각형(simple polygon)의 삼각화 알고리즘의 다른 구현이다. 이전 구현과 같은 Ear-Clipping 알고리즘을 써서 구현한 것이다. 폴리곤의 각 꼭짓점을 돌면서 현재 꼭짓점의 직전과 직후의 꼭짓점들로 구성한 삼각형이 Ear인가를 체크해서 맞으면 이 삼각형을 추가하고, 현 꼭짓점은 주어진 폴리곤에서 제거를 한 후 나머지에 대해 테스트를 계속한다. $N$ 개의 꼭짓점을 가지는 단순 폴리곤의 경우에 $N-2$ 개의 삼각형으로 분해가 된다.
주어진 폴리곤이 시계방향으로 정렬된 경우는 반시계 방향으로 바꾸어서 정렬시킨 후 작업을 한다. 구현된 알고리즘의 복잡도는 ${\cal O}(N^2)$이다.
last update: 2012.10.03;
static
bool leftSide(const CPoint &P, const CPoint &A, const CPoint &B) {
// P가 선분(A,B)의 왼쪽 또는 위에 있는가?
return ((B.x-A.x)*(P.y-A.y) - (B.y-A.y)*(P.x-A.x)) >= 0;
}
static
bool insideTriangle(const CPoint &P, const CPoint &A, const CPoint &B, const CPoint &C) {
// 반시계방향으로 정렬된 삼각형(A,B,C)의 안에 점 P가 있는가?(경계포함);
if (!leftSide(P, A, B)) return false;
if (!leftSide(P, B, C)) return false;
if (!leftSide(P, C, A)) return false;
return true;
};
static
int polygonArea2(const std::vector<CPoint> &pts) {
double area2 = 0;
for (int p = pts.size() - 1, q = 0; q < pts.size(); p = q++)
area2 += pts[p].x * pts[q].y - pts[p].y * pts[q].x;
return area2;
}
static
bool earTest(int a, int b, int c,
int nv, const std::vector<int> &V,
const std::vector<CPoint> &pts)
{
// Check the triangle {V[a], V[b], V[c]} is an Ear;
const CPoint &A = pts[V[a]];
const CPoint &B = pts[V[b]];
const CPoint &C = pts[V[c]];
// colinear(A->B->C) or concave인가를 체크; 반시계 방향으로 정렬된 입력점이므로 왼편에
// 있으면 concave이므로 ear가 안된다;
if (leftSide(B, A, C)) return false;
// 이미 만들어진 삼각형의 꼭지점이 포함되는지 체크한다;
for (int k = 0; k < nv; ++k) {
if ((k == a) || (k == b) || (k == c)) continue;
if (insideTriangle(pts[V[k]], A, B, C)) return false;
}
return true;
};
/* Polygon should be simple(ccw or cw);*/
std::vector<Triple> polygonTriangulation(const std::vector<CPoint>& pts) {
if (pts.size() < 3) return std::vector<Triple> (); //null_vec;
std::vector<int> V(pts.size()); // contains vertex indices;
// 주어진 단순폴리곤을 ccw-ordering;
if (polygonArea2(pts) > 0)
for (int v = pts.size(); v-- > 0;) V[v] = v;
else
for (int v = pts.size(), k = 0; v-- > 0;) V[v] = k++;
std::vector<Triple> triples;
triples.reserve(pts.size() - 2);
// (pts.size()-2)개 꼭지점을 차례로 제거한다.
// 꼭지점이 제거될때마다 삼각형이 하나씩 생긴다;
int nv = pts.size();
int err_detect = 2 * nv; /* error detection */
for (int v = nv - 1; nv > 2; ) {
// 단순폴리곤이 아닌 경우에 무한루프를 도는 것을 방지;
if (0 >= err_detect--) {
TRACE("ERROR::non-simple polygon\n");
return std::vector<Triple> ();
}
// <u,v,w>는 연속하는 세 꼭지점의 인덱스임;
int u = v % nv;
v = (u + 1) % nv;
int w = (v + 1) % nv;
if (earTest(u, v, w, nv, V, pts)) {
triples.push_back(Triple(V[u], V[v], V[w]));
// 꼭지점 V[v]를 제거함;
for (int k = v, j = k + 1; j < nv;) V[k++] = V[j++];
/* resest error detection counter */
err_detect = 2 * (--nv);
}
}
return triples;
}
Polygon의 내부를 겹치지 않게 삼각형으로 분할하는 것을 polygon의 삼각화라고 한다. $N$개의 꼭짓점이 있는 simple polygon의 경우 처음 세 꼭짓점이 만드는 삼각형이 생성된 후 꼭짓점을 하나씩 추가할 때마다 삼각형도 하나씩 늘어난다. 따라서 simple polygon 내부에는 $N-2$개의 서로 겹치지 않은 삼각형으로 분할이 되고 polygon의 경계와 겹치지 않는 $N-3$개의 내부 경계 라인을(diagonal)을 가지게 된다.
분할 삼각형의 외접원 반지름에 차이가 많이 나는 경우 삼각형들의 면적에 편차가 심하게 된다. 좀 더 균일한 면적을 가지는 삼각형으로 분할하기 위해서는 인접하는 두 삼각형으로 만들어지는 사변형에서 현재 대각 변에 교차하는 대각 변으로 재분할을 시도해 본다. 새로이 만들어지는 두 삼각형을 검사하여 이전보다 면적이 균일하면 이 분할을 사용하면 된다. 물론 삼각형의 외접원의 반경이 작을 때는 거의 균일한 분할을 얻을 수 있다. Polygon Optimizing을 참고하기 바란다. *simple polygon이 아닌 경우는 아래의 링크를 참조하기 바란다. 1). O(n log(n)) 복잡도를 갖는 polygon triangulation algorithm; ==> 'Fast Polygon Triangulation based on Seidel's Algorithm'. ==>http://www.cs.unc.edu/~dm/CODE/GEM/chapter.html 2). Poly2Tri 홈페이지: ==> Fast and Robust Simple Polygon Triangulation With/Without Holes by Sweep Line Algorithm ==>http://www.mema.ucl.ac.be/~wu/Poly2Tri/poly2tri.html
평면에서 보로노이 다이어그램을 계산하는 알고리즘 중에 1980년대에 Fortune에 의해서 발견이 된 Sweep line 알고리즘에 대해서 살펴보겠다. 이 알고리즘은 보로노이 다이어그램 계산문제의 optimal 한 답 중에 하나이다(~O(n * log(n)). Fortune이 자신의 알고리즘을 직접 C로 구현한 것이 있는데 이것은 웹상에서 쉽게 구할 수 있다. 그러나 알고리즘의 얼개는 간단하나 구현에 들어간 자료구조 등을 간단히 않아서 코드를 이해하는 것이 쉽지 않았다(이에 비하면 incremental delaunay triangulation은 매우 단순하다). 몇 번을 간단한 자료구조를 이용해서 재 작성을 해보려고 했는데 잘 되지 않았다. 인터넷에서 구할 수 있는 다양한 구현 소스를 분석해 보았으나 쉬운 것을 찾기가 힘들었다(triangle.c에서도 구현되어 있고, 자바로도 구현이 된 것을 찾을 수 있다.) 원래의 목적은 Fortune의 코드가 메모리 처리를 잘못하여 메모리 leak이 발생하여서 이 문제를 해결하려고 시도한 것인데, 메모리 문제는 자바의 경우에 처럼 가비지 콜렉터를 만들어서 처리할 수는 있다. 이 알고리즘의 구현에 들어가는 기본적인 자료구조는 priority queue와 balanced binary tree이다. 이들의 기본적인 구현은 이미 잘 알려져 있으므로 이것들을 직접적으로 알고리즘에 적용하는 문제만 남는다. balanced binary tree로 만들어진 자료구조를 쓰면 자료를 찾는 시간이 O(log(n))의 시간 복잡도를 주지만, 이것이 알고리즘을 직관적으로 보는 것을 방해하므로 여기서는 이중 연결 리스트를 이용하도록 한다. 전체적으로 알고리즘은 O(n^2)의 복잡도를 가진다.
1. priority queue 구성: 주어진 입력점들을 가지고 구성한다. 알고리즘 중간에 세 점으로 원이 구성이 되는 circle event가 생성이 되는데 이것은 따로 event 큐를 구성한다. event의 우선순위는 x 좌표의 순서에 의해서 구성하고 동일한 x 좌표값에 대해서는 y값이 작은 순서대로 구성한다(compare 구조체에서 point와 event에 대한 비교 연산자를 정의한다).
2. site event 처리: sweep line(x=const)은 전체 평면을 반으로 나누는 역할을 한다. 알고리즘은 sweep 라인이 쓸고 지나간 영역에서만 관심 영역이다. sweep라인과 주어진 점은 하나의 포물선을 정의할 수 있다(주어짐 점이 포물선의 초점을 구성함). 따라서 sweep 라인이 진행하면 sweep line의 왼편의 입력점들은 각각 하나의 포물선을 형성하고, 이 포물선들의 구간 중에서 sweep라인에 가장 가까운 부분은 포물선 arc로 연결인 되는 하나의 beach line을 형성한다. sweep라인이 새로운 입력점을 지나면 여기서 새로운 포물선 arc가 생기는데, 이것은 이미 만들어진 beach라인을 구성하는 포물선 arc 중 하나를 둘로 가르고, 그 자신이 새로운 beach 라인의 일부를 구성하게 된다; sweep라인이 전진하면서 비치라인을 구성하는 arc가 다른 arc에 파 묻혀서 없어지는 경우가 생긴다. 이 경우가 circle event가 만들어지는 시점이다.
.......................
* BEFORE AFTER * * new point new arc * | | * __ | _____ __ | _____ * \ \|/ | \ \|/ arc a * \ ` | \ ` __|__ * \ X | p------X __|__<-- arc c * a | arc a a | | * / | / arc a * / | / | * /__ __|__ /__ __|__ * __/ \ | __/ \ | * \ | \ | * \ | \ | * b | arc b b | arc b * / | / | * / | / | * __/ __|__ __/ __|__
3. circle event 처리; sweep라인이 각각의 site에 도달하면 새로운 포물선의 arc가 비치라인에 추가가 되는데 sweep 라인으로부터 멀어진 arc들은 어느 시점에서 없어져야 한다. 이것은 인접하는 세 개의 arc들의 교점이 현재의 sweepline 위치에서 하나로 만들어져서 중앙의 arc가 없어지는 경우로 이것을 circle event라고 한다. 이때 만들어지는 교점은 보로노이 에지의 꼭짓점을 구성하게 된다. circle event에서는 나머지 두 개의 남은 arc의 교점이 trace 하는 새로운 에지가 생기게 된다.
* __ * \ * a \ * \ * __ | * ` / * \ / * b X <-- Arc b is overtaken at point X (this is a circle center) * / \ * __, \ * \ * \ * c | * / * / * /
* * BEFORE AFTER * * \ arc a \ * \ \ a.s0 * \ <-- a.s0 \ | * \ \ \|/ * \ \ ` * \ \ * arc b X <-- termination point .--------------------<-- new segment * / / * / / . * / <-- c.s1 / /|\ * / / | * / arc c / c.s1 * /
//포물선 arc를 정의하기 위한 클래스;
struct arc {
point p; //focus of parabola;
arc *prev, *next; //double linked-list;
event *e;
////////////////////////////////////////////////////
seg *s0 ; //edge of voronoi starting from break point;
seg *s1; //edge of voronoi starting from break point;
arc(point _p, arc *_a=0, arc *_b=0)
: p(_p), prev(_a), next(_b), e(0), s0(0), s1(0) {}
};
// voronoi edge를 정의하는 segment class;
struct seg {
point start, end; //defines segment;
bool done;
int ref ; //referece count in order to distingush the double edges;
seg(point _p, int _ref=0)
: start(_p), end(0,0), done(false), ref(_ref)
{ output.push_back(this); } //garbage collector;
void finish(point p) { if (done) return; end = p; done = true; }
};
// 메인 wrapper;
int fortune_main(std::vector<POINT>& pts, //voronoi points;
std::vector<std::pair<CPoint,CPoint>>& edges) //voronoi edges;
{
if (pts.empty()) return 0;
// point priority_queue구성; // degeneray를 막기 위해 작은 random #을 추가하는 것이 좋다;
for (int i = 0; i < pts.size(); i++)
points.push(point(pts[i].x. pts[i].y));
// body of algorithm ;
while (!points.empty()) {
if (!events.empty() && events.top()->x <= points.top().x)
process_event() ; //circle event ;
else
process_site() ; //site event;
}
//for remaining circle events if any;
while (!events.empty())
process_event() ;
//2. prepare output edges and clean memory;
}
arc * root=NULL; //global variable;
void process_site() {
point P= points.top(); points.pop(); //because popping return void in STL;
if (!root) { // root of double linked list for arcs(global);
root = new arc(P) ; return ;
}
//
for (arc *a = root; a; a = a->next) {
point Q;
if (intersect(P, a, &Q)) {//arc a와 점P에서의 포물선(degenerate 된)이 만나는 점 Q;
duplicate_arc(P, a) ; //교차하는 arc를 둘로 만든다
//(다음 arc가 있고, 만약에 이것과도 교차하면 circle event이으로 제외)
insert_new_arc(P, a, a->next); //새 arc를 중간에 삽입한다;
// 새 arc의 에지 세팅 ::교차점을 출발점으로 하는 두개의 반직선을 만든다:
// 도착점은 circle event에 대부분 결정되고, 나머지는 후처리 과정에서 결정)
a = a->next ;
a->prev->s1 = a->s0 = new seg(Q, ref_count) ;
a->next->s0 = a->s1 = new seg(Q, ref_count++) ; //쌍으로 생성되는 반직선을 identify하기 위해서 동일번호를 부여함.
//새 site 추가로 비치라인을 구성하는 arc의 초점들과 circle event를 만들 수 있는가 체크;
check_circle_event(a, P.x) ;
check_circle_event(a->prev, P.x) ;
check_circle_event(a->next, P.x) ;
}
} // for-;
};
//
void process_event(){
event *e = events.top(); events.pop();
if (e->valid) {
// start a new edge.(single edges)
seg *s = new seg(e->p, ref_count++);
// remove the associated arc from the front. and attach a new segment;
arc *a = e->a;
remove_arc(a, s);
// finish the edges before and after a==>new voronoi vertex;
if (a->s0) a->s0->finish(e->p);
if (a->s1) a->s1->finish(e->p);
// recheck circle events on either side of p:
if (a->prev) check_circle_event(a->prev, e->x);
if (a->next) check_circle_event(a->next, e->x);
delete a;
}
delete e;
};
// Look for a new circle event for arc a.
void check_circle_event(arc *a, double x0/*=current sweep-line*/) {
// Invalidate any old event.
if (a->e && a->e->x != x0) a->e->valid = false;
a->e = NULL;
if (!a->prev || !a->next) return;
point &A = a->prev->p ;
point &B = a->p ;
point &C = a->next->p ;
double maxx;
point center;
//collinear가 아니고, 최대값이 현재의 site event위치보다도 큰 경우에;
//점 A,B,C에 의해서 형성이 되는 원 중심(center)과 최우측x(maxx) 좌표;
if (circle(A, B, C, &maxx, ¢er) && maxx > x0) {
// create new event.
a->e = new event(maxx, center, a);
events.push(a->e);
}
}
나머지 함수들은 모두 단순한 구현이므로 생략한다;
보로노이 에지는 seg collector에서 각각의 segment를 끄집어내어서 그리면 된다. 그러나 각각의 segment는 보로노이 에지를 전체를 커버하는 것이 아니다. site event인 경우에는 항상 듀얼로 반직선을 생성하는데. 이 두 개의 반직선이 하나의 보로노이 에지를 정의한다(따라서 ref를 참조하면 온전한 하나의 에지를 찾을 수 있다). circle event의 경우에는 에지가 듀얼로 생성하지 않았으므로 이 경우에는 하나의 segment가 그 에지를 표현한다. 에지 액세스를 쉽게 하기 위해서는 모든 에지를 듀얼 구조로 만들어서 사용할 수 있다.
n 개의 점들의 보로노이 다이어그램은 얼마나 많은 꼭짓점과 에지를 가질까?
이것은 오일러 공식 V-E+F=2를 사용하면 된다. 보로노이 다이어그램은 경우 바깥쪽의 에지들은 무한대로 연결이 되어 있다. 이것을 무한대에서 가상의 꼭짓점을 가정하고 그것에 연결이 되어 있다고 생각하면 된다. 따라서 V --> V+1로 계산을 해야 한다. 오일러 공식 : V+1-E + F= 2; 여기서 F=n임을 알 수 있다. (n개의 점들이 하나의 face상에 놓여 있음)그리고 하나의 에지가 두 개의 꼭짓점을 연결하므로 각각의 꼭짓점에서 나간 에지의 합(=deg(v))을 계산하면 2*E 개 임을 알 수 있다. sum deg(v) = 2 * E;그런데 각각의 꼭짓점에서는 적어도 3개 이상의 에지와 연결이 되어 있으므로 위식의 좌변은 (V+1) * 3 <= 2 * E;를 준다. 따라서 오일러 공식과 이 부등식을 연관시키면 V <= 2*n - 5; E <= 3*n - 6;임을 알 수 있다. 즉 필요한 메모리의 양은 입력점의 수의 선형적으로 비례한다.
따라서 전체 이벤트의 개수는 site 이벤트와 꼭짓점에 해당하는 circle event 만큼이 있으므로 O(3*n) 정도이다. 각각의 event에 대해서 arc노드를 검색해야 하므로 O(n) 번 탐색을 해야 한다(balanced binary tree의 경우에는 log(n)). 따라서 알고리즘의 복잡도는 O(n^2 ( n*log(n))이다./**** http://blog.naver.com/helloktk/80041603288*/
** 첨부된 실행파일로 알고리즘을 테스트해 볼 수 있다(중복된 입력점은 제거하였고, 세 점이 일직선상에 놓인 것을 방지하기 위해서 작은 랜덤 값을 첨가하였다)
