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

update src

Browse Source
This commit is contained in:
张壹
2021-09-16 23:09:12 +08:00
parent f5e1e808e1
commit d6984a3270
9 changed files with 1010370 additions and 393 deletions
Download Patch File Download Diff File Expand all files Collapse all files

168
tin.h
View File

@@ -131,6 +131,7 @@ struct dem_point
{
double x, y; // position of the DEM location
double elev; // elevation at the DEM location
double err;
triangle *host; // host triangle of the DEM location
std::vector<triangle*> circum_host; // triangles which circumcircles include the location
@@ -139,10 +140,60 @@ struct dem_point
void set(double inx, double iny, double inelev)
{
x = inx; y = iny; elev = inelev; host = nullptr;
circum_host.clear();
return;
}
};
/**
* @brief Utility function of the heap_sort function.
*
* @param a Input vector
* @param[in] i vector's index i
* @param[in] n vector's index n
*/
void update_heap(std::vector<dem_point*> &a, int i, int n)
{
int iMax = i, iLeft = 2 * i + 1, iRight = 2 * (i + 1);
if (iLeft < n && a[iMax]->err > a[iLeft]->err)
{
iMax = iLeft;
}
if (iRight < n && a[iMax]->err > a[iRight]->err)
{
iMax = iRight;
}
if (iMax != i)
{
dem_point *tmp = a[iMax]; a[iMax] = a[i]; a[i] = tmp;
update_heap(a, iMax, n);
}
return;
}
/**
* @brief Heap sort of the dem_point vector in a descending order with respect to the error values
*
* @param a Input vector
*/
void heap_sort(std::vector<dem_point*> &a)
{
int n = a.size();
for (int i = (n - 1) / 2; i >= 0; i--)
{
update_heap(a, i, n);
}
dem_point *tmp;
for (int i = n - 1; i > 0; --i)
{
tmp = a[i]; a[i] = a[0]; a[0] = tmp;
update_heap(a, 0, i);
}
return;
}
// End DEM definition
/**
@@ -177,33 +228,46 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
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;
dem_point *tmp_dem;
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]);
dem_grid[j + i*xnum] = new dem_point(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
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);
d_iter = dem_grid.begin();
tmp_dem = *d_iter; delete tmp_dem;
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
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*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
d_iter = dem_grid.begin() + (xnum - 2);
tmp_dem = *d_iter; delete tmp_dem;
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-1) - 2); dem_grid.erase(d_iter);
d_iter = dem_grid.begin() + (xnum*ynum - 3);
tmp_dem = *d_iter; delete tmp_dem;
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);
tmp_dem = *d_iter; delete tmp_dem;
dem_grid.erase(d_iter);
triangle *tmp_tri = nullptr;
std::vector<triangle*> cnst_tri, new_tri;
@@ -226,9 +290,9 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
{
for (int t = 0; t < out_tris.size(); ++t)
{
if (out_tris[t]->bound_location(dem_grid[i].x, dem_grid[i].y))
if (out_tris[t]->bound_location(dem_grid[i]->x, dem_grid[i]->y))
{
dem_grid[i].host = out_tris[t];
dem_grid[i]->host = out_tris[t];
break; // already found, no need to search more
}
}
@@ -240,52 +304,45 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
{
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);
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]);
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;
// loop all DEM data to find the location with maximal error
for (int i = 0; i < dem_grid.size(); ++i)
{
dem_grid[i]->err = fabs(dem_grid[i]->host->interpolate(dem_grid[i]->x, dem_grid[i]->y) - dem_grid[i]->elev);
}
heap_sort(dem_grid);
bool removed;
edge tmp_edge;
std::vector<edge> cnst_edge;
std::vector<edge>::iterator e_iter;
do // quit til the threshold is meet
while (dem_grid[0]->err >= maxi_err) // 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);
err_records->push_back(dem_grid[0]->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());
tmp_vert = new vertex2dc(dem_grid[0]->x, dem_grid[0]->y, dem_grid[0]->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)
for (int i = 0; i < dem_grid[0]->circum_host.size(); ++i)
{
cnst_tri.push_back(dem_grid[now_maxi_id].circum_host[i]);
cnst_tri.push_back(dem_grid[0]->circum_host[i]);
}
for (int i = 0; i < cnst_tri.size(); ++i)
@@ -303,9 +360,10 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
}
// 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);
dem_grid[0]->circum_host.clear();
d_iter = dem_grid.begin();
tmp_dem = *d_iter; delete tmp_dem;
dem_grid.erase(d_iter);
// loop to remove duplicate edges
cnst_edge.clear();
@@ -351,11 +409,11 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
{
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(); )
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);
t_iter = dem_grid[i]->circum_host.erase(t_iter);
break; // no need to search more
}
else t_iter++;
@@ -368,9 +426,10 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
{
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))
if (new_tri[n]->bound_location(dem_grid[i]->x, dem_grid[i]->y))
{
dem_grid[i].host = new_tri[n];
dem_grid[i]->host = new_tri[n];
dem_grid[i]->err = fabs(new_tri[n]->interpolate(dem_grid[i]->x, dem_grid[i]->y) - dem_grid[i]->elev);
break; // already found, no need to search more
}
}
@@ -381,11 +440,11 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
{
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);
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]);
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
}
}
@@ -398,7 +457,20 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
delete tmp_tri; tmp_tri = nullptr;
}
} while (now_maxi_err >= maxi_err);
heap_sort(dem_grid);
}
if (err_records != nullptr)
{
err_records->push_back(dem_grid[0]->err);
}
// destroy remaining DEM data
for (int i = 0; i < dem_grid.size(); ++i)
{
tmp_dem = dem_grid[i];
delete tmp_dem; tmp_dem = nullptr;
}
return;
}