complete tin
This commit is contained in:
parent
f705e487c7
commit
85f2c3bf1d
2
demo.cpp
2
demo.cpp
@ -1,4 +1,4 @@
|
|||||||
#include "tin_backup.h"
|
#include "tin.h"
|
||||||
#include "iostream"
|
#include "iostream"
|
||||||
#include "fstream"
|
#include "fstream"
|
||||||
#include "iomanip"
|
#include "iomanip"
|
||||||
|
|||||||
785
tin.h
785
tin.h
@ -13,8 +13,6 @@
|
|||||||
#include "vector"
|
#include "vector"
|
||||||
#include "algorithm"
|
#include "algorithm"
|
||||||
|
|
||||||
#include "iostream"
|
|
||||||
|
|
||||||
#define ZERO 1e-5
|
#define ZERO 1e-5
|
||||||
|
|
||||||
// Start vertex definition
|
// Start vertex definition
|
||||||
@ -51,56 +49,77 @@ bool is_collinear(vertex2dc *a_ptr, vertex2dc *b_ptr, vertex2dc *c_ptr) // Test
|
|||||||
}
|
}
|
||||||
return false;
|
return false;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
void circumcircle(vertex2dc *v0, vertex2dc *v1, vertex2dc *v2, double &cx, double &cy, double &cr) // calculate the circumcircle from three points
|
||||||
|
{
|
||||||
|
double s = 0.5 / ((v1->x - v0->x) * (v2->y - v0->y) - (v1->y - v0->y) * (v2->x - v0->x));
|
||||||
|
double m = v1->x*v1->x - v0->x*v0->x + v1->y*v1->y - v0->y*v0->y;
|
||||||
|
double u = v2->x*v2->x - v0->x*v0->x + v2->y*v2->y - v0->y*v0->y;
|
||||||
|
|
||||||
|
cx = ((v2->y - v0->y)*m + (v0->y - v1->y)*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
|
||||||
|
return;
|
||||||
|
}
|
||||||
// End vertex definition
|
// End vertex definition
|
||||||
|
|
||||||
// Start edge definition
|
// Start DEM definition
|
||||||
struct edge
|
struct triangle;
|
||||||
{
|
|
||||||
vertex2dc *vert[2]; // vertex of the edge
|
|
||||||
|
|
||||||
edge() {vert[0] = vert[1] = nullptr;}
|
struct dem_point
|
||||||
edge(vertex2dc *v0ptr, vertex2dc *v1ptr) {set(v0ptr, v1ptr);}
|
{
|
||||||
void set(vertex2dc *v0ptr, vertex2dc *v1ptr)
|
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;
|
||||||
|
|
||||||
|
dem_point() : x(NAN), y(NAN), elev(NAN), err(0.0), host(nullptr) {}
|
||||||
|
dem_point(double inx, double iny, double inelev) {set(inx, iny, inelev);}
|
||||||
|
void set(double inx, double iny, double inelev)
|
||||||
{
|
{
|
||||||
vert[0] = v0ptr; vert[1] = v1ptr;
|
x = inx; y = iny; elev = inelev; err = 0.0; host = nullptr;
|
||||||
return;
|
return;
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
bool operator ==(const edge &a, const edge &b) // overload the == operator for edge type
|
bool compare_dem_point(dem_point *a, dem_point *b) // determination function for std::sort
|
||||||
{
|
{
|
||||||
if((a.vert[0] == b.vert[0] && a.vert[1] == b.vert[1]) ||
|
if (a->err > b->err) return true;
|
||||||
(a.vert[0] == b.vert[1] && a.vert[1] == b.vert[0]))
|
|
||||||
{
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
return false;
|
return false;
|
||||||
}
|
}
|
||||||
// End edge definition
|
// End DEM definition
|
||||||
|
|
||||||
// Start triangle definition
|
|
||||||
struct dem_point;
|
|
||||||
|
|
||||||
|
/* Start triangle definition
|
||||||
|
* v2
|
||||||
|
* /\
|
||||||
|
* / \
|
||||||
|
* n2 / \ n1
|
||||||
|
* / \
|
||||||
|
* /------------\
|
||||||
|
* v0 n0 v1
|
||||||
|
*/
|
||||||
struct triangle
|
struct triangle
|
||||||
{
|
{
|
||||||
|
int id;
|
||||||
vertex2dc *vert[3]; // vertex of the triangle
|
vertex2dc *vert[3]; // vertex of the triangle
|
||||||
|
triangle *neigh[3]; // neighbors of the triangle
|
||||||
double cx, cy; // center of the triangle's circumcircle
|
double cx, cy; // center of the triangle's circumcircle
|
||||||
double cr; // radius of the circumcircle
|
double cr; // radius of the circumcircle
|
||||||
std::vector<dem_point*> circum_dem;
|
std::vector<dem_point*> hosted_dem;
|
||||||
|
|
||||||
triangle() {vert[0] = vert[1] = vert[2] = nullptr;}
|
triangle() {vert[0] = vert[1] = vert[2] = nullptr; neigh[0] = neigh[1] = neigh[2] = nullptr;}
|
||||||
triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) {set(v0ptr, v1ptr, v2ptr);}
|
triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) {set(v0ptr, v1ptr, v2ptr);}
|
||||||
void set(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr)
|
void set(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);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
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));
|
void set_neighbor(triangle *n0ptr, triangle *n1ptr, triangle *n2ptr)
|
||||||
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;
|
neigh[0] = n0ptr; neigh[1] = n1ptr; neigh[2] = n2ptr;
|
||||||
|
|
||||||
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;
|
return;
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -130,32 +149,288 @@ struct triangle
|
|||||||
return (a1*vert[0]->elev + a2*vert[1]->elev + a3*vert[2]->elev)/(a1 + a2 + a3);
|
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
|
* @brief Flip neighboring triangles and their neighbors
|
||||||
|
*
|
||||||
|
* original
|
||||||
|
*
|
||||||
|
* /\
|
||||||
|
* / \
|
||||||
|
* / \
|
||||||
|
* / t \
|
||||||
|
* t_id-------\ t_id (0, 1 or 2)
|
||||||
|
* \--------/
|
||||||
|
* \ /
|
||||||
|
* \ n /
|
||||||
|
* \ /
|
||||||
|
* \/
|
||||||
|
* n_id (0, 1 or 2)
|
||||||
|
*
|
||||||
|
* flipped
|
||||||
|
*
|
||||||
|
* /|\
|
||||||
|
* / | \
|
||||||
|
* / | \
|
||||||
|
* / | \
|
||||||
|
* t_id | \ t_id (0, 1 or 2)
|
||||||
|
* \ t | n /
|
||||||
|
* \ | /
|
||||||
|
* \ | /
|
||||||
|
* \ | /
|
||||||
|
* \|/
|
||||||
|
* n_id (0, 1 or 2)
|
||||||
|
*
|
||||||
|
* @param t target triangle
|
||||||
|
* @param n neighboring triangle
|
||||||
|
* @param t_vid reference index of the target triangle
|
||||||
|
* @param n_vid reference index of the neighboring triangle
|
||||||
|
*/
|
||||||
|
void flip_neighboring_triangles(triangle *t, triangle *n, int t_id, int n_id)
|
||||||
{
|
{
|
||||||
double x, y; // position of the DEM location
|
t->vert[(t_id+1)%3] = n->vert[n_id]; // flip t
|
||||||
double elev; // elevation at the DEM location
|
circumcircle(t->vert[0], t->vert[1], t->vert[2], t->cx, t->cy, t->cr); // update circumcircle
|
||||||
double err; // error of the TIN with respect to the elevation
|
|
||||||
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) {}
|
n->vert[(n_id+2)%3] = t->vert[(t_id+2)%3]; // flip n
|
||||||
dem_point(double inx, double iny, double inelev) {set(inx, iny, inelev);}
|
circumcircle(n->vert[0], n->vert[1], n->vert[2], n->cx, n->cy, n->cr); // update circumcircle
|
||||||
void set(double inx, double iny, double inelev)
|
|
||||||
|
// set side neighbors
|
||||||
|
t->neigh[t_id] = n->neigh[(n_id+2)%3];
|
||||||
|
n->neigh[(n_id+1)%3] = t->neigh[(t_id+1)%3];
|
||||||
|
|
||||||
|
// set opposite neighbors
|
||||||
|
t->neigh[(t_id+1)%3] = n;
|
||||||
|
n->neigh[(n_id+2)%3] = t;
|
||||||
|
|
||||||
|
// set oppsite neighbors
|
||||||
|
if (t->neigh[t_id] != nullptr)
|
||||||
{
|
{
|
||||||
x = inx; y = iny; elev = inelev; host = nullptr;
|
for (int i = 0; i < 3; i++)
|
||||||
return;
|
{
|
||||||
|
if (t->neigh[t_id]->neigh[i] == n)
|
||||||
|
{
|
||||||
|
t->neigh[t_id]->neigh[i] = t;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
if (n->neigh[(n_id+1)%3] != nullptr)
|
||||||
|
{
|
||||||
|
for (int i = 0; i < 3; i++)
|
||||||
|
{
|
||||||
|
if (n->neigh[(n_id+1)%3]->neigh[i] == t)
|
||||||
|
{
|
||||||
|
n->neigh[(n_id+1)%3]->neigh[i] = n;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
}
|
}
|
||||||
};
|
|
||||||
|
|
||||||
bool compare_dem_point(dem_point *a, dem_point *b)
|
// move hosted DEM points
|
||||||
{
|
dem_point *tmp_dem;
|
||||||
if (a->err > b->err) return true;
|
std::vector<dem_point*>::iterator d_iter;
|
||||||
return false;
|
for (d_iter = t->hosted_dem.begin(); d_iter != t->hosted_dem.end(); )
|
||||||
|
{
|
||||||
|
tmp_dem = *d_iter;
|
||||||
|
if (n->bound_location(tmp_dem->x, tmp_dem->y))
|
||||||
|
{
|
||||||
|
tmp_dem->host = n;
|
||||||
|
n->hosted_dem.push_back(tmp_dem);
|
||||||
|
d_iter = t->hosted_dem.erase(d_iter);
|
||||||
|
}
|
||||||
|
else d_iter++;
|
||||||
|
}
|
||||||
|
|
||||||
|
for (d_iter = n->hosted_dem.begin(); d_iter != n->hosted_dem.end(); )
|
||||||
|
{
|
||||||
|
tmp_dem = *d_iter;
|
||||||
|
if (t->bound_location(tmp_dem->x, tmp_dem->y))
|
||||||
|
{
|
||||||
|
tmp_dem->host = t;
|
||||||
|
t->hosted_dem.push_back(tmp_dem);
|
||||||
|
d_iter = n->hosted_dem.erase(d_iter);
|
||||||
|
}
|
||||||
|
else d_iter++;
|
||||||
|
}
|
||||||
|
|
||||||
|
// update errors for hosted DEM data
|
||||||
|
for (int i = 0; i < n->hosted_dem.size(); i++)
|
||||||
|
{
|
||||||
|
tmp_dem = n->hosted_dem[i];
|
||||||
|
tmp_dem->err = fabs(n->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
|
||||||
|
}
|
||||||
|
|
||||||
|
for (int i = 0; i < t->hosted_dem.size(); i++)
|
||||||
|
{
|
||||||
|
tmp_dem = t->hosted_dem[i];
|
||||||
|
tmp_dem->err = fabs(t->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
|
||||||
|
}
|
||||||
|
|
||||||
|
// Sort maximal errors for triangles t and n
|
||||||
|
std::sort(t->hosted_dem.begin(), t->hosted_dem.end(), compare_dem_point);
|
||||||
|
std::sort(n->hosted_dem.begin(), n->hosted_dem.end(), compare_dem_point);
|
||||||
|
return;
|
||||||
}
|
}
|
||||||
// End DEM definition
|
|
||||||
|
/**
|
||||||
|
* @brief Make sure that the input triangle meets the empty circumcircle condition
|
||||||
|
*
|
||||||
|
* @param t Input triangle
|
||||||
|
*/
|
||||||
|
void make_delaunay(triangle *t)
|
||||||
|
{
|
||||||
|
double dist;
|
||||||
|
vertex2dc *n_vert;
|
||||||
|
triangle *n_tri;
|
||||||
|
dem_point *tmp_dem;
|
||||||
|
for (int n = 0; n < 3; n++)
|
||||||
|
{
|
||||||
|
if (t->neigh[n] != nullptr) // must has neighbor on this side
|
||||||
|
{
|
||||||
|
n_tri = t->neigh[n];
|
||||||
|
for (int v = 0; v < 3; v++)
|
||||||
|
{
|
||||||
|
n_vert = n_tri->vert[v];
|
||||||
|
if (n_vert != t->vert[n] && n_vert != t->vert[(n+1)%3]) // find the opposite vertex
|
||||||
|
{
|
||||||
|
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
|
||||||
|
{
|
||||||
|
flip_neighboring_triangles(t, n_tri, n, v);
|
||||||
|
|
||||||
|
// Make sure the triangles also meet the empty circumcircle condition after flipping
|
||||||
|
make_delaunay(t);
|
||||||
|
make_delaunay(n_tri);
|
||||||
|
return; // Neighborhood changed. The current loop is not valid any more.
|
||||||
|
}
|
||||||
|
break; // no need to search more
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
triangle *split_triangle(vertex2dc *v, triangle *t, triangle *new_t[4])
|
||||||
|
{
|
||||||
|
vertex2dc *tmp_vert;
|
||||||
|
triangle *tmp_tri;
|
||||||
|
|
||||||
|
new_t[0] = new_t[1] = new_t[2] = new_t[3] = nullptr;
|
||||||
|
|
||||||
|
// Check for collinear
|
||||||
|
for (int i = 0; i < 3; i++)
|
||||||
|
{
|
||||||
|
if (is_collinear(t->vert[i], t->vert[(i+1)%3], v)) // the new vertex is on edge
|
||||||
|
{
|
||||||
|
if (t->neigh[i] == nullptr) // no neighboring triangle. create two new triangles
|
||||||
|
{
|
||||||
|
tmp_tri = new triangle(t->vert[i], v, t->vert[(i+2)%3]); new_t[0] = tmp_tri;
|
||||||
|
tmp_tri = new triangle(t->vert[(i+2)%3], v, t->vert[(i+1)%3]); new_t[1] = tmp_tri;
|
||||||
|
new_t[0]->set_neighbor(nullptr, new_t[1], t->neigh[(i+2)%3]);
|
||||||
|
new_t[1]->set_neighbor(new_t[0], nullptr, t->neigh[(i+1)%3]);
|
||||||
|
|
||||||
|
for (int n = 0; n < 2; n++)
|
||||||
|
{
|
||||||
|
if (new_t[n]->neigh[2] != nullptr)
|
||||||
|
{
|
||||||
|
for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle
|
||||||
|
{
|
||||||
|
if (new_t[n]->neigh[2]->neigh[k] == t)
|
||||||
|
{
|
||||||
|
new_t[n]->neigh[2]->neigh[k] = new_t[n];
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return nullptr;
|
||||||
|
}
|
||||||
|
|
||||||
|
// has a neighboring triangle. create four new triangles
|
||||||
|
for (int k = 0; k < 3; k++)
|
||||||
|
{
|
||||||
|
tmp_vert = t->neigh[i]->vert[k];
|
||||||
|
if (tmp_vert != t->vert[i] && tmp_vert != t->vert[(i+1)%3])
|
||||||
|
{
|
||||||
|
tmp_tri = new triangle(t->vert[i], v, t->vert[(i+2)%3]); new_t[0] = tmp_tri;
|
||||||
|
tmp_tri = new triangle(t->vert[(i+2)%3], v, t->vert[(i+1)%3]); new_t[1] = tmp_tri;
|
||||||
|
tmp_tri = new triangle(tmp_vert, v, t->vert[i]); new_t[2] = tmp_tri;
|
||||||
|
tmp_tri = new triangle(t->vert[(i+1)%3], v, tmp_vert); new_t[3] = tmp_tri;
|
||||||
|
new_t[0]->set_neighbor(new_t[2], new_t[1], t->neigh[(i+2)%3]);
|
||||||
|
new_t[1]->set_neighbor(new_t[0], new_t[3], t->neigh[(i+1)%3]);
|
||||||
|
new_t[2]->set_neighbor(new_t[3], new_t[0], t->neigh[i]->neigh[(k+2)%3]);
|
||||||
|
new_t[3]->set_neighbor(new_t[1], new_t[2], t->neigh[i]->neigh[k]);
|
||||||
|
|
||||||
|
for (int n = 0; n < 2; n++)
|
||||||
|
{
|
||||||
|
if (new_t[n]->neigh[2] != nullptr)
|
||||||
|
{
|
||||||
|
for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle
|
||||||
|
{
|
||||||
|
if (new_t[n]->neigh[2]->neigh[k] == t)
|
||||||
|
{
|
||||||
|
new_t[n]->neigh[2]->neigh[k] = new_t[n];
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
for (int n = 2; n < 4; n++)
|
||||||
|
{
|
||||||
|
if (new_t[n]->neigh[2] != nullptr)
|
||||||
|
{
|
||||||
|
for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle
|
||||||
|
{
|
||||||
|
if (new_t[n]->neigh[2]->neigh[k] == t->neigh[i])
|
||||||
|
{
|
||||||
|
new_t[n]->neigh[2]->neigh[k] = new_t[n];
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return t->neigh[i]; // Return the neighboring tiangle to be deleted
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// The new vertex is inside the triangle. create three new triangles
|
||||||
|
for (int n = 0; n < 3; ++n)
|
||||||
|
{
|
||||||
|
tmp_tri = new triangle(t->vert[n], t->vert[(n+1)%3], v);
|
||||||
|
new_t[n] = tmp_tri;
|
||||||
|
}
|
||||||
|
|
||||||
|
// sort neighbors for new triangles
|
||||||
|
for (int n = 0; n < 3; ++n)
|
||||||
|
{
|
||||||
|
if (t->neigh[n] == nullptr)
|
||||||
|
{
|
||||||
|
new_t[n]->set_neighbor(nullptr, new_t[(n+1)%3], new_t[(n+2)%3]);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
new_t[n]->set_neighbor(t->neigh[n], new_t[(n+1)%3], new_t[(n+2)%3]);
|
||||||
|
for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle
|
||||||
|
{
|
||||||
|
if (t->neigh[n]->neigh[k] == t)
|
||||||
|
{
|
||||||
|
t->neigh[n]->neigh[k] = new_t[n];
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return nullptr;
|
||||||
|
}
|
||||||
|
// End triangle definition
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* @brief Generate the TIN from the DEM grid
|
* @brief Generate the TIN from the DEM grid
|
||||||
@ -189,274 +464,252 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
|
|||||||
if (dem.size() != xnum*ynum) return;
|
if (dem.size() != xnum*ynum) return;
|
||||||
|
|
||||||
// Prepare the DEM points
|
// Prepare the DEM points
|
||||||
dem_point *tmp_dem = nullptr;;
|
dem_point *tmp_dem = nullptr;
|
||||||
std::vector<dem_point*> dem_grid(xnum*ynum);
|
std::vector<dem_point*> dem_tri;
|
||||||
std::vector<dem_point*>::iterator d_iter;
|
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] = new dem_point(xmin + dx*j, ymin + dy*i, dem[j + i*xnum]);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
vertex2dc *tmp_vert = nullptr;
|
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[0], out_verts.size()); // lower left corner
|
||||||
out_verts.push_back(tmp_vert);
|
out_verts.push_back(tmp_vert);
|
||||||
|
|
||||||
d_iter = dem_grid.begin();
|
tmp_vert = new vertex2dc(xmax, ymin, dem[xnum-1], out_verts.size()); // lower right corner
|
||||||
tmp_dem = *d_iter; delete tmp_dem;
|
|
||||||
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);
|
out_verts.push_back(tmp_vert);
|
||||||
|
|
||||||
d_iter = dem_grid.begin() + (xnum - 2);
|
tmp_vert = new vertex2dc(xmax, ymax, dem[xnum*ynum-1], out_verts.size()); // upper right corner
|
||||||
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);
|
out_verts.push_back(tmp_vert);
|
||||||
|
|
||||||
d_iter = dem_grid.begin() + (xnum*ynum - 3);
|
tmp_vert = new vertex2dc(xmin, ymax, dem[xnum*(ynum-1)], out_verts.size()); // upper left corner
|
||||||
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);
|
out_verts.push_back(tmp_vert);
|
||||||
|
|
||||||
d_iter = dem_grid.begin() + (xnum*(ynum-1) - 2);
|
triangle *old_tri = nullptr, *tmp_tri = nullptr;
|
||||||
tmp_dem = *d_iter; delete tmp_dem;
|
triangle *cnst_tri[4];
|
||||||
dem_grid.erase(d_iter);
|
|
||||||
|
|
||||||
triangle *tmp_tri = nullptr;
|
|
||||||
std::vector<triangle*> cnst_tri, new_tri;
|
|
||||||
std::vector<triangle*>::iterator t_iter;
|
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
|
||||||
{
|
{
|
||||||
tmp_tri = new triangle(out_verts[0], out_verts[1], out_verts[2]); // order the vertex anti-clock wise
|
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);
|
out_tris.push_back(tmp_tri); tmp_tri = nullptr;
|
||||||
}
|
}
|
||||||
|
|
||||||
if (!is_collinear(out_verts[0], out_verts[2], out_verts[3]))
|
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
|
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);
|
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);
|
||||||
|
|
||||||
// Find host triangle for all DEM locations
|
// Find host triangle for all DEM locations
|
||||||
for (int i = 0; i < dem_grid.size(); ++i)
|
int tmp_id;
|
||||||
|
for (int i = 0; i < ynum; ++i)
|
||||||
{
|
{
|
||||||
for (int t = 0; t < out_tris.size(); ++t)
|
for (int j = 0; j < xnum; ++j)
|
||||||
{
|
{
|
||||||
if (out_tris[t]->bound_location(dem_grid[i]->x, dem_grid[i]->y))
|
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
|
||||||
{
|
{
|
||||||
dem_grid[i]->host = out_tris[t];
|
tmp_dem = new dem_point(xmin + dx*j, ymin + dy*i, dem[j + i*xnum]);
|
||||||
break; // already found, no need to search more
|
for (int t = 0; t < out_tris.size(); ++t)
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// 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]);
|
|
||||||
out_tris[t]->circum_dem.push_back(dem_grid[i]);
|
|
||||||
// no beak here. There might be more than one triangle's circumcircle includes the DEM location
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// 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);
|
|
||||||
}
|
|
||||||
|
|
||||||
// Sort dem_grid in the desceding order with respect to the error
|
|
||||||
std::sort(dem_grid.begin(), dem_grid.end(), compare_dem_point);
|
|
||||||
|
|
||||||
bool removed;
|
|
||||||
edge tmp_edge;
|
|
||||||
std::vector<edge> cnst_edge;
|
|
||||||
std::vector<edge>::iterator e_iter;
|
|
||||||
|
|
||||||
while (dem_grid[0]->err >= maxi_err) // quit til the threshold is meet
|
|
||||||
{
|
|
||||||
if (err_records != nullptr)
|
|
||||||
{
|
|
||||||
err_records->push_back(dem_grid[0]->err);
|
|
||||||
}
|
|
||||||
|
|
||||||
// create a new vertex
|
|
||||||
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[0]->circum_host.size(); ++i)
|
|
||||||
{
|
|
||||||
cnst_tri.push_back(dem_grid[0]->circum_host[i]);
|
|
||||||
}
|
|
||||||
|
|
||||||
for (int c = 0; c < cnst_tri.size(); ++c)
|
|
||||||
{
|
|
||||||
for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
|
|
||||||
{
|
|
||||||
tmp_tri = *t_iter;
|
|
||||||
if (cnst_tri[c] == tmp_tri)
|
|
||||||
{
|
{
|
||||||
t_iter = out_tris.erase(t_iter);
|
if (out_tris[t]->bound_location(tmp_dem->x, tmp_dem->y))
|
||||||
break; // no need to search more
|
|
||||||
}
|
|
||||||
else t_iter++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// remove cnst_tri from its circumed DEM's circum triangle list
|
|
||||||
for (int c = 0; c < cnst_tri.size(); ++c)
|
|
||||||
{
|
|
||||||
for (int i = 0; i < cnst_tri[c]->circum_dem.size(); ++i)
|
|
||||||
{
|
|
||||||
tmp_dem = cnst_tri[c]->circum_dem[i];
|
|
||||||
for (t_iter = tmp_dem->circum_host.begin(); t_iter != tmp_dem->circum_host.end(); )
|
|
||||||
{
|
|
||||||
if (cnst_tri[c] == *t_iter)
|
|
||||||
{
|
{
|
||||||
t_iter = tmp_dem->circum_host.erase(t_iter);
|
tmp_dem->host = out_tris[t];
|
||||||
break;
|
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);
|
||||||
else t_iter++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// remove dem_grid[0] from its circumed triangle's circum DEM list
|
|
||||||
for (int c = 0; c < cnst_tri.size(); ++c)
|
|
||||||
{
|
|
||||||
for (d_iter = cnst_tri[c]->circum_dem.begin(); d_iter != cnst_tri[c]->circum_dem.end(); )
|
|
||||||
{
|
|
||||||
if (dem_grid[0] == *d_iter)
|
|
||||||
{
|
|
||||||
d_iter = cnst_tri[c]->circum_dem.erase(d_iter);
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
else d_iter++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// clear host and circumcircle triangles for the used DEM location
|
|
||||||
d_iter = dem_grid.begin();
|
|
||||||
tmp_dem = *d_iter; tmp_dem->circum_host.clear(); delete tmp_dem;
|
|
||||||
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);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// loop all DEM data to update host triangles
|
|
||||||
for (int c = 0; c < cnst_tri.size(); ++c)
|
|
||||||
{
|
|
||||||
for (int i = 0; i < cnst_tri[c]->circum_dem.size(); ++i)
|
|
||||||
{
|
|
||||||
tmp_dem = cnst_tri[c]->circum_dem[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(tmp_dem->x, tmp_dem->y))
|
|
||||||
{
|
|
||||||
tmp_dem->host = new_tri[n];
|
|
||||||
tmp_dem->err = fabs(new_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
|
|
||||||
break; // already found, no need to search more
|
break; // already found, no need to search more
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
}
|
||||||
|
|
||||||
// Find circum_host triangles for all DEM locations
|
// Sort hosted_dem in the desceding order with respect to the error. Add maximal zeros to dem_tri
|
||||||
// cnst_tri's circum area doesn't over cover new_tri's circum area
|
for (int t = 0; t < out_tris.size(); ++t)
|
||||||
for (int i = 0; i < dem_grid.size(); ++i)
|
{
|
||||||
|
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]);
|
||||||
|
}
|
||||||
|
|
||||||
|
// 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
|
||||||
|
{
|
||||||
|
if (err_records != nullptr)
|
||||||
{
|
{
|
||||||
for (int n = 0; n < new_tri.size(); ++n) // search in newly created triangles to find new circumcircle triangles
|
err_records->push_back(dem_tri[0]->err);
|
||||||
|
}
|
||||||
|
|
||||||
|
// find the triangle that includes dem_tri[0] 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)
|
||||||
{
|
{
|
||||||
dist = (new_tri[n]->cx - dem_grid[i]->x) * (new_tri[n]->cx - dem_grid[i]->x)
|
t_iter = out_tris.erase(t_iter);
|
||||||
+ (new_tri[n]->cy - dem_grid[i]->y) * (new_tri[n]->cy - dem_grid[i]->y);
|
break;
|
||||||
if ((dist - new_tri[n]->cr) <= ZERO) // Points on the circumcircle are also included
|
}
|
||||||
|
else t_iter++;
|
||||||
|
}
|
||||||
|
|
||||||
|
// remove dem_tri[0] 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)
|
||||||
|
{
|
||||||
|
d_iter = old_tri->hosted_dem.erase(d_iter);
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
else d_iter++;
|
||||||
|
}
|
||||||
|
|
||||||
|
// create a new vertex
|
||||||
|
tmp_vert = new vertex2dc(dem_tri[0]->x, dem_tri[0]->y, dem_tri[0]->elev, out_verts.size());
|
||||||
|
out_verts.push_back(tmp_vert);
|
||||||
|
|
||||||
|
// Delete dem_tri[0]
|
||||||
|
tmp_dem = dem_tri[0]; delete tmp_dem;
|
||||||
|
|
||||||
|
// 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)
|
||||||
{
|
{
|
||||||
new_tri[n]->circum_dem.push_back(dem_grid[i]);
|
t_iter = out_tris.erase(t_iter);
|
||||||
dem_grid[i]->circum_host.push_back(new_tri[n]);
|
break;
|
||||||
// no beak here. There might be more than one triangle's circumcircle includes the DEM location
|
|
||||||
}
|
}
|
||||||
|
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 and sort dem_tri
|
||||||
|
dem_tri.clear(); dem_tri.reserve(out_tris.size());
|
||||||
|
for (int t = 0; t < out_tris.size(); t++)
|
||||||
|
{
|
||||||
|
if (!out_tris[t]->hosted_dem.empty())
|
||||||
|
{
|
||||||
|
dem_tri.push_back(out_tris[t]->hosted_dem[0]);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
std::sort(dem_tri.begin(), dem_tri.end(), compare_dem_point);
|
||||||
// destroy memories used by cnst_edge
|
|
||||||
for (int c = 0; c < cnst_tri.size(); ++c)
|
|
||||||
{
|
|
||||||
tmp_tri = cnst_tri[c];
|
|
||||||
tmp_tri->circum_dem.clear();
|
|
||||||
delete tmp_tri; tmp_tri = nullptr;
|
|
||||||
}
|
|
||||||
|
|
||||||
// Sort dem_grid in the desceding order with respect to the error
|
|
||||||
std::sort(dem_grid.begin(), dem_grid.end(), compare_dem_point);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
if (err_records != nullptr)
|
if (err_records != nullptr)
|
||||||
{
|
{
|
||||||
err_records->push_back(dem_grid[0]->err);
|
err_records->push_back(dem_tri[0]->err);
|
||||||
}
|
}
|
||||||
|
|
||||||
// destroy remaining DEM data
|
// assign triangles index
|
||||||
for (int i = 0; i < dem_grid.size(); ++i)
|
for (int i = 0; i < out_tris.size(); i++)
|
||||||
{
|
{
|
||||||
tmp_dem = dem_grid[i];
|
out_tris[i]->id = i;
|
||||||
delete tmp_dem; tmp_dem = nullptr;
|
// 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;
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
return;
|
return;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|||||||
791
tin_backup.h
791
tin_backup.h
@ -13,8 +13,6 @@
|
|||||||
#include "vector"
|
#include "vector"
|
||||||
#include "algorithm"
|
#include "algorithm"
|
||||||
|
|
||||||
#include "iostream"
|
|
||||||
|
|
||||||
#define ZERO 1e-5
|
#define ZERO 1e-5
|
||||||
|
|
||||||
// Start vertex definition
|
// Start vertex definition
|
||||||
@ -51,77 +49,56 @@ bool is_collinear(vertex2dc *a_ptr, vertex2dc *b_ptr, vertex2dc *c_ptr) // Test
|
|||||||
}
|
}
|
||||||
return false;
|
return false;
|
||||||
}
|
}
|
||||||
|
|
||||||
void circumcircle(vertex2dc *v0, vertex2dc *v1, vertex2dc *v2, double &cx, double &cy, double &cr) // calculate the circumcircle from three points
|
|
||||||
{
|
|
||||||
double s = 0.5 / ((v1->x - v0->x) * (v2->y - v0->y) - (v1->y - v0->y) * (v2->x - v0->x));
|
|
||||||
double m = v1->x*v1->x - v0->x*v0->x + v1->y*v1->y - v0->y*v0->y;
|
|
||||||
double u = v2->x*v2->x - v0->x*v0->x + v2->y*v2->y - v0->y*v0->y;
|
|
||||||
|
|
||||||
cx = ((v2->y - v0->y)*m + (v0->y - v1->y)*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
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
// End vertex definition
|
// End vertex definition
|
||||||
|
|
||||||
// Start DEM definition
|
// Start edge definition
|
||||||
struct triangle;
|
struct edge
|
||||||
|
|
||||||
struct dem_point
|
|
||||||
{
|
{
|
||||||
double x, y; // position of the DEM location
|
vertex2dc *vert[2]; // vertex of the edge
|
||||||
double elev; // elevation at the DEM location
|
|
||||||
double err; // error of the TIN with respect to the elevation
|
|
||||||
triangle *host;
|
|
||||||
|
|
||||||
dem_point() : x(NAN), y(NAN), elev(NAN), err(0.0), host(nullptr) {}
|
edge() {vert[0] = vert[1] = nullptr;}
|
||||||
dem_point(double inx, double iny, double inelev) {set(inx, iny, inelev);}
|
edge(vertex2dc *v0ptr, vertex2dc *v1ptr) {set(v0ptr, v1ptr);}
|
||||||
void set(double inx, double iny, double inelev)
|
void set(vertex2dc *v0ptr, vertex2dc *v1ptr)
|
||||||
{
|
{
|
||||||
x = inx; y = iny; elev = inelev; err = 0.0; host = nullptr;
|
vert[0] = v0ptr; vert[1] = v1ptr;
|
||||||
return;
|
return;
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
bool compare_dem_point(dem_point *a, dem_point *b) // determination function for std::sort
|
bool operator ==(const edge &a, const edge &b) // overload the == operator for edge type
|
||||||
{
|
{
|
||||||
if (a->err > b->err) return true;
|
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;
|
return false;
|
||||||
}
|
}
|
||||||
// End DEM definition
|
// End edge definition
|
||||||
|
|
||||||
|
// Start triangle definition
|
||||||
|
struct dem_point;
|
||||||
|
|
||||||
/* Start triangle definition
|
|
||||||
* v2
|
|
||||||
* /\
|
|
||||||
* / \
|
|
||||||
* n2 / \ n1
|
|
||||||
* / \
|
|
||||||
* /------------\
|
|
||||||
* v0 n0 v1
|
|
||||||
*/
|
|
||||||
struct triangle
|
struct triangle
|
||||||
{
|
{
|
||||||
int id;
|
|
||||||
vertex2dc *vert[3]; // vertex of the triangle
|
vertex2dc *vert[3]; // vertex of the triangle
|
||||||
triangle *neigh[3]; // neighbors of the triangle
|
|
||||||
double cx, cy; // center of the triangle's circumcircle
|
double cx, cy; // center of the triangle's circumcircle
|
||||||
double cr; // radius of the circumcircle
|
double cr; // radius of the circumcircle
|
||||||
std::vector<dem_point*> hosted_dem;
|
std::vector<dem_point*> circum_dem;
|
||||||
|
|
||||||
triangle() {vert[0] = vert[1] = vert[2] = nullptr; neigh[0] = neigh[1] = neigh[2] = nullptr;}
|
triangle() {vert[0] = vert[1] = vert[2] = nullptr;}
|
||||||
triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) {set(v0ptr, v1ptr, v2ptr);}
|
triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) {set(v0ptr, v1ptr, v2ptr);}
|
||||||
void set(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr)
|
void set(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);
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
void set_neighbor(triangle *n0ptr, triangle *n1ptr, triangle *n2ptr)
|
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;
|
||||||
neigh[0] = n0ptr; neigh[1] = n1ptr; neigh[2] = n2ptr;
|
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;
|
return;
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -151,285 +128,33 @@ struct triangle
|
|||||||
return (a1*vert[0]->elev + a2*vert[1]->elev + a3*vert[2]->elev)/(a1 + a2 + a3);
|
return (a1*vert[0]->elev + a2*vert[1]->elev + a3*vert[2]->elev)/(a1 + a2 + a3);
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
/**
|
|
||||||
* @brief Flip neighboring triangles and their neighbors
|
|
||||||
*
|
|
||||||
* original
|
|
||||||
*
|
|
||||||
* /\
|
|
||||||
* / \
|
|
||||||
* / \
|
|
||||||
* / t \
|
|
||||||
* t_id-------\ t_id (0, 1 or 2)
|
|
||||||
* \--------/
|
|
||||||
* \ /
|
|
||||||
* \ n /
|
|
||||||
* \ /
|
|
||||||
* \/
|
|
||||||
* n_id (0, 1 or 2)
|
|
||||||
*
|
|
||||||
* flipped
|
|
||||||
*
|
|
||||||
* /|\
|
|
||||||
* / | \
|
|
||||||
* / | \
|
|
||||||
* / | \
|
|
||||||
* t_id | \ t_id (0, 1 or 2)
|
|
||||||
* \ t | n /
|
|
||||||
* \ | /
|
|
||||||
* \ | /
|
|
||||||
* \ | /
|
|
||||||
* \|/
|
|
||||||
* n_id (0, 1 or 2)
|
|
||||||
*
|
|
||||||
* @param t target triangle
|
|
||||||
* @param n neighboring triangle
|
|
||||||
* @param t_vid reference index of the target triangle
|
|
||||||
* @param n_vid reference index of the neighboring triangle
|
|
||||||
*/
|
|
||||||
void flip_neighboring_triangles(triangle *t, triangle *n, int t_id, int n_id)
|
|
||||||
{
|
|
||||||
t->vert[(t_id+1)%3] = n->vert[n_id]; // flip t
|
|
||||||
circumcircle(t->vert[0], t->vert[1], t->vert[2], t->cx, t->cy, t->cr); // update circumcircle
|
|
||||||
|
|
||||||
n->vert[(n_id+2)%3] = t->vert[(t_id+2)%3]; // flip n
|
|
||||||
circumcircle(n->vert[0], n->vert[1], n->vert[2], n->cx, n->cy, n->cr); // update circumcircle
|
|
||||||
|
|
||||||
// set side neighbors
|
|
||||||
t->neigh[t_id] = n->neigh[(n_id+2)%3];
|
|
||||||
n->neigh[(n_id+1)%3] = t->neigh[(t_id+1)%3];
|
|
||||||
|
|
||||||
// set opposite neighbors
|
|
||||||
t->neigh[(t_id+1)%3] = n;
|
|
||||||
n->neigh[(n_id+2)%3] = t;
|
|
||||||
|
|
||||||
// set oppsite neighbors
|
|
||||||
if (t->neigh[t_id] != nullptr)
|
|
||||||
{
|
|
||||||
for (int i = 0; i < 3; i++)
|
|
||||||
{
|
|
||||||
if (t->neigh[t_id]->neigh[i] == n)
|
|
||||||
{
|
|
||||||
t->neigh[t_id]->neigh[i] = t;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
if (n->neigh[(n_id+1)%3] != nullptr)
|
|
||||||
{
|
|
||||||
for (int i = 0; i < 3; i++)
|
|
||||||
{
|
|
||||||
if (n->neigh[(n_id+1)%3]->neigh[i] == t)
|
|
||||||
{
|
|
||||||
n->neigh[(n_id+1)%3]->neigh[i] = n;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// move hosted DEM points
|
|
||||||
dem_point *tmp_dem;
|
|
||||||
std::vector<dem_point*>::iterator d_iter;
|
|
||||||
for (d_iter = t->hosted_dem.begin(); d_iter != t->hosted_dem.end(); )
|
|
||||||
{
|
|
||||||
tmp_dem = *d_iter;
|
|
||||||
if (n->bound_location(tmp_dem->x, tmp_dem->y))
|
|
||||||
{
|
|
||||||
tmp_dem->host = n;
|
|
||||||
n->hosted_dem.push_back(tmp_dem);
|
|
||||||
d_iter = t->hosted_dem.erase(d_iter);
|
|
||||||
}
|
|
||||||
else d_iter++;
|
|
||||||
}
|
|
||||||
|
|
||||||
for (d_iter = n->hosted_dem.begin(); d_iter != n->hosted_dem.end(); )
|
|
||||||
{
|
|
||||||
tmp_dem = *d_iter;
|
|
||||||
if (t->bound_location(tmp_dem->x, tmp_dem->y))
|
|
||||||
{
|
|
||||||
tmp_dem->host = t;
|
|
||||||
t->hosted_dem.push_back(tmp_dem);
|
|
||||||
d_iter = n->hosted_dem.erase(d_iter);
|
|
||||||
}
|
|
||||||
else d_iter++;
|
|
||||||
}
|
|
||||||
|
|
||||||
for (int i = 0; i < n->hosted_dem.size(); i++)
|
|
||||||
{
|
|
||||||
tmp_dem = n->hosted_dem[i];
|
|
||||||
tmp_dem->err = fabs(n->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
|
|
||||||
}
|
|
||||||
|
|
||||||
for (int i = 0; i < t->hosted_dem.size(); i++)
|
|
||||||
{
|
|
||||||
tmp_dem = t->hosted_dem[i];
|
|
||||||
tmp_dem->err = fabs(t->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
|
|
||||||
}
|
|
||||||
|
|
||||||
std::sort(t->hosted_dem.begin(), t->hosted_dem.end(), compare_dem_point);
|
|
||||||
std::sort(n->hosted_dem.begin(), n->hosted_dem.end(), compare_dem_point);
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @brief Make sure that the input triangle meets the empty circumcircle condition
|
|
||||||
*
|
|
||||||
* @param t Input triangle
|
|
||||||
*/
|
|
||||||
void make_delaunay(triangle *t)
|
|
||||||
{
|
|
||||||
double dist;
|
|
||||||
vertex2dc *n_vert;
|
|
||||||
triangle *n_tri;
|
|
||||||
dem_point *tmp_dem;
|
|
||||||
for (int n = 0; n < 3; n++)
|
|
||||||
{
|
|
||||||
if (t->neigh[n] != nullptr) // must has neighbor on this side
|
|
||||||
{
|
|
||||||
n_tri = t->neigh[n];
|
|
||||||
for (int v = 0; v < 3; v++)
|
|
||||||
{
|
|
||||||
n_vert = n_tri->vert[v];
|
|
||||||
if (n_vert != t->vert[n] && n_vert != t->vert[(n+1)%3]) // find the opposite vertex
|
|
||||||
{
|
|
||||||
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
|
|
||||||
{
|
|
||||||
flip_neighboring_triangles(t, n_tri, n, v);
|
|
||||||
|
|
||||||
// Make sure the triangles also meet the empty circumcircle condition after flipping
|
|
||||||
make_delaunay(t);
|
|
||||||
make_delaunay(n_tri);
|
|
||||||
return; // Neighborhood changed. The current loop is not valid any more.
|
|
||||||
}
|
|
||||||
break; // no need to search more
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
triangle *split_triangle(vertex2dc *v, triangle *t, triangle *new_t[4])
|
|
||||||
{
|
|
||||||
vertex2dc *tmp_vert;
|
|
||||||
triangle *tmp_tri;
|
|
||||||
|
|
||||||
new_t[0] = new_t[1] = new_t[2] = new_t[3] = nullptr;
|
|
||||||
|
|
||||||
// Check for collinear
|
|
||||||
for (int i = 0; i < 3; i++)
|
|
||||||
{
|
|
||||||
if (is_collinear(t->vert[i], t->vert[(i+1)%3], v))
|
|
||||||
{
|
|
||||||
if (t->neigh[i] == nullptr)
|
|
||||||
{
|
|
||||||
tmp_tri = new triangle(t->vert[i], v, t->vert[(i+2)%3]); new_t[0] = tmp_tri;
|
|
||||||
tmp_tri = new triangle(t->vert[(i+2)%3], v, t->vert[(i+1)%3]); new_t[1] = tmp_tri;
|
|
||||||
new_t[0]->set_neighbor(nullptr, new_t[1], t->neigh[(i+2)%3]);
|
|
||||||
new_t[1]->set_neighbor(new_t[0], nullptr, t->neigh[(i+1)%3]);
|
|
||||||
|
|
||||||
for (int n = 0; n < 2; n++)
|
|
||||||
{
|
|
||||||
if (new_t[n]->neigh[2] != nullptr)
|
|
||||||
{
|
|
||||||
for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle
|
|
||||||
{
|
|
||||||
if (new_t[n]->neigh[2]->neigh[k] == t)
|
|
||||||
{
|
|
||||||
new_t[n]->neigh[2]->neigh[k] = new_t[n];
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
return nullptr;
|
|
||||||
}
|
|
||||||
|
|
||||||
for (int k = 0; k < 3; k++)
|
|
||||||
{
|
|
||||||
tmp_vert = t->neigh[i]->vert[k];
|
|
||||||
if (tmp_vert != t->vert[i] && tmp_vert != t->vert[(i+1)%3])
|
|
||||||
{
|
|
||||||
tmp_tri = new triangle(t->vert[i], v, t->vert[(i+2)%3]); new_t[0] = tmp_tri;
|
|
||||||
tmp_tri = new triangle(t->vert[(i+2)%3], v, t->vert[(i+1)%3]); new_t[1] = tmp_tri;
|
|
||||||
tmp_tri = new triangle(tmp_vert, v, t->vert[i]); new_t[2] = tmp_tri;
|
|
||||||
tmp_tri = new triangle(t->vert[(i+1)%3], v, tmp_vert); new_t[3] = tmp_tri;
|
|
||||||
new_t[0]->set_neighbor(new_t[2], new_t[1], t->neigh[(i+2)%3]);
|
|
||||||
new_t[1]->set_neighbor(new_t[0], new_t[3], t->neigh[(i+1)%3]);
|
|
||||||
new_t[2]->set_neighbor(new_t[3], new_t[0], t->neigh[i]->neigh[(k+2)%3]);
|
|
||||||
new_t[3]->set_neighbor(new_t[1], new_t[2], t->neigh[i]->neigh[k]);
|
|
||||||
|
|
||||||
for (int n = 0; n < 2; n++)
|
|
||||||
{
|
|
||||||
if (new_t[n]->neigh[2] != nullptr)
|
|
||||||
{
|
|
||||||
for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle
|
|
||||||
{
|
|
||||||
if (new_t[n]->neigh[2]->neigh[k] == t)
|
|
||||||
{
|
|
||||||
new_t[n]->neigh[2]->neigh[k] = new_t[n];
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
for (int n = 2; n < 4; n++)
|
|
||||||
{
|
|
||||||
if (new_t[n]->neigh[2] != nullptr)
|
|
||||||
{
|
|
||||||
for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle
|
|
||||||
{
|
|
||||||
if (new_t[n]->neigh[2]->neigh[k] == t->neigh[i])
|
|
||||||
{
|
|
||||||
new_t[n]->neigh[2]->neigh[k] = new_t[n];
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
return t->neigh[i];;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
for (int n = 0; n < 3; ++n)
|
|
||||||
{
|
|
||||||
tmp_tri = new triangle(t->vert[n], t->vert[(n+1)%3], v);
|
|
||||||
new_t[n] = tmp_tri;
|
|
||||||
}
|
|
||||||
|
|
||||||
// sort neighbors for new triangles
|
|
||||||
for (int n = 0; n < 3; ++n)
|
|
||||||
{
|
|
||||||
if (t->neigh[n] == nullptr)
|
|
||||||
{
|
|
||||||
new_t[n]->set_neighbor(nullptr, new_t[(n+1)%3], new_t[(n+2)%3]);
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
new_t[n]->set_neighbor(t->neigh[n], new_t[(n+1)%3], new_t[(n+2)%3]);
|
|
||||||
for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle
|
|
||||||
{
|
|
||||||
if (t->neigh[n]->neigh[k] == t)
|
|
||||||
{
|
|
||||||
t->neigh[n]->neigh[k] = new_t[n];
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
return nullptr;
|
|
||||||
}
|
|
||||||
// End triangle definition
|
// End triangle definition
|
||||||
|
|
||||||
|
// Start DEM definition
|
||||||
|
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; // 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;
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
bool compare_dem_point(dem_point *a, dem_point *b)
|
||||||
|
{
|
||||||
|
if (a->err > b->err) return true;
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
// End DEM definition
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* @brief Generate the TIN from the DEM grid
|
* @brief Generate the TIN from the DEM grid
|
||||||
*
|
*
|
||||||
@ -462,252 +187,274 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
|
|||||||
if (dem.size() != xnum*ynum) return;
|
if (dem.size() != xnum*ynum) return;
|
||||||
|
|
||||||
// Prepare the DEM points
|
// Prepare the DEM points
|
||||||
dem_point *tmp_dem = nullptr;
|
dem_point *tmp_dem = nullptr;;
|
||||||
std::vector<dem_point*> dem_tri;
|
std::vector<dem_point*> dem_grid(xnum*ynum);
|
||||||
std::vector<dem_point*>::iterator d_iter;
|
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] = new dem_point(xmin + dx*j, ymin + dy*i, dem[j + i*xnum]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
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_grid[0]->elev, out_verts.size()); // lower left corner
|
||||||
out_verts.push_back(tmp_vert);
|
out_verts.push_back(tmp_vert);
|
||||||
|
|
||||||
tmp_vert = new vertex2dc(xmax, ymin, dem[xnum-1], out_verts.size()); // lower right corner
|
d_iter = dem_grid.begin();
|
||||||
|
tmp_dem = *d_iter; delete tmp_dem;
|
||||||
|
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);
|
out_verts.push_back(tmp_vert);
|
||||||
|
|
||||||
tmp_vert = new vertex2dc(xmax, ymax, dem[xnum*ynum-1], out_verts.size()); // upper right corner
|
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);
|
out_verts.push_back(tmp_vert);
|
||||||
|
|
||||||
tmp_vert = new vertex2dc(xmin, ymax, dem[xnum*(ynum-1)], out_verts.size()); // upper left corner
|
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);
|
out_verts.push_back(tmp_vert);
|
||||||
|
|
||||||
triangle *old_tri = nullptr, *tmp_tri = nullptr;
|
d_iter = dem_grid.begin() + (xnum*(ynum-1) - 2);
|
||||||
triangle *cnst_tri[4];
|
tmp_dem = *d_iter; delete tmp_dem;
|
||||||
|
dem_grid.erase(d_iter);
|
||||||
|
|
||||||
|
triangle *tmp_tri = nullptr;
|
||||||
|
std::vector<triangle*> cnst_tri, new_tri;
|
||||||
std::vector<triangle*>::iterator t_iter;
|
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
|
||||||
{
|
{
|
||||||
tmp_tri = new triangle(out_verts[0], out_verts[1], out_verts[2]); // order the vertex anti-clock wise
|
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 = nullptr;
|
out_tris.push_back(tmp_tri);
|
||||||
}
|
}
|
||||||
|
|
||||||
if (!is_collinear(out_verts[0], out_verts[2], out_verts[3]))
|
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
|
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); tmp_tri = nullptr;
|
out_tris.push_back(tmp_tri);
|
||||||
}
|
}
|
||||||
|
|
||||||
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);
|
|
||||||
|
|
||||||
// Find host triangle for all DEM locations
|
// Find host triangle for all DEM locations
|
||||||
int tmp_id;
|
for (int i = 0; i < dem_grid.size(); ++i)
|
||||||
for (int i = 0; i < ynum; ++i)
|
|
||||||
{
|
{
|
||||||
for (int j = 0; j < xnum; ++j)
|
for (int t = 0; t < out_tris.size(); ++t)
|
||||||
{
|
{
|
||||||
tmp_id = j + i*xnum;
|
if (out_tris[t]->bound_location(dem_grid[i]->x, dem_grid[i]->y))
|
||||||
if (tmp_id != 0 && tmp_id != (xnum-1) && tmp_id != (xnum*ynum-1) && tmp_id != (xnum*(ynum-1)))
|
|
||||||
{
|
{
|
||||||
tmp_dem = new dem_point(xmin + dx*j, ymin + dy*i, dem[j + i*xnum]);
|
dem_grid[i]->host = out_tris[t];
|
||||||
for (int t = 0; t < out_tris.size(); ++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]);
|
||||||
|
out_tris[t]->circum_dem.push_back(dem_grid[i]);
|
||||||
|
// no beak here. There might be more than one triangle's circumcircle includes the DEM location
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// 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);
|
||||||
|
}
|
||||||
|
|
||||||
|
// Sort dem_grid in the desceding order with respect to the error
|
||||||
|
std::sort(dem_grid.begin(), dem_grid.end(), compare_dem_point);
|
||||||
|
|
||||||
|
bool removed;
|
||||||
|
edge tmp_edge;
|
||||||
|
std::vector<edge> cnst_edge;
|
||||||
|
std::vector<edge>::iterator e_iter;
|
||||||
|
|
||||||
|
while (dem_grid[0]->err >= maxi_err) // quit til the threshold is meet
|
||||||
|
{
|
||||||
|
if (err_records != nullptr)
|
||||||
|
{
|
||||||
|
err_records->push_back(dem_grid[0]->err);
|
||||||
|
}
|
||||||
|
|
||||||
|
// create a new vertex
|
||||||
|
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[0]->circum_host.size(); ++i)
|
||||||
|
{
|
||||||
|
cnst_tri.push_back(dem_grid[0]->circum_host[i]);
|
||||||
|
}
|
||||||
|
|
||||||
|
for (int c = 0; c < cnst_tri.size(); ++c)
|
||||||
|
{
|
||||||
|
for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
|
||||||
|
{
|
||||||
|
tmp_tri = *t_iter;
|
||||||
|
if (cnst_tri[c] == tmp_tri)
|
||||||
{
|
{
|
||||||
if (out_tris[t]->bound_location(tmp_dem->x, tmp_dem->y))
|
t_iter = out_tris.erase(t_iter);
|
||||||
|
break; // no need to search more
|
||||||
|
}
|
||||||
|
else t_iter++;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// remove cnst_tri from its circumed DEM's circum triangle list
|
||||||
|
for (int c = 0; c < cnst_tri.size(); ++c)
|
||||||
|
{
|
||||||
|
for (int i = 0; i < cnst_tri[c]->circum_dem.size(); ++i)
|
||||||
|
{
|
||||||
|
tmp_dem = cnst_tri[c]->circum_dem[i];
|
||||||
|
for (t_iter = tmp_dem->circum_host.begin(); t_iter != tmp_dem->circum_host.end(); )
|
||||||
|
{
|
||||||
|
if (cnst_tri[c] == *t_iter)
|
||||||
{
|
{
|
||||||
tmp_dem->host = out_tris[t];
|
t_iter = tmp_dem->circum_host.erase(t_iter);
|
||||||
tmp_dem->err = fabs(out_tris[t]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
|
break;
|
||||||
out_tris[t]->hosted_dem.push_back(tmp_dem);
|
}
|
||||||
|
else t_iter++;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// remove dem_grid[0] from its circumed triangle's circum DEM list
|
||||||
|
for (int c = 0; c < cnst_tri.size(); ++c)
|
||||||
|
{
|
||||||
|
for (d_iter = cnst_tri[c]->circum_dem.begin(); d_iter != cnst_tri[c]->circum_dem.end(); )
|
||||||
|
{
|
||||||
|
if (dem_grid[0] == *d_iter)
|
||||||
|
{
|
||||||
|
d_iter = cnst_tri[c]->circum_dem.erase(d_iter);
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
else d_iter++;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// clear host and circumcircle triangles for the used DEM location
|
||||||
|
d_iter = dem_grid.begin();
|
||||||
|
tmp_dem = *d_iter; tmp_dem->circum_host.clear(); delete tmp_dem;
|
||||||
|
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);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// loop all DEM data to update host triangles
|
||||||
|
for (int c = 0; c < cnst_tri.size(); ++c)
|
||||||
|
{
|
||||||
|
for (int i = 0; i < cnst_tri[c]->circum_dem.size(); ++i)
|
||||||
|
{
|
||||||
|
tmp_dem = cnst_tri[c]->circum_dem[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(tmp_dem->x, tmp_dem->y))
|
||||||
|
{
|
||||||
|
tmp_dem->host = new_tri[n];
|
||||||
|
tmp_dem->err = fabs(new_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
|
||||||
break; // already found, no need to search more
|
break; // already found, no need to search more
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
|
||||||
|
|
||||||
// Sort hosted_dem in the desceding order with respect to the error. Add maximal zeros to dem_tri
|
// Find circum_host triangles for all DEM locations
|
||||||
for (int t = 0; t < out_tris.size(); ++t)
|
// cnst_tri's circum area doesn't over cover new_tri's circum area
|
||||||
{
|
for (int i = 0; i < dem_grid.size(); ++i)
|
||||||
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]);
|
|
||||||
}
|
|
||||||
|
|
||||||
// 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
|
|
||||||
{
|
|
||||||
if (err_records != nullptr)
|
|
||||||
{
|
{
|
||||||
err_records->push_back(dem_tri[0]->err);
|
for (int n = 0; n < new_tri.size(); ++n) // search in newly created triangles to find new circumcircle triangles
|
||||||
}
|
|
||||||
|
|
||||||
// find the triangle that includes dem_tri[0] 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)
|
|
||||||
{
|
{
|
||||||
t_iter = out_tris.erase(t_iter);
|
dist = (new_tri[n]->cx - dem_grid[i]->x) * (new_tri[n]->cx - dem_grid[i]->x)
|
||||||
break;
|
+ (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
|
||||||
else t_iter++;
|
|
||||||
}
|
|
||||||
|
|
||||||
// remove dem_tri[0] 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)
|
|
||||||
{
|
|
||||||
d_iter = old_tri->hosted_dem.erase(d_iter);
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
else d_iter++;
|
|
||||||
}
|
|
||||||
|
|
||||||
// create a new vertex
|
|
||||||
tmp_vert = new vertex2dc(dem_tri[0]->x, dem_tri[0]->y, dem_tri[0]->elev, out_verts.size());
|
|
||||||
out_verts.push_back(tmp_vert);
|
|
||||||
|
|
||||||
// Delete dem_tri[0]
|
|
||||||
tmp_dem = dem_tri[0]; delete tmp_dem;
|
|
||||||
|
|
||||||
// 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);
|
new_tri[n]->circum_dem.push_back(dem_grid[i]);
|
||||||
break;
|
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
|
||||||
}
|
}
|
||||||
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 and sort dem_tri
|
|
||||||
dem_tri.clear(); dem_tri.reserve(out_tris.size());
|
|
||||||
for (int t = 0; t < out_tris.size(); t++)
|
|
||||||
{
|
|
||||||
if (!out_tris[t]->hosted_dem.empty())
|
|
||||||
{
|
|
||||||
dem_tri.push_back(out_tris[t]->hosted_dem[0]);
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
std::sort(dem_tri.begin(), dem_tri.end(), compare_dem_point);
|
|
||||||
|
// destroy memories used by cnst_edge
|
||||||
|
for (int c = 0; c < cnst_tri.size(); ++c)
|
||||||
|
{
|
||||||
|
tmp_tri = cnst_tri[c];
|
||||||
|
tmp_tri->circum_dem.clear();
|
||||||
|
delete tmp_tri; tmp_tri = nullptr;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Sort dem_grid in the desceding order with respect to the error
|
||||||
|
std::sort(dem_grid.begin(), dem_grid.end(), compare_dem_point);
|
||||||
}
|
}
|
||||||
|
|
||||||
if (err_records != nullptr)
|
if (err_records != nullptr)
|
||||||
{
|
{
|
||||||
err_records->push_back(dem_tri[0]->err);
|
err_records->push_back(dem_grid[0]->err);
|
||||||
}
|
}
|
||||||
|
|
||||||
// assign triangles index
|
// destroy remaining DEM data
|
||||||
for (int i = 0; i < out_tris.size(); i++)
|
for (int i = 0; i < dem_grid.size(); ++i)
|
||||||
{
|
{
|
||||||
out_tris[i]->id = i;
|
tmp_dem = dem_grid[i];
|
||||||
// destroy remaining DEM data
|
delete tmp_dem; tmp_dem = nullptr;
|
||||||
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;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
return;
|
return;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user