update src
This commit is contained in:
parent
83660da4d1
commit
209fca546e
61
demo.cpp
61
demo.cpp
@ -1,8 +1,67 @@
|
||||
#include "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("topo.txt");
|
||||
for (int i = 0; i < 10201; ++i)
|
||||
{
|
||||
infile >> topo[i];
|
||||
}
|
||||
infile.close();
|
||||
|
||||
std::vector<vertex2dc*> tin_vert;
|
||||
std::vector<triangle*> tin_ele;
|
||||
dem2tin(topo, 0, 1000, 0, 1000, 10, 10, tin_vert, tin_ele, 3.0);
|
||||
|
||||
// Write a Gmsh's .msh file
|
||||
std::ofstream outfile("topo_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();
|
||||
|
||||
// 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;
|
||||
}
|
||||
188
tin.h
188
tin.h
@ -12,6 +12,8 @@
|
||||
#include "cmath"
|
||||
#include "vector"
|
||||
|
||||
#include "iostream"
|
||||
|
||||
#define ZERO 1e-5
|
||||
|
||||
// Start vertex definition
|
||||
@ -19,12 +21,13 @@ 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), 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)
|
||||
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; id = inid;
|
||||
x = inx; y = iny; elev = inelev; id = inid;
|
||||
return;
|
||||
}
|
||||
};
|
||||
@ -86,81 +89,136 @@ struct triangle
|
||||
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)
|
||||
{
|
||||
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)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
double interpolate(double inx, double iny)
|
||||
{
|
||||
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
|
||||
|
||||
/**
|
||||
* @brief 2D Delaunay triangulation of some given points using the Bowyer-Watson algorithm.
|
||||
* @brief Generate the TIN from the DEM grid
|
||||
*
|
||||
* @param in_verts Input vertexes. Defined by the user.
|
||||
* @param out_tris Output triangles. Compute by the function.
|
||||
* @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-2
|
||||
*/
|
||||
void triangulation(std::vector<vertex2dc> &in_verts, std::vector<triangle> &out_tris)
|
||||
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-2)
|
||||
{
|
||||
if (!out_tris.empty()) out_tris.clear();
|
||||
if (in_verts.size() < 3) return;
|
||||
if (!out_verts.empty()) out_verts.clear();
|
||||
if (!out_tris.empty()) out_tris.clear();
|
||||
|
||||
// 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);
|
||||
}
|
||||
if (xmin >= xmax || ymin >= ymax || (xmin + dx) > xmax || (ymin + dy) > ymax) return;
|
||||
|
||||
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
|
||||
int xnum = round((xmax - xmin)/dx) + 1;
|
||||
int ynum = round((ymax - ymin)/dy) + 1;
|
||||
|
||||
if (dem.size() != xnum*ynum) return;
|
||||
|
||||
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(xmin, ymin, dem[0], out_verts.size()); // lower left corner
|
||||
out_verts.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(xmax, ymin, dem[xnum-1], out_verts.size()); // lower right corner
|
||||
out_verts.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(xmax, ymax, dem[xnum*ynum-1], out_verts.size()); // upper right corner
|
||||
out_verts.push_back(tmp_vert);
|
||||
|
||||
tmp_vert = new vertex2dc(midx - maxi_s, midy + maxi_s); // upper left corner
|
||||
assit_vert.push_back(tmp_vert);
|
||||
tmp_vert = new vertex2dc(xmin, ymax, dem[xnum*(ynum-1)], out_verts.size()); // upper left corner
|
||||
out_verts.push_back(tmp_vert);
|
||||
|
||||
triangle *tmp_tri = nullptr;
|
||||
std::vector<triangle*> exist_tri, cnst_tri;
|
||||
std::vector<triangle*> 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(out_verts[0], out_verts[1], out_verts[2]); // order the vertex anti-clock wise
|
||||
out_tris.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);
|
||||
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);
|
||||
|
||||
int now_maxi_id;
|
||||
double now_x, now_y, now_err;
|
||||
double now_maxi_err;
|
||||
|
||||
// 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)
|
||||
|
||||
do // quit til the threshold is meet
|
||||
{
|
||||
// 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
|
||||
// loop all DEM data to find the location with maximal error
|
||||
// this part is very time consuming. We will fix it later
|
||||
now_maxi_err = -1.0;
|
||||
for (int i = 0; i < xnum*ynum; ++i)
|
||||
{
|
||||
now_x = (i%xnum)*dx + xmin;
|
||||
now_y = (i/xnum)*dy + ymin;
|
||||
for (int e = 0; e < out_tris.size(); ++e)
|
||||
{
|
||||
if (out_tris[e]->bound_location(now_x, now_y))
|
||||
{
|
||||
now_err = fabs(out_tris[e]->interpolate(now_x, now_y) - dem[i]);
|
||||
if (now_err > now_maxi_err)
|
||||
{
|
||||
now_maxi_err = now_err;
|
||||
now_maxi_id = i;
|
||||
}
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// create a new vertex
|
||||
now_x = (now_maxi_id%xnum)*dx + xmin;
|
||||
now_y = (now_maxi_id/xnum)*dy + ymin;
|
||||
tmp_vert = new vertex2dc(now_x, now_y, dem[now_maxi_id], out_verts.size());
|
||||
out_verts.push_back(tmp_vert);
|
||||
|
||||
// determine triangles that include the point and add the triangle to the cnst_tri and remove the triangle from out_tris
|
||||
// this is also a 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(); )
|
||||
for (t_iter = out_tris.begin(); t_iter != out_tris.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);
|
||||
|
||||
dist = (tmp_tri->cx - now_x) * (tmp_tri->cx - now_x) + (tmp_tri->cy - now_y) * (tmp_tri->cy - now_y);
|
||||
if ((dist - tmp_tri->cr) <= ZERO) // Points on the circumcircle are also included
|
||||
{
|
||||
t_iter = exist_tri.erase(t_iter);
|
||||
t_iter = out_tris.erase(t_iter);
|
||||
cnst_tri.push_back(tmp_tri);
|
||||
}
|
||||
else t_iter++;
|
||||
@ -193,11 +251,11 @@ void triangulation(std::vector<vertex2dc> &in_verts, std::vector<triangle> &out_
|
||||
}
|
||||
}
|
||||
|
||||
// construct new triangles and add to exist_tri
|
||||
// construct new triangles and add to out_tris
|
||||
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);
|
||||
tmp_tri = new triangle(cnst_edge[c].vert[0], cnst_edge[c].vert[1], out_verts.back()); // order the vertex anti-clock wise
|
||||
out_tris.push_back(tmp_tri);
|
||||
}
|
||||
|
||||
// destroy memories used by cnst_edge
|
||||
@ -206,36 +264,8 @@ void triangulation(std::vector<vertex2dc> &in_verts, std::vector<triangle> &out_
|
||||
tmp_tri = cnst_tri[c];
|
||||
delete tmp_tri; tmp_tri = nullptr;
|
||||
}
|
||||
}
|
||||
} while (now_maxi_err >= maxi_err);
|
||||
|
||||
// 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;
|
||||
}
|
||||
|
||||
|
||||
2917
topo_TIN.msh
Normal file
2917
topo_TIN.msh
Normal file
File diff suppressed because it is too large
Load Diff
Loading…
Reference in New Issue
Block a user