458 lines
12 KiB
C++
458 lines
12 KiB
C++
/**
|
|
* @defgroup DELAUNAY
|
|
*
|
|
* @brief An implementation of the 2D Delaunay triangulation using the flip algorithm.
|
|
*
|
|
* @author Yi Zhang
|
|
* @date 2021-09-19
|
|
*/
|
|
|
|
#ifndef _FLIP_2D_DELAUNAY_H
|
|
#define _FLIP_2D_DELAUNAY_H
|
|
#include "cmath"
|
|
#include "vector"
|
|
|
|
#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;
|
|
}
|
|
|
|
void circumcircle(vertex2dc *v0, vertex2dc *v1, vertex2dc *v2, double &cx, double &cy, double &cr) // calculate the circumcircle from three points
|
|
{
|
|
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;
|
|
|
|
cx = ((v2->y - v0->y)*m + (v0->y - v1->y)*u)*s;
|
|
cy = ((v0->x - v2->x)*m + (v1->x - v0->x)*u)*s;
|
|
cr = (v0->x - cx)*(v0->x - cx) + (v0->y - cy)*(v0->y - cy); // not need to calculate the squared root here
|
|
return;
|
|
}
|
|
// End vertex definition
|
|
|
|
// Start triangle definition
|
|
struct triangle
|
|
{
|
|
int id; // index 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 cr; // radius of the circumcircle
|
|
|
|
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);}
|
|
void set(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr)
|
|
{
|
|
vert[0] = v0ptr; vert[1] = v1ptr; vert[2] = v2ptr;
|
|
neigh[0] = neigh[1] = neigh[2] = nullptr;
|
|
circumcircle(vert[0], vert[1], vert[2], cx, cy, cr);
|
|
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 meets 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) < -1.0*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
|
|
|
|
/**
|
|
* @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)
|
|
{
|
|
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*> box_vert;
|
|
|
|
tmp_vert = new vertex2dc(midx - maxi_s, midy - maxi_s); // lower left corner
|
|
box_vert.push_back(tmp_vert);
|
|
|
|
tmp_vert = new vertex2dc(midx + maxi_s, midy - maxi_s); // lower right corner
|
|
box_vert.push_back(tmp_vert);
|
|
|
|
tmp_vert = new vertex2dc(midx + maxi_s, midy + maxi_s); // upper right corner
|
|
box_vert.push_back(tmp_vert);
|
|
|
|
tmp_vert = new vertex2dc(midx - maxi_s, midy + maxi_s); // upper left corner
|
|
box_vert.push_back(tmp_vert);
|
|
|
|
triangle *old_tri = nullptr, *tmp_tri = nullptr;
|
|
triangle *cnst_tri[3];
|
|
std::vector<triangle*>::iterator t_iter;
|
|
|
|
tmp_tri = new triangle(box_vert[0], box_vert[1], box_vert[2]); // order the vertex anti-clock wise
|
|
out_tris.push_back(tmp_tri);
|
|
|
|
tmp_tri = new triangle(box_vert[0], box_vert[2], box_vert[3]); // order the vertex anti-clock wise
|
|
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
|
|
for (int i = 0; i < in_verts.size(); ++i)
|
|
{
|
|
// determine the triangle that includes the new vertex and remove it from out_tris
|
|
for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
|
|
{
|
|
old_tri = *t_iter;
|
|
if (old_tri->bound_location(in_verts[i].x, in_verts[i].y))
|
|
{
|
|
t_iter = out_tris.erase(t_iter);
|
|
break;
|
|
}
|
|
else t_iter++;
|
|
}
|
|
|
|
// build three new triangles
|
|
for (int n = 0; n < 3; ++n)
|
|
{
|
|
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)
|
|
{
|
|
cnst_tri[n]->set_neighbor(nullptr, cnst_tri[(n+1)%3], cnst_tri[(n+2)%3]);
|
|
}
|
|
else
|
|
{
|
|
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 (old_tri->neigh[n]->neigh[k] == old_tri)
|
|
{
|
|
old_tri->neigh[n]->neigh[k] = cnst_tri[n];
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// delete the old triangle
|
|
delete old_tri; old_tri = nullptr;
|
|
|
|
// Make sure cnst_tri meet the empty circumcircle condition
|
|
for (int n = 0; n < 3; ++n)
|
|
{
|
|
make_delaunay(cnst_tri[n]);
|
|
}
|
|
}
|
|
|
|
// remove any triangles has an box vertex from out_tris
|
|
for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
|
|
{
|
|
tmp_tri = *t_iter;
|
|
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])
|
|
{
|
|
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
|
|
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;
|
|
}
|
|
else t_iter++;
|
|
}
|
|
|
|
// assign triangles index
|
|
for (int i = 0; i < out_tris.size(); i++)
|
|
{
|
|
out_tris[i]->id = i;
|
|
}
|
|
|
|
// destroy memories located for box_vert
|
|
for (int i = 0; i < 4; ++i)
|
|
{
|
|
delete box_vert[i]; box_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 false;
|
|
|
|
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 // _FLIP_2D_DELAUNAY_H
|