This time, I had more luck in my attempt not to reinvent the wheel. I stumbled upon some old John W. Ratcliff's code which seemed quite promising. I tested it to find that the algorithm failed on the 'hard' cases that interested me, such as shown below.
I managed to modify the source code, so that it uses a more robust (at least in my case) indexed approach to describe the polygon. Allowing to disambiguate the vertices on the hole junctions segments. I post the modified source code hereafter, with all my thanks to John W. Ratcliff for posting his useful code snippets.
It seems that this code also does a good job at triangulating polygons (with holes this time) but is more heavyweight and less directly usable.
// COTD Entry submitted by John W. Ratcliff [jratcliff@verant.com]
// ** THIS IS A CODE SNIPPET WHICH WILL EFFICIEINTLY TRIANGULATE ANY
// ** POLYGON/CONTOUR (without holes) AS A STATIC CLASS. THIS SNIPPET
// ** IS COMPRISED OF 3 FILES, TRIANGULATE.H, THE HEADER FILE FOR THE
// ** TRIANGULATE BASE CLASS, TRIANGULATE.CPP, THE IMPLEMENTATION OF
// ** THE TRIANGULATE BASE CLASS, AND TEST.CPP, A SMALL TEST PROGRAM
// ** DEMONSTRATING THE USAGE OF THE TRIANGULATOR. THE TRIANGULATE
// ** BASE CLASS ALSO PROVIDES TWO USEFUL HELPER METHODS, ONE WHICH
// ** COMPUTES THE AREA OF A POLYGON, AND ANOTHER WHICH DOES AN EFFICENT
// ** POINT IN A TRIANGLE TEST.
// ** SUBMITTED BY JOHN W. RATCLIFF (jratcliff@verant.com) July 22, 2000
// ** SOME CHANGES BY ROTOGLUP ON March 28, 2008 TO ADDED INDEXED
// ** POLYGON SUPPORT, THUS ALLOWING TO HANDLE MORE COMPLEX POLYGONS
// ** (STILL WITHOUT REAL HOLES)
/**********************************************************************/
/************ HEADER FILE FOR TRIANGULATE.H ***************************/
/**********************************************************************/
#ifndef TRIANGULATE_H
#define TRIANGULATE_H
/*****************************************************************/
/** Static class to triangulate any contour/polygon efficiently **/
/** You should replace Vector2d with whatever your own Vector **/
/** class might be. Does not support polygons with holes. **/
/** Uses STL vectors to represent a dynamic array of vertices. **/
/** This code snippet was submitted to FlipCode.com by **/
/** John W. Ratcliff (jratcliff@verant.com) on July 22, 2000 **/
/** I did not write the original code/algorithm for this **/
/** this triangulator, in fact, I can't even remember where I **/
/** found it in the first place. However, I did rework it into **/
/** the following black-box static class so you can make easy **/
/** use of it in your own code. Simply replace Vector2d with **/
/** whatever your own Vector implementation might be. **/
/*****************************************************************/
#include <vector> // Include STL vector class.
class Vector2d
{
public:
Vector2d(float x,float y)
{
Set(x,y);
};
float GetX(void) const { return mX; };
float GetY(void) const { return mY; };
void Set(float x,float y)
{
mX = x;
mY = y;
};
private:
float mX;
float mY;
};
// Typedef an STL vector of vertices which are used to represent
// a polygon/contour and a series of triangles.
typedef std::vector< Vector2d > Vector2dVector;
typedef std::vector< int > IndexVector;
class Triangulate
{
public:
// triangulate a contour/polygon, places results in STL vector
// as series of triangles.
static bool Process(const Vector2dVector &contour, Vector2dVector& result);
static bool Process(const Vector2dVector &contour, const IndexVector& indices, IndexVector& result );
// compute area of a contour/polygon
static float Area(const Vector2dVector &contour, int const* V, int numPoints);
// decide if point Px/Py is inside triangle defined by
// (Ax,Ay) (Bx,By) (Cx,Cy)
static bool InsideTriangle(float Ax, float Ay,
float Bx, float By,
float Cx, float Cy,
float Px, float Py);
private:
static bool Snip(const Vector2dVector &contour,int u,int v,int w,int n,int const *V);
};
#endif
/**************************************************************************/
/*** END OF HEADER FILE TRIANGULATE.H BEGINNING OF CODE TRIANGULATE.CPP ***/
/**************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include "triangulate.h"
static const float EPSILON=0.0000000001f;
float Triangulate::Area(const Vector2dVector &contour, int const* V, int numPoints)
{
int n = numPoints;
float A=0.0f;
for(int _p=n-1,_q=0; _q<n; _p=_q++)
{
int p = V[_p];
int q = V[_q];
A += contour[p].GetX()*contour[q].GetY() - contour[q].GetX()*contour[p].GetY();
}
return A*0.5f;
}
/*
InsideTriangle decides if a point P is Inside of the triangle
defined by A, B, C.
*/
bool Triangulate::InsideTriangle(float Ax, float Ay,
float Bx, float By,
float Cx, float Cy,
float Px, float Py)
{
float ax, ay, bx, by, cx, cy, apx, apy, bpx, bpy, cpx, cpy;
float cCROSSap, bCROSScp, aCROSSbp;
ax = Cx - Bx; ay = Cy - By;
bx = Ax - Cx; by = Ay - Cy;
cx = Bx - Ax; cy = By - Ay;
apx= Px - Ax; apy= Py - Ay;
bpx= Px - Bx; bpy= Py - By;
cpx= Px - Cx; cpy= Py - Cy;
aCROSSbp = ax*bpy - ay*bpx;
cCROSSap = cx*apy - cy*apx;
bCROSScp = bx*cpy - by*cpx;
return ((aCROSSbp >= 0.0f) && (bCROSScp >= 0.0f) && (cCROSSap >= 0.0f));
};
bool Triangulate::Snip(const Vector2dVector &contour,int u,int v,int w,int n,int const *V)
{
int p;
float Ax, Ay, Bx, By, Cx, Cy, Px, Py;
int vu = V[u];
int vv = V[v];
int vw = V[w];
Ax = contour[vu].GetX();
Ay = contour[vu].GetY();
Bx = contour[vv].GetX();
By = contour[vv].GetY();
Cx = contour[vw].GetX();
Cy = contour[vw].GetY();
if ( EPSILON > (((Bx-Ax)*(Cy-Ay)) - ((By-Ay)*(Cx-Ax))) ) return false;
for (p=0;p<n;p++)
{
int vp = V[p];
if( (vp == vu) || (vp == vv) || (vp == vw) ) continue;
Px = contour[vp].GetX();
Py = contour[vp].GetY();
if (InsideTriangle(Ax,Ay,Bx,By,Cx,Cy,Px,Py)) return false;
}
return true;
}
bool Triangulate::Process(const Vector2dVector &contour, const IndexVector& indices, IndexVector& result )
{
/* allocate and initialize list of Vertices in polygon */
int const n = indices.size();
if ( n < 3 ) return false;
int* V = new int[n];
/* we want a counter-clockwise polygon in V */
if ( 0.0f < Area(contour, &indices[0], n) )
{
std::copy(indices.begin(), indices.end(), V);
}
else
{
std::copy(indices.rbegin(), indices.rend(), V);
}
int nv = n;
/* remove nv-2 Vertices, creating 1 triangle every time */
int count = 2*nv; /* error detection */
for(int m=0, v=nv-1; nv>2; )
{
/* if we loop, it is probably a non-simple polygon */
if (0 >= (count--))
{
//** Triangulate: ERROR - probable bad polygon!
return false;
}
/* three consecutive vertices in current polygon, <u,v,w> */
int u = v ; if (nv <= u) u = 0; /* previous */
v = u+1; if (nv <= v) v = 0; /* new v */
int w = v+1; if (nv <= w) w = 0; /* next */
if ( Snip(contour,u,v,w,nv,&V[0]) )
{
int a,b,c,s,t;
/* true names of the vertices */
a = V[u]; b = V[v]; c = V[w];
/* output counter-clockwise Triangle */
result.push_back( a );
result.push_back( b );
result.push_back( c );
m++;
/* remove v from remaining polygon */
for(s=v,t=v+1;t<nv;s++,t++) V[s] = V[t]; nv--;
/* resest error detection counter */
count = 2*nv;
}
}
delete V;
return true;
}
bool Triangulate::Process(const Vector2dVector &contour, Vector2dVector& result )
{
IndexVector indices;
IndexVector indices_result;
int const n = contour.size();
for (int i=0; i<n; i++) indices.push_back(i);
bool success = Process(contour, indices, indices_result);
if (success)
{
int const t = indices_result.size();
for (int i=0; i<t; i++)
{
result.push_back( contour[indices_result[i]] );
}
}
return success;
}
/************************************************************************/
/*** END OF CODE SECTION TRIANGULATE.CPP BEGINNING OF TEST.CPP A SMALL **/
/*** TEST APPLICATION TO DEMONSTRATE THE USAGE OF THE TRIANGULATOR **/
/************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include "triangulate.h"
void main(int argc,char **argv)
{
// Small test application demonstrating the usage of the triangulate
// class.
{
printf("** Raw positions usage\n");
// Create a pretty complicated little contour by pushing them onto
// an stl vector.
Vector2dVector a;
a.push_back( Vector2d(0,6));
a.push_back( Vector2d(0,0));
a.push_back( Vector2d(3,0));
a.push_back( Vector2d(4,1));
a.push_back( Vector2d(6,1));
a.push_back( Vector2d(8,0));
a.push_back( Vector2d(12,0));
a.push_back( Vector2d(13,2));
a.push_back( Vector2d(8,2));
a.push_back( Vector2d(8,4));
a.push_back( Vector2d(11,4));
a.push_back( Vector2d(11,6));
a.push_back( Vector2d(6,6));
a.push_back( Vector2d(4,3));
a.push_back( Vector2d(2,6));
// allocate an STL vector to hold the answer.
Vector2dVector result;
// Invoke the triangulator to triangulate this polygon.
Triangulate::Process(a, result);
// print out the results.
int tcount = result.size()/3;
for (int i=0; i<tcount; i++)
{
const Vector2d &p1 = result[i*3+0];
const Vector2d &p2 = result[i*3+1];
const Vector2d &p3 = result[i*3+2];
printf("Triangle %d => (%0.0f,%0.0f) (%0.0f,%0.0f) (%0.0f,%0.0f)\n",i+1,p1.GetX(),p1.GetY(),p2.GetX(),p2.GetY(),p3.GetX(),p3.GetY());
}
}
{
printf("\n** Indexed positions usage\n");
// Create a pretty complicated little contour by pushing them onto
// an stl vector, that would fail without indices to desambiguate vertices positions
Vector2dVector a;
/*0*/a.push_back( Vector2d(-1,-1));
/*1*/a.push_back( Vector2d(-1,1));
/*2*/a.push_back( Vector2d(1,1));
/*3*/a.push_back( Vector2d(1,-1));
/*4*/a.push_back( Vector2d(0,-1));
/*5*/a.push_back( Vector2d(0,0));
/*6*/a.push_back( Vector2d(-0.25f,0.25f));
/*7*/a.push_back( Vector2d(0,0.5f));
/*8*/a.push_back( Vector2d(0.25f,0.25f));
/*9*/a.push_back( Vector2d(0,0)); /* unused in indexed mode */
/*10*/a.push_back( Vector2d(0,-1)); /* unused in indexed mode */
IndexVector indices;
indices.push_back(0);
indices.push_back(1);
indices.push_back(2);
indices.push_back(3);
indices.push_back(4);
indices.push_back(5);
indices.push_back(6);
indices.push_back(7);
indices.push_back(8);
indices.push_back(5);
indices.push_back(4);
// allocate an STL vector to hold the answer.
Vector2dVector result;
// Invoke the triangulator to triangulate this polygon.
bool success = Triangulate::Process(a, result);
if (!success)
{
printf("1. Failed _without_ indices...\n");
}
else
{
printf("1. Succeeded _without_ indices...\n");
// print out the results.
int tcount = result.size()/3;
for (int i=0; i<tcount; i++)
{
const Vector2d &p1 = result[i*3+0];
const Vector2d &p2 = result[i*3+1];
const Vector2d &p3 = result[i*3+2];
printf("Triangle %d => (%0.0f,%0.0f) (%0.0f,%0.0f) (%0.0f,%0.0f)\n",i+1,p1.GetX(),p1.GetY(),p2.GetX(),p2.GetY(),p3.GetX(),p3.GetY());
}
}
IndexVector result_indices;
// Invoke the triangulator to triangulate this polygon.
bool success_with_indices = Triangulate::Process(a, indices, result_indices);
if (!success_with_indices)
{
printf("2. Failed _with_ indices...\n");
}
else
{
printf("2. Succeeded _with_ indices...\n");
// print out the results.
int tcount = result_indices.size()/3;
for (int i=0; i<tcount; i++)
{
const Vector2d &p1 = a[ result_indices[i*3+0] ];
const Vector2d &p2 = a[ result_indices[i*3+1] ];
const Vector2d &p3 = a[ result_indices[i*3+2] ];
printf("Triangle %d => (%0.0f,%0.0f) (%0.0f,%0.0f) (%0.0f,%0.0f)\n",i+1,p1.GetX(),p1.GetY(),p2.GetX(),p2.GetY(),p3.GetX(),p3.GetY());
}
}
}
}
Posts
No comments:
Post a Comment