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

tmp update

Browse Source
This commit is contained in:
张壹
2021-10-15 16:53:12 +08:00
parent e3812057b7
commit 10bb022a07
20 changed files with 44068 additions and 122 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
src/CMakeLists.txt
View File

@@ -66,3 +66,6 @@ macro(add_demo name)
endmacro()
add_demo(demo)
add_demo(demo2)
add_demo(demo3)
add_demo(demo4)

2
src/demo/demo.cpp
View File

@@ -17,7 +17,7 @@ int main(int argc, char const *argv[])
std::vector<double> err_records;
std::vector<vertex2dc*> tin_vert;
std::vector<triangle*> tin_ele;
dem2tin(topo, 0, 1000, 0, 1000, 10, 10, tin_vert, tin_ele, 1.0, &err_records);
dem2tin(topo, 0, 1000, 0, 1000, 10, 10, tin_vert, tin_ele, 1.0, nullptr, &err_records);
// Write a log file
std::ofstream logfile("data/topo_TIN.log");

94
src/demo/demo2.cpp Normal file
View File

@@ -0,0 +1,94 @@
#include "../lib/tin.h"
#include "iostream"
#include "fstream"
#include "iomanip"
int main(int argc, char const *argv[])
{
// read dem grid
std::vector<dem_point> topo(8000);
std::ifstream infile("data/topo_rnd");
for (int i = 0; i < 8000; ++i)
{
infile >> topo[i].x >> topo[i].y >> topo[i].elev;
}
infile.close();
std::vector<double> err_records;
std::vector<vertex2dc*> tin_vert;
std::vector<triangle*> tin_ele;
rnd2tin(topo, tin_vert, tin_ele, 1.0, nullptr, &err_records);
// Write a log file
std::ofstream logfile("data/topo_rnd_TIN.log");
logfile << "# Insertion Maxi-Error\n";
for (int i = 0; i < err_records.size(); ++i)
{
logfile << i+1 << " " << err_records[i] << std::endl;
}
logfile.close();
// Write a Gmsh's .msh file
std::ofstream outfile("data/topo_rnd_TIN.msh");
outfile << "$MeshFormat" << std::endl << "2.2 0 8" << std::endl << "$EndMeshFormat "<<std::endl;
outfile << "$Nodes" << std::endl << tin_vert.size() << std::endl;
for (int i = 0; i < tin_vert.size(); i++)
{
outfile << tin_vert[i]->id + 1 << " " << std::setprecision(16)
<< tin_vert[i]->x << " " << tin_vert[i]->y << " " << tin_vert[i]->elev << std::endl;
}
outfile<<"$EndNodes"<<std::endl;
outfile << "$Elements" << std::endl << tin_ele.size() <<std::endl;
for (int i = 0; i < tin_ele.size(); i++)
{
outfile << i + 1 << " 2 0";
for (int j = 0; j < 3; j++)
{
outfile << " " << tin_ele[i]->vert[j]->id + 1;
}
outfile << std::endl;
}
outfile << "$EndElements"<< std::endl;
outfile<<"$NodeData"<<std::endl;
outfile<<1<<std::endl
<<"\"Topography (m)\"" <<std::endl
<< 1 <<std::endl<< 0.0 <<std::endl
<< 3 <<std::endl<< 0<<std::endl
<< 1 <<std::endl<< tin_vert.size() <<std::endl;
for (int i = 0; i < tin_vert.size(); i++)
{
outfile << tin_vert[i]->id + 1 << " " << std::setprecision(16) << tin_vert[i]->elev << std::endl;
}
outfile << "$EndNodeData" << std::endl;
outfile.close();
// write a neighbor file
outfile.open("data/topo_rnd_TIN.neigh");
outfile << tin_ele.size() << std::endl;
for (int i = 0; i < tin_ele.size(); i++)
{
outfile << i + 1;
for (int j = 0; j < 3; j++)
{
if (tin_ele[i]->neigh[j] != nullptr)
{
outfile << " " << tin_ele[i]->neigh[j]->id + 1;
}
else outfile << " -1";
}
outfile << std::endl;
}
outfile.close();
// Destroy memories allocated by the dem2tin function
for (int i = 0; i < tin_vert.size(); ++i)
{
delete tin_vert[i];
}
for (int i = 0; i < tin_ele.size(); ++i)
{
delete tin_ele[i];
}
return 0;
}

105
src/demo/demo3.cpp Normal file
View File

@@ -0,0 +1,105 @@
#include "../lib/tin.h"
#include "iostream"
#include "fstream"
#include "iomanip"
int main(int argc, char const *argv[])
{
// read dem grid
std::vector<double> topo(10201);
std::ifstream infile("data/topo");
for (int i = 0; i < 10201; ++i)
{
infile >> topo[i];
}
infile.close();
// Set outline polygon
std::vector<vertex2dc> valid_area(8);
valid_area[0].set(0, 500, 0, 0);
valid_area[1].set(58, 365, 0, 1);
valid_area[2].set(314, 158, 0, 2);
valid_area[3].set(681, 22, 0, 3);
valid_area[4].set(942, 105, 0, 4);
valid_area[5].set(1000, 360, 0, 5);
valid_area[6].set(1000, 1000, 0, 6);
valid_area[7].set(0, 1000, 0, 7);
std::vector<double> err_records;
std::vector<vertex2dc*> tin_vert;
std::vector<triangle*> tin_ele;
dem2tin(topo, 0, 1000, 0, 1000, 10, 10, tin_vert, tin_ele, 1.0, &valid_area, &err_records);
// Write a log file
std::ofstream logfile("data/topo_TIN2.log");
logfile << "# Insertion Maxi-Error\n";
for (int i = 0; i < err_records.size(); ++i)
{
logfile << i+1 << " " << err_records[i] << std::endl;
}
logfile.close();
// Write a Gmsh's .msh file
std::ofstream outfile("data/topo_TIN2.msh");
outfile << "$MeshFormat" << std::endl << "2.2 0 8" << std::endl << "$EndMeshFormat "<<std::endl;
outfile << "$Nodes" << std::endl << tin_vert.size() << std::endl;
for (int i = 0; i < tin_vert.size(); i++)
{
outfile << tin_vert[i]->id + 1 << " " << std::setprecision(16)
<< tin_vert[i]->x << " " << tin_vert[i]->y << " " << tin_vert[i]->elev << std::endl;
}
outfile<<"$EndNodes"<<std::endl;
outfile << "$Elements" << std::endl << tin_ele.size() <<std::endl;
for (int i = 0; i < tin_ele.size(); i++)
{
outfile << i + 1 << " 2 0";
for (int j = 0; j < 3; j++)
{
outfile << " " << tin_ele[i]->vert[j]->id + 1;
}
outfile << std::endl;
}
outfile << "$EndElements"<< std::endl;
outfile<<"$NodeData"<<std::endl;
outfile<<1<<std::endl
<<"\"Topography (m)\"" <<std::endl
<< 1 <<std::endl<< 0.0 <<std::endl
<< 3 <<std::endl<< 0<<std::endl
<< 1 <<std::endl<< tin_vert.size() <<std::endl;
for (int i = 0; i < tin_vert.size(); i++)
{
outfile << tin_vert[i]->id + 1 << " " << std::setprecision(16) << tin_vert[i]->elev << std::endl;
}
outfile << "$EndNodeData" << std::endl;
outfile.close();
// write a neighbor file
outfile.open("data/topo_TIN2.neigh");
outfile << tin_ele.size() << std::endl;
for (int i = 0; i < tin_ele.size(); i++)
{
outfile << i + 1;
for (int j = 0; j < 3; j++)
{
if (tin_ele[i]->neigh[j] != nullptr)
{
outfile << " " << tin_ele[i]->neigh[j]->id + 1;
}
else outfile << " -1";
}
outfile << std::endl;
}
outfile.close();
// Destroy memories allocated by the dem2tin function
for (int i = 0; i < tin_vert.size(); ++i)
{
delete tin_vert[i];
}
for (int i = 0; i < tin_ele.size(); ++i)
{
delete tin_ele[i];
}
return 0;
}

105
src/demo/demo4.cpp Normal file
View File

@@ -0,0 +1,105 @@
#include "../lib/tin.h"
#include "iostream"
#include "fstream"
#include "iomanip"
int main(int argc, char const *argv[])
{
// read dem grid
std::vector<dem_point> topo(8000);
std::ifstream infile("data/topo_rnd");
for (int i = 0; i < 8000; ++i)
{
infile >> topo[i].x >> topo[i].y >> topo[i].elev;
}
infile.close();
// Set outline polygon
std::vector<vertex2dc> valid_area(8);
valid_area[0].set(0, 500, 0, 0);
valid_area[1].set(58, 365, 0, 1);
valid_area[2].set(314, 158, 0, 2);
valid_area[3].set(681, 22, 0, 3);
valid_area[4].set(942, 105, 0, 4);
valid_area[5].set(1000, 360, 0, 5);
valid_area[6].set(1000, 1000, 0, 6);
valid_area[7].set(0, 1000, 0, 7);
std::vector<double> err_records;
std::vector<vertex2dc*> tin_vert;
std::vector<triangle*> tin_ele;
rnd2tin(topo, tin_vert, tin_ele, 1.0, &valid_area, &err_records);
// Write a log file
std::ofstream logfile("data/topo_rnd_TIN2.log");
logfile << "# Insertion Maxi-Error\n";
for (int i = 0; i < err_records.size(); ++i)
{
logfile << i+1 << " " << err_records[i] << std::endl;
}
logfile.close();
// Write a Gmsh's .msh file
std::ofstream outfile("data/topo_rnd_TIN2.msh");
outfile << "$MeshFormat" << std::endl << "2.2 0 8" << std::endl << "$EndMeshFormat "<<std::endl;
outfile << "$Nodes" << std::endl << tin_vert.size() << std::endl;
for (int i = 0; i < tin_vert.size(); i++)
{
outfile << tin_vert[i]->id + 1 << " " << std::setprecision(16)
<< tin_vert[i]->x << " " << tin_vert[i]->y << " " << tin_vert[i]->elev << std::endl;
}
outfile<<"$EndNodes"<<std::endl;
outfile << "$Elements" << std::endl << tin_ele.size() <<std::endl;
for (int i = 0; i < tin_ele.size(); i++)
{
outfile << i + 1 << " 2 0";
for (int j = 0; j < 3; j++)
{
outfile << " " << tin_ele[i]->vert[j]->id + 1;
}
outfile << std::endl;
}
outfile << "$EndElements"<< std::endl;
outfile<<"$NodeData"<<std::endl;
outfile<<1<<std::endl
<<"\"Topography (m)\"" <<std::endl
<< 1 <<std::endl<< 0.0 <<std::endl
<< 3 <<std::endl<< 0<<std::endl
<< 1 <<std::endl<< tin_vert.size() <<std::endl;
for (int i = 0; i < tin_vert.size(); i++)
{
outfile << tin_vert[i]->id + 1 << " " << std::setprecision(16) << tin_vert[i]->elev << std::endl;
}
outfile << "$EndNodeData" << std::endl;
outfile.close();
// write a neighbor file
outfile.open("data/topo_rnd_TIN2.neigh");
outfile << tin_ele.size() << std::endl;
for (int i = 0; i < tin_ele.size(); i++)
{
outfile << i + 1;
for (int j = 0; j < 3; j++)
{
if (tin_ele[i]->neigh[j] != nullptr)
{
outfile << " " << tin_ele[i]->neigh[j]->id + 1;
}
else outfile << " -1";
}
outfile << std::endl;
}
outfile.close();
// Destroy memories allocated by the dem2tin function
for (int i = 0; i < tin_vert.size(); ++i)
{
delete tin_vert[i];
}
for (int i = 0; i < tin_ele.size(); ++i)
{
delete tin_ele[i];
}
return 0;
}

521
src/lib/tin.cpp
View File

@@ -156,7 +156,7 @@ double triangle::interpolate(double inx, double iny)
* / \
* / \
* / t \
* t_id-------\ t_id (0, 1 or 2)
* t_id-------\ t_id+1 (0, 1 or 2)
* \--------/
* \ /
* \ n /
@@ -294,7 +294,9 @@ void make_delaunay(triangle *t)
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) // A very restrict condition. The testing point must be really inside the circumcircle
// We have to use a very restrict condition here. The testing point must be really inside the circumcircle.
// Otherwise, this recursive function may fall into endless callings till the segment fails.
if ((dist - t->cr) < -1.0*ZERO)
{
flip_neighboring_triangles(t, n_tri, n, v);
@@ -429,6 +431,74 @@ triangle *split_triangle(vertex2dc *v, triangle *t, triangle *new_t[4])
}
// End triangle definition
/**
* @brief Test if the input triangle is inside of the given polygon
*
* @param tri_p Pointer of a test triangle
* @param poly_vert Vertexes of a polygon
* @return true The test triangle is inside of the polygon
* @return false The test triangle is outside of the polygon
*/
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 (poly_vert->size() < 3)
{
return false;
}
else
{
int pnum = poly_vert->size();
//确定外接矩形
double xmin = poly_vert->at(0).x, ymin = poly_vert->at(0).y;
double xmax = poly_vert->at(0).x, ymax = poly_vert->at(0).y;
for (int j = 0; j < pnum; ++j)
{
xmin = std::min(xmin, poly_vert->at(j).x);
xmax = std::max(xmax, poly_vert->at(j).x);
ymin = std::min(ymin, poly_vert->at(j).y);
ymax = std::max(ymax, poly_vert->at(j).y);
}
int count;
double tmp_x;
vertex2dc *one_p;
for (int t = 0; t < 3; ++t)
{
one_p = tri_p->vert[t];
// the testing point is outside of the surrounding box, return false
if (one_p->x < xmin || one_p->x > xmax || one_p->y < ymin || one_p->y > ymax)
{
return false;
}
else
{
count = 0;
for (int i = 0; i < pnum; ++i)
{
if ((one_p->y >= poly_vert->at(i).y && one_p->y < poly_vert->at((i+1)%pnum).y) ||
(one_p->y <= poly_vert->at(i).y && one_p->y > poly_vert->at((i+1)%pnum).y))
{
tmp_x = (poly_vert->at((i+1)%pnum).x - poly_vert->at(i).x)
* (one_p->y - poly_vert->at(i).y)/(poly_vert->at((i+1)%pnum).y - poly_vert->at(i).y)
+ poly_vert->at(i).x;
if (one_p->x <= tmp_x)
count++;
}
}
if (pow(-1, count) > 0) return false;
}
}
// all vertexes are inside the polygon
return true;
}
}
/**
* @brief Generate the TIN from the DEM grid
*
@@ -445,7 +515,8 @@ triangle *split_triangle(vertex2dc *v, triangle *t, triangle *new_t[4])
* @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, std::vector<double> *err_records)
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)
{
if (!out_verts.empty()) out_verts.clear();
if (!out_tris.empty()) out_tris.clear();
@@ -461,7 +532,7 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
// Prepare the DEM points
dem_point *tmp_dem = nullptr;
std::vector<dem_point*> dem_tri;
dem_point *maxi_err_dem = nullptr;
std::vector<dem_point*>::iterator d_iter;
vertex2dc *tmp_vert = nullptr;
@@ -522,28 +593,30 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
}
}
// Sort hosted_dem in the desceding order with respect to the error. Add maximal zeros to dem_tri
// 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;
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);
dem_tri.push_back(out_tris[t]->hosted_dem[0]);
if (out_tris[t]->hosted_dem[0]->err > maxi_err_tmp)
{
maxi_err_tmp = out_tris[t]->hosted_dem[0]->err;
maxi_err_dem = out_tris[t]->hosted_dem[0];
}
}
// Sort dem_tri
std::sort(dem_tri.begin(), dem_tri.end(), compare_dem_point);
while (dem_tri[0]->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)
{
err_records->push_back(dem_tri[0]->err);
err_records->push_back(maxi_err_dem->err);
}
// find the triangle that includes dem_tri[0] and remove it from out_tris
// find the triangle that includes maxi_err_dem and remove it from out_tris
for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
{
old_tri = *t_iter;
if (old_tri == dem_tri[0]->host)
if (old_tri == maxi_err_dem->host)
{
t_iter = out_tris.erase(t_iter);
break;
@@ -551,10 +624,10 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
else t_iter++;
}
// remove dem_tri[0] from its host triangle's hosted DEM list
// remove maxi_err_dem from its host triangle's hosted DEM list
for (d_iter = old_tri->hosted_dem.begin(); d_iter != old_tri->hosted_dem.end(); )
{
if (dem_tri[0] == *d_iter)
if (maxi_err_dem == *d_iter)
{
d_iter = old_tri->hosted_dem.erase(d_iter);
break;
@@ -563,11 +636,11 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
}
// create a new vertex
tmp_vert = new vertex2dc(dem_tri[0]->x, dem_tri[0]->y, dem_tri[0]->elev, out_verts.size());
tmp_vert = new vertex2dc(maxi_err_dem->x, maxi_err_dem->y, maxi_err_dem->elev, out_verts.size());
out_verts.push_back(tmp_vert);
// Delete dem_tri[0]
tmp_dem = dem_tri[0]; delete tmp_dem;
// Delete maxi_err_dem
delete maxi_err_dem; maxi_err_dem = nullptr;
// build new triangles
tmp_tri = split_triangle(tmp_vert, old_tri, cnst_tri);
@@ -677,22 +750,52 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
}
}
// get maximal errors from out_tris and sort dem_tri
dem_tri.clear(); dem_tri.reserve(out_tris.size());
// get maximal errors from out_tris
maxi_err_tmp = -1.0;
for (int t = 0; t < out_tris.size(); t++)
{
if (!out_tris[t]->hosted_dem.empty())
if (!out_tris[t]->hosted_dem.empty() && out_tris[t]->hosted_dem[0]->err > maxi_err_tmp)
{
dem_tri.push_back(out_tris[t]->hosted_dem[0]);
maxi_err_tmp = out_tris[t]->hosted_dem[0]->err;
maxi_err_dem = out_tris[t]->hosted_dem[0];
}
}
std::sort(dem_tri.begin(), dem_tri.end(), compare_dem_point);
}
if (err_records != nullptr)
{
err_records->push_back(dem_tri[0]->err);
err_records->push_back(maxi_err_dem->err);
}
// Cut outline if there is one
if (outline_poly != nullptr)
{
for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
{
tmp_tri = *t_iter;
if (!triangle_inside_polygon(tmp_tri, outline_poly))
{
// set neighbors
for (int i = 0; i < 3; ++i)
{
if (tmp_tri->neigh[i] != nullptr)
{
for (int j = 0; j < 3; ++j)
{
if (tmp_tri->neigh[i]->neigh[j] == tmp_tri)
{
tmp_tri->neigh[i]->neigh[j] = 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
@@ -709,18 +812,378 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
return;
}
/**
* @brief Generate the TIN from random DEM points
*
* @param dem Input DEM points
* @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 dem Input DEM points
* @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] outline_poly If this pointer is not NULL, Cut triangle outside the polygon.
* @param[in] outline_poly If this pointer is not NULL, Cut triangle outside of or intersected with the polygon.
* @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,
double maxi_err, std::vector<dem_point> *outline_poly, std::vector<double> *err_records)
double maxi_err, std::vector<vertex2dc> *outline_poly, std::vector<double> *err_records)
{
if (dem.size() < 3) return;
if (maxi_err <= 0.0) return;
if (!out_verts.empty()) out_verts.clear();
if (!out_tris.empty()) out_tris.clear();
if (err_records != nullptr && !err_records->empty()) err_records->clear();
// locate the surrounding box and initiate the staring two triangles
double xmin = dem[0].x, xmax = dem[0].x;
double ymin = dem[0].y, ymax = dem[0].y;
for (int i = 0; i < dem.size(); ++i)
{
xmin = std::min(xmin, dem[i].x);
xmax = std::max(xmax, dem[i].x);
ymin = std::min(ymin, dem[i].y);
ymax = std::max(ymax, dem[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, 0.0, 0); // lower left corner
box_vert.push_back(tmp_vert);
tmp_vert = new vertex2dc(midx + maxi_s, midy - maxi_s, 0.0, 1); // lower right corner
box_vert.push_back(tmp_vert);
tmp_vert = new vertex2dc(midx + maxi_s, midy + maxi_s, 0.0, 2); // upper right corner
box_vert.push_back(tmp_vert);
tmp_vert = new vertex2dc(midx - maxi_s, midy + maxi_s, 0.0, 3); // upper left corner
box_vert.push_back(tmp_vert);
triangle *old_tri = nullptr, *tmp_tri = nullptr;
triangle *cnst_tri[4];
std::vector<triangle*>::iterator t_iter;
if (!is_collinear(box_vert[0], box_vert[1], box_vert[2])) // Do not create triangle if the vertexes are collinear
{
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 = nullptr;
}
if (!is_collinear(box_vert[0], box_vert[2], box_vert[3]))
{
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); tmp_tri = nullptr;
}
if (out_tris.size() != 2) return;
out_tris[0]->set_neighbor(nullptr, nullptr, out_tris[1]);
out_tris[1]->set_neighbor(out_tris[0], nullptr, nullptr);
// Prepare the DEM points
dem_point *tmp_dem = nullptr;
dem_point *maxi_err_dem = nullptr;
std::vector<dem_point*>::iterator d_iter;
// Find host triangle for all DEM locations
for (int d = 0; d < dem.size(); ++d)
{
tmp_dem = new dem_point(dem[d].x, dem[d].y, dem[d].elev);
for (int t = 0; t < out_tris.size(); ++t)
{
if (out_tris[t]->bound_location(tmp_dem->x, tmp_dem->y))
{
tmp_dem->host = out_tris[t];
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);
break; // already found, no need to search more
}
}
}
// 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;
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);
if (out_tris[t]->hosted_dem[0]->err > maxi_err_tmp)
{
maxi_err_tmp = out_tris[t]->hosted_dem[0]->err;
maxi_err_dem = out_tris[t]->hosted_dem[0];
}
}
while (maxi_err_dem->err >= maxi_err) // quit til the threshold is meet
{
if (err_records != nullptr)
{
err_records->push_back(maxi_err_dem->err);
}
// find the triangle that includes maxi_err_dem and remove it from out_tris
for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
{
old_tri = *t_iter;
if (old_tri == maxi_err_dem->host)
{
t_iter = out_tris.erase(t_iter);
break;
}
else t_iter++;
}
// remove maxi_err_dem from its host triangle's hosted DEM list
for (d_iter = old_tri->hosted_dem.begin(); d_iter != old_tri->hosted_dem.end(); )
{
if (maxi_err_dem == *d_iter)
{
d_iter = old_tri->hosted_dem.erase(d_iter);
break;
}
else d_iter++;
}
// create a new vertex
tmp_vert = new vertex2dc(maxi_err_dem->x, maxi_err_dem->y, maxi_err_dem->elev, out_verts.size());
out_verts.push_back(tmp_vert);
// Delete maxi_err_dem
delete maxi_err_dem; maxi_err_dem = nullptr;
// build new triangles
tmp_tri = split_triangle(tmp_vert, old_tri, cnst_tri);
for (int n = 0; n < 4; ++n)
{
if (cnst_tri[n] != nullptr)
{
out_tris.push_back(cnst_tri[n]);
}
}
if (tmp_tri != nullptr)
{
for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
{
if (tmp_tri == *t_iter)
{
t_iter = out_tris.erase(t_iter);
break;
}
else t_iter++;
}
// build hosted dem for the new triangles
for (int d = 0; d < old_tri->hosted_dem.size(); d++)
{
tmp_dem = old_tri->hosted_dem[d];
for (int n = 0; n < 4; n++)
{
if (cnst_tri[n] != nullptr && cnst_tri[n]->bound_location(tmp_dem->x, tmp_dem->y))
{
tmp_dem->host = cnst_tri[n];
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);
break;
}
}
}
for (int d = 0; d < tmp_tri->hosted_dem.size(); d++)
{
tmp_dem = tmp_tri->hosted_dem[d];
for (int n = 0; n < 4; n++)
{
if (cnst_tri[n] != nullptr && cnst_tri[n]->bound_location(tmp_dem->x, tmp_dem->y))
{
tmp_dem->host = cnst_tri[n];
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);
break;
}
}
}
for (int n = 0; n < 4; n++)
{
if (cnst_tri[n] != nullptr)
{
std::sort(cnst_tri[n]->hosted_dem.begin(), cnst_tri[n]->hosted_dem.end(), compare_dem_point);
}
}
// delete the old triangle
old_tri->hosted_dem.clear();
delete old_tri; old_tri = nullptr;
tmp_tri->hosted_dem.clear();
delete tmp_tri; tmp_tri = nullptr;
}
else
{
// build hosted dem for the new triangles
for (int d = 0; d < old_tri->hosted_dem.size(); d++)
{
tmp_dem = old_tri->hosted_dem[d];
for (int n = 0; n < 4; n++)
{
if (cnst_tri[n] != nullptr && cnst_tri[n]->bound_location(tmp_dem->x, tmp_dem->y))
{
tmp_dem->host = cnst_tri[n];
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);
break;
}
}
}
for (int n = 0; n < 4; n++)
{
if (cnst_tri[n] != nullptr)
{
std::sort(cnst_tri[n]->hosted_dem.begin(), cnst_tri[n]->hosted_dem.end(), compare_dem_point);
}
}
// delete the old triangle
old_tri->hosted_dem.clear();
delete old_tri; old_tri = nullptr;
}
// Make sure cnst_tri meet the empty circumcircle condition
for (int n = 0; n < 4; ++n)
{
if (cnst_tri[n] != nullptr)
{
make_delaunay(cnst_tri[n]);
}
}
// get maximal errors from out_tris
maxi_err_tmp = -1.0;
for (int t = 0; t < out_tris.size(); t++)
{
if (!out_tris[t]->hosted_dem.empty() && out_tris[t]->hosted_dem[0]->err > maxi_err_tmp)
{
maxi_err_tmp = out_tris[t]->hosted_dem[0]->err;
maxi_err_dem = out_tris[t]->hosted_dem[0];
}
}
}
if (err_records != nullptr)
{
err_records->push_back(maxi_err_dem->err);
}
// 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++;
}
// Cut outline if there is one
if (outline_poly != nullptr)
{
for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
{
tmp_tri = *t_iter;
if (!triangle_inside_polygon(tmp_tri, outline_poly))
{
// set neighbors
for (int i = 0; i < 3; ++i)
{
if (tmp_tri->neigh[i] != nullptr)
{
for (int j = 0; j < 3; ++j)
{
if (tmp_tri->neigh[i]->neigh[j] == tmp_tri)
{
tmp_tri->neigh[i]->neigh[j] = 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 remaining DEM data
for (int d = 0; d < out_tris[i]->hosted_dem.size(); d++)
{
tmp_dem = out_tris[i]->hosted_dem[d];
delete tmp_dem; tmp_dem = nullptr;
}
}
// destroy memories located for box_vert
for (int i = 0; i < 4; ++i)
{
delete box_vert[i]; box_vert[i] = nullptr;
}
return;
}

27
src/lib/tin.h
View File

@@ -51,7 +51,7 @@ struct vertex2dc
};
/**
* @brief Compare two vertexes
* @brief Compare two vertexes. Index of the vertexes are not used for comparsion
*
* @param a vertex a
* @param b vertex b
@@ -92,7 +92,7 @@ struct dem_point
double x, y; // position of the DEM location
double elev; // elevation at the DEM location
double err; // error of the TIN with respect to the elevation
triangle *host;
triangle *host; // pointer of the triangle that the dem_point falls inside of
dem_point();
@@ -185,7 +185,17 @@ struct triangle
// End triangle definition
/**
* @brief Generate the TIN from the DEM grid
* @brief Test if the input triangle is inside of the given polygon
*
* @param tri_p Pointer of a test triangle
* @param poly_vert Vertexes of a polygon
* @return true The test triangle is inside of the polygon
* @return false The test triangle is outside of the polygon
*/
bool triangle_inside_polygon(triangle *tri_p, std::vector<vertex2dc> *poly_vert);
/**
* @brief Generate the TIN from a dense 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
@@ -197,24 +207,25 @@ struct triangle
* @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] outline_poly If this pointer is not NULL, Cut triangle outside of or intersected with the polygon.
* @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);
double maxi_err = 1e-0, std::vector<vertex2dc> *outline_poly = nullptr, std::vector<double> *err_records = nullptr);
/**
* @brief Generate the TIN from random DEM points
*
* @param dem Input DEM points
* @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 dem Input DEM points
* @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] outline_poly If this pointer is not NULL, Cut triangle outside the polygon.
* @param[in] outline_poly If this pointer is not NULL, Cut triangle outside of or intersected with the polygon.
* @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, double maxi_err = 1e-0, std::vector<dem_point> *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);
#endif // _TIN_DELAUNAY_H