2015-11-08 145 views
1

我正試圖在C++中實現Delaunay三角剖分。目前它正在工作,但我沒有得到正確數量的三角形。 我用4點的方格式來嘗試它:(0,0),(1,0),(0,1),(1,1)。Delaunay三角剖分:太多的三角形

這是我使用的算法:

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 

我沒有我dea爲什麼有一個與(0,0)和(1,1)的三角形。 我需要一個外部的眼睛來審查代碼,並找出發生了什麼問題。

所有來源都在我的Github repo上。隨意分叉它,並PR您的代碼。

謝謝!

+0

歡迎so so!我馬上發現的是,第一個,有缺陷的三角形的集合和預期結果的集合都是均勻的。 –

+0

我真的不明白你在說什麼(英語不是我的第一語言)。你在談論我處理超級三角形或程序輸出的方式? – Bl4ckb0ne

+0

輸出。只需畫一個正方形,標出點,然後用鉛筆沿着點到點的路徑「走」。 –

回答

0

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

這是我用於我的C++版本。但是我發現C版有點難以閱讀。 – Bl4ckb0ne

0

我沒有用調試器去,但是從所得到的三角形似乎這是一個精度/不確定性問題的樣子。

當你正在對一個正方形進行三角測量時,有兩種方法可以將它分成三角形,並且從Delaunay準則(外接圓的中心位於三角形的邊界上)都是正確的。

因此,如果您獨立評估每個三角形,有時甚至可能會得到4個三角形(取決於實施情況)。

通常在這種情況下,我建議將算法構建爲一系列不能產生矛盾答案的問題。在這種情況下,問題是「基於邊緣(1,0) - (1,1)」到哪個點變爲三角形。但通常這需要對算法進行重大更改。

快速修復通常包括爲比較和額外檢查(如不相交的三角形)添加一些公差。但通常它只會讓問題更爲罕見。

0

很可能你並沒有刪除所有的雙邊,特別是不是來自同一個三角形的邊,而是隻有其他順序的頂點。正確的功能來自@cMinor的答案。