406 lines
13 KiB
C++
406 lines
13 KiB
C++
/**
|
||
* @defgroup TIN
|
||
*
|
||
* @brief Generation of a Triangular Irregular Network (TIN) from a dense DEM grid.
|
||
*
|
||
* @author Yi Zhang
|
||
* @date 2021-09-16
|
||
*/
|
||
|
||
#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
|
||
double elev; // elevation at the vertex
|
||
|
||
vertex2dc() : x(NAN), y(NAN), elev(NAN), id(0) {}
|
||
vertex2dc(double inx, double iny, double inelev, unsigned int inid = 0) {set(inx, iny, inelev, inid);}
|
||
void set(double inx, double iny, double inelev, unsigned int inid = 0)
|
||
{
|
||
x = inx; y = iny; elev = inelev; 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;
|
||
}
|
||
|
||
bool is_collinear(vertex2dc *a_ptr, vertex2dc *b_ptr, vertex2dc *c_ptr) // Test if the three points are on the same line
|
||
{
|
||
// ｜(y3−y1)(x2−x1)−(y2−y1)(x3−x1)｜
|
||
if (fabs((c_ptr->y - a_ptr->y)*(b_ptr->x - a_ptr->x) - (b_ptr->y - a_ptr->y)*(c_ptr->x - a_ptr->x)) <= 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;
|
||
}
|
||
|
||
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;
|
||
}
|
||
|
||
double interpolate(double inx, double iny) // Interpolate the elevation of the given location inside the triangle
|
||
{
|
||
double a1 = 0.5 * ((vert[1]->x - inx)*(vert[2]->y - iny) - (vert[1]->y - iny)*(vert[2]->x - inx));
|
||
double a2 = 0.5 * ((vert[2]->x - inx)*(vert[0]->y - iny) - (vert[2]->y - iny)*(vert[0]->x - inx));
|
||
double a3 = 0.5 * ((vert[0]->x - inx)*(vert[1]->y - iny) - (vert[0]->y - iny)*(vert[1]->x - inx));
|
||
return (a1*vert[0]->elev + a2*vert[1]->elev + a3*vert[2]->elev)/(a1 + a2 + a3);
|
||
}
|
||
};
|
||
// End triangle definition
|
||
|
||
// Start DEM definition
|
||
struct dem_point
|
||
{
|
||
double x, y; // position of the DEM location
|
||
double elev; // elevation at the DEM location
|
||
triangle *host; // host triangle of the DEM location
|
||
std::vector<triangle*> circum_host; // triangles which circumcircles include the location
|
||
|
||
dem_point() : x(NAN), y(NAN), elev(NAN), host(nullptr) {}
|
||
dem_point(double inx, double iny, double inelev) {set(inx, iny, inelev);}
|
||
void set(double inx, double iny, double inelev)
|
||
{
|
||
x = inx; y = iny; elev = inelev; host = nullptr;
|
||
circum_host.clear();
|
||
return;
|
||
}
|
||
};
|
||
// End DEM definition
|
||
|
||
/**
|
||
* @brief Generate the TIN from the DEM grid
|
||
*
|
||
* @param[in] dem Input DEM grid (Ordered from lower left corner to the upper right corner)
|
||
* @param[in] xmin The minimal coordinate of the DEM grid on the x-axis
|
||
* @param[in] xmax The maximal coordinate of the DEM grid on the x-axis
|
||
* @param[in] ymin The minimal coordinate of the DEM grid on the y-axis
|
||
* @param[in] ymax The maximal coordinate of the DEM grid on the y-axis
|
||
* @param[in] dx Data spacing of the DEM grid on the x-axis
|
||
* @param[in] dy Data spacing of the DEM grid on the y-axis
|
||
* @param out_verts The output vector of vertex's pointers. The user need to destroy the memories allocated by the function before destroy the vector
|
||
* @param out_tris The output vector of triangle's pointers. The user need to destroy the memories allocated by the function before destroy the vector
|
||
* @param[in] maxi_err Threshold to quit the algorithm. The default is 1e-0
|
||
* @param[in] err_records If this pointer is not NULL, record maximal error values after each insertion of vertex.
|
||
*/
|
||
void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ymin, double ymax,
|
||
double dx, double dy, std::vector<vertex2dc*> &out_verts, std::vector<triangle*> &out_tris,
|
||
double maxi_err = 1e-0, std::vector<double> *err_records = nullptr)
|
||
{
|
||
if (!out_verts.empty()) out_verts.clear();
|
||
if (!out_tris.empty()) out_tris.clear();
|
||
if (err_records != nullptr && !err_records->empty()) err_records->clear();
|
||
|
||
if (dx <= 0.0 || dy <= 0.0 || maxi_err <= 0.0) return;
|
||
if (xmin >= xmax || ymin >= ymax || (xmin + dx) > xmax || (ymin + dy) > ymax) return;
|
||
|
||
int xnum = round((xmax - xmin)/dx) + 1;
|
||
int ynum = round((ymax - ymin)/dy) + 1;
|
||
|
||
if (dem.size() != xnum*ynum) return;
|
||
|
||
// Prepare the DEM points
|
||
std::vector<dem_point> dem_grid(xnum*ynum);
|
||
std::vector<dem_point>::iterator d_iter;
|
||
for (int i = 0; i < ynum; ++i)
|
||
{
|
||
for (int j = 0; j < xnum; ++j)
|
||
{
|
||
dem_grid[j + i*xnum].set(xmin + dx*j, ymin + dy*i, dem[j + i*xnum]);
|
||
}
|
||
}
|
||
|
||
vertex2dc *tmp_vert = nullptr;
|
||
|
||
tmp_vert = new vertex2dc(xmin, ymin, dem_grid[0].elev, out_verts.size()); // lower left corner
|
||
out_verts.push_back(tmp_vert);
|
||
d_iter = dem_grid.begin(); dem_grid.erase(d_iter);
|
||
|
||
tmp_vert = new vertex2dc(xmax, ymin, dem_grid[xnum-2].elev, out_verts.size()); // lower right corner. Note the first location is already erased
|
||
out_verts.push_back(tmp_vert);
|
||
d_iter = dem_grid.begin() + (xnum - 2); dem_grid.erase(d_iter);
|
||
|
||
tmp_vert = new vertex2dc(xmax, ymax, dem_grid[xnum*ynum-3].elev, out_verts.size()); // upper right corner. Note the first two locations are already erased
|
||
out_verts.push_back(tmp_vert);
|
||
d_iter = dem_grid.begin() + (xnum*ynum - 3); dem_grid.erase(d_iter);
|
||
|
||
tmp_vert = new vertex2dc(xmin, ymax, dem_grid[xnum*(ynum-1) - 2].elev, out_verts.size()); // upper left corner. Note the first two locations are already erased
|
||
out_verts.push_back(tmp_vert);
|
||
d_iter = dem_grid.begin() + (xnum*(ynum-1) - 2); dem_grid.erase(d_iter);
|
||
|
||
triangle *tmp_tri = nullptr;
|
||
std::vector<triangle*> cnst_tri, new_tri;
|
||
std::vector<triangle*>::iterator t_iter;
|
||
|
||
if (!is_collinear(out_verts[0], out_verts[1], out_verts[2])) // Do not create triangle if the vertexes are collinear
|
||
{
|
||
tmp_tri = new triangle(out_verts[0], out_verts[1], out_verts[2]); // order the vertex anti-clock wise
|
||
out_tris.push_back(tmp_tri);
|
||
}
|
||
|
||
if (!is_collinear(out_verts[0], out_verts[2], out_verts[3]))
|
||
{
|
||
tmp_tri = new triangle(out_verts[0], out_verts[2], out_verts[3]); // order the vertex anti-clock wise
|
||
out_tris.push_back(tmp_tri);
|
||
}
|
||
|
||
// Find host triangle for all DEM locations
|
||
for (int i = 0; i < dem_grid.size(); ++i)
|
||
{
|
||
for (int t = 0; t < out_tris.size(); ++t)
|
||
{
|
||
if (out_tris[t]->bound_location(dem_grid[i].x, dem_grid[i].y))
|
||
{
|
||
dem_grid[i].host = out_tris[t];
|
||
break; // already found, no need to search more
|
||
}
|
||
}
|
||
}
|
||
|
||
// Find circum_host triangles for all DEM locations
|
||
double dist;
|
||
for (int i = 0; i < dem_grid.size(); ++i)
|
||
{
|
||
for (int t = 0; t < out_tris.size(); ++t)
|
||
{
|
||
dist = (out_tris[t]->cx - dem_grid[i].x) * (out_tris[t]->cx - dem_grid[i].x)
|
||
+ (out_tris[t]->cy - dem_grid[i].y) * (out_tris[t]->cy - dem_grid[i].y);
|
||
if ((dist - out_tris[t]->cr) <= ZERO) // Points on the circumcircle are also included
|
||
{
|
||
dem_grid[i].circum_host.push_back(out_tris[t]);
|
||
// no beak here. There might be more than one triangle's circumcircle includes the DEM location
|
||
}
|
||
}
|
||
}
|
||
|
||
int now_maxi_id;
|
||
double now_err, now_maxi_err;
|
||
|
||
bool removed;
|
||
edge tmp_edge;
|
||
std::vector<edge> cnst_edge;
|
||
std::vector<edge>::iterator e_iter;
|
||
|
||
do // quit til the threshold is meet
|
||
{
|
||
// loop all DEM data to find the location with maximal error
|
||
now_maxi_err = -1.0;
|
||
for (int i = 0; i < dem_grid.size(); ++i)
|
||
{
|
||
now_err = fabs(dem_grid[i].host->interpolate(dem_grid[i].x, dem_grid[i].y) - dem_grid[i].elev);
|
||
if (now_err > now_maxi_err)
|
||
{
|
||
now_maxi_err = now_err;
|
||
now_maxi_id = i;
|
||
}
|
||
}
|
||
|
||
if (err_records != nullptr)
|
||
{
|
||
err_records->push_back(now_maxi_err);
|
||
}
|
||
|
||
// create a new vertex
|
||
tmp_vert = new vertex2dc(dem_grid[now_maxi_id].x, dem_grid[now_maxi_id].y, dem_grid[now_maxi_id].elev, out_verts.size());
|
||
out_verts.push_back(tmp_vert);
|
||
|
||
// Move triangles which circumcircles include the new vertex to the cnst_tri and remove it from out_tris
|
||
cnst_tri.clear();
|
||
for (int i = 0; i < dem_grid[now_maxi_id].circum_host.size(); ++i)
|
||
{
|
||
cnst_tri.push_back(dem_grid[now_maxi_id].circum_host[i]);
|
||
}
|
||
|
||
for (int i = 0; i < cnst_tri.size(); ++i)
|
||
{
|
||
for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
|
||
{
|
||
tmp_tri = *t_iter;
|
||
if (cnst_tri[i] == tmp_tri)
|
||
{
|
||
t_iter = out_tris.erase(t_iter);
|
||
break; // no need to search more
|
||
}
|
||
else t_iter++;
|
||
}
|
||
}
|
||
|
||
// clear host and circumcircle triangles for the used DEM location
|
||
dem_grid[now_maxi_id].host = nullptr;
|
||
dem_grid[now_maxi_id].circum_host.clear();
|
||
d_iter = dem_grid.begin() + now_maxi_id; dem_grid.erase(d_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 out_tris
|
||
new_tri.clear();
|
||
for (int c = 0; c < cnst_edge.size(); ++c)
|
||
{
|
||
if (!is_collinear(cnst_edge[c].vert[0], cnst_edge[c].vert[1], tmp_vert)) // Do not create triangle if the vertexes are collinear
|
||
{
|
||
tmp_tri = new triangle(cnst_edge[c].vert[0], cnst_edge[c].vert[1], tmp_vert); // order the vertex anti-clock wise
|
||
out_tris.push_back(tmp_tri);
|
||
new_tri.push_back(tmp_tri);
|
||
}
|
||
}
|
||
|
||
// purge circumcircle triangles for all DEM data
|
||
for (int c = 0; c < cnst_tri.size(); ++c)
|
||
{
|
||
for (int i = 0; i < dem_grid.size(); ++i)
|
||
{
|
||
for (t_iter = dem_grid[i].circum_host.begin(); t_iter != dem_grid[i].circum_host.end(); )
|
||
{
|
||
if (cnst_tri[c] == *t_iter)
|
||
{
|
||
t_iter = dem_grid[i].circum_host.erase(t_iter);
|
||
break; // no need to search more
|
||
}
|
||
else t_iter++;
|
||
}
|
||
}
|
||
}
|
||
|
||
// loop all DEM data to update host and circumcircle triangles
|
||
for (int i = 0; i < dem_grid.size(); ++i)
|
||
{
|
||
for (int n = 0; n < new_tri.size(); ++n) // search in newly created triangles to find new host
|
||
{
|
||
if (new_tri[n]->bound_location(dem_grid[i].x, dem_grid[i].y))
|
||
{
|
||
dem_grid[i].host = new_tri[n];
|
||
break; // already found, no need to search more
|
||
}
|
||
}
|
||
}
|
||
|
||
// Find circum_host triangles for all DEM locations
|
||
for (int i = 0; i < dem_grid.size(); ++i)
|
||
{
|
||
for (int n = 0; n < new_tri.size(); ++n) // search in newly created triangles to find new circumcircle triangles
|
||
{
|
||
dist = (new_tri[n]->cx - dem_grid[i].x) * (new_tri[n]->cx - dem_grid[i].x)
|
||
+ (new_tri[n]->cy - dem_grid[i].y) * (new_tri[n]->cy - dem_grid[i].y);
|
||
if ((dist - new_tri[n]->cr) <= ZERO) // Points on the circumcircle are also included
|
||
{
|
||
dem_grid[i].circum_host.push_back(new_tri[n]);
|
||
// no beak here. There might be more than one triangle's circumcircle includes the DEM location
|
||
}
|
||
}
|
||
}
|
||
|
||
// 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;
|
||
}
|
||
|
||
} while (now_maxi_err >= maxi_err);
|
||
|
||
return;
|
||
}
|
||
|
||
#endif // _TIN_DELAUNAY_H
|