我正在尝试用 C++ 实现 Delaunay 三角剖分。目前它正在工作，但我没有得到正确数量的三角形。 我用方形图案中的 4 个点进行尝试:(0,0)、(1,0)、(0,1)、(1,1)。
这是我使用的算法:
std::vector<Triangle> Delaunay::triangulate(std::vector<Vec2f> &vertices) {
// Determinate the super triangle
float minX = vertices[0].getX();
float minY = vertices[0].getY();
float maxX = minX;
float maxY = minY;
for(std::size_t i = 0; i < vertices.size(); ++i) {
if (vertices[i].getX() < minX) minX = vertices[i].getX();
if (vertices[i].getY() < minY) minY = vertices[i].getY();
if (vertices[i].getX() > maxX) maxX = vertices[i].getX();
if (vertices[i].getY() > maxY) maxY = vertices[i].getY();
}
float dx = maxX - minX;
float dy = maxY - minY;
float deltaMax = std::max(dx, dy);
float midx = (minX + maxX) / 2.f;
float midy = (minY + maxY) / 2.f;
Vec2f p1(midx - 20 * deltaMax, midy - deltaMax);
Vec2f p2(midx, midy + 20 * deltaMax);
Vec2f p3(midx + 20 * deltaMax, midy - deltaMax);
// Add the super triangle vertices to the end of the vertex list
vertices.push_back(p1);
vertices.push_back(p2);
vertices.push_back(p3);
// Add the super triangle to the triangle list
std::vector<Triangle> triangleList = {Triangle(p1, p2, p3)};
// For each point in the vertex list
for(auto point = begin(vertices); point != end(vertices); point++)
{
// Initialize the edges buffer
std::vector<Edge> edgesBuff;
// For each triangles currently in the triangle list
for(auto triangle = begin(triangleList); triangle != end(triangleList);)
{
if(triangle->inCircumCircle(*point))
{
Edge tmp[3] = {triangle->getE1(), triangle->getE2(), triangle->getE3()};
edgesBuff.insert(end(edgesBuff), tmp, tmp + 3);
triangle = triangleList.erase(triangle);
}
else
{
triangle++;
}
}
// Delete all doubly specified edges from the edge buffer
// Black magic by https://github.com/MechaRage
auto ite = begin(edgesBuff), last = end(edgesBuff);
while(ite != last) {
// Search for at least one duplicate of the current element
auto twin = std::find(ite + 1, last, *ite);
if(twin != last)
// If one is found, push them all to the end.
last = std::partition(ite, last, [&ite](auto const &o){ return !(o == *ite); });
else
++ite;
}
// Remove all the duplicates, which have been shoved past "last".
edgesBuff.erase(last, end(edgesBuff));
// Add the triangle to the list
for(auto edge = begin(edgesBuff); edge != end(edgesBuff); edge++)
triangleList.push_back(Triangle(edge->getP1(), edge->getP2(), *point));
}
// Remove any triangles from the triangle list that use the supertriangle vertices
triangleList.erase(std::remove_if(begin(triangleList), end(triangleList), [p1, p2, p3](auto t){
return t.containsVertex(p1) || t.containsVertex(p2) || t.containsVertex(p3);
}), end(triangleList));
return triangleList;
这是我得到的:
Triangle:
Point x: 1 y: 0
Point x: 0 y: 0
Point x: 1 y: 1
Triangle:
Point x: 1 y: 0
Point x: 1 y: 1
Point x: 0 y: 1
Triangle:
Point x: 0 y: 0
Point x: 1 y: 1
Point x: 0 y: 1
虽然这是正确的输出:
Triangle:
Point x: 1 y: 0
Point x: 0 y: 0
Point x: 0 y: 1
Triangle:
Point x: 1 y: 0
Point x: 1 y: 1
Point x: 0 y: 1
我不知道为什么会有 (0, 0) 和 (1, 1) 的三角形。 我需要一个局外人来检查代码并找出问题所在。
所有资源都在我的 Github repo 上.随意 fork 它并公开您的代码。
谢谢!
最佳答案
Paul Bourke 的 Delaunay 三角剖分算法的实现情况如何。看看 Triangulate() 我已经多次使用这个源，没有任何提示
#include <iostream>
#include <stdlib.h> // for C qsort
#include <cmath>
#include <time.h> // for random
const int MaxVertices = 500;
const int MaxTriangles = 1000;
//const int n_MaxPoints = 10; // for the test programm
const double EPSILON = 0.000001;
struct ITRIANGLE{
int p1, p2, p3;
};
struct IEDGE{
int p1, p2;
};
struct XYZ{
double x, y, z;
};
int XYZCompare(const void *v1, const void *v2);
int Triangulate(int nv, XYZ pxyz[], ITRIANGLE v[], int &ntri);
int CircumCircle(double, double, double, double, double, double, double, double, double&, double&, double&);
using namespace std;
////////////////////////////////////////////////////////////////////////
// CircumCircle() :
// Return true if a point (xp,yp) is inside the circumcircle made up
// of the points (x1,y1), (x2,y2), (x3,y3)
// The circumcircle centre is returned in (xc,yc) and the radius r
// Note : A point on the edge is inside the circumcircle
////////////////////////////////////////////////////////////////////////
int CircumCircle(double xp, double yp, double x1, double y1, double x2,
double y2, double x3, double y3, double &xc, double &yc, double &r){
double m1, m2, mx1, mx2, my1, my2;
double dx, dy, rsqr, drsqr;
/* Check for coincident points */
if(abs(y1 - y2) < EPSILON && abs(y2 - y3) < EPSILON)
return(false);
if(abs(y2-y1) < EPSILON){
m2 = - (x3 - x2) / (y3 - y2);
mx2 = (x2 + x3) / 2.0;
my2 = (y2 + y3) / 2.0;
xc = (x2 + x1) / 2.0;
yc = m2 * (xc - mx2) + my2;
}else if(abs(y3 - y2) < EPSILON){
m1 = - (x2 - x1) / (y2 - y1);
mx1 = (x1 + x2) / 2.0;
my1 = (y1 + y2) / 2.0;
xc = (x3 + x2) / 2.0;
yc = m1 * (xc - mx1) + my1;
}else{
m1 = - (x2 - x1) / (y2 - y1);
m2 = - (x3 - x2) / (y3 - y2);
mx1 = (x1 + x2) / 2.0;
mx2 = (x2 + x3) / 2.0;
my1 = (y1 + y2) / 2.0;
my2 = (y2 + y3) / 2.0;
xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
yc = m1 * (xc - mx1) + my1;
}
dx = x2 - xc;
dy = y2 - yc;
rsqr = dx * dx + dy * dy;
r = sqrt(rsqr);
dx = xp - xc;
dy = yp - yc;
drsqr = dx * dx + dy * dy;
return((drsqr <= rsqr) ? true : false);
}
///////////////////////////////////////////////////////////////////////////////
// Triangulate() :
// Triangulation subroutine
// Takes as input NV vertices in array pxyz
// Returned is a list of ntri triangular faces in the array v
// These triangles are arranged in a consistent clockwise order.
// The triangle array 'v' should be malloced to 3 * nv
// The vertex array pxyz must be big enough to hold 3 more points
// The vertex array must be sorted in increasing x values say
//
// qsort(p,nv,sizeof(XYZ),XYZCompare);
///////////////////////////////////////////////////////////////////////////////
int Triangulate(int nv, XYZ pxyz[], ITRIANGLE v[], int &ntri){
int *complete = NULL;
IEDGE *edges = NULL;
IEDGE *p_EdgeTemp;
int nedge = 0;
int trimax, emax = 200;
int status = 0;
int inside;
int i, j, k;
double xp, yp, x1, y1, x2, y2, x3, y3, xc, yc, r;
double xmin, xmax, ymin, ymax, xmid, ymid;
double dx, dy, dmax;
/* Allocate memory for the completeness list, flag for each triangle */
trimax = 4 * nv;
complete = new int[trimax];
/* Allocate memory for the edge list */
edges = new IEDGE[emax];
/*
Find the maximum and minimum vertex bounds.
This is to allow calculation of the bounding triangle
*/
xmin = pxyz[0].x;
ymin = pxyz[0].y;
xmax = xmin;
ymax = ymin;
for(i = 1; i < nv; i++){
if (pxyz[i].x < xmin) xmin = pxyz[i].x;
if (pxyz[i].x > xmax) xmax = pxyz[i].x;
if (pxyz[i].y < ymin) ymin = pxyz[i].y;
if (pxyz[i].y > ymax) ymax = pxyz[i].y;
}
dx = xmax - xmin;
dy = ymax - ymin;
dmax = (dx > dy) ? dx : dy;
xmid = (xmax + xmin) / 2.0;
ymid = (ymax + ymin) / 2.0;
/*
Set up the supertriangle
his is a triangle which encompasses all the sample points.
The supertriangle coordinates are added to the end of the
vertex list. The supertriangle is the first triangle in
the triangle list.
*/
pxyz[nv+0].x = xmid - 20 * dmax;
pxyz[nv+0].y = ymid - dmax;
pxyz[nv+1].x = xmid;
pxyz[nv+1].y = ymid + 20 * dmax;
pxyz[nv+2].x = xmid + 20 * dmax;
pxyz[nv+2].y = ymid - dmax;
v[0].p1 = nv;
v[0].p2 = nv+1;
v[0].p3 = nv+2;
complete[0] = false;
ntri = 1;
/*
Include each point one at a time into the existing mesh
*/
for(i = 0; i < nv; i++){
xp = pxyz[i].x;
yp = pxyz[i].y;
nedge = 0;
/*
Set up the edge buffer.
If the point (xp,yp) lies inside the circumcircle then the
three edges of that triangle are added to the edge buffer
and that triangle is removed.
*/
for(j = 0; j < ntri; j++){
if(complete[j])
continue;
x1 = pxyz[v[j].p1].x;
y1 = pxyz[v[j].p1].y;
x2 = pxyz[v[j].p2].x;
y2 = pxyz[v[j].p2].y;
x3 = pxyz[v[j].p3].x;
y3 = pxyz[v[j].p3].y;
inside = CircumCircle(xp, yp, x1, y1, x2, y2, x3, y3, xc, yc, r);
if (xc + r < xp)
// Suggested
// if (xc + r + EPSILON < xp)
complete[j] = true;
if(inside){
/* Check that we haven't exceeded the edge list size */
if(nedge + 3 >= emax){
emax += 100;
p_EdgeTemp = new IEDGE[emax];
for (int i = 0; i < nedge; i++) { // Fix by John Bowman
p_EdgeTemp[i] = edges[i];
}
delete []edges;
edges = p_EdgeTemp;
}
edges[nedge+0].p1 = v[j].p1;
edges[nedge+0].p2 = v[j].p2;
edges[nedge+1].p1 = v[j].p2;
edges[nedge+1].p2 = v[j].p3;
edges[nedge+2].p1 = v[j].p3;
edges[nedge+2].p2 = v[j].p1;
nedge += 3;
v[j] = v[ntri-1];
complete[j] = complete[ntri-1];
ntri--;
j--;
}
}
/*
Tag multiple edges
Note: if all triangles are specified anticlockwise then all
interior edges are opposite pointing in direction.
*/
for(j = 0; j < nedge - 1; j++){
for(k = j + 1; k < nedge; k++){
if((edges[j].p1 == edges[k].p2) && (edges[j].p2 == edges[k].p1)){
edges[j].p1 = -1;
edges[j].p2 = -1;
edges[k].p1 = -1;
edges[k].p2 = -1;
}
/* Shouldn't need the following, see note above */
if((edges[j].p1 == edges[k].p1) && (edges[j].p2 == edges[k].p2)){
edges[j].p1 = -1;
edges[j].p2 = -1;
edges[k].p1 = -1;
edges[k].p2 = -1;
}
}
}
/*
Form new triangles for the current point
Skipping over any tagged edges.
All edges are arranged in clockwise order.
*/
for(j = 0; j < nedge; j++) {
if(edges[j].p1 < 0 || edges[j].p2 < 0)
continue;
v[ntri].p1 = edges[j].p1;
v[ntri].p2 = edges[j].p2;
v[ntri].p3 = i;
complete[ntri] = false;
ntri++;
}
}
/*
Remove triangles with supertriangle vertices
These are triangles which have a vertex number greater than nv
*/
for(i = 0; i < ntri; i++) {
if(v[i].p1 >= nv || v[i].p2 >= nv || v[i].p3 >= nv) {
v[i] = v[ntri-1];
ntri--;
i--;
}
}
delete[] edges;
delete[] complete;
return 0;
}
int XYZCompare(const void *v1, const void *v2){
XYZ *p1, *p2;
p1 = (XYZ*)v1;
p2 = (XYZ*)v2;
if(p1->x < p2->x)
return(-1);
else if(p1->x > p2->x)
return(1);
else
return(0);
}
关于c++ - Delaunay 三角剖分 : too many triangles，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/33590159/