【问题标题】:Delaunay triangulation : too many trianglesDelaunay三角剖分:三角形太多
【发布时间】:2016-02-08 23:30:19
【问题描述】:

我正在尝试在 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 它并 PR 你的代码。

谢谢!

【问题讨论】:

  • 欢迎您!我立即发现的是,第一组有缺陷的三角形的缠绕和预期结果的缠绕都不是均匀的。
  • 我不太明白你在说什么(英语不是我的第一语言)。你说的是我处理超三角形的方式,还是程序的输出?
  • 输出。只需画一个正方形,标记点,然后用铅笔从一个点到另一个点“走”。
  • 我正在考虑使用 SFML 显示输入,但我真的想在做之前把它做好。但我明天会做。
  • 对于两个三角形,一张纸就足够了。不过“明天”是个好关键词,那我再仔细看看你的德劳内。

标签: c++ algorithm delaunay


【解决方案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);
}

【讨论】:

  • 那是我用于我的 C++ 版本的那个。但我觉得 C 版本有点难读。
【解决方案2】:

我没有使用调试器,但从生成的三角形来看,这似乎是一个准确性/模糊性问题。

当您对一个正方形进行三角剖分时,有两种方法可以将其拆分为三角形,并且都符合 Delaunay 标准(外接圆中心位于三角形的边界上)。

因此,如果您独立评估每个三角形,有时甚至可能得到 4 个三角形(取决于实现)。

通常在这种情况下,我建议将算法构建为一系列不能产生矛盾答案的问题。在这种情况下,问题是“基于边 (1,0)-(1,1) 的三角形到哪一点”。但这通常需要对算法进行重大更改。

快速修复通常涉及为比较和额外检查添加一些容差(例如不相交的三角形)。但通常它只会让问题变得更少。

【讨论】:

    【解决方案3】:

    很可能您没有删除所有的双边,尤其是没有删除来自相同三角形的边,但顶点仅以另一个顺序排列。正确的功能在@cMinor 的答案中。

    【讨论】:

      猜你喜欢
      • 2014-01-21
      • 2015-01-07
      • 2012-03-04
      • 2011-07-15
      • 1970-01-01
      • 2019-04-18
      • 2020-10-17
      • 2019-09-28
      • 2021-05-14
      相关资源
      最近更新 更多