242 lines
7.4 KiB
C++
242 lines
7.4 KiB
C++
/**
|
|
* @defgroup TIN
|
|
*
|
|
* @brief Generation of a Triangular Irregular Network (TIN) from a dense DEM grid
|
|
*
|
|
* @author Yi Zhang
|
|
* @date 2021-09-15
|
|
*/
|
|
|
|
#ifndef _TIN_DELAUNAY_H
|
|
#define _TIN_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;
|
|
}
|
|
// 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)
|
|
{
|
|
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); // lower left corner
|
|
assit_vert.push_back(tmp_vert);
|
|
|
|
tmp_vert = new vertex2dc(midx + maxi_s, midy - maxi_s); // lower right corner
|
|
assit_vert.push_back(tmp_vert);
|
|
|
|
tmp_vert = new vertex2dc(midx + maxi_s, midy + maxi_s); // upper right corner
|
|
assit_vert.push_back(tmp_vert);
|
|
|
|
tmp_vert = new vertex2dc(midx - maxi_s, midy + maxi_s); // 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)
|
|
{
|
|
// 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;
|
|
}
|
|
}
|
|
|
|
// 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;
|
|
}
|
|
|
|
#endif // _TIN_DELAUNAY_H
|