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

add regions to grd2tin and rnd2tin

Browse Source
This commit is contained in:
张壹 2022-03-07 14:47:04 +08:00
parent f112ca41ee
commit 1d3dd9ccf9
7 changed files with 12464 additions and 14692 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3474
data/topo_TIN.log
View File

File diff suppressed because it is too large Load Diff

15128
data/topo_TIN.msh
View File

File diff suppressed because it is too large Load Diff

8180
data/topo_TIN.neigh
View File

File diff suppressed because it is too large Load Diff

12
src/demo/demo.cpp
View File

@ -14,10 +14,20 @@ int main(int argc, char const *argv[])
} }
infile.close(); infile.close();
std::vector<vertex2dc> box(4);
box[0].id = 0; box[0].x = 500; box[0].y = 250;
box[1].id = 1; box[1].x = 750; box[1].y = 250;
box[2].id = 2; box[2].x = 750; box[2].y = 500;
box[3].id = 3; box[3].x = 500; box[3].y = 500;
std::vector<region> box_region(1);
box_region[0].set(box, 5.0);
std::vector<double> err_records; std::vector<double> err_records;
std::vector<vertex2dc*> tin_vert; std::vector<vertex2dc*> tin_vert;
std::vector<triangle*> tin_ele; std::vector<triangle*> tin_ele;
grd2tin(topo, 0, 1000, 0, 1000, 10, 10, tin_vert, tin_ele, 1.0, nullptr, &err_records); //grd2tin(topo, 0, 1000, 0, 1000, 10, 10, tin_vert, tin_ele, 1.0, nullptr, &err_records, nullptr);
grd2tin(topo, 0, 1000, 0, 1000, 10, 10, tin_vert, tin_ele, 1.0, nullptr, &err_records, &box_region);
// Write a log file // Write a log file
std::ofstream logfile("data/topo_TIN.log"); std::ofstream logfile("data/topo_TIN.log");

320
src/lib/tin.cpp
View File

@ -13,6 +13,10 @@
*/ */
#include "tin.h" #include "tin.h"
#include "cmath"
#include "algorithm"
#include "exception"
#include "stdexcept"
// Start vertex definition // Start vertex definition
vertex2dc::vertex2dc() vertex2dc::vertex2dc()
@ -21,12 +25,12 @@ vertex2dc::vertex2dc()
id = 0; id = 0;
} }
vertex2dc::vertex2dc(double inx, double iny, double inelev, unsigned int inid) vertex2dc::vertex2dc(double inx, double iny, double inelev, size_t inid)
{ {
set(inx, iny, inelev, inid); set(inx, iny, inelev, inid);
} }
void vertex2dc::set(double inx, double iny, double inelev, unsigned int inid) void vertex2dc::set(double inx, double iny, double inelev, size_t inid)
{ {
x = inx; y = iny; elev = inelev; id = inid; x = inx; y = iny; elev = inelev; id = inid;
return; return;
@ -42,6 +46,15 @@ bool operator ==(const vertex2dc &a, const vertex2dc &b)
return false; return false;
} }
bool duplicated_vertex(const vertex2dc &a, const vertex2dc &b)
{
if (a == b && a.id == b.id)
{
return true;
}
return false;
}
// Test if the three points are on the same line // Test if the three points are on the same line
bool is_collinear(vertex2dc *a_ptr, vertex2dc *b_ptr, vertex2dc *c_ptr) bool is_collinear(vertex2dc *a_ptr, vertex2dc *b_ptr, vertex2dc *c_ptr)
{ {
@ -62,16 +75,47 @@ void circumcircle(vertex2dc *v0, vertex2dc *v1, vertex2dc *v2, double &cx, doubl
cx = ((v2->y - v0->y)*m + (v0->y - v1->y)*u)*s; cx = ((v2->y - v0->y)*m + (v0->y - v1->y)*u)*s;
cy = ((v0->x - v2->x)*m + (v1->x - v0->x)*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 // not need to calculate the squared root here
cr = (v0->x - cx)*(v0->x - cx) + (v0->y - cy)*(v0->y - cy);
return; return;
} }
// End vertex definition // End vertex definition
// Start region definition
region::region()
{
maxi_err = 0.0;
}
region::region(const std::vector<vertex2dc> &outline_vert, double err)
{
set(outline_vert, err);
}
void region::set(const std::vector<vertex2dc> &outline_vert, double err)
{
if (err <= 0.0)
{
throw std::invalid_argument("Invalid maximal allowed error.");
}
if (outline_vert.size() < 3)
{
throw std::runtime_error("The input polygon must have at least 3 vertexes.");
}
maxi_err = err;
poly_vert = outline_vert;
return;
}
// End region definition
// Start DEM definition // Start DEM definition
dem_point::dem_point() dem_point::dem_point()
{ {
x = y = elev = NAN; x = y = elev = NAN;
err = 0.0; err = 0.0;
allowed_err = -1.0;
host = nullptr; host = nullptr;
} }
@ -82,11 +126,14 @@ dem_point::dem_point(double inx, double iny, double inelev)
void dem_point::set(double inx, double iny, double inelev) void dem_point::set(double inx, double iny, double inelev)
{ {
x = inx; y = iny; elev = inelev; err = 0.0; host = nullptr; x = inx; y = iny; elev = inelev;
err = 0.0;
allowed_err = -1.0;
host = nullptr;
return; return;
} }
bool compare_dem_point(dem_point *a, dem_point *b) // determination function for std::sort bool compare_dem_error(dem_point *a, dem_point *b) // determination function for std::sort
{ {
if (a->err > b->err) return true; if (a->err > b->err) return true;
return false; return false;
@ -102,13 +149,13 @@ triangle::triangle()
triangle::triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) triangle::triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr)
{ {
set(v0ptr, v1ptr, v2ptr); set_vertex(v0ptr, v1ptr, v2ptr);
neigh[0] = neigh[1] = neigh[2] = nullptr;
} }
void triangle::set(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) void triangle::set_vertex(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;
circumcircle(vert[0], vert[1], vert[2], cx, cy, cr); circumcircle(vert[0], vert[1], vert[2], cx, cy, cr);
return; return;
} }
@ -256,17 +303,27 @@ void flip_neighboring_triangles(triangle *t, triangle *n, int t_id, int n_id)
{ {
tmp_dem = n->hosted_dem[i]; tmp_dem = n->hosted_dem[i];
tmp_dem->err = fabs(n->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev); tmp_dem->err = fabs(n->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
// If the dem point belongs to a region, then set err to 0 if satisfy the region's threshold
if (tmp_dem->allowed_err > 0.0 && tmp_dem->err <= tmp_dem->allowed_err)
{
tmp_dem->err = 0.0;
}
} }
for (int i = 0; i < t->hosted_dem.size(); i++) for (int i = 0; i < t->hosted_dem.size(); i++)
{ {
tmp_dem = t->hosted_dem[i]; tmp_dem = t->hosted_dem[i];
tmp_dem->err = fabs(t->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev); tmp_dem->err = fabs(t->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
// If the dem point belongs to a region, then set err to 0 if satisfy the region's threshold
if (tmp_dem->allowed_err > 0.0 && tmp_dem->err <= tmp_dem->allowed_err)
{
tmp_dem->err = 0.0;
}
} }
// Sort maximal errors for triangles t and n // Sort maximal errors for triangles t and n
std::sort(t->hosted_dem.begin(), t->hosted_dem.end(), compare_dem_point); std::sort(t->hosted_dem.begin(), t->hosted_dem.end(), compare_dem_error);
std::sort(n->hosted_dem.begin(), n->hosted_dem.end(), compare_dem_point); std::sort(n->hosted_dem.begin(), n->hosted_dem.end(), compare_dem_error);
return; return;
} }
@ -396,7 +453,7 @@ triangle *split_triangle(vertex2dc *v, triangle *t, triangle *new_t[4])
break; break;
} }
} }
return t->neigh[i]; // Return the neighboring tiangle to be deleted return t->neigh[i]; // Return the neighboring triangle to be deleted
} }
} }
@ -444,7 +501,7 @@ bool triangle_inside_polygon(triangle *tri_p, std::vector<vertex2dc> *poly_vert)
// If any vertexes of the input triangle is outside of the polygon, return false. Otherwise, return true // If any vertexes of the input triangle is outside of the polygon, return false. Otherwise, return true
if (poly_vert->size() < 3) if (poly_vert->size() < 3)
{ {
return true; throw std::runtime_error("The input polygon must have at least 3 vertexes.");
} }
else else
{ {
@ -500,7 +557,7 @@ bool triangle_inside_polygon(triangle *tri_p, std::vector<vertex2dc> *poly_vert)
} }
/** /**
* @brief Test if the input locaiton is inside of the given polygon * @brief Test if the input location is inside of the given polygon
* *
* @param x Input x coordinate * @param x Input x coordinate
* @param y Input y coordinate * @param y Input y coordinate
@ -513,7 +570,7 @@ bool location_inside_polygon(double x, double y, std::vector<vertex2dc> *poly_ve
// If any vertexes of the input triangle is outside of the polygon, return false. Otherwise, return true // If any vertexes of the input triangle is outside of the polygon, return false. Otherwise, return true
if (poly_vert->size() < 3) if (poly_vert->size() < 3)
{ {
return true; throw std::runtime_error("The input polygon must have at least 3 vertexes.");
} }
else else
{ {
@ -566,7 +623,7 @@ bool location_inside_polygon(double x, double y, std::vector<vertex2dc> *poly_ve
/** /**
* @brief Generate the TIN from the DEM grid * @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] dem Input DEM grid (Ordered from lower left corner to the upper right corner in a row-major fashion)
* @param[in] xmin The minimal coordinate of the DEM grid on the x-axis * @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] 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] ymin The minimal coordinate of the DEM grid on the y-axis
@ -580,24 +637,48 @@ bool location_inside_polygon(double x, double y, std::vector<vertex2dc> *poly_ve
*/ */
void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ymin, double ymax, double dx, void grd2tin(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, double dy, std::vector<vertex2dc*> &out_verts, std::vector<triangle*> &out_tris, double maxi_err,
std::vector<vertex2dc> *outline_poly, std::vector<double> *err_records) std::vector<vertex2dc> *outline_poly, std::vector<double> *err_records,
std::vector<region> *regions, std::vector<double> *dem_error)
{ {
if (!out_verts.empty()) out_verts.clear(); if (!out_verts.empty()) out_verts.clear();
if (!out_tris.empty()) out_tris.clear(); if (!out_tris.empty()) out_tris.clear();
if (err_records != nullptr && !err_records->empty()) err_records->clear(); if (err_records != nullptr && !err_records->empty()) err_records->clear();
if (dx <= 0.0 || dy <= 0.0 || maxi_err <= 0.0) return; if (maxi_err <= 0.0)
if (xmin >= xmax || ymin >= ymax || (xmin + dx) > xmax || (ymin + dy) > ymax) return; {
throw std::invalid_argument("Invalid maximal allowed error.");
}
if (dx <= 0.0 || dy <= 0.0 || xmin >= xmax || ymin >= ymax || (xmin + dx) > xmax || (ymin + dy) > ymax)
{
throw std::invalid_argument("Invalid DEM parameters.");
}
int xnum = round((xmax - xmin)/dx) + 1; int xnum = round((xmax - xmin)/dx) + 1;
int ynum = round((ymax - ymin)/dy) + 1; int ynum = round((ymax - ymin)/dy) + 1;
if (dem.size() != xnum*ynum) return; if (dem.size() != xnum*ynum)
{
throw std::invalid_argument("Invalid size of the input DEM values w.r.t. the DEM parameters.");
}
// Prepare the DEM points if (dem_error != nullptr)
dem_point *tmp_dem = nullptr; {
dem_point *maxi_err_dem = nullptr; if (dem_error->size() != xnum*ynum)
std::vector<dem_point*>::iterator d_iter; {
throw std::invalid_argument("Invalid size of the input DEM errors w.r.t. the DEM parameters.");
}
for (size_t i = 0; i < xnum*ynum; i++)
{
if (dem_error->at(i) <= 0.0)
{
throw std::runtime_error("Invalid input DEM error value.");
}
}
}
// Prepare the initial vertexes and triangles from corner DEM points
vertex2dc *tmp_vert = nullptr; vertex2dc *tmp_vert = nullptr;
tmp_vert = new vertex2dc(xmin, ymin, dem[0], out_verts.size()); // lower left corner tmp_vert = new vertex2dc(xmin, ymin, dem[0], out_verts.size()); // lower left corner
@ -612,9 +693,7 @@ void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
tmp_vert = new vertex2dc(xmin, ymax, dem[xnum*(ynum-1)], out_verts.size()); // upper left corner tmp_vert = new vertex2dc(xmin, ymax, dem[xnum*(ynum-1)], out_verts.size()); // upper left corner
out_verts.push_back(tmp_vert); out_verts.push_back(tmp_vert);
triangle *old_tri = nullptr, *tmp_tri = nullptr; triangle *tmp_tri = nullptr;
triangle *cnst_tri[4];
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 if (!is_collinear(out_verts[0], out_verts[1], out_verts[2])) // Do not create triangle if the vertexes are collinear
{ {
@ -628,13 +707,18 @@ void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
out_tris.push_back(tmp_tri); tmp_tri = nullptr; out_tris.push_back(tmp_tri); tmp_tri = nullptr;
} }
if (out_tris.size() != 2) return; if (out_tris.size() != 2)
{
throw std::runtime_error("Fail to construct the initial triangles.");
}
out_tris[0]->set_neighbor(nullptr, nullptr, out_tris[1]); out_tris[0]->set_neighbor(nullptr, nullptr, out_tris[1]);
out_tris[1]->set_neighbor(out_tris[0], nullptr, nullptr); out_tris[1]->set_neighbor(out_tris[0], nullptr, nullptr);
// Find host triangle for all DEM locations // Find host triangle for all DEM locations
int tmp_id; int tmp_id;
dem_point *tmp_dem = nullptr;
if (outline_poly == nullptr) if (outline_poly == nullptr)
{ {
for (int i = 0; i < ynum; ++i) for (int i = 0; i < ynum; ++i)
@ -642,9 +726,11 @@ void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
for (int j = 0; j < xnum; ++j) for (int j = 0; j < xnum; ++j)
{ {
tmp_id = j + i*xnum; tmp_id = j + i*xnum;
if (tmp_id != 0 && tmp_id != (xnum-1) && tmp_id != (xnum*ynum-1) && tmp_id != (xnum*(ynum-1))) // the four corners are already used if (tmp_id != 0 && tmp_id != (xnum-1) &&
tmp_id != (xnum*ynum-1) && tmp_id != (xnum*(ynum-1))) // the four corners are already used
{ {
tmp_dem = new dem_point(xmin + dx*j, ymin + dy*i, dem[j + i*xnum]); tmp_dem = new dem_point(xmin + dx*j, ymin + dy*i, dem[j + i*xnum]);
for (int t = 0; t < out_tris.size(); ++t) for (int t = 0; t < out_tris.size(); ++t)
{ {
if (out_tris[t]->bound_location(tmp_dem->x, tmp_dem->y)) if (out_tris[t]->bound_location(tmp_dem->x, tmp_dem->y))
@ -652,6 +738,29 @@ void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
tmp_dem->host = out_tris[t]; tmp_dem->host = out_tris[t];
tmp_dem->err = fabs(out_tris[t]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev); tmp_dem->err = fabs(out_tris[t]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
out_tris[t]->hosted_dem.push_back(tmp_dem); out_tris[t]->hosted_dem.push_back(tmp_dem);
// Find region
if (regions != nullptr)
{
for (size_t r = 0; r < regions->size(); r++)
{
if (location_inside_polygon(tmp_dem->x, tmp_dem->y, &(regions->at(r).poly_vert)))
{
tmp_dem->allowed_err = regions->at(r).maxi_err;
break;
}
}
}
if (dem_error != nullptr)
{
tmp_dem->allowed_err = dem_error->at(tmp_id);
}
// If the dem point belongs to a region, then set err to 0 if satisfy the region's threshold
if (tmp_dem->allowed_err > 0.0 && tmp_dem->err <= tmp_dem->allowed_err)
{
tmp_dem->err = 0.0;
}
break; // already found, no need to search more break; // already found, no need to search more
} }
} }
@ -666,11 +775,14 @@ void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
for (int j = 0; j < xnum; ++j) for (int j = 0; j < xnum; ++j)
{ {
tmp_id = j + i*xnum; tmp_id = j + i*xnum;
// the four corners are already used
// We also check to see if the candidate dem point is inside the outline polygon
if (tmp_id != 0 && tmp_id != (xnum-1) && if (tmp_id != 0 && tmp_id != (xnum-1) &&
tmp_id != (xnum*ynum-1) && tmp_id != (xnum*(ynum-1)) && tmp_id != (xnum*ynum-1) && tmp_id != (xnum*(ynum-1)) &&
location_inside_polygon(xmin+dx*j, ymin+dy*i, outline_poly)) // the four corners are already used location_inside_polygon(xmin+dx*j, ymin+dy*i, outline_poly))
{ {
tmp_dem = new dem_point(xmin + dx*j, ymin + dy*i, dem[j + i*xnum]); tmp_dem = new dem_point(xmin + dx*j, ymin + dy*i, dem[j + i*xnum]);
for (int t = 0; t < out_tris.size(); ++t) for (int t = 0; t < out_tris.size(); ++t)
{ {
if (out_tris[t]->bound_location(tmp_dem->x, tmp_dem->y)) if (out_tris[t]->bound_location(tmp_dem->x, tmp_dem->y))
@ -678,6 +790,29 @@ void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
tmp_dem->host = out_tris[t]; tmp_dem->host = out_tris[t];
tmp_dem->err = fabs(out_tris[t]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev); tmp_dem->err = fabs(out_tris[t]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
out_tris[t]->hosted_dem.push_back(tmp_dem); out_tris[t]->hosted_dem.push_back(tmp_dem);
// Find region
if (regions != nullptr)
{
for (size_t r = 0; r < regions->size(); r++)
{
if (location_inside_polygon(tmp_dem->x, tmp_dem->y, &(regions->at(r).poly_vert)))
{
tmp_dem->allowed_err = regions->at(r).maxi_err;
break;
}
}
}
if (dem_error != nullptr)
{
tmp_dem->allowed_err = dem_error->at(tmp_id);
}
// If the dem point belongs to a region, then set err to 0 if satisfy the region's threshold
if (tmp_dem->allowed_err > 0.0 && tmp_dem->err <= tmp_dem->allowed_err)
{
tmp_dem->err = 0.0;
}
break; // already found, no need to search more break; // already found, no need to search more
} }
} }
@ -688,9 +823,10 @@ void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
// Sort hosted_dem in the desceding order with respect to the error. Get pointer of the dem_point with maximal error // Sort hosted_dem in the desceding order with respect to the error. Get pointer of the dem_point with maximal error
double maxi_err_tmp = -1.0; double maxi_err_tmp = -1.0;
dem_point *maxi_err_dem = nullptr;
for (int t = 0; t < out_tris.size(); ++t) for (int t = 0; t < out_tris.size(); ++t)
{ {
std::sort(out_tris[t]->hosted_dem.begin(), out_tris[t]->hosted_dem.end(), compare_dem_point); std::sort(out_tris[t]->hosted_dem.begin(), out_tris[t]->hosted_dem.end(), compare_dem_error);
if (out_tris[t]->hosted_dem[0]->err > maxi_err_tmp) if (out_tris[t]->hosted_dem[0]->err > maxi_err_tmp)
{ {
maxi_err_tmp = out_tris[t]->hosted_dem[0]->err; maxi_err_tmp = out_tris[t]->hosted_dem[0]->err;
@ -698,6 +834,9 @@ void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
} }
} }
triangle *old_tri = nullptr, *cnst_tri[4];
std::vector<dem_point*>::iterator d_iter;
std::vector<triangle*>::iterator t_iter;
while (maxi_err_dem->err >= maxi_err) // quit til the threshold is meet while (maxi_err_dem->err >= maxi_err) // quit til the threshold is meet
{ {
if (err_records != nullptr) if (err_records != nullptr)
@ -768,6 +907,11 @@ void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
tmp_dem->host = cnst_tri[n]; tmp_dem->host = cnst_tri[n];
tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev); tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
cnst_tri[n]->hosted_dem.push_back(tmp_dem); cnst_tri[n]->hosted_dem.push_back(tmp_dem);
// If the dem point belongs to a region, then set err to 0 if satisfy the region's threshold
if (tmp_dem->allowed_err > 0.0 && tmp_dem->err <= tmp_dem->allowed_err)
{
tmp_dem->err = 0.0;
}
break; break;
} }
} }
@ -783,6 +927,11 @@ void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
tmp_dem->host = cnst_tri[n]; tmp_dem->host = cnst_tri[n];
tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev); tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
cnst_tri[n]->hosted_dem.push_back(tmp_dem); cnst_tri[n]->hosted_dem.push_back(tmp_dem);
// If the dem point belongs to a region, then set err to 0 if satisfy the region's threshold
if (tmp_dem->allowed_err > 0.0 && tmp_dem->err <= tmp_dem->allowed_err)
{
tmp_dem->err = 0.0;
}
break; break;
} }
} }
@ -792,7 +941,7 @@ void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
{ {
if (cnst_tri[n] != nullptr) if (cnst_tri[n] != nullptr)
{ {
std::sort(cnst_tri[n]->hosted_dem.begin(), cnst_tri[n]->hosted_dem.end(), compare_dem_point); std::sort(cnst_tri[n]->hosted_dem.begin(), cnst_tri[n]->hosted_dem.end(), compare_dem_error);
} }
} }
@ -816,6 +965,11 @@ void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
tmp_dem->host = cnst_tri[n]; tmp_dem->host = cnst_tri[n];
tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev); tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
cnst_tri[n]->hosted_dem.push_back(tmp_dem); cnst_tri[n]->hosted_dem.push_back(tmp_dem);
// If the dem point belongs to a region, then set err to 0 if satisfy the region's threshold
if (tmp_dem->allowed_err > 0.0 && tmp_dem->err <= tmp_dem->allowed_err)
{
tmp_dem->err = 0.0;
}
break; break;
} }
} }
@ -825,7 +979,7 @@ void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
{ {
if (cnst_tri[n] != nullptr) if (cnst_tri[n] != nullptr)
{ {
std::sort(cnst_tri[n]->hosted_dem.begin(), cnst_tri[n]->hosted_dem.end(), compare_dem_point); std::sort(cnst_tri[n]->hosted_dem.begin(), cnst_tri[n]->hosted_dem.end(), compare_dem_error);
} }
} }
@ -917,15 +1071,39 @@ void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
* @param[in] err_records If this pointer is not NULL, record maximal error values after each insertion of vertex. * @param[in] err_records If this pointer is not NULL, record maximal error values after each insertion of vertex.
*/ */
void rnd2tin(const std::vector<dem_point> &dem, std::vector<vertex2dc*> &out_verts, std::vector<triangle*> &out_tris, void rnd2tin(const std::vector<dem_point> &dem, std::vector<vertex2dc*> &out_verts, std::vector<triangle*> &out_tris,
double maxi_err, std::vector<vertex2dc> *outline_poly, std::vector<double> *err_records) double maxi_err, std::vector<vertex2dc> *outline_poly, std::vector<double> *err_records,
std::vector<region> *regions, std::vector<double> *dem_error)
{ {
if (dem.size() < 3) return; if (dem.size() < 3)
if (maxi_err <= 0.0) return; {
throw std::runtime_error("The input DEM must have at least 3 points.");
}
if (maxi_err <= 0.0)
{
throw std::invalid_argument("Invalid maximal allowed error.");
}
if (!out_verts.empty()) out_verts.clear(); if (!out_verts.empty()) out_verts.clear();
if (!out_tris.empty()) out_tris.clear(); if (!out_tris.empty()) out_tris.clear();
if (err_records != nullptr && !err_records->empty()) err_records->clear(); if (err_records != nullptr && !err_records->empty()) err_records->clear();
if (dem_error != nullptr)
{
if (dem_error->size() != dem.size())
{
throw std::invalid_argument("Invalid size of the input DEM errors w.r.t. the DEM parameters.");
}
for (size_t i = 0; i < dem.size(); i++)
{
if (dem_error->at(i) <= 0.0)
{
throw std::runtime_error("Invalid input DEM error value.");
}
}
}
// locate the surrounding box and initiate the staring two triangles // locate the surrounding box and initiate the staring two triangles
double xmin = dem[0].x, xmax = dem[0].x; double xmin = dem[0].x, xmax = dem[0].x;
double ymin = dem[0].y, ymax = dem[0].y; double ymin = dem[0].y, ymax = dem[0].y;
@ -972,7 +1150,10 @@ void rnd2tin(const std::vector<dem_point> &dem, std::vector<vertex2dc*> &out_ver
out_tris.push_back(tmp_tri); tmp_tri = nullptr; out_tris.push_back(tmp_tri); tmp_tri = nullptr;
} }
if (out_tris.size() != 2) return; if (out_tris.size() != 2)
{
throw std::runtime_error("Fail to construct the initial triangles.");
}
out_tris[0]->set_neighbor(nullptr, nullptr, out_tris[1]); out_tris[0]->set_neighbor(nullptr, nullptr, out_tris[1]);
out_tris[1]->set_neighbor(out_tris[0], nullptr, nullptr); out_tris[1]->set_neighbor(out_tris[0], nullptr, nullptr);
@ -995,6 +1176,29 @@ void rnd2tin(const std::vector<dem_point> &dem, std::vector<vertex2dc*> &out_ver
tmp_dem->host = out_tris[t]; tmp_dem->host = out_tris[t];
tmp_dem->err = fabs(out_tris[t]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev); tmp_dem->err = fabs(out_tris[t]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
out_tris[t]->hosted_dem.push_back(tmp_dem); out_tris[t]->hosted_dem.push_back(tmp_dem);
// Find region
if (regions != nullptr)
{
for (size_t r = 0; r < regions->size(); r++)
{
if (location_inside_polygon(tmp_dem->x, tmp_dem->y, &(regions->at(r).poly_vert)))
{
tmp_dem->allowed_err = regions->at(r).maxi_err;
break;
}
}
}
if (dem_error != nullptr)
{
tmp_dem->allowed_err = dem_error->at(d);
}
// If the dem point belongs to a region, then set err to 0 if satisfy the region's threshold
if (tmp_dem->allowed_err > 0.0 && tmp_dem->err <= tmp_dem->allowed_err)
{
tmp_dem->err = 0.0;
}
break; // already found, no need to search more break; // already found, no need to search more
} }
} }
@ -1014,6 +1218,29 @@ void rnd2tin(const std::vector<dem_point> &dem, std::vector<vertex2dc*> &out_ver
tmp_dem->host = out_tris[t]; tmp_dem->host = out_tris[t];
tmp_dem->err = fabs(out_tris[t]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev); tmp_dem->err = fabs(out_tris[t]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
out_tris[t]->hosted_dem.push_back(tmp_dem); out_tris[t]->hosted_dem.push_back(tmp_dem);
// Find region
if (regions != nullptr)
{
for (size_t r = 0; r < regions->size(); r++)
{
if (location_inside_polygon(tmp_dem->x, tmp_dem->y, &(regions->at(r).poly_vert)))
{
tmp_dem->allowed_err = regions->at(r).maxi_err;
break;
}
}
}
if (dem_error != nullptr)
{
tmp_dem->allowed_err = dem_error->at(d);
}
// If the dem point belongs to a region, then set err to 0 if satisfy the region's threshold
if (tmp_dem->allowed_err > 0.0 && tmp_dem->err <= tmp_dem->allowed_err)
{
tmp_dem->err = 0.0;
}
break; // already found, no need to search more break; // already found, no need to search more
} }
} }
@ -1025,7 +1252,7 @@ void rnd2tin(const std::vector<dem_point> &dem, std::vector<vertex2dc*> &out_ver
double maxi_err_tmp = -1.0; double maxi_err_tmp = -1.0;
for (int t = 0; t < out_tris.size(); ++t) for (int t = 0; t < out_tris.size(); ++t)
{ {
std::sort(out_tris[t]->hosted_dem.begin(), out_tris[t]->hosted_dem.end(), compare_dem_point); std::sort(out_tris[t]->hosted_dem.begin(), out_tris[t]->hosted_dem.end(), compare_dem_error);
if (out_tris[t]->hosted_dem[0]->err > maxi_err_tmp) if (out_tris[t]->hosted_dem[0]->err > maxi_err_tmp)
{ {
maxi_err_tmp = out_tris[t]->hosted_dem[0]->err; maxi_err_tmp = out_tris[t]->hosted_dem[0]->err;
@ -1103,6 +1330,11 @@ void rnd2tin(const std::vector<dem_point> &dem, std::vector<vertex2dc*> &out_ver
tmp_dem->host = cnst_tri[n]; tmp_dem->host = cnst_tri[n];
tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev); tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
cnst_tri[n]->hosted_dem.push_back(tmp_dem); cnst_tri[n]->hosted_dem.push_back(tmp_dem);
// If the dem point belongs to a region, then set err to 0 if satisfy the region's threshold
if (tmp_dem->allowed_err > 0.0 && tmp_dem->err <= tmp_dem->allowed_err)
{
tmp_dem->err = 0.0;
}
break; break;
} }
} }
@ -1118,6 +1350,11 @@ void rnd2tin(const std::vector<dem_point> &dem, std::vector<vertex2dc*> &out_ver
tmp_dem->host = cnst_tri[n]; tmp_dem->host = cnst_tri[n];
tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev); tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
cnst_tri[n]->hosted_dem.push_back(tmp_dem); cnst_tri[n]->hosted_dem.push_back(tmp_dem);
// If the dem point belongs to a region, then set err to 0 if satisfy the region's threshold
if (tmp_dem->allowed_err > 0.0 && tmp_dem->err <= tmp_dem->allowed_err)
{
tmp_dem->err = 0.0;
}
break; break;
} }
} }
@ -1127,7 +1364,7 @@ void rnd2tin(const std::vector<dem_point> &dem, std::vector<vertex2dc*> &out_ver
{ {
if (cnst_tri[n] != nullptr) if (cnst_tri[n] != nullptr)
{ {
std::sort(cnst_tri[n]->hosted_dem.begin(), cnst_tri[n]->hosted_dem.end(), compare_dem_point); std::sort(cnst_tri[n]->hosted_dem.begin(), cnst_tri[n]->hosted_dem.end(), compare_dem_error);
} }
} }
@ -1151,6 +1388,11 @@ void rnd2tin(const std::vector<dem_point> &dem, std::vector<vertex2dc*> &out_ver
tmp_dem->host = cnst_tri[n]; tmp_dem->host = cnst_tri[n];
tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev); tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
cnst_tri[n]->hosted_dem.push_back(tmp_dem); cnst_tri[n]->hosted_dem.push_back(tmp_dem);
// If the dem point belongs to a region, then set err to 0 if satisfy the region's threshold
if (tmp_dem->allowed_err > 0.0 && tmp_dem->err <= tmp_dem->allowed_err)
{
tmp_dem->err = 0.0;
}
break; break;
} }
} }
@ -1160,7 +1402,7 @@ void rnd2tin(const std::vector<dem_point> &dem, std::vector<vertex2dc*> &out_ver
{ {
if (cnst_tri[n] != nullptr) if (cnst_tri[n] != nullptr)
{ {
std::sort(cnst_tri[n]->hosted_dem.begin(), cnst_tri[n]->hosted_dem.end(), compare_dem_point); std::sort(cnst_tri[n]->hosted_dem.begin(), cnst_tri[n]->hosted_dem.end(), compare_dem_error);
} }
} }

41
src/lib/tin.h
View File

@ -14,16 +14,14 @@
#ifndef _TIN_DELAUNAY_H #ifndef _TIN_DELAUNAY_H
#define _TIN_DELAUNAY_H #define _TIN_DELAUNAY_H
#include "cmath"
#include "vector" #include "vector"
#include "algorithm"
#define ZERO 1e-5 #define ZERO 1e-5
// Start vertex definition // Start vertex definition
struct vertex2dc struct vertex2dc
{ {
unsigned int id; // index of the vertex size_t id; // index of the vertex
double x, y; // position of the vertex double x, y; // position of the vertex
double elev; // elevation at the vertex double elev; // elevation at the vertex
@ -37,7 +35,7 @@ struct vertex2dc
* @param inelev initial elevation * @param inelev initial elevation
* @param inid initial index * @param inid initial index
*/ */
vertex2dc(double inx, double iny, double inelev, unsigned int inid); vertex2dc(double inx, double iny, double inelev, size_t inid);
/** /**
* @brief Set a vertex2dc object * @brief Set a vertex2dc object
@ -47,7 +45,7 @@ struct vertex2dc
* @param inelev initial elevation * @param inelev initial elevation
* @param inid initial index * @param inid initial index
*/ */
void set(double inx, double iny, double inelev, unsigned int inid); void set(double inx, double iny, double inelev, size_t inid);
}; };
/** /**
@ -60,6 +58,16 @@ struct vertex2dc
*/ */
bool operator ==(const vertex2dc &a, const vertex2dc &b); bool operator ==(const vertex2dc &a, const vertex2dc &b);
/**
* @brief Compare two vertexes to see if they are duplicated.
*
* @param a vertex a
* @param b vertex b
* @return true the two vertexes are at the same location and a same index
* @return false the two vertexes are not at the same location and a same index
*/
bool duplicated_vertex(const vertex2dc &a, const vertex2dc &b);
/** /**
* @brief Test if the three points are on the same line * @brief Test if the three points are on the same line
* *
@ -84,6 +92,20 @@ bool is_collinear(vertex2dc *a_ptr, vertex2dc *b_ptr, vertex2dc *c_ptr);
void circumcircle(vertex2dc *v0, vertex2dc *v1, vertex2dc *v2, double &cx, double &cy, double &cr); void circumcircle(vertex2dc *v0, vertex2dc *v1, vertex2dc *v2, double &cx, double &cy, double &cr);
// End vertex definition // End vertex definition
// Start region definition
struct region
{
double maxi_err;
std::vector<vertex2dc> poly_vert;
region();
region(const std::vector<vertex2dc> &outline_vert, double err);
void set(const std::vector<vertex2dc> &outline_vert, double err);
};
// End region definition
// Start DEM definition // Start DEM definition
struct triangle; struct triangle;
@ -92,6 +114,7 @@ struct dem_point
double x, y; // position of the DEM location double x, y; // position of the DEM location
double elev; // elevation at the DEM location double elev; // elevation at the DEM location
double err; // error of the TIN with respect to the elevation double err; // error of the TIN with respect to the elevation
double allowed_err; // point specified threshold of maximal allowed error
triangle *host; // pointer of the triangle that the dem_point falls inside of triangle *host; // pointer of the triangle that the dem_point falls inside of
dem_point(); dem_point();
@ -152,7 +175,7 @@ struct triangle
* @param v1ptr pointer of the vertex 1 * @param v1ptr pointer of the vertex 1
* @param v2ptr pointer of the vertex 2 * @param v2ptr pointer of the vertex 2
*/ */
void set(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr); void set_vertex(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr);
/** /**
* @brief Set neighbors of a triangle object * @brief Set neighbors of a triangle object
@ -202,7 +225,8 @@ struct triangle
*/ */
void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ymin, double ymax, void grd2tin(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 dx, double dy, std::vector<vertex2dc*> &out_verts, std::vector<triangle*> &out_tris,
double maxi_err = 1e-0, std::vector<vertex2dc> *outline_poly = nullptr, std::vector<double> *err_records = nullptr); double maxi_err = 1e-0, std::vector<vertex2dc> *outline_poly = nullptr, std::vector<double> *err_records = nullptr,
std::vector<region> *regions = nullptr, std::vector<double> *dem_error = nullptr);
/** /**
* @brief Generate the TIN from random DEM points * @brief Generate the TIN from random DEM points
@ -216,6 +240,7 @@ void grd2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
*/ */
void rnd2tin(const std::vector<dem_point> &dem, std::vector<vertex2dc*> &out_verts, void rnd2tin(const std::vector<dem_point> &dem, std::vector<vertex2dc*> &out_verts,
std::vector<triangle*> &out_tris, double maxi_err = 1e-0, std::vector<vertex2dc> *outline_poly = nullptr, std::vector<triangle*> &out_tris, double maxi_err = 1e-0, std::vector<vertex2dc> *outline_poly = nullptr,
std::vector<double> *err_records = nullptr); std::vector<double> *err_records = nullptr, std::vector<region> *regions = nullptr,
std::vector<double> *dem_error = nullptr);
#endif // _TIN_DELAUNAY_H #endif // _TIN_DELAUNAY_H

1
src/tool/grd2tin.cpp
View File

@ -5,6 +5,7 @@
#include "fstream" #include "fstream"
#include "sstream" #include "sstream"
#include "string" #include "string"
#include "cmath"
#include "getopt.h" #include "getopt.h"
void display_help(std::string exe_name) void display_help(std::string exe_name)