// 头文件#pragma once#include <set> #include <algorithm> #include <math.h>using namespace std;#ifndef _GDIPLUS_H// I designed this with GDI+ in mind. However, this particular code doesn't // use GDI+ at all, only some of it's variable types. // These definitions are substitutes for those of GDI+. typedef float REAL; struct PointF { PointF() : X(0), Y(0) {propZ=0;} PointF(const PointF& p) : X(p.X), Y(p.Y) {propZ=p.propZ;} PointF(REAL x, REAL y, REAL z) : X(x), Y(y) {propZ=z;} PointF operator+(const PointF& p) const { return PointF(X + p.X, Y + p.Y, propZ + p.propZ); } PointF operator-(const PointF& p) const { return PointF(X - p.X, Y - p.Y, propZ - p.propZ); } REAL X; REAL Y; REAL propZ; // 用于等值线的点值,add by zhaoyang. };const REAL REAL_EPSILON = 1.192092896e-07F; // = 2^-23; I've no idea why this is a good value, but GDI+ has it.#endif // _GDIPLUS_H/////////////////// // vertexclass vertex { public: vertex() : m_Pnt(0.0F, 0.0F, 0.0F) {} vertex(const vertex& v) : m_Pnt(v.m_Pnt) {} vertex(const PointF& pnt) : m_Pnt(pnt) {} vertex(REAL x, REAL y, REAL z) : m_Pnt(x, y, z) {} vertex(int x, int y, int z) : m_Pnt((REAL) x, (REAL) y, (REAL) z) {} bool operator<(const vertex& v) const { if (m_Pnt.X == v.m_Pnt.X) return m_Pnt.Y < v.m_Pnt.Y; return m_Pnt.X < v.m_Pnt.X; } bool operator==(const vertex& v) const { return m_Pnt.X == v.m_Pnt.X && m_Pnt.Y == v.m_Pnt.Y; }
REAL GetX() const { return m_Pnt.X; } REAL GetY() const { return m_Pnt.Y; } REAL GetZ() const { return m_Pnt.propZ;}; void SetX(REAL x) { m_Pnt.X = x; } void SetY(REAL y) { m_Pnt.Y = y; } void SetZ(REAL z) { m_Pnt.propZ = z; } const PointF& GetPoint() const { return m_Pnt; } protected: PointF m_Pnt; };typedef set<vertex> vertexSet; typedef set<vertex>::iterator vIterator; typedef set<vertex>::const_iterator cvIterator;/////////////////// // triangleclass triangle { public: triangle(const triangle& tri) : m_Center(tri.m_Center) , m_R(tri.m_R) , m_R2(tri.m_R2) { for (int i = 0; i < 3; i++) m_Vertices[i] = tri.m_Vertices[i]; } triangle(const vertex * p0, const vertex * p1, const vertex * p2) { m_Vertices[0] = p0; m_Vertices[1] = p1; m_Vertices[2] = p2; SetCircumCircle(); } triangle(const vertex * pV) { for (int i = 0; i < 3; i++) m_Vertices[i] = pV++; SetCircumCircle(); } bool operator<(const triangle& tri) const { if (m_Center.X == tri.m_Center.X) return m_Center.Y < tri.m_Center.Y; return m_Center.X < tri.m_Center.X; } const vertex * GetVertex(int i) const { ASSERT(i >= 0); ASSERT(i < 3); return m_Vertices[i]; } bool IsLeftOf(cvIterator itVertex) const { // returns true if * itVertex is to the right of the triangle's circumcircle return itVertex->GetPoint().X > (m_Center.X + m_R); } bool CCEncompasses(cvIterator itVertex) const { // Returns true if * itVertex is in the triangle's circumcircle. // A vertex exactly on the circle is also considered to be in the circle. // These next few lines look like optimisation, however in practice they are not. // They even seem to slow down the algorithm somewhat. // Therefore, I've commented them out. // First check boundary box. // REAL x = itVertex->GetPoint().X; // // if (x > (m_Center.X + m_R)) return false; // if (x < (m_Center.X - m_R)) return false; // // REAL y = itVertex->GetPoint().Y; // // if (y > (m_Center.Y + m_R)) return false; // if (y < (m_Center.Y - m_R)) return false; PointF dist = itVertex->GetPoint() - m_Center; // the distance between v and the circle center REAL dist2 = dist.X * dist.X + dist.Y * dist.Y; // squared return dist2 <= m_R2; // compare with squared radius } protected: const vertex * m_Vertices[3]; // the three triangle vertices PointF m_Center; // center of circumcircle REAL m_R; // radius of circumcircle REAL m_R2; // radius of circumcircle, squared void SetCircumCircle(); };// Changed in verion 1.1: collect triangles in a multiset. // In version 1.0, I used a set, preventing the creation of multiple // triangles with identical center points. Therefore, more than three // co-circular vertices yielded incorrect results. Thanks to Roger Labbe. typedef multiset<triangle> triangleSet; typedef multiset<triangle>::iterator tIterator; typedef multiset<triangle>::const_iterator ctIterator;/////////////////// // edgeclass edge { public: edge(const edge& e) : m_pV0(e.m_pV0), m_pV1(e.m_pV1) {} edge(const vertex * pV0, const vertex * pV1) : m_pV0(pV0), m_pV1(pV1) { } bool operator<(const edge& e) const { if (m_pV0 == e.m_pV0) return * m_pV1 < * e.m_pV1; return * m_pV0 < * e.m_pV0; } const vertex * m_pV0; const vertex * m_pV1; };typedef set<edge> edgeSet; typedef set<edge>::iterator edgeIterator; typedef set<edge>::const_iterator cedgeIterator;/////////////////// // Delaunayclass Delaunay { public: // Calculate the Delaunay triangulation for the given set of vertices. void Triangulate(const vertexSet& vertices, triangleSet& output); // Put the edges of the triangles in an edgeSet, eliminating double edges. // This comes in useful for drawing the triangulation. void TrianglesToEdges(const triangleSet& triangles, edgeSet& edges); protected: void HandleEdge(const vertex * p0, const vertex * p1, edgeSet& edges); };
//CPP #include "StdAfx.h" #include "Delaunay.h"const REAL sqrt3 = 1.732050808F;void triangle::SetCircumCircle() { REAL x0 = m_Vertices[0]->GetX(); REAL y0 = m_Vertices[0]->GetY(); REAL x1 = m_Vertices[1]->GetX(); REAL y1 = m_Vertices[1]->GetY(); REAL x2 = m_Vertices[2]->GetX(); REAL y2 = m_Vertices[2]->GetY(); REAL y10 = y1 - y0; REAL y21 = y2 - y1; bool b21zero = y21 > -REAL_EPSILON && y21 < REAL_EPSILON; if (y10 > -REAL_EPSILON && y10 < REAL_EPSILON) { if (b21zero) // All three vertices are on one horizontal line. { if (x1 > x0) { if (x2 > x1) x1 = x2; } else { if (x2 < x0) x0 = x2; } m_Center.X = (x0 + x1) * .5F; m_Center.Y = y0; } else // m_Vertices[0] and m_Vertices[1] are on one horizontal line. { REAL m1 = - (x2 - x1) / y21; REAL mx1 = (x1 + x2) * .5F; REAL my1 = (y1 + y2) * .5F; m_Center.X = (x0 + x1) * .5F; m_Center.Y = m1 * (m_Center.X - mx1) + my1; } } else if (b21zero) // m_Vertices[1] and m_Vertices[2] are on one horizontal line. { REAL m0 = - (x1 - x0) / y10; REAL mx0 = (x0 + x1) * .5F; REAL my0 = (y0 + y1) * .5F; m_Center.X = (x1 + x2) * .5F; m_Center.Y = m0 * (m_Center.X - mx0) + my0; } else // 'Common' cases, no multiple vertices are on one horizontal line. { REAL m0 = - (x1 - x0) / y10; REAL m1 = - (x2 - x1) / y21; REAL mx0 = (x0 + x1) * .5F; REAL my0 = (y0 + y1) * .5F; REAL mx1 = (x1 + x2) * .5F; REAL my1 = (y1 + y2) * .5F; m_Center.X = (m0 * mx0 - m1 * mx1 + my1 - my0) / (m0 - m1); m_Center.Y = m0 * (m_Center.X - mx0) + my0; } REAL dx = x0 - m_Center.X; REAL dy = y0 - m_Center.Y; m_R2 = dx * dx + dy * dy; // the radius of the circumcircle, squared m_R = (REAL) sqrt(m_R2); // the proper radius // Version 1.1: make m_R2 slightly higher to ensure that all edges // of co-circular vertices will be caught. // Note that this is a compromise. In fact, the algorithm isn't really // suited for very many co-circular vertices. m_R2 *= 1.000001f; }// Function object to check whether a triangle has one of the vertices in SuperTriangle. // operator() returns true if it does. class triangleHasVertex { public: triangleHasVertex(const vertex SuperTriangle[3]) : m_pSuperTriangle(SuperTriangle) {} bool operator()(const triangle& tri) const { for (int i = 0; i < 3; i++) { const vertex * p = tri.GetVertex(i); if (p >= m_pSuperTriangle && p < (m_pSuperTriangle + 3)) return true; } return false; } protected: const vertex * m_pSuperTriangle; };// Function object to check whether a triangle is 'completed', i.e. doesn't need to be checked // again in the algorithm, i.e. it won't be changed anymore. // Therefore it can be removed from the workset. // A triangle is completed if the circumcircle is completely to the left of the current vertex. // If a triangle is completed, it will be inserted in the output set, unless one or more of it's vertices // belong to the 'super triangle'. class triangleIsCompleted { public: triangleIsCompleted(cvIterator itVertex, triangleSet& output, const vertex SuperTriangle[3]) : m_itVertex(itVertex) , m_Output(output) , m_pSuperTriangle(SuperTriangle) {} bool operator()(const triangle& tri) const { bool b = tri.IsLeftOf(m_itVertex); if (b) { triangleHasVertex thv(m_pSuperTriangle); if (! thv(tri)) m_Output.insert(tri); } return b; }protected: cvIterator m_itVertex; triangleSet& m_Output; const vertex * m_pSuperTriangle; };
// Function object to check whether vertex is in circumcircle of triangle. // operator() returns true if it does. // The edges of a 'hot' triangle are stored in the edgeSet edges. class vertexIsInCircumCircle { public: vertexIsInCircumCircle(cvIterator itVertex, edgeSet& edges) : m_itVertex(itVertex), m_Edges(edges) {} bool operator()(const triangle& tri) const { bool b = tri.CCEncompasses(m_itVertex); if (b) { HandleEdge(tri.GetVertex(0), tri.GetVertex(1)); HandleEdge(tri.GetVertex(1), tri.GetVertex(2)); HandleEdge(tri.GetVertex(2), tri.GetVertex(0)); } return b; } protected: void HandleEdge(const vertex * p0, const vertex * p1) const { const vertex * pVertex0(NULL); const vertex * pVertex1(NULL); // Create a normalized edge, in which the smallest vertex comes first. if (* p0 < * p1) { pVertex0 = p0; pVertex1 = p1; } else { pVertex0 = p1; pVertex1 = p0; } edge e(pVertex0, pVertex1); // Check if this edge is already in the buffer edgeIterator found = m_Edges.find(e); if (found == m_Edges.end()) m_Edges.insert(e); // no, it isn't, so insert else m_Edges.erase(found); // yes, it is, so erase it to eliminate double edges } cvIterator m_itVertex; edgeSet& m_Edges; };void Delaunay::Triangulate(const vertexSet& vertices, triangleSet& output) { if (vertices.size() < 3) return; // nothing to handle // Determine the bounding box. cvIterator itVertex = vertices.begin(); REAL xMin = itVertex->GetX(); REAL yMin = itVertex->GetY(); REAL xMax = xMin; REAL yMax = yMin; ++itVertex; // If we're here, we know that vertices is not empty. for (; itVertex != vertices.end(); itVertex++) { xMax = itVertex->GetX(); // Vertices are sorted along the x-axis, so the last one stored will be the biggest. REAL y = itVertex->GetY(); if (y < yMin) yMin = y; if (y > yMax) yMax = y; } REAL dx = xMax - xMin; REAL dy = yMax - yMin; // Make the bounding box slightly bigger, just to feel safe. REAL ddx = dx * 0.01F; REAL ddy = dy * 0.01F; xMin -= ddx; xMax += ddx; dx += 2 * ddx; yMin -= ddy; yMax += ddy; dy += 2 * ddy; // Create a 'super triangle', encompassing all the vertices. We choose an equilateral triangle with horizontal base. // We could have made the 'super triangle' simply very big. However, the algorithm is quite sensitive to // rounding errors, so it's better to make the 'super triangle' just big enough, like we do here. vertex vSuper[3]; vSuper[0] = vertex(xMin - dy * sqrt3 / 3.0F, yMin,0.0F); // Simple highschool geometry, believe me. vSuper[1] = vertex(xMax + dy * sqrt3 / 3.0F, yMin,0.0F); vSuper[2] = vertex((xMin + xMax) * 0.5F, yMax + dx * sqrt3 * 0.5F,0.0F); triangleSet workset; workset.insert(triangle(vSuper)); for (itVertex = vertices.begin(); itVertex != vertices.end(); itVertex++) { // First, remove all 'completed' triangles from the workset. // A triangle is 'completed' if its circumcircle is entirely to the left of the current vertex. // Vertices are sorted in x-direction (the set container does this automagically). // Unless they are part of the 'super triangle', copy the 'completed' triangles to the output. // The algorithm also works without this step, but it is an important optimalization for bigger numbers of vertices. // It makes the algorithm about five times faster for 2000 vertices, and for 10000 vertices, // it's thirty times faster. For smaller numbers, the difference is negligible. tIterator itEnd = remove_if(workset.begin(), workset.end(), triangleIsCompleted(itVertex, output, vSuper)); edgeSet edges; // A triangle is 'hot' if the current vertex v is inside the circumcircle. // Remove all hot triangles, but keep their edges. itEnd = remove_if(workset.begin(), itEnd, vertexIsInCircumCircle(itVertex, edges)); workset.erase(itEnd, workset.end()); // remove_if doesn't actually remove; we have to do this explicitly. // Create new triangles from the edges and the current vertex. for (edgeIterator it = edges.begin(); it != edges.end(); it++) workset.insert(triangle(it->m_pV0, it->m_pV1, & (* itVertex))); } // Finally, remove all the triangles belonging to the 'super triangle' and move the remaining // triangles tot the output; remove_copy_if lets us do that in one go. tIterator where = output.begin(); remove_copy_if(workset.begin(), workset.end(), inserter(output, where), triangleHasVertex(vSuper)); }void Delaunay::TrianglesToEdges(const triangleSet& triangles, edgeSet& edges) { for (ctIterator it = triangles.begin(); it != triangles.end(); ++it) { HandleEdge(it->GetVertex(0), it->GetVertex(1), edges); HandleEdge(it->GetVertex(1), it->GetVertex(2), edges); HandleEdge(it->GetVertex(2), it->GetVertex(0), edges); } }void Delaunay::HandleEdge(const vertex * p0, const vertex * p1, edgeSet& edges) { const vertex * pV0(NULL); const vertex * pV1(NULL); if (* p0 < * p1) { pV0 = p0; pV1 = p1; } else { pV0 = p1; pV1 = p0; } // Insert a normalized edge. If it's already in edges, insertion will fail, // thus leaving only unique edges. edges.insert(edge(pV0, pV1)); }
http://www.codeguru.com/Cpp/data/mfc_database/misc/article.php/c8901/
有我回复过的
// 头文件#pragma once#include <set>
#include <algorithm>
#include <math.h>using namespace std;#ifndef _GDIPLUS_H// I designed this with GDI+ in mind. However, this particular code doesn't
// use GDI+ at all, only some of it's variable types.
// These definitions are substitutes for those of GDI+.
typedef float REAL;
struct PointF
{
PointF() : X(0), Y(0) {propZ=0;}
PointF(const PointF& p) : X(p.X), Y(p.Y) {propZ=p.propZ;}
PointF(REAL x, REAL y, REAL z) : X(x), Y(y) {propZ=z;}
PointF operator+(const PointF& p) const { return PointF(X + p.X, Y + p.Y, propZ + p.propZ); }
PointF operator-(const PointF& p) const { return PointF(X - p.X, Y - p.Y, propZ - p.propZ); }
REAL X;
REAL Y;
REAL propZ; // 用于等值线的点值,add by zhaoyang.
};const REAL REAL_EPSILON = 1.192092896e-07F; // = 2^-23; I've no idea why this is a good value, but GDI+ has it.#endif // _GDIPLUS_H///////////////////
// vertexclass vertex
{
public:
vertex() : m_Pnt(0.0F, 0.0F, 0.0F) {}
vertex(const vertex& v) : m_Pnt(v.m_Pnt) {}
vertex(const PointF& pnt) : m_Pnt(pnt) {}
vertex(REAL x, REAL y, REAL z) : m_Pnt(x, y, z) {}
vertex(int x, int y, int z) : m_Pnt((REAL) x, (REAL) y, (REAL) z) {} bool operator<(const vertex& v) const
{
if (m_Pnt.X == v.m_Pnt.X) return m_Pnt.Y < v.m_Pnt.Y;
return m_Pnt.X < v.m_Pnt.X;
} bool operator==(const vertex& v) const
{
return m_Pnt.X == v.m_Pnt.X && m_Pnt.Y == v.m_Pnt.Y;
}
REAL GetX() const { return m_Pnt.X; }
REAL GetY() const { return m_Pnt.Y; }
REAL GetZ() const { return m_Pnt.propZ;}; void SetX(REAL x) { m_Pnt.X = x; }
void SetY(REAL y) { m_Pnt.Y = y; }
void SetZ(REAL z) { m_Pnt.propZ = z; } const PointF& GetPoint() const { return m_Pnt; }
protected:
PointF m_Pnt;
};typedef set<vertex> vertexSet;
typedef set<vertex>::iterator vIterator;
typedef set<vertex>::const_iterator cvIterator;///////////////////
// triangleclass triangle
{
public:
triangle(const triangle& tri)
: m_Center(tri.m_Center)
, m_R(tri.m_R)
, m_R2(tri.m_R2)
{
for (int i = 0; i < 3; i++) m_Vertices[i] = tri.m_Vertices[i];
}
triangle(const vertex * p0, const vertex * p1, const vertex * p2)
{
m_Vertices[0] = p0;
m_Vertices[1] = p1;
m_Vertices[2] = p2;
SetCircumCircle();
}
triangle(const vertex * pV)
{
for (int i = 0; i < 3; i++) m_Vertices[i] = pV++;
SetCircumCircle();
} bool operator<(const triangle& tri) const
{
if (m_Center.X == tri.m_Center.X) return m_Center.Y < tri.m_Center.Y;
return m_Center.X < tri.m_Center.X;
} const vertex * GetVertex(int i) const
{
ASSERT(i >= 0);
ASSERT(i < 3);
return m_Vertices[i];
} bool IsLeftOf(cvIterator itVertex) const
{
// returns true if * itVertex is to the right of the triangle's circumcircle
return itVertex->GetPoint().X > (m_Center.X + m_R);
} bool CCEncompasses(cvIterator itVertex) const
{
// Returns true if * itVertex is in the triangle's circumcircle.
// A vertex exactly on the circle is also considered to be in the circle. // These next few lines look like optimisation, however in practice they are not.
// They even seem to slow down the algorithm somewhat.
// Therefore, I've commented them out. // First check boundary box.
// REAL x = itVertex->GetPoint().X;
//
// if (x > (m_Center.X + m_R)) return false;
// if (x < (m_Center.X - m_R)) return false;
//
// REAL y = itVertex->GetPoint().Y;
//
// if (y > (m_Center.Y + m_R)) return false;
// if (y < (m_Center.Y - m_R)) return false; PointF dist = itVertex->GetPoint() - m_Center; // the distance between v and the circle center
REAL dist2 = dist.X * dist.X + dist.Y * dist.Y; // squared
return dist2 <= m_R2; // compare with squared radius
}
protected:
const vertex * m_Vertices[3]; // the three triangle vertices
PointF m_Center; // center of circumcircle
REAL m_R; // radius of circumcircle
REAL m_R2; // radius of circumcircle, squared void SetCircumCircle();
};// Changed in verion 1.1: collect triangles in a multiset.
// In version 1.0, I used a set, preventing the creation of multiple
// triangles with identical center points. Therefore, more than three
// co-circular vertices yielded incorrect results. Thanks to Roger Labbe.
typedef multiset<triangle> triangleSet;
typedef multiset<triangle>::iterator tIterator;
typedef multiset<triangle>::const_iterator ctIterator;///////////////////
// edgeclass edge
{
public:
edge(const edge& e) : m_pV0(e.m_pV0), m_pV1(e.m_pV1) {}
edge(const vertex * pV0, const vertex * pV1)
: m_pV0(pV0), m_pV1(pV1)
{
} bool operator<(const edge& e) const
{
if (m_pV0 == e.m_pV0) return * m_pV1 < * e.m_pV1;
return * m_pV0 < * e.m_pV0;
} const vertex * m_pV0;
const vertex * m_pV1;
};typedef set<edge> edgeSet;
typedef set<edge>::iterator edgeIterator;
typedef set<edge>::const_iterator cedgeIterator;///////////////////
// Delaunayclass Delaunay
{
public:
// Calculate the Delaunay triangulation for the given set of vertices.
void Triangulate(const vertexSet& vertices, triangleSet& output); // Put the edges of the triangles in an edgeSet, eliminating double edges.
// This comes in useful for drawing the triangulation.
void TrianglesToEdges(const triangleSet& triangles, edgeSet& edges);
protected:
void HandleEdge(const vertex * p0, const vertex * p1, edgeSet& edges);
};
#include "StdAfx.h"
#include "Delaunay.h"const REAL sqrt3 = 1.732050808F;void triangle::SetCircumCircle()
{
REAL x0 = m_Vertices[0]->GetX();
REAL y0 = m_Vertices[0]->GetY(); REAL x1 = m_Vertices[1]->GetX();
REAL y1 = m_Vertices[1]->GetY(); REAL x2 = m_Vertices[2]->GetX();
REAL y2 = m_Vertices[2]->GetY(); REAL y10 = y1 - y0;
REAL y21 = y2 - y1; bool b21zero = y21 > -REAL_EPSILON && y21 < REAL_EPSILON; if (y10 > -REAL_EPSILON && y10 < REAL_EPSILON)
{
if (b21zero) // All three vertices are on one horizontal line.
{
if (x1 > x0)
{
if (x2 > x1) x1 = x2;
}
else
{
if (x2 < x0) x0 = x2;
}
m_Center.X = (x0 + x1) * .5F;
m_Center.Y = y0;
}
else // m_Vertices[0] and m_Vertices[1] are on one horizontal line.
{
REAL m1 = - (x2 - x1) / y21; REAL mx1 = (x1 + x2) * .5F;
REAL my1 = (y1 + y2) * .5F; m_Center.X = (x0 + x1) * .5F;
m_Center.Y = m1 * (m_Center.X - mx1) + my1;
}
}
else if (b21zero) // m_Vertices[1] and m_Vertices[2] are on one horizontal line.
{
REAL m0 = - (x1 - x0) / y10; REAL mx0 = (x0 + x1) * .5F;
REAL my0 = (y0 + y1) * .5F; m_Center.X = (x1 + x2) * .5F;
m_Center.Y = m0 * (m_Center.X - mx0) + my0;
}
else // 'Common' cases, no multiple vertices are on one horizontal line.
{
REAL m0 = - (x1 - x0) / y10;
REAL m1 = - (x2 - x1) / y21; REAL mx0 = (x0 + x1) * .5F;
REAL my0 = (y0 + y1) * .5F; REAL mx1 = (x1 + x2) * .5F;
REAL my1 = (y1 + y2) * .5F; m_Center.X = (m0 * mx0 - m1 * mx1 + my1 - my0) / (m0 - m1);
m_Center.Y = m0 * (m_Center.X - mx0) + my0;
} REAL dx = x0 - m_Center.X;
REAL dy = y0 - m_Center.Y; m_R2 = dx * dx + dy * dy; // the radius of the circumcircle, squared
m_R = (REAL) sqrt(m_R2); // the proper radius // Version 1.1: make m_R2 slightly higher to ensure that all edges
// of co-circular vertices will be caught.
// Note that this is a compromise. In fact, the algorithm isn't really
// suited for very many co-circular vertices.
m_R2 *= 1.000001f;
}// Function object to check whether a triangle has one of the vertices in SuperTriangle.
// operator() returns true if it does.
class triangleHasVertex
{
public:
triangleHasVertex(const vertex SuperTriangle[3]) : m_pSuperTriangle(SuperTriangle) {}
bool operator()(const triangle& tri) const
{
for (int i = 0; i < 3; i++)
{
const vertex * p = tri.GetVertex(i);
if (p >= m_pSuperTriangle && p < (m_pSuperTriangle + 3)) return true;
}
return false;
}
protected:
const vertex * m_pSuperTriangle;
};// Function object to check whether a triangle is 'completed', i.e. doesn't need to be checked
// again in the algorithm, i.e. it won't be changed anymore.
// Therefore it can be removed from the workset.
// A triangle is completed if the circumcircle is completely to the left of the current vertex.
// If a triangle is completed, it will be inserted in the output set, unless one or more of it's vertices
// belong to the 'super triangle'.
class triangleIsCompleted
{
public:
triangleIsCompleted(cvIterator itVertex, triangleSet& output, const vertex SuperTriangle[3])
: m_itVertex(itVertex)
, m_Output(output)
, m_pSuperTriangle(SuperTriangle)
{}
bool operator()(const triangle& tri) const
{
bool b = tri.IsLeftOf(m_itVertex); if (b)
{
triangleHasVertex thv(m_pSuperTriangle);
if (! thv(tri)) m_Output.insert(tri);
}
return b;
}protected:
cvIterator m_itVertex;
triangleSet& m_Output;
const vertex * m_pSuperTriangle;
};
// operator() returns true if it does.
// The edges of a 'hot' triangle are stored in the edgeSet edges.
class vertexIsInCircumCircle
{
public:
vertexIsInCircumCircle(cvIterator itVertex, edgeSet& edges) : m_itVertex(itVertex), m_Edges(edges) {}
bool operator()(const triangle& tri) const
{
bool b = tri.CCEncompasses(m_itVertex); if (b)
{
HandleEdge(tri.GetVertex(0), tri.GetVertex(1));
HandleEdge(tri.GetVertex(1), tri.GetVertex(2));
HandleEdge(tri.GetVertex(2), tri.GetVertex(0));
}
return b;
}
protected:
void HandleEdge(const vertex * p0, const vertex * p1) const
{
const vertex * pVertex0(NULL);
const vertex * pVertex1(NULL); // Create a normalized edge, in which the smallest vertex comes first.
if (* p0 < * p1)
{
pVertex0 = p0;
pVertex1 = p1;
}
else
{
pVertex0 = p1;
pVertex1 = p0;
} edge e(pVertex0, pVertex1); // Check if this edge is already in the buffer
edgeIterator found = m_Edges.find(e); if (found == m_Edges.end()) m_Edges.insert(e); // no, it isn't, so insert
else m_Edges.erase(found); // yes, it is, so erase it to eliminate double edges
} cvIterator m_itVertex;
edgeSet& m_Edges;
};void Delaunay::Triangulate(const vertexSet& vertices, triangleSet& output)
{
if (vertices.size() < 3) return; // nothing to handle // Determine the bounding box.
cvIterator itVertex = vertices.begin(); REAL xMin = itVertex->GetX();
REAL yMin = itVertex->GetY();
REAL xMax = xMin;
REAL yMax = yMin; ++itVertex; // If we're here, we know that vertices is not empty.
for (; itVertex != vertices.end(); itVertex++)
{
xMax = itVertex->GetX(); // Vertices are sorted along the x-axis, so the last one stored will be the biggest.
REAL y = itVertex->GetY();
if (y < yMin) yMin = y;
if (y > yMax) yMax = y;
} REAL dx = xMax - xMin;
REAL dy = yMax - yMin; // Make the bounding box slightly bigger, just to feel safe.
REAL ddx = dx * 0.01F;
REAL ddy = dy * 0.01F; xMin -= ddx;
xMax += ddx;
dx += 2 * ddx; yMin -= ddy;
yMax += ddy;
dy += 2 * ddy; // Create a 'super triangle', encompassing all the vertices. We choose an equilateral triangle with horizontal base.
// We could have made the 'super triangle' simply very big. However, the algorithm is quite sensitive to
// rounding errors, so it's better to make the 'super triangle' just big enough, like we do here.
vertex vSuper[3]; vSuper[0] = vertex(xMin - dy * sqrt3 / 3.0F, yMin,0.0F); // Simple highschool geometry, believe me.
vSuper[1] = vertex(xMax + dy * sqrt3 / 3.0F, yMin,0.0F);
vSuper[2] = vertex((xMin + xMax) * 0.5F, yMax + dx * sqrt3 * 0.5F,0.0F); triangleSet workset;
workset.insert(triangle(vSuper)); for (itVertex = vertices.begin(); itVertex != vertices.end(); itVertex++)
{
// First, remove all 'completed' triangles from the workset.
// A triangle is 'completed' if its circumcircle is entirely to the left of the current vertex.
// Vertices are sorted in x-direction (the set container does this automagically).
// Unless they are part of the 'super triangle', copy the 'completed' triangles to the output.
// The algorithm also works without this step, but it is an important optimalization for bigger numbers of vertices.
// It makes the algorithm about five times faster for 2000 vertices, and for 10000 vertices,
// it's thirty times faster. For smaller numbers, the difference is negligible.
tIterator itEnd = remove_if(workset.begin(), workset.end(), triangleIsCompleted(itVertex, output, vSuper)); edgeSet edges; // A triangle is 'hot' if the current vertex v is inside the circumcircle.
// Remove all hot triangles, but keep their edges.
itEnd = remove_if(workset.begin(), itEnd, vertexIsInCircumCircle(itVertex, edges));
workset.erase(itEnd, workset.end()); // remove_if doesn't actually remove; we have to do this explicitly. // Create new triangles from the edges and the current vertex.
for (edgeIterator it = edges.begin(); it != edges.end(); it++)
workset.insert(triangle(it->m_pV0, it->m_pV1, & (* itVertex)));
} // Finally, remove all the triangles belonging to the 'super triangle' and move the remaining
// triangles tot the output; remove_copy_if lets us do that in one go.
tIterator where = output.begin();
remove_copy_if(workset.begin(), workset.end(), inserter(output, where), triangleHasVertex(vSuper));
}void Delaunay::TrianglesToEdges(const triangleSet& triangles, edgeSet& edges)
{
for (ctIterator it = triangles.begin(); it != triangles.end(); ++it)
{
HandleEdge(it->GetVertex(0), it->GetVertex(1), edges);
HandleEdge(it->GetVertex(1), it->GetVertex(2), edges);
HandleEdge(it->GetVertex(2), it->GetVertex(0), edges);
}
}void Delaunay::HandleEdge(const vertex * p0, const vertex * p1, edgeSet& edges)
{
const vertex * pV0(NULL);
const vertex * pV1(NULL); if (* p0 < * p1)
{
pV0 = p0;
pV1 = p1;
}
else
{
pV0 = p1;
pV1 = p0;
} // Insert a normalized edge. If it's already in edges, insertion will fail,
// thus leaving only unique edges.
edges.insert(edge(pV0, pV1));
}