tmp update

demo.cpp

@@ -1,4 +1,4 @@
-#include "tin.h"
+#include "tin_backup.h"
 #include "iostream"
 #include "fstream"
 #include "iomanip"

tin.h

@@ -289,7 +289,8 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
 	std::vector<edge> cnst_edge;
 	std::vector<edge>::iterator e_iter;
 
-	while (dem_grid[0]->err >= maxi_err) // quit til the threshold is meet
+	//while (dem_grid[0]->err >= maxi_err) // quit til the threshold is meet
+	for (int q = 0; q < 28; q++)
 	{
 		if (err_records != nullptr)
 		{

tin_backup.h

@@ -0,0 +1,722 @@
+/**
+ * @defgroup   TIN
+ *
+ * @brief      Generation of a Triangular Irregular Network (TIN) from a dense DEM grid.
+ *
+ * @author     Yi Zhang
+ * @date       2021-09-16
+ */
+
+#ifndef _TIN_DELAUNAY_H
+#define _TIN_DELAUNAY_H
+#include "cmath"
+#include "vector"
+#include "algorithm"
+
+#include "iostream"
+
+#define ZERO 1e-5
+
+// Start vertex definition
+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), elev(NAN), id(0) {}
+	vertex2dc(double inx, double iny, double inelev, unsigned int inid) {set(inx, iny, inelev, inid);}
+	void set(double inx, double iny, double inelev, unsigned int inid)
+	{
+		x = inx; y = iny; elev = inelev; id = inid;
+		return;
+	}
+};
+
+bool operator ==(const vertex2dc &a, const vertex2dc &b) // overload the == operator for vertex2dc type
+{
+	if(fabs(a.x - b.x) <= ZERO && fabs(a.y - b.y) <= ZERO)
+	{
+		return true;
+	}
+	return false;
+}
+
+bool is_collinear(vertex2dc *a_ptr, vertex2dc *b_ptr, vertex2dc *c_ptr) // Test if the three points are on the same line
+{
+	// ｜(y3−y1)(x2−x1)−(y2−y1)(x3−x1)｜
+	if (fabs((c_ptr->y - a_ptr->y)*(b_ptr->x - a_ptr->x) - (b_ptr->y - a_ptr->y)*(c_ptr->x - a_ptr->x)) <= ZERO)
+	{
+		return true;
+	}
+	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
+
+// Start DEM definition
+struct triangle;
+
+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;
+
+	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)
+	{
+		x = inx; y = iny; elev = inelev; err = 0.0; 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
+
+// Start triangle definition
+struct triangle
+{
+	int id;
+	vertex2dc *vert[3]; // vertex of the triangle
+	triangle *neigh[3]; // neighbors of the triangle
+	double cx, cy; // center of the triangle's circumcircle
+	double cr; // radius of the circumcircle
+	std::vector<dem_point*> hosted_dem;
+
+	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);}
+	void set(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *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)
+	{
+		neigh[0] = n0ptr; neigh[1] = n1ptr; neigh[2] = n2ptr;
+		return;
+	}
+
+	bool bound_location(double inx, double iny) // Test if the location is inside the triangle
+	{
+		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) // This condition includes points on the triangle's edge
+			{
+				return false;
+			}
+		}
+		return true;
+	}
+
+	double interpolate(double inx, double iny) // Interpolate the elevation of the given location inside the triangle
+	{
+		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);
+	}
+};
+
+/**
+ * @brief    Flip neighboring triangles and their neighbors
+ * 
+ * original
+ * 
+ *        /\
+ *       /  \
+ *      /    \
+ *     /   t  \
+ *  t_id-------\  t_id (0, 1 or 2)
+ *    \--------/
+ *     \      /
+ *      \  n /
+ *       \  /
+ *        \/
+ *        n_id (0, 1 or 2)
+ * 
+ * fliped
+ * 
+ *        /|\
+ *       / | \
+ *      /  |  \
+ *     /   |   \
+ *  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);
+	}	
+	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) < 0.0) // need to be flipped
+					{
+						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]);
+
+				if (new_t[0]->neigh[2] != nullptr)
+				{
+					for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle 
+					{
+						if (new_t[0]->neigh[2]->neigh[k] == t)
+						{
+							new_t[0]->neigh[2]->neigh[k] = new_t[0];
+							break;
+						}
+					}
+				}
+				
+				if (new_t[1]->neigh[2] != nullptr)
+				{
+					for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle 
+					{
+						if (new_t[1]->neigh[2]->neigh[k] == t)
+						{
+							new_t[1]->neigh[2]->neigh[k] = new_t[1];
+							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]);
+
+					if (new_t[0]->neigh[2] != nullptr)
+					{
+						for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle 
+						{
+							if (new_t[0]->neigh[2]->neigh[k] == t)
+							{
+								new_t[0]->neigh[2]->neigh[k] = new_t[0];
+								break;
+							}
+						}
+					}
+					
+					if (new_t[1]->neigh[2] != nullptr)
+					{
+						for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle 
+						{
+							if (new_t[1]->neigh[2]->neigh[k] == t)
+							{
+								new_t[1]->neigh[2]->neigh[k] = new_t[1];
+								break;
+							}
+						}
+					}
+					
+					if (new_t[2]->neigh[2] != nullptr)
+					{
+						for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle 
+						{
+							if (new_t[2]->neigh[2]->neigh[k] == t->neigh[i])
+							{
+								new_t[2]->neigh[2]->neigh[k] = new_t[2];
+								break;
+							}
+						}
+					}
+					
+					if (new_t[3]->neigh[2] != nullptr)
+					{
+						for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle 
+						{
+							if (new_t[3]->neigh[2]->neigh[k] == t->neigh[i])
+							{
+								new_t[3]->neigh[2]->neigh[k] = new_t[3];
+								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
+
+/**
+ * @brief      Generate the TIN from the DEM grid
+ *
+ * @param[in]  dem          Input DEM grid (Ordered from lower left corner to the upper right corner)
+ * @param[in]  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-0
+ * @param[in]  err_records  If this pointer is not NULL, record maximal error values after each insertion of vertex.
+ */
+void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ymin, double ymax, 
+	double dx, double dy, std::vector<vertex2dc*> &out_verts, std::vector<triangle*> &out_tris, 
+	double maxi_err = 1e-0, std::vector<double> *err_records = nullptr)
+{
+	if (!out_verts.empty()) out_verts.clear();
+	if (!out_tris.empty())  out_tris.clear();
+	if (err_records != nullptr && !err_records->empty()) err_records->clear();
+
+	if (dx <= 0.0 || dy <= 0.0 || maxi_err <= 0.0) return;
+	if (xmin >= xmax || ymin >= ymax || (xmin + dx) > xmax || (ymin + dy) > ymax) return;
+
+	int xnum = round((xmax - xmin)/dx) + 1;
+	int ynum = round((ymax - ymin)/dy) + 1;
+
+	if (dem.size() != xnum*ynum) return;
+
+	// Prepare the DEM points
+	dem_point *tmp_dem = nullptr;
+	std::vector<dem_point*> cnst_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] = 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
+	out_verts.push_back(tmp_vert);
+
+	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);
+
+	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 - 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 *old_tri = nullptr, *tmp_tri = nullptr;
+	triangle *cnst_tri[4];
+	std::vector<triangle*>::iterator t_iter;
+
+	if (!is_collinear(out_verts[0], out_verts[1], out_verts[2])) // Do not create triangle if the vertexes are collinear
+	{
+		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;
+	}
+
+	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
+		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
+	for (int i = 0; i < dem_grid.size(); ++i)
+	{
+		for (int t = 0; t < out_tris.size(); ++t)
+		{
+			if (out_tris[t]->bound_location(dem_grid[i]->x, dem_grid[i]->y))
+			{
+				dem_grid[i]->host = out_tris[t];
+				out_tris[t]->hosted_dem.push_back(dem_grid[i]);
+				break; // already found, no need to search more
+			}
+		}
+	}
+
+	// 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);
+
+	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);
+		}
+
+		// find the triangle that includes dem_grid[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_grid[0]->host)
+			{
+				t_iter = out_tris.erase(t_iter);
+				break;
+			}
+			else t_iter++;
+		}
+
+		// remove dem_grid[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_grid[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_grid[0]->x, dem_grid[0]->y, dem_grid[0]->elev, out_verts.size());
+		out_verts.push_back(tmp_vert);
+
+		// Remove dem_grid[0] from the list and delete it
+		d_iter = dem_grid.begin();
+		tmp_dem = *d_iter; delete tmp_dem;
+		dem_grid.erase(d_iter);
+
+		// build new triangles
+		tmp_tri = split_triangle(tmp_vert, old_tri, cnst_tri);
+		for (int n = 0; n < 4; ++n)
+		{
+			if (cnst_tri[n] != nullptr)
+			{
+				out_tris.push_back(cnst_tri[n]);
+			}
+		}
+
+		if (tmp_tri != nullptr)
+		{
+			for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
+			{
+				if (tmp_tri == *t_iter)
+				{
+					t_iter = out_tris.erase(t_iter);
+					break;
+				}
+				else t_iter++;
+			}
+
+			// build hosted dem for the new triangles
+			for (int d = 0; d < old_tri->hosted_dem.size(); d++)
+			{
+				tmp_dem = old_tri->hosted_dem[d];
+				for (int n = 0; n < 4; n++)
+				{
+					if (cnst_tri[n] != nullptr && cnst_tri[n]->bound_location(tmp_dem->x, tmp_dem->y))
+					{
+						tmp_dem->host = cnst_tri[n];
+						tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
+						cnst_tri[n]->hosted_dem.push_back(tmp_dem);
+						break;
+					}
+				}
+			}
+
+			for (int d = 0; d < tmp_tri->hosted_dem.size(); d++)
+			{
+				tmp_dem = tmp_tri->hosted_dem[d];
+				for (int n = 0; n < 4; n++)
+				{
+					if (cnst_tri[n] != nullptr && cnst_tri[n]->bound_location(tmp_dem->x, tmp_dem->y))
+					{
+						tmp_dem->host = cnst_tri[n];
+						tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
+						cnst_tri[n]->hosted_dem.push_back(tmp_dem);
+						break;
+					}
+				}
+			}
+
+			// 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;
+					}
+				}
+			}
+
+			// 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]);
+			}
+		}
+
+		// 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)
+	{
+		err_records->push_back(dem_grid[0]->err);
+	}
+
+	// assign triangles index
+	for (int i = 0; i < out_tris.size(); i++)
+	{
+		out_tris[i]->id = i;
+	}
+
+	// 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;
+}
+
+#endif // _TIN_DELAUNAY_H