c++ - Delaunay 三角剖分 : too many triangles

Question

我正在尝试用 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 它并公开您的代码。

谢谢!

Answer 1

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);
}