Delaunay triangulation in 2D

들로네 삼각분할(delaunay triangulation). 삼각분할(triangulation) 기법 중 하나로 보로노이 다이어그램(voronoi diagram)과 듀얼이미지인 그래프를 그리게된다.

 

들로네 삼각분할은 다음 조건을 만족하는 삼각형들의 집합으로 설명할 수 있다.

* simplex로 이루어져있다

* 모든 simplex의 가장 작은 내각이 최대값을 가진다.

* simplex로 이루어지는 circumcircle (3d 라면 circumsphere, 4d 라면...)의 내부에

   simplex 포인트 외의 다른 포인트가 존재하지 않는다.

 

[그림출처: http://en.wikipedia.org/wiki/Delaunay_triangulation]

 

위 그림으로 쉽게 이해할 수 있을것이다.

 

2차원 들로네 삼각분할한 그래프에서 삼각형 외접원의 중심을 이으면 보로노이 다이어그램이 된다고한다. 그래서 듀얼이미지..

 

  

[그림출처(역시 위키..): http://en.wikipedia.org/wiki/Delaunay_triangulation]

 

들로네 삼각분할은 incremental 알고리즘과 divide and conquer 알고리즘으로 구분할 수 있는데, 난 incremental algorithm을 작성해보았다. 둘 다 최대수행비용은 비슷한듯?

 

OpenGL로 삼각형을 표시하게 했고, 핵심 클래스는 DelaunayTriangulation 을 보면 된다.

 

/* Write by alleysark */

#include <iostream>
#include <vector>
#include <list>
#include <GL/glut.h>
using namespace std;

#pragma comment(lib, "opengl32")
#pragma comment(lib, "glu32")
#pragma comment(lib, "glut32")

#define INIT_SUPERTRI_SIZE 1.0f

#define sqre(x) (x*x)

class Vertex{
public:
 float x, y;

 Vertex(){ x=y=0.0f; }
 Vertex(const Vertex& v){
  x=v.x; y=v.y;
 }
 Vertex(float _x, float _y){
  x=_x; y=_y;
 }

 bool operator==(const Vertex& v){
  return (x==v.x && y==v.y);

 }
};

class Triangle{
public:
 int vt[3];

 Triangle(){ memset(vt, 0, sizeof(vt)); }
 Triangle(const Triangle& t){ memcpy(vt, t.vt, sizeof(vt)); }
 Triangle(int v1, int v2, int v3){
  vt[0]=v1; vt[1]=v2; vt[2]=v3;
 }
};

class DelaunayTriangulation{
private:
 vector<Vertex> mVertices;
 vector<Triangle> mTriangles;

public:
 DelaunayTriangulation(){
  mVertices.push_back( Vertex(0, INIT_SUPERTRI_SIZE) );
  mVertices.push_back( Vertex(-INIT_SUPERTRI_SIZE, -INIT_SUPERTRI_SIZE) );
  mVertices.push_back( Vertex(INIT_SUPERTRI_SIZE, -INIT_SUPERTRI_SIZE) );

  mTriangles.push_back( Triangle(0,1,2) );
 }

 void PushVertex(float x, float y){
  Vertex vt(x,y);
  mVertices.push_back( vt );
  int iVtx = mVertices.size()-1;

  int iTri = SearchCoverTriangle(vt);
  while( iTri==-1 ){
   for(int i=0; i<3; i++){
    mVertices[i].x*2.0f;
    mVertices[i].y*2.0f;
   }
   iTri = SearchCoverTriangle(vt);
  }

  Triangle mainTri = mTriangles[iTri];
  mTriangles.erase(mTriangles.begin()+iTri);

  for(int i=0; i<3; i++){
   recTriangulation(mainTri.vt[i], mainTri.vt[(i+1)%3], iVtx);
  }

 }

 // print to consol
 void PrintTriangles(){
  int nTri = mTriangles.size();
  for(int i=0; i<nTri; i++){
   for(int j=0; j<3; j++){
    cout<<"("<<mVertices[mTriangles[i].vt[j]].x<<","<<mVertices[mTriangles[i].vt[j]].y<<")  ";
   }
   cout<<endl;
  }
 }

 void RenderVertices(){
  unsigned vtxCount = mVertices.size();

  glBegin(GL_POINTS);
  glColor3f(1, 0, 0);
  for(unsigned i=0; i<vtxCount; i++){
   glVertex2f(mVertices[i].x, mVertices[i].y);
  }
  glEnd();
 }

 void RenderTriangles(){
  unsigned triCount = mTriangles.size();

  glBegin(GL_TRIANGLES);
  for(unsigned i=0; i<triCount; i++){
   glColor3f((float)(i*2)/(float)triCount, 1.0f/(float)i, (float)i/(float)triCount);
   for(int k=0; k<3; k++)
    glVertex2f(mVertices[mTriangles[i].vt[k]].x, mVertices[mTriangles[i].vt[k]].y);
  }
  glEnd();
 }

private:
 void recTriangulation(int edgeVtx1, int edgeVtx2, int iVtx){
  int iAdjTri, adjEdgeOrder;
  iAdjTri = SearchAdjecentTriangle(edgeVtx2, edgeVtx1, &adjEdgeOrder);
  if( iAdjTri!=-1 ){
   if( InCircumcircle(mTriangles[iAdjTri], mVertices[iVtx]) ){
    Triangle adjTri = mTriangles[iAdjTri];
    mTriangles.erase(mTriangles.begin()+iAdjTri);

    recTriangulation(adjTri.vt[adjEdgeOrder], adjTri.vt[(2+adjEdgeOrder)%3], iVtx);
    recTriangulation(adjTri.vt[(1+adjEdgeOrder)%3], adjTri.vt[(2+adjEdgeOrder)%3], iVtx);
    return;
   }
  }

  if( TriangleAreaX2(mVertices[iVtx], mVertices[edgeVtx1], mVertices[edgeVtx2]) >= 0 )
   mTriangles.push_back( Triangle(iVtx, edgeVtx1, edgeVtx2) );
  else
   mTriangles.push_back( Triangle(iVtx, edgeVtx2, edgeVtx1) );
 }

 // > 0 : ccw, < 0 : cw,  == 0 : on same line
 float TriangleAreaX2(const Vertex& v1, const Vertex& v2, const Vertex& v3){
  return (v2.x - v1.x)*(v3.y - v1.y) - (v2.y - v1.y)*(v3.x - v1.x);
 }

 bool IsIn(const Triangle& tri, const Vertex& vtx){
  if( TriangleAreaX2(vtx, mVertices[tri.vt[0]], mVertices[tri.vt[1]]) < 0 )
   return false;
  if( TriangleAreaX2(vtx, mVertices[tri.vt[1]], mVertices[tri.vt[2]]) < 0 )
   return false;
  if( TriangleAreaX2(vtx, mVertices[tri.vt[2]], mVertices[tri.vt[0]]) < 0 )
   return false;

  return true;
 }

 int SearchCoverTriangle(const Vertex& v){
  for(int i=mTriangles.size()-1; i>=0; i--){
   if( IsIn( mTriangles[i], v) )
    return i;
  }
  return -1;
 }

 // return triangle index. pOutEdgeOrder is important.
 // triangle index set and edge order is; 0-1: 0, 1-2: 1, 2-0: 2
 int SearchAdjecentTriangle(int iv1, int iv2, int* pOutEdgeOrder){
  for(int i=mTriangles.size()-1; i>=0; i--){
   for(int j=0; j<3; j++){
    if( mTriangles[i].vt[j]==iv1 && mTriangles[i].vt[(j+1)%3]==iv2 ){
     if( pOutEdgeOrder ) *pOutEdgeOrder = j;
     return i;
    }
    if( mTriangles[i].vt[j]==iv2 && mTriangles[i].vt[(j+1)%3]==iv1 ){
     if( pOutEdgeOrder ) *pOutEdgeOrder = j;
     return i;
    }
   }
  }
  return -1;
 }

 bool InCircumcircle(const Triangle& tri, const Vertex& vtx){
  return (sqre(mVertices[tri.vt[0]].x) + sqre(mVertices[tri.vt[0]].y)) * TriangleAreaX2(mVertices[tri.vt[1]], mVertices[tri.vt[2]], vtx) -
   (sqre(mVertices[tri.vt[1]].x) + sqre(mVertices[tri.vt[1]].y)) * TriangleAreaX2(mVertices[tri.vt[0]], mVertices[tri.vt[2]], vtx) +
   (sqre(mVertices[tri.vt[2]].x) + sqre(mVertices[tri.vt[2]].y)) * TriangleAreaX2(mVertices[tri.vt[0]], mVertices[tri.vt[1]], vtx) -
   (sqre(vtx.x) + sqre(vtx.y)) * TriangleAreaX2(mVertices[tri.vt[0]], mVertices[tri.vt[1]], mVertices[tri.vt[2]]) > 0;
 }
};

// global variables
DelaunayTriangulation gDelaunayTri;

// gl methods

void RenderScene(){
 glClear(GL_COLOR_BUFFER_BIT);
 glMatrixMode(GL_MODELVIEW);
 glLoadIdentity();

 gDelaunayTri.RenderVertices();

 gDelaunayTri.RenderTriangles();
 
 glFlush();
}

void Resize(GLsizei width, GLsizei height){
 GLfloat nRange = 100;

 if( height == 0 )
  height = 1;

 glViewport(0, 0, width, height);

 glMatrixMode(GL_PROJECTION);
 glLoadIdentity();

 if( width <= height )
  glOrtho(-nRange, nRange, -nRange*height/width, nRange*height/width, 1, 500);
 else
  glOrtho(-nRange*width/height, nRange*width/height, -nRange, nRange, 1, 500);
}

void MouseHandler(int button, int state, int x, int y){
 if( button == GLUT_LEFT_BUTTON && state == GLUT_DOWN ){
  float fx = (float)(x - 300)/300.0f;
  float fy = (float)(300 - y)/300.0f;

  gDelaunayTri.PushVertex(fx,fy);

  gDelaunayTri.PrintTriangles();
  cout<<"---------------------------"<<endl;
 }

 glutPostRedisplay();
}

int main(){
 glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
 
 glutInitWindowSize(600, 600);
 glutCreateWindow("Delaunay Triangulation in 2D");
 
 glutMouseFunc(MouseHandler);
 glutDisplayFunc(RenderScene);

 // set render scene
 glClearColor(1,1,1,1);
 glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
 glutMainLoop();
 return 0;
}

 

super triangle의 크기를 재설정하는 부분이 조금 이상해서 처음 삼각형 밖에 점을 찍으면 프로그램이 오작동한다. 귀찮으므로 패스...

 

들로네 삼각분할을 재밌게 설명한 글이 있다. http://darkpgmr.tistory.com/96 

이게 어디에 쓰일 수 있는지 볼수있다.