zhangyiss/delaunay2d
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

complete the flip algorithm

Browse Source
This commit is contained in:
张壹
2021-09-19 22:19:08 +08:00
parent 2a69d72aec
commit 30abc3ede5
9 changed files with 955 additions and 146 deletions
Download Patch File Download Diff File Expand all files Collapse all files

19
.vscode/settings.json vendored Normal file
View File

@@ -0,0 +1,19 @@
{
"files.associations": {
"__hash_table": "cpp",
"__split_buffer": "cpp",
"array": "cpp",
"bitset": "cpp",
"deque": "cpp",
"initializer_list": "cpp",
"iterator": "cpp",
"queue": "cpp",
"random": "cpp",
"stack": "cpp",
"string": "cpp",
"string_view": "cpp",
"unordered_map": "cpp",
"utility": "cpp",
"vector": "cpp"
}
}

379
delaunay.h
View File

@@ -1,14 +1,14 @@
/** /**
* @defgroup DELAUNAY * @defgroup DELAUNAY
* *
* @brief An implementation of the 2D Delaunay triangulation using the Bowyer-Watson algorithm. * @brief An implementation of the 2D Delaunay triangulation using the flip algorithm.
* *
* @author Yi Zhang * @author Yi Zhang
* @date 2021-09-12 * @date 2021-09-19
*/ */
#ifndef _BW_2D_DELAUNAY_H #ifndef _FLIP_2D_DELAUNAY_H
#define _BW_2D_DELAUNAY_H #define _FLIP_2D_DELAUNAY_H
#include "cmath" #include "cmath"
#include "vector" #include "vector"
@@ -37,56 +37,181 @@ bool operator ==(const vertex2dc &a, const vertex2dc &b) // overload the == oper
} }
return false; return false;
} }
// End vertex definition
// Start edge definition void circumcircle(vertex2dc *v0, vertex2dc *v1, vertex2dc *v2, double &cx, double &cy, double &cr) // calculate the circumcircle from three points
struct edge
{ {
vertex2dc *vert[2]; // vertex of the edge double s = 0.5 / ((v1->x - v0->x) * (v2->y - v0->y) - (v1->y - v0->y) * (v2->x - v0->x));
double m = v1->x*v1->x - v0->x*v0->x + v1->y*v1->y - v0->y*v0->y;
double u = v2->x*v2->x - v0->x*v0->x + v2->y*v2->y - v0->y*v0->y;
edge() {vert[0] = vert[1] = nullptr;} cx = ((v2->y - v0->y)*m + (v0->y - v1->y)*u)*s;
edge(vertex2dc *v0ptr, vertex2dc *v1ptr) {set(v0ptr, v1ptr);} cy = ((v0->x - v2->x)*m + (v1->x - v0->x)*u)*s;
void set(vertex2dc *v0ptr, vertex2dc *v1ptr) cr = (v0->x - cx)*(v0->x - cx) + (v0->y - cy)*(v0->y - cy); // not need to calculate the squared root here
{ return;
vert[0] = v0ptr; vert[1] = v1ptr;
return;
}
};
bool operator ==(const edge &a, const edge &b) // overload the == operator for edge type
{
if((a.vert[0] == b.vert[0] && a.vert[1] == b.vert[1]) ||
(a.vert[0] == b.vert[1] && a.vert[1] == b.vert[0]))
{
return true;
}
return false;
} }
// End edge definition // End vertex definition
// Start triangle definition // Start triangle definition
struct triangle struct triangle
{ {
int id; // index of the triangle
vertex2dc *vert[3]; // vertex of the triangle vertex2dc *vert[3]; // vertex of the triangle
triangle *neigh[3]; // neighbors of the triangle
double cx, cy; // center of the triangle's circumcircle double cx, cy; // center of the triangle's circumcircle
double cr; // radius of the circumcircle double cr; // radius of the circumcircle
triangle() {vert[0] = vert[1] = vert[2] = nullptr;} triangle() {vert[0] = vert[1] = vert[2] = nullptr; neigh[0] = neigh[1] = neigh[2] = nullptr;}
triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) {set(v0ptr, v1ptr, v2ptr);} triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) {set(v0ptr, v1ptr, v2ptr);}
void set(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) void set(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr)
{ {
vert[0] = v0ptr; vert[1] = v1ptr; vert[2] = v2ptr; vert[0] = v0ptr; vert[1] = v1ptr; vert[2] = v2ptr;
neigh[0] = neigh[1] = neigh[2] = nullptr;
double s = 0.5 / ((vert[1]->x - vert[0]->x) * (vert[2]->y - vert[0]->y) - (vert[1]->y - vert[0]->y) * (vert[2]->x - vert[0]->x)); circumcircle(vert[0], vert[1], vert[2], cx, cy, cr);
double m = vert[1]->x * vert[1]->x - vert[0]->x * vert[0]->x + vert[1]->y * vert[1]->y - vert[0]->y * vert[0]->y;
double u = vert[2]->x * vert[2]->x - vert[0]->x * vert[0]->x + vert[2]->y * vert[2]->y - vert[0]->y * vert[0]->y;
cx = ((vert[2]->y - vert[0]->y) * m + (vert[0]->y - vert[1]->y) * u) * s;
cy = ((vert[0]->x - vert[2]->x) * m + (vert[1]->x - vert[0]->x) * u) * s;
cr = (vert[0]->x - cx) * (vert[0]->x - cx) + (vert[0]->y - cy) * (vert[0]->y - cy); // not need to sqrt() here
return; return;
} }
void set_neighbor(triangle *n0ptr, triangle *n1ptr, triangle *n2ptr)
{
neigh[0] = n0ptr; neigh[1] = n1ptr; neigh[2] = n2ptr;
return;
}
bool bound_location(double inx, double iny) // Test if the location is inside the triangle
{
double l1x, l1y, l2x, l2y;
for (int i = 0; i < 3; i++)
{
l1x = vert[(i+1)%3]->x - vert[i]->x;
l1y = vert[(i+1)%3]->y - vert[i]->y;
l2x = inx - vert[i]->x;
l2y = iny - vert[i]->y;
if ((l1x*l2y - l1y*l2x) < 0) // This condition includes points on the triangle's edge
{
return false;
}
}
return true;
}
}; };
/**
* @brief Flip neighboring triangles and their neighbors
*
* original
*
* /\
* / \
* / \
* / t \
* t_id-------\ t_id (0, 1 or 2)
* \--------/
* \ /
* \ n /
* \ /
* \/
* n_id (0, 1 or 2)
*
* fliped
*
* /|\
* / | \
* / | \
* / | \
* t_id | \ t_id (0, 1 or 2)
* \ t | n /
* \ | /
* \ | /
* \ | /
* \|/
* n_id (0, 1 or 2)
*
* @param t target triangle
* @param n neighboring triangle
* @param t_vid reference index of the target triangle
* @param n_vid reference index of the neighboring triangle
*/
void flip_neighboring_triangles(triangle *t, triangle *n, int t_id, int n_id)
{
t->vert[(t_id+1)%3] = n->vert[n_id]; // flip t
circumcircle(t->vert[0], t->vert[1], t->vert[2], t->cx, t->cy, t->cr); // update circumcircle
n->vert[(n_id+2)%3] = t->vert[(t_id+2)%3]; // flip n
circumcircle(n->vert[0], n->vert[1], n->vert[2], n->cx, n->cy, n->cr); // update circumcircle
// set side neighbors
t->neigh[t_id] = n->neigh[(n_id+2)%3];
n->neigh[(n_id+1)%3] = t->neigh[(t_id+1)%3];
// set opposite neighbors
t->neigh[(t_id+1)%3] = n;
n->neigh[(n_id+2)%3] = t;
// set oppsite neighbors
if (t->neigh[t_id] != nullptr)
{
for (int i = 0; i < 3; i++)
{
if (t->neigh[t_id]->neigh[i] == n)
{
t->neigh[t_id]->neigh[i] = t;
break;
}
}
}
if (n->neigh[(n_id+1)%3] != nullptr)
{
for (int i = 0; i < 3; i++)
{
if (n->neigh[(n_id+1)%3]->neigh[i] == t)
{
n->neigh[(n_id+1)%3]->neigh[i] = n;
break;
}
}
}
return;
}
/**
* @brief Make sure that the input triangle meet the empty circumcircle condition
*
* @param t Input triangle
*/
void make_delaunay(triangle *t)
{
double dist;
vertex2dc *n_vert;
triangle *n_tri;
for (int n = 0; n < 3; ++n)
{
if (t->neigh[n] != nullptr) // must has neighbor on this side
{
n_tri = t->neigh[n];
for (int v = 0; v < 3; ++v)
{
n_vert = n_tri->vert[v];
if (n_vert != t->vert[n] && n_vert != t->vert[(n+1)%3]) // find the opposite vertex
{
dist = (t->cx - n_vert->x) * (t->cx - n_vert->x) +
(t->cy - n_vert->y) * (t->cy - n_vert->y);
if ((dist - t->cr) <= ZERO) // need to be flipped
{
flip_neighboring_triangles(t, n_tri, n, v);
// Make sure the triangles also meet the empty circumcircle condition after flipping
make_delaunay(t);
make_delaunay(n_tri);
return; // Neighborhood changed. The current loop is not valid any more.
}
break; // no need to search more
}
}
}
}
return;
}
// End triangle definition // End triangle definition
/** /**
@@ -95,7 +220,7 @@ struct triangle
* @param in_verts Input vertexes. Defined by the user. * @param in_verts Input vertexes. Defined by the user.
* @param out_tris Output triangles. Compute by the function. * @param out_tris Output triangles. Compute by the function.
*/ */
void triangulation(std::vector<vertex2dc> &in_verts, std::vector<triangle> &out_tris) void triangulation(std::vector<vertex2dc> &in_verts, std::vector<triangle*> &out_tris)
{ {
if (!out_tris.empty()) out_tris.clear(); if (!out_tris.empty()) out_tris.clear();
if (in_verts.size() < 3) return; if (in_verts.size() < 3) return;
@@ -116,125 +241,155 @@ void triangulation(std::vector<vertex2dc> &in_verts, std::vector<triangle> &out_
double maxi_s = std::max(xmax - xmin, ymax - ymin); // use an four times bigger rectangle to include all points double maxi_s = std::max(xmax - xmin, ymax - ymin); // use an four times bigger rectangle to include all points
vertex2dc *tmp_vert = nullptr; vertex2dc *tmp_vert = nullptr;
std::vector<vertex2dc*> assit_vert; std::vector<vertex2dc*> box_vert;
tmp_vert = new vertex2dc(midx - maxi_s, midy - maxi_s); // lower left corner tmp_vert = new vertex2dc(midx - maxi_s, midy - maxi_s); // lower left corner
assit_vert.push_back(tmp_vert); box_vert.push_back(tmp_vert);
tmp_vert = new vertex2dc(midx + maxi_s, midy - maxi_s); // lower right corner tmp_vert = new vertex2dc(midx + maxi_s, midy - maxi_s); // lower right corner
assit_vert.push_back(tmp_vert); box_vert.push_back(tmp_vert);
tmp_vert = new vertex2dc(midx + maxi_s, midy + maxi_s); // upper right corner tmp_vert = new vertex2dc(midx + maxi_s, midy + maxi_s); // upper right corner
assit_vert.push_back(tmp_vert); box_vert.push_back(tmp_vert);
tmp_vert = new vertex2dc(midx - maxi_s, midy + maxi_s); // upper left corner tmp_vert = new vertex2dc(midx - maxi_s, midy + maxi_s); // upper left corner
assit_vert.push_back(tmp_vert); box_vert.push_back(tmp_vert);
triangle *tmp_tri = nullptr; triangle *old_tri = nullptr, *tmp_tri = nullptr;
std::vector<triangle*> exist_tri, cnst_tri; triangle *cnst_tri[3];
std::vector<triangle*>::iterator t_iter; std::vector<triangle*>::iterator t_iter;
tmp_tri = new triangle(assit_vert[0], assit_vert[1], assit_vert[2]); // order the vertex anti-clock wise tmp_tri = new triangle(box_vert[0], box_vert[1], box_vert[2]); // order the vertex anti-clock wise
exist_tri.push_back(tmp_tri); out_tris.push_back(tmp_tri);
tmp_tri = new triangle(assit_vert[0], assit_vert[2], assit_vert[3]); // order the vertex anti-clock wise tmp_tri = new triangle(box_vert[0], box_vert[2], box_vert[3]); // order the vertex anti-clock wise
exist_tri.push_back(tmp_tri); out_tris.push_back(tmp_tri);
out_tris[0]->set_neighbor(nullptr, nullptr, out_tris[1]);
out_tris[1]->set_neighbor(out_tris[0], nullptr, nullptr);
// loop all input vertice // loop all input vertice
bool removed;
double dist;
edge tmp_edge;
std::vector<edge> cnst_edge;
std::vector<edge>::iterator e_iter;
for (int i = 0; i < in_verts.size(); ++i) for (int i = 0; i < in_verts.size(); ++i)
{ {
// determine triangles that include the point and add the triangle to the cnst_tri and remove the triangle from exist_tri // determine the triangle that includes the new vertex and remove it from out_tris
// this is the part that could take a lot of time if we are working with a large amount of points. We will fix this later for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
cnst_tri.clear();
for (t_iter = exist_tri.begin(); t_iter != exist_tri.end(); )
{ {
tmp_tri = *t_iter; old_tri = *t_iter;
if (old_tri->bound_location(in_verts[i].x, in_verts[i].y))
dist = (tmp_tri->cx - in_verts[i].x) * (tmp_tri->cx - in_verts[i].x) +
(tmp_tri->cy - in_verts[i].y) * (tmp_tri->cy - in_verts[i].y);
if ((dist - tmp_tri->cr) <= ZERO) // Points on the circumcircle are also included
{ {
t_iter = exist_tri.erase(t_iter); t_iter = out_tris.erase(t_iter);
cnst_tri.push_back(tmp_tri); break;
} }
else t_iter++; else t_iter++;
} }
// loop to remove duplicate edges // build three new triangles
cnst_edge.clear(); for (int n = 0; n < 3; ++n)
for (int c = 0; c < cnst_tri.size(); ++c)
{ {
for (int e = 0; e < 3; ++e) tmp_tri = new triangle(old_tri->vert[n], old_tri->vert[(n+1)%3], &in_verts[i]);
cnst_tri[n] = tmp_tri;
out_tris.push_back(tmp_tri);
}
// sort neighbors
for (int n = 0; n < 3; ++n)
{
if (old_tri->neigh[n] == nullptr)
{ {
tmp_edge.set(cnst_tri[c]->vert[e], cnst_tri[c]->vert[(e+1)%3]); cnst_tri[n]->set_neighbor(nullptr, cnst_tri[(n+1)%3], cnst_tri[(n+2)%3]);
}
removed = false; else
for (e_iter = cnst_edge.begin(); e_iter != cnst_edge.end(); ) {
cnst_tri[n]->set_neighbor(old_tri->neigh[n], cnst_tri[(n+1)%3], cnst_tri[(n+2)%3]);
for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle
{ {
if (tmp_edge == *e_iter) // duplicate edge, remove from cnst_edge if (old_tri->neigh[n]->neigh[k] == old_tri)
{ {
e_iter = cnst_edge.erase(e_iter); old_tri->neigh[n]->neigh[k] = cnst_tri[n];
removed = true; break;
break; // no need to search more
} }
else e_iter++;
}
if (!removed) // not a duplicate edge, add to the cnst_edge
{
cnst_edge.push_back(tmp_edge);
} }
} }
} }
// construct new triangles and add to exist_tri // delete the old triangle
for (int c = 0; c < cnst_edge.size(); ++c) delete old_tri; old_tri = nullptr;
{
tmp_tri = new triangle(cnst_edge[c].vert[0], cnst_edge[c].vert[1], &in_verts[i]); // order the vertex anti-clock wise
exist_tri.push_back(tmp_tri);
}
// destroy memories used by cnst_edge // Make sure cnst_tri meet the empty circumcircle condition
for (int c = 0; c < cnst_tri.size(); ++c) for (int n = 0; n < 3; ++n)
{ {
tmp_tri = cnst_tri[c]; make_delaunay(cnst_tri[n]);
delete tmp_tri; tmp_tri = nullptr;
} }
} }
// remove any triangles has an assistant vertex from exist_tri // remove any triangles has an box vertex from out_tris
for (t_iter = exist_tri.begin(); t_iter != exist_tri.end(); ) for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
{ {
tmp_tri = *t_iter; tmp_tri = *t_iter;
if (tmp_tri->vert[0] == assit_vert[0] || tmp_tri->vert[0] == assit_vert[1] || tmp_tri->vert[0] == assit_vert[2] || tmp_tri->vert[0] == assit_vert[3] || if (tmp_tri->vert[0] == box_vert[0] || tmp_tri->vert[0] == box_vert[1] || tmp_tri->vert[0] == box_vert[2] || tmp_tri->vert[0] == box_vert[3])
tmp_tri->vert[1] == assit_vert[0] || tmp_tri->vert[1] == assit_vert[1] || tmp_tri->vert[1] == assit_vert[2] || tmp_tri->vert[1] == assit_vert[3] ||
tmp_tri->vert[2] == assit_vert[0] || tmp_tri->vert[2] == assit_vert[1] || tmp_tri->vert[2] == assit_vert[2] || tmp_tri->vert[2] == assit_vert[3])
{ {
if (tmp_tri->neigh[1] != nullptr)
{
for (int k = 0; k < 3; ++k)
{
if (tmp_tri->neigh[1]->neigh[k] == tmp_tri)
{
tmp_tri->neigh[1]->neigh[k] = nullptr;
break;
}
}
}
// destroy the memories located and remove from the vector // destroy the memories located and remove from the vector
t_iter = exist_tri.erase(t_iter); t_iter = out_tris.erase(t_iter);
delete tmp_tri; tmp_tri = nullptr;
}
else if (tmp_tri->vert[1] == box_vert[0] || tmp_tri->vert[1] == box_vert[1] || tmp_tri->vert[1] == box_vert[2] || tmp_tri->vert[1] == box_vert[3])
{
if (tmp_tri->neigh[2] != nullptr)
{
for (int k = 0; k < 3; ++k)
{
if (tmp_tri->neigh[2]->neigh[k] == tmp_tri)
{
tmp_tri->neigh[2]->neigh[k] = nullptr;
break;
}
}
}
// destroy the memories located and remove from the vector
t_iter = out_tris.erase(t_iter);
delete tmp_tri; tmp_tri = nullptr;
}
else if (tmp_tri->vert[2] == box_vert[0] || tmp_tri->vert[2] == box_vert[1] || tmp_tri->vert[2] == box_vert[2] || tmp_tri->vert[2] == box_vert[3])
{
if (tmp_tri->neigh[0] != nullptr)
{
for (int k = 0; k < 3; ++k)
{
if (tmp_tri->neigh[0]->neigh[k] == tmp_tri)
{
tmp_tri->neigh[0]->neigh[k] = nullptr;
break;
}
}
}
// destroy the memories located and remove from the vector
t_iter = out_tris.erase(t_iter);
delete tmp_tri; tmp_tri = nullptr; delete tmp_tri; tmp_tri = nullptr;
} }
else t_iter++; else t_iter++;
} }
// copy exist_tri to out_tris and destroy memories located // assign triangles index
out_tris.resize(exist_tri.size()); for (int i = 0; i < out_tris.size(); i++)
for (int i = 0; i < exist_tri.size(); ++i)
{ {
out_tris[i].set(exist_tri[i]->vert[0], exist_tri[i]->vert[1], exist_tri[i]->vert[2]); out_tris[i]->id = i;
delete exist_tri[i]; exist_tri[i] = nullptr;
} }
// destroy memories located for assit_vert // destroy memories located for box_vert
for (int i = 0; i < 4; ++i) for (int i = 0; i < 4; ++i)
{ {
delete assit_vert[i]; assit_vert[i] = nullptr; delete box_vert[i]; box_vert[i] = nullptr;
} }
return; return;
} }
@@ -271,9 +426,9 @@ bool duplicated_vertex(const std::vector<vertex2dc> &in_verts)
* *
* @return If the triangulation is fully delaunay * @return If the triangulation is fully delaunay
*/ */
bool fully_delaunay(const std::vector<triangle> &in_tris, const std::vector<vertex2dc> &in_verts) bool fully_delaunay(const std::vector<triangle*> &in_tris, const std::vector<vertex2dc> &in_verts)
{ {
if (in_tris.empty()) return true; if (in_tris.empty()) return false;
int count; int count;
double dist; double dist;
@@ -282,10 +437,10 @@ bool fully_delaunay(const std::vector<triangle> &in_tris, const std::vector<vert
count = 0; count = 0;
for (int j = 0; j < in_verts.size(); ++j) for (int j = 0; j < in_verts.size(); ++j)
{ {
dist = (in_tris[i].cx - in_verts[j].x) * (in_tris[i].cx - in_verts[j].x) + dist = (in_tris[i]->cx - in_verts[j].x) * (in_tris[i]->cx - in_verts[j].x) +
(in_tris[i].cy - in_verts[j].y) * (in_tris[i].cy - in_verts[j].y); (in_tris[i]->cy - in_verts[j].y) * (in_tris[i]->cy - in_verts[j].y);
if ((dist - in_tris[i].cr) <= ZERO) if ((dist - in_tris[i]->cr) <= ZERO)
{ {
count++; count++;
} }
@@ -300,4 +455,4 @@ bool fully_delaunay(const std::vector<triangle> &in_tris, const std::vector<vert
return true; return true;
} }
#endif // _BW_2D_DELAUNAY_H #endif // _FLIP_2D_DELAUNAY_H

337
delaunay_backup.h Normal file
View File

@@ -0,0 +1,337 @@
/**
* @defgroup DELAUNAY
*
* @brief An implementation of the 2D Delaunay triangulation using the Bowyer-Watson algorithm.
*
* @author Yi Zhang
* @date 2021-09-12
*/
#ifndef _BW_2D_DELAUNAY_H
#define _BW_2D_DELAUNAY_H
#include "cmath"
#include "vector"
#include "iostream"
#include "fstream"
#include "iomanip"
#define ZERO 1e-5
// Start vertex definition
struct vertex2dc
{
unsigned int id; // index of the vertex
double x, y; // position of the vertex
vertex2dc() : x(NAN), y(NAN), id(0) {}
vertex2dc(double inx, double iny, unsigned int inid = 0) {set(inx, iny, inid);}
void set(double inx, double iny, unsigned int inid = 0)
{
x = inx; y = iny; id = inid;
return;
}
};
bool operator ==(const vertex2dc &a, const vertex2dc &b) // overload the == operator for vertex2dc type
{
if(fabs(a.x - b.x) <= ZERO && fabs(a.y - b.y) <= ZERO)
{
return true;
}
return false;
}
// End vertex definition
// Start edge definition
struct edge
{
vertex2dc *vert[2]; // vertex of the edge
edge() {vert[0] = vert[1] = nullptr;}
edge(vertex2dc *v0ptr, vertex2dc *v1ptr) {set(v0ptr, v1ptr);}
void set(vertex2dc *v0ptr, vertex2dc *v1ptr)
{
vert[0] = v0ptr; vert[1] = v1ptr;
return;
}
};
bool operator ==(const edge &a, const edge &b) // overload the == operator for edge type
{
if((a.vert[0] == b.vert[0] && a.vert[1] == b.vert[1]) ||
(a.vert[0] == b.vert[1] && a.vert[1] == b.vert[0]))
{
return true;
}
return false;
}
// End edge definition
// Start triangle definition
struct triangle
{
vertex2dc *vert[3]; // vertex of the triangle
double cx, cy; // center of the triangle's circumcircle
double cr; // radius of the circumcircle
triangle() {vert[0] = vert[1] = vert[2] = nullptr;}
triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) {set(v0ptr, v1ptr, v2ptr);}
void set(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr)
{
vert[0] = v0ptr; vert[1] = v1ptr; vert[2] = v2ptr;
double s = 0.5 / ((vert[1]->x - vert[0]->x) * (vert[2]->y - vert[0]->y) - (vert[1]->y - vert[0]->y) * (vert[2]->x - vert[0]->x));
double m = vert[1]->x * vert[1]->x - vert[0]->x * vert[0]->x + vert[1]->y * vert[1]->y - vert[0]->y * vert[0]->y;
double u = vert[2]->x * vert[2]->x - vert[0]->x * vert[0]->x + vert[2]->y * vert[2]->y - vert[0]->y * vert[0]->y;
cx = ((vert[2]->y - vert[0]->y) * m + (vert[0]->y - vert[1]->y) * u) * s;
cy = ((vert[0]->x - vert[2]->x) * m + (vert[1]->x - vert[0]->x) * u) * s;
cr = (vert[0]->x - cx) * (vert[0]->x - cx) + (vert[0]->y - cy) * (vert[0]->y - cy); // not need to sqrt() here
return;
}
};
// End triangle definition
/**
* @brief 2D Delaunay triangulation of some given points using the Bowyer-Watson algorithm.
*
* @param in_verts Input vertexes. Defined by the user.
* @param out_tris Output triangles. Compute by the function.
*/
void triangulation(std::vector<vertex2dc> &in_verts, std::vector<triangle> &out_tris, int valid_size)
{
if (!out_tris.empty()) out_tris.clear();
if (in_verts.size() < 3) return;
// locate the surrounding box and initiate the staring two triangles
double xmin = in_verts[0].x, xmax = in_verts[0].x;
double ymin = in_verts[0].y, ymax = in_verts[0].y;
for (int i = 0; i < in_verts.size(); ++i)
{
xmin = std::min(xmin, in_verts[i].x);
xmax = std::max(xmax, in_verts[i].x);
ymin = std::min(ymin, in_verts[i].y);
ymax = std::max(ymax, in_verts[i].y);
}
double midx = 0.5*(xmin + xmax);
double midy = 0.5*(ymin + ymax);
double maxi_s = std::max(xmax - xmin, ymax - ymin); // use an four times bigger rectangle to include all points
vertex2dc *tmp_vert = nullptr;
std::vector<vertex2dc*> assit_vert;
tmp_vert = new vertex2dc(midx - maxi_s, midy - maxi_s, 0); // lower left corner
assit_vert.push_back(tmp_vert);
tmp_vert = new vertex2dc(midx + maxi_s, midy - maxi_s, 1); // lower right corner
assit_vert.push_back(tmp_vert);
tmp_vert = new vertex2dc(midx + maxi_s, midy + maxi_s, 2); // upper right corner
assit_vert.push_back(tmp_vert);
tmp_vert = new vertex2dc(midx - maxi_s, midy + maxi_s, 3); // upper left corner
assit_vert.push_back(tmp_vert);
triangle *tmp_tri = nullptr;
std::vector<triangle*> exist_tri, cnst_tri;
std::vector<triangle*>::iterator t_iter;
tmp_tri = new triangle(assit_vert[0], assit_vert[1], assit_vert[2]); // order the vertex anti-clock wise
exist_tri.push_back(tmp_tri);
tmp_tri = new triangle(assit_vert[0], assit_vert[2], assit_vert[3]); // order the vertex anti-clock wise
exist_tri.push_back(tmp_tri);
// loop all input vertice
bool removed;
double dist;
edge tmp_edge;
std::vector<edge> cnst_edge;
std::vector<edge>::iterator e_iter;
//for (int i = 0; i < in_verts.size(); ++i)
for (int i = 0; i < valid_size; ++i)
{
// determine triangles that include the point and add the triangle to the cnst_tri and remove the triangle from exist_tri
// this is the part that could take a lot of time if we are working with a large amount of points. We will fix this later
cnst_tri.clear();
for (t_iter = exist_tri.begin(); t_iter != exist_tri.end(); )
{
tmp_tri = *t_iter;
dist = (tmp_tri->cx - in_verts[i].x) * (tmp_tri->cx - in_verts[i].x) +
(tmp_tri->cy - in_verts[i].y) * (tmp_tri->cy - in_verts[i].y);
if ((dist - tmp_tri->cr) <= ZERO) // Points on the circumcircle are also included
{
t_iter = exist_tri.erase(t_iter);
cnst_tri.push_back(tmp_tri);
}
else t_iter++;
}
// loop to remove duplicate edges
cnst_edge.clear();
for (int c = 0; c < cnst_tri.size(); ++c)
{
for (int e = 0; e < 3; ++e)
{
tmp_edge.set(cnst_tri[c]->vert[e], cnst_tri[c]->vert[(e+1)%3]);
removed = false;
for (e_iter = cnst_edge.begin(); e_iter != cnst_edge.end(); )
{
if (tmp_edge == *e_iter) // duplicate edge, remove from cnst_edge
{
e_iter = cnst_edge.erase(e_iter);
removed = true;
break; // no need to search more
}
else e_iter++;
}
if (!removed) // not a duplicate edge, add to the cnst_edge
{
cnst_edge.push_back(tmp_edge);
}
}
}
// construct new triangles and add to exist_tri
for (int c = 0; c < cnst_edge.size(); ++c)
{
tmp_tri = new triangle(cnst_edge[c].vert[0], cnst_edge[c].vert[1], &in_verts[i]); // order the vertex anti-clock wise
exist_tri.push_back(tmp_tri);
}
// destroy memories used by cnst_edge
for (int c = 0; c < cnst_tri.size(); ++c)
{
tmp_tri = cnst_tri[c];
delete tmp_tri; tmp_tri = nullptr;
}
}
// Start testing code
std::ofstream outfile("test_backup.msh");
outfile << "$MeshFormat" << std::endl << "2.2 0 8" << std::endl << "$EndMeshFormat "<<std::endl;
outfile << "$Nodes" << std::endl << valid_size + 4 << std::endl;
for (int i = 0; i < 4; i++)
{
outfile << i + 1 << " " << std::setprecision(16)
<< assit_vert[i]->x << " " << assit_vert[i]->y << " 0.0" << std::endl;
}
for (int i = 0; i < valid_size; i++)
{
outfile << i + 5 << " " << std::setprecision(16)
<< in_verts[i].x << " " << in_verts[i].y << " 0.0" << std::endl;
}
outfile<<"$EndNodes"<<std::endl;
outfile << "$Elements" << std::endl << exist_tri.size() <<std::endl;
for (int i = 0; i < exist_tri.size(); i++)
{
outfile << i + 1 << " 2 0";
for (int j = 0; j < 3; j++)
{
outfile << " " << exist_tri[i]->vert[j]->id + 1;
}
outfile << std::endl;
}
outfile << "$EndElements"<< std::endl;
outfile.close();
// End testing code
// remove any triangles has an assistant vertex from exist_tri
for (t_iter = exist_tri.begin(); t_iter != exist_tri.end(); )
{
tmp_tri = *t_iter;
if (tmp_tri->vert[0] == assit_vert[0] || tmp_tri->vert[0] == assit_vert[1] || tmp_tri->vert[0] == assit_vert[2] || tmp_tri->vert[0] == assit_vert[3] ||
tmp_tri->vert[1] == assit_vert[0] || tmp_tri->vert[1] == assit_vert[1] || tmp_tri->vert[1] == assit_vert[2] || tmp_tri->vert[1] == assit_vert[3] ||
tmp_tri->vert[2] == assit_vert[0] || tmp_tri->vert[2] == assit_vert[1] || tmp_tri->vert[2] == assit_vert[2] || tmp_tri->vert[2] == assit_vert[3])
{
// destroy the memories located and remove from the vector
t_iter = exist_tri.erase(t_iter);
delete tmp_tri; tmp_tri = nullptr;
}
else t_iter++;
}
// copy exist_tri to out_tris and destroy memories located
out_tris.resize(exist_tri.size());
for (int i = 0; i < exist_tri.size(); ++i)
{
out_tris[i].set(exist_tri[i]->vert[0], exist_tri[i]->vert[1], exist_tri[i]->vert[2]);
delete exist_tri[i]; exist_tri[i] = nullptr;
}
// destroy memories located for assit_vert
for (int i = 0; i < 4; ++i)
{
delete assit_vert[i]; assit_vert[i] = nullptr;
}
return;
}
/**
* @brief Check for duplicated vertex
*
* @param[in] in_verts Input vertexes
*
* @return If there is duplicated vertex
*/
bool duplicated_vertex(const std::vector<vertex2dc> &in_verts)
{
if (in_verts.empty()) return false;
for (int i = 0; i < in_verts.size()-1; ++i)
{
for (int j = i+1; j < in_verts.size(); ++j)
{
if (in_verts[i] == in_verts[j] && in_verts[i].id != in_verts[j].id)
{
return true;
}
}
}
return false;
}
/**
* @brief Check to see if the triangulation is fully delaunay
*
* @param[in] in_tris Input triangles
* @param[in] in_verts Input vertexes
*
* @return If the triangulation is fully delaunay
*/
bool fully_delaunay(const std::vector<triangle> &in_tris, const std::vector<vertex2dc> &in_verts)
{
if (in_tris.empty()) return true;
int count;
double dist;
for (int i = 0; i < in_tris.size(); ++i)
{
count = 0;
for (int j = 0; j < in_verts.size(); ++j)
{
dist = (in_tris[i].cx - in_verts[j].x) * (in_tris[i].cx - in_verts[j].x) +
(in_tris[i].cy - in_verts[j].y) * (in_tris[i].cy - in_verts[j].y);
if ((dist - in_tris[i].cr) <= ZERO)
{
count++;
}
}
if (count > 3)
{
return false;
}
}
return true;
}
#endif // _BW_2D_DELAUNAY_H

58
demo.cpp
View File

@@ -5,6 +5,7 @@
int main(int argc, char const *argv[]) int main(int argc, char const *argv[])
{ {
std::vector<vertex2dc> points(21); std::vector<vertex2dc> points(21);
points[0].set(-0.8, -0.8, 0); points[0].set(-0.8, -0.8, 0);
points[1].set(0.4, -1.2, 1); points[1].set(0.4, -1.2, 1);
@@ -28,13 +29,38 @@ int main(int argc, char const *argv[])
points[19].set(3.5, 1.8, 19); points[19].set(3.5, 1.8, 19);
points[20].set(3.6, 3.1, 20); points[20].set(3.6, 3.1, 20);
/*
std::vector<vertex2dc> points(21);
points[0].set(-0.8, -0.8, 4);
points[1].set(0.4, -1.2, 5);
points[2].set(1.2, -0.9, 6);
points[3].set(1.6, 0.1, 7);
points[4].set(2.5, 0.5, 8);
points[5].set(4.1, 0.7, 9);
points[6].set(5.7, 1.8, 10);
points[7].set(5.1, 3.4, 11);
points[8].set(2.5, 4.4, 12);
points[9].set(1.2, 3.7, 13);
points[10].set(-1.2, 3.9, 14);
points[11].set(-3.2, 5.1, 15);
points[12].set(-4.3, 2.9, 16);
points[13].set(-3.1, 0.7, 17);
points[14].set(-1.3, 0.6, 18);
points[15].set(-2.1, 2.9, 19);
points[16].set(0.6, 1.2, 20);
points[17].set(0.1, 2.4, 21);
points[18].set(2.4, 2.8, 22);
points[19].set(3.5, 1.8, 23);
points[20].set(3.6, 3.1, 24);
*/
if (duplicated_vertex(points)) if (duplicated_vertex(points))
{ {
std::cerr << "Duplicated vertexes detected.\n"; std::cerr << "Duplicated vertice detected.\n";
return 0; return 0;
} }
std::vector<triangle> elements; std::vector<triangle*> elements;
triangulation(points, elements); triangulation(points, elements);
if (fully_delaunay(elements, points)) if (fully_delaunay(elements, points))
@@ -59,13 +85,35 @@ int main(int argc, char const *argv[])
outfile << i + 1 << " 2 0"; outfile << i + 1 << " 2 0";
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
{ {
outfile << " " << elements[i].vert[j]->id + 1; outfile << " " << elements[i]->vert[j]->id + 1;
} }
outfile << std::endl; outfile << std::endl;
} }
outfile << "$EndElements"<< std::endl; outfile << "$EndElements"<< std::endl;
outfile.close(); outfile.close();
// write a neighbor file
outfile.open("demo.neigh");
outfile << elements.size() << std::endl;
for (int i = 0; i < elements.size(); i++)
{
outfile << i + 1;
for (int j = 0; j < 3; j++)
{
if (elements[i]->neigh[j] != nullptr)
{
outfile << " " << elements[i]->neigh[j]->id + 1;
}
else outfile << " -1";
}
outfile << std::endl;
}
outfile.close();
// destroy allocated memories
for (int i = 0; i < elements.size(); i++)
{
delete elements[i];
}
return 0; return 0;
} }

58
demo.msh
View File

@@ -27,34 +27,34 @@ $Nodes
$EndNodes $EndNodes
$Elements $Elements
30 30
1 2 0 2 3 4 1 2 0 12 13 16
2 2 0 4 3 5 2 2 0 17 5 19
3 2 0 5 3 6 3 2 0 3 4 2
4 2 0 10 9 11 4 2 0 17 2 4
5 2 0 11 9 12 5 2 0 4 3 5
6 2 0 14 1 15 6 2 0 7 8 20
7 2 0 11 12 16 7 2 0 5 3 6
8 2 0 12 13 16 8 2 0 18 16 15
9 2 0 13 14 16 9 2 0 10 9 11
10 2 0 14 15 16 10 2 0 11 9 12
11 2 0 1 2 17 11 2 0 1 2 17
12 2 0 2 4 17 12 2 0 17 18 15
13 2 0 4 5 17 13 2 0 14 1 15
14 2 0 15 1 17 14 2 0 13 14 16
15 2 0 10 11 18 15 2 0 14 15 16
16 2 0 11 16 18 16 2 0 11 12 16
17 2 0 16 15 18 17 2 0 15 1 17
18 2 0 15 17 18 18 2 0 4 5 17
19 2 0 9 10 19 19 2 0 10 11 18
20 2 0 17 5 19 20 2 0 11 16 18
21 2 0 10 18 19 21 2 0 17 19 18
22 2 0 18 17 19 22 2 0 10 18 19
23 2 0 6 7 20 23 2 0 5 6 20
24 2 0 7 8 20 24 2 0 9 10 19
25 2 0 5 6 20 25 2 0 19 5 20
26 2 0 19 5 20 26 2 0 6 7 20
27 2 0 8 9 21 27 2 0 8 21 20
28 2 0 9 19 21 28 2 0 19 20 21
29 2 0 19 20 21 29 2 0 8 9 21
30 2 0 20 8 21 30 2 0 9 19 21
$EndElements $EndElements

31
demo.neigh Normal file
View File

@@ -0,0 +1,31 @@
30
1 -1 14 16
2 18 25 21
3 5 4 -1
4 11 3 18
5 3 7 18
6 -1 27 26
7 5 -1 23
8 20 15 12
9 24 10 19
10 9 -1 16
11 -1 4 17
12 21 8 17
13 -1 17 15
14 -1 15 1
15 13 8 14
16 10 1 20
17 13 11 12
18 5 2 4
19 9 20 22
20 16 8 19
21 2 22 12
22 19 21 24
23 7 26 25
24 9 22 30
25 2 23 28
26 -1 6 23
27 29 28 6
28 25 27 30
29 -1 30 27
30 24 28 29

96
demo_backup.cpp Normal file
View File

@@ -0,0 +1,96 @@
#include "delaunay_backup.h"
#include "iostream"
#include "fstream"
#include "iomanip"
int main(int argc, char const *argv[])
{
/*
std::vector<vertex2dc> points(21);
points[0].set(-0.8, -0.8, 0);
points[1].set(0.4, -1.2, 1);
points[2].set(1.2, -0.9, 2);
points[3].set(1.6, 0.1, 3);
points[4].set(2.5, 0.5, 4);
points[5].set(4.1, 0.7, 5);
points[6].set(5.7, 1.8, 6);
points[7].set(5.1, 3.4, 7);
points[8].set(2.5, 4.4, 8);
points[9].set(1.2, 3.7, 9);
points[10].set(-1.2, 3.9, 10);
points[11].set(-3.2, 5.1, 11);
points[12].set(-4.3, 2.9, 12);
points[13].set(-3.1, 0.7, 13);
points[14].set(-1.3, 0.6, 14);
points[15].set(-2.1, 2.9, 15);
points[16].set(0.6, 1.2, 16);
points[17].set(0.1, 2.4, 17);
points[18].set(2.4, 2.8, 18);
points[19].set(3.5, 1.8, 19);
points[20].set(3.6, 3.1, 20);
*/
std::vector<vertex2dc> points(21);
points[0].set(-0.8, -0.8, 4);
points[1].set(0.4, -1.2, 5);
points[2].set(1.2, -0.9, 6);
points[3].set(1.6, 0.1, 7);
points[4].set(2.5, 0.5, 8);
points[5].set(4.1, 0.7, 9);
points[6].set(5.7, 1.8, 10);
points[7].set(5.1, 3.4, 11);
points[8].set(2.5, 4.4, 12);
points[9].set(1.2, 3.7, 13);
points[10].set(-1.2, 3.9, 14);
points[11].set(-3.2, 5.1, 15);
points[12].set(-4.3, 2.9, 16);
points[13].set(-3.1, 0.7, 17);
points[14].set(-1.3, 0.6, 18);
points[15].set(-2.1, 2.9, 19);
points[16].set(0.6, 1.2, 20);
points[17].set(0.1, 2.4, 21);
points[18].set(2.4, 2.8, 22);
points[19].set(3.5, 1.8, 23);
points[20].set(3.6, 3.1, 24);
if (duplicated_vertex(points))
{
std::cerr << "Duplicated vertexes detected.\n";
return 0;
}
std::vector<triangle> elements;
triangulation(points, elements, atoi(argv[1]));
if (fully_delaunay(elements, points))
{
std::clog << "The triangulation is fully delaunay.\n";
}
else std::clog << "The triangulation is not fully delaunay\n";
// Write a Gmsh's .msh file
std::ofstream outfile("demo.msh");
outfile << "$MeshFormat" << std::endl << "2.2 0 8" << std::endl << "$EndMeshFormat "<<std::endl;
outfile << "$Nodes" << std::endl << points.size() << std::endl;
for (int i = 0; i < points.size(); i++)
{
outfile << points[i].id + 1 << " " << std::setprecision(16)
<< points[i].x << " " << points[i].y << " 0.0" << std::endl;
}
outfile<<"$EndNodes"<<std::endl;
outfile << "$Elements" << std::endl << elements.size() <<std::endl;
for (int i = 0; i < elements.size(); i++)
{
outfile << i + 1 << " 2 0";
for (int j = 0; j < 3; j++)
{
outfile << " " << elements[i].vert[j]->id + 1;
}
outfile << std::endl;
}
outfile << "$EndElements"<< std::endl;
outfile.close();
return 0;
}

78
test.msh Normal file
View File

@@ -0,0 +1,78 @@
$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
25
1 -9.300000000000001 -8.050000000000001 0.0
2 10.7 -8.050000000000001 0.0
3 10.7 11.95 0.0
4 -9.300000000000001 11.95 0.0
5 -0.8 -0.8 0.0
6 0.4 -1.2 0.0
7 1.2 -0.9 0.0
8 1.6 0.1 0.0
9 2.5 0.5 0.0
10 4.1 0.7 0.0
11 5.7 1.8 0.0
12 5.1 3.4 0.0
13 2.5 4.4 0.0
14 1.2 3.7 0.0
15 -1.2 3.9 0.0
16 -3.2 5.1 0.0
17 -4.3 2.9 0.0
18 -3.1 0.7 0.0
19 -1.3 0.6 0.0
20 -2.1 2.9 0.0
21 0.6 1.2 0.0
22 0.1 2.4 0.0
23 2.4 2.8 0.0
24 3.5 1.8 0.0
25 3.6 3.1 0.0
$EndNodes
$Elements
44
1 2 0 1 2 6
2 2 0 7 8 6
3 2 0 5 1 6
4 2 0 6 2 7
5 2 0 5 6 8
6 2 0 7 2 10
7 2 0 8 7 9
8 2 0 11 12 10
9 2 0 9 7 10
10 2 0 2 3 11
11 2 0 10 2 11
12 2 0 3 4 13
13 2 0 11 3 12
14 2 0 12 3 13
15 2 0 8 9 23
16 2 0 13 4 16
17 2 0 14 13 15
18 2 0 17 18 16
19 2 0 15 13 16
20 2 0 4 1 17
21 2 0 16 4 17
22 2 0 1 5 18
23 2 0 5 19 18
24 2 0 17 1 18
25 2 0 20 18 19
26 2 0 5 8 21
27 2 0 15 16 20
28 2 0 16 18 20
29 2 0 19 21 20
30 2 0 19 5 21
31 2 0 8 14 21
32 2 0 22 20 21
33 2 0 21 14 22
34 2 0 14 15 22
35 2 0 15 20 22
36 2 0 9 10 24
37 2 0 13 14 23
38 2 0 14 8 23
39 2 0 23 9 24
40 2 0 10 12 24
41 2 0 25 23 24
42 2 0 13 23 25
43 2 0 24 12 25
44 2 0 12 13 25
$EndElements

45
test_backup.msh Normal file
View File

@@ -0,0 +1,45 @@
$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
14
1 -9.300000000000001 -8.050000000000001 0.0
2 10.7 -8.050000000000001 0.0
3 10.7 11.95 0.0
4 -9.300000000000001 11.95 0.0
5 -0.8 -0.8 0.0
6 0.4 -1.2 0.0
7 1.2 -0.9 0.0
8 1.6 0.1 0.0
9 2.5 0.5 0.0
10 4.1 0.7 0.0
11 5.7 1.8 0.0
12 5.1 3.4 0.0
13 2.5 4.4 0.0
14 1.2 3.7 0.0
$EndNodes
$Elements
22
1 2 0 4 1 5
2 2 0 1 2 6
3 2 0 5 1 6
4 2 0 6 2 7
5 2 0 5 6 8
6 2 0 6 7 8
7 2 0 8 7 9
8 2 0 7 2 10
9 2 0 9 7 10
10 2 0 2 3 11
11 2 0 10 2 11
12 2 0 11 3 12
13 2 0 10 11 12
14 2 0 3 4 13
15 2 0 12 3 13
16 2 0 9 10 13
17 2 0 10 12 13
18 2 0 4 5 14
19 2 0 13 4 14
20 2 0 5 8 14
21 2 0 8 9 14
22 2 0 9 13 14
$EndElements