diff --git a/demo.cpp b/demo.cpp
index 8738bfb..e85f11c 100644
--- a/demo.cpp
+++ b/demo.cpp
@@ -1,4 +1,4 @@
-#include "tin_backup.h"
+#include "tin.h"
 #include "iostream"
 #include "fstream"
 #include "iomanip"
diff --git a/tin.h b/tin.h
index e563bd3..3b207dc 100644
--- a/tin.h
+++ b/tin.h
@@ -13,8 +13,6 @@
 #include "vector"
 #include "algorithm"
 
-#include "iostream"
-
 #define ZERO 1e-5
 
 // Start vertex definition
@@ -51,56 +49,77 @@ bool is_collinear(vertex2dc *a_ptr, vertex2dc *b_ptr, vertex2dc *c_ptr) // Test
 	}
 	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 edge definition
-struct edge
-{
-	vertex2dc *vert[2]; // vertex of the edge
+// Start DEM definition
+struct triangle;
 
-	edge() {vert[0] = vert[1] = nullptr;}
-	edge(vertex2dc *v0ptr, vertex2dc *v1ptr) {set(v0ptr, v1ptr);}
-	void set(vertex2dc *v0ptr, vertex2dc *v1ptr)
+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)
 	{
-		vert[0] = v0ptr; vert[1] = v1ptr;
+		x = inx; y = iny; elev = inelev; err = 0.0; host = nullptr;
 		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]) || 
-		(a.vert[0] == b.vert[1] && a.vert[1] == b.vert[0]))
-	{
-		return true;
-	}
+	if (a->err > b->err) return true;
 	return false;
 }
-// End edge definition
-
-// Start triangle definition
-struct dem_point;
+// End DEM definition
 
+/* Start triangle definition
+ *            v2
+ *            /\
+ *          /   \
+ *     n2 /      \ n1
+ *      /         \
+ *    /------------\
+ *   v0     n0      v1
+ */
 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*> 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);}
 	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;
+	}
 
-		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;
-		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
+	void set_neighbor(triangle *n0ptr, triangle *n1ptr, triangle *n2ptr)
+	{
+		neigh[0] = n0ptr; neigh[1] = n1ptr; neigh[2] = n2ptr;
 		return;
 	}
 
@@ -130,32 +149,288 @@ struct triangle
 		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
-	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
+	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
 
-	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)
+	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)
 	{
-		x = inx; y = iny; elev = inelev; host = nullptr;
-		return;
+		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;
+			}
+		}
 	}
-};
 
-bool compare_dem_point(dem_point *a, dem_point *b)
-{
-	if (a->err > b->err) return true;
-	return false;
+	// 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++;
+	}
+
+	// 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
@@ -189,274 +464,252 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
 	if (dem.size() != xnum*ynum) return;
 
 	// Prepare the DEM points
-	dem_point *tmp_dem = nullptr;;
-	std::vector<dem_point*> dem_grid(xnum*ynum);
+	dem_point *tmp_dem = nullptr;
+	std::vector<dem_point*> dem_tri;
 	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
+	tmp_vert = new vertex2dc(xmin, ymin, dem[0], 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
+	tmp_vert = new vertex2dc(xmax, ymin, dem[xnum-1], out_verts.size()); // lower right corner
 	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
+	tmp_vert = new vertex2dc(xmax, ymax, dem[xnum*ynum-1], out_verts.size()); // upper right corner
 	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
+	tmp_vert = new vertex2dc(xmin, ymax, dem[xnum*(ynum-1)], out_verts.size()); // upper left corner
 	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;
+	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);
+		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);
+		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)
+	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];
-				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)
+				tmp_dem = new dem_point(xmin + dx*j, ymin + dy*i, dem[j + i*xnum]);
+				for (int t = 0; t < out_tris.size(); ++t)
 				{
-					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)
+					if (out_tris[t]->bound_location(tmp_dem->x, tmp_dem->y))
 					{
-						t_iter = tmp_dem->circum_host.erase(t_iter);
-						break;
-					}
-					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);
+						tmp_dem->host = out_tris[t];
+						tmp_dem->err = fabs(out_tris[t]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
+						out_tris[t]->hosted_dem.push_back(tmp_dem);
 						break; // already found, no need to search more
 					}
 				}
 			}
 		}
+	}
 
-		// Find circum_host triangles for all DEM locations
-		// cnst_tri's circum area doesn't over cover new_tri's circum area
-		for (int i = 0; i < dem_grid.size(); ++i)
+	// Sort hosted_dem in the desceding order with respect to the error. Add maximal zeros to dem_tri
+	for (int t = 0; t < out_tris.size(); ++t)
+	{
+		std::sort(out_tris[t]->hosted_dem.begin(), out_tris[t]->hosted_dem.end(), compare_dem_point);
+		dem_tri.push_back(out_tris[t]->hosted_dem[0]);
+	}
+
+	// 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) 
-					+ (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
+				t_iter = out_tris.erase(t_iter);
+				break;
+			}
+			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]);
-					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
+					t_iter = out_tris.erase(t_iter);
+					break;
 				}
+				else t_iter++;
+			}
+
+			// build hosted dem for the new triangles
+			for (int d = 0; d < old_tri->hosted_dem.size(); d++)
+			{
+				tmp_dem = old_tri->hosted_dem[d];
+				for (int n = 0; n < 4; n++)
+				{
+					if (cnst_tri[n] != nullptr && cnst_tri[n]->bound_location(tmp_dem->x, tmp_dem->y))
+					{
+						tmp_dem->host = cnst_tri[n];
+						tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
+						cnst_tri[n]->hosted_dem.push_back(tmp_dem);
+						break;
+					}
+				}
+			}
+
+			for (int d = 0; d < tmp_tri->hosted_dem.size(); d++)
+			{
+				tmp_dem = tmp_tri->hosted_dem[d];
+				for (int n = 0; n < 4; n++)
+				{
+					if (cnst_tri[n] != nullptr && cnst_tri[n]->bound_location(tmp_dem->x, tmp_dem->y))
+					{
+						tmp_dem->host = cnst_tri[n];
+						tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
+						cnst_tri[n]->hosted_dem.push_back(tmp_dem);
+						break;
+					}
+				}
+			}
+
+			for (int n = 0; n < 4; n++)
+			{
+				if (cnst_tri[n] != nullptr)
+				{
+					std::sort(cnst_tri[n]->hosted_dem.begin(), cnst_tri[n]->hosted_dem.end(), compare_dem_point);
+				}
+			}
+
+			// delete the old triangle
+			old_tri->hosted_dem.clear();
+			delete old_tri; old_tri = nullptr;
+
+			tmp_tri->hosted_dem.clear();
+			delete tmp_tri; tmp_tri = nullptr;
+		}
+		else
+		{
+			// build hosted dem for the new triangles
+			for (int d = 0; d < old_tri->hosted_dem.size(); d++)
+			{
+				tmp_dem = old_tri->hosted_dem[d];
+				for (int n = 0; n < 4; n++)
+				{
+					if (cnst_tri[n] != nullptr && cnst_tri[n]->bound_location(tmp_dem->x, tmp_dem->y))
+					{
+						tmp_dem->host = cnst_tri[n];
+						tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
+						cnst_tri[n]->hosted_dem.push_back(tmp_dem);
+						break;
+					}
+				}
+			}
+
+			for (int n = 0; n < 4; n++)
+			{
+				if (cnst_tri[n] != nullptr)
+				{
+					std::sort(cnst_tri[n]->hosted_dem.begin(), cnst_tri[n]->hosted_dem.end(), compare_dem_point);
+				}
+			}
+
+			// delete the old triangle
+			old_tri->hosted_dem.clear();
+			delete old_tri; old_tri = nullptr;
+		}
+
+		// Make sure cnst_tri meet the empty circumcircle condition
+		for (int n = 0; n < 4; ++n)
+		{
+			if (cnst_tri[n] != nullptr)
+			{
+				make_delaunay(cnst_tri[n]);
+			}
+		}
+
+		// get maximal errors from out_tris 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]);
 			}
 		}
 		
-
-		// 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);
+		std::sort(dem_tri.begin(), dem_tri.end(), compare_dem_point);
 	}
 
 	if (err_records != nullptr)
 	{
-		err_records->push_back(dem_grid[0]->err);
+		err_records->push_back(dem_tri[0]->err);
 	}
 
-	// destroy remaining DEM data
-	for (int i = 0; i < dem_grid.size(); ++i)
+	// assign triangles index
+	for (int i = 0; i < out_tris.size(); i++)
 	{
-		tmp_dem = dem_grid[i];
-		delete tmp_dem; tmp_dem = nullptr;
+		out_tris[i]->id = i;
+		// destroy remaining DEM data
+		for (int d = 0; d < out_tris[i]->hosted_dem.size(); d++)
+		{
+			tmp_dem = out_tris[i]->hosted_dem[d];
+			delete tmp_dem; tmp_dem = nullptr;
+		}
 	}
-
 	return;
 }
 
diff --git a/tin_backup.h b/tin_backup.h
index 263ff39..d13f24c 100644
--- a/tin_backup.h
+++ b/tin_backup.h
@@ -13,8 +13,6 @@
 #include "vector"
 #include "algorithm"
 
-#include "iostream"
-
 #define ZERO 1e-5
 
 // Start vertex definition
@@ -51,77 +49,56 @@ bool is_collinear(vertex2dc *a_ptr, vertex2dc *b_ptr, vertex2dc *c_ptr) // Test
 	}
 	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
+// Start edge definition
+struct edge
 {
-	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;
+	vertex2dc *vert[2]; // vertex of the edge
 
-	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)
+	edge() {vert[0] = vert[1] = nullptr;}
+	edge(vertex2dc *v0ptr, vertex2dc *v1ptr) {set(v0ptr, v1ptr);}
+	void set(vertex2dc *v0ptr, vertex2dc *v1ptr)
 	{
-		x = inx; y = iny; elev = inelev; err = 0.0; host = nullptr;
+		vert[0] = v0ptr; vert[1] = v1ptr;
 		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;
 }
-// End DEM definition
+// End edge definition
+
+// Start triangle definition
+struct dem_point;
 
-/* Start triangle definition
- *            v2
- *            /\
- *          /   \
- *     n2 /      \ n1
- *      /         \
- *    /------------\
- *   v0     n0      v1
- */
 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;
+	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);}
 	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;
+		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;
+		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;
 	}
 
@@ -151,285 +128,33 @@ struct triangle
 		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
 
+// 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
  *
@@ -462,252 +187,274 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
 	if (dem.size() != xnum*ynum) return;
 
 	// Prepare the DEM points
-	dem_point *tmp_dem = nullptr;
-	std::vector<dem_point*> dem_tri;
+	dem_point *tmp_dem = nullptr;;
+	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[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);
 
-	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);
 
-	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);
 
-	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);
 
-	triangle *old_tri = nullptr, *tmp_tri = nullptr;
-	triangle *cnst_tri[4];
+	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;
 	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;
+		out_tris.push_back(tmp_tri);
 	}
 
 	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;
+		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
-	int tmp_id;
-	for (int i = 0; i < ynum; ++i)
+	for (int i = 0; i < dem_grid.size(); ++i)
 	{
-		for (int j = 0; j < xnum; ++j)
+		for (int t = 0; t < out_tris.size(); ++t)
 		{
-			tmp_id = j + i*xnum;
-			if (tmp_id != 0 && tmp_id != (xnum-1) && tmp_id != (xnum*ynum-1) && tmp_id != (xnum*(ynum-1)))
+			if (out_tris[t]->bound_location(dem_grid[i]->x, dem_grid[i]->y))
 			{
-				tmp_dem = new dem_point(xmin + dx*j, ymin + dy*i, dem[j + i*xnum]);
-				for (int t = 0; t < out_tris.size(); ++t)
+				dem_grid[i]->host = out_tris[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];
-						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);
+						t_iter = tmp_dem->circum_host.erase(t_iter);
+						break;
+					}
+					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
 					}
 				}
 			}
 		}
-	}
 
-	// Sort hosted_dem in the desceding order with respect to the error. Add maximal zeros to dem_tri
-	for (int t = 0; t < out_tris.size(); ++t)
-	{
-		std::sort(out_tris[t]->hosted_dem.begin(), out_tris[t]->hosted_dem.end(), compare_dem_point);
-		dem_tri.push_back(out_tris[t]->hosted_dem[0]);
-	}
-
-	// 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)
+		// Find circum_host triangles for all DEM locations
+		// cnst_tri's circum area doesn't over cover new_tri's circum area
+		for (int i = 0; i < dem_grid.size(); ++i)
 		{
-			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)
+			for (int n = 0; n < new_tri.size(); ++n) // search in newly created triangles to find new circumcircle triangles
 			{
-				t_iter = out_tris.erase(t_iter);
-				break;
-			}
-			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)
+				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
 				{
-					t_iter = out_tris.erase(t_iter);
-					break;
+					new_tri[n]->circum_dem.push_back(dem_grid[i]);
+					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)
 	{
-		err_records->push_back(dem_tri[0]->err);
+		err_records->push_back(dem_grid[0]->err);
 	}
 
-	// assign triangles index
-	for (int i = 0; i < out_tris.size(); i++)
+	// destroy remaining DEM data
+	for (int i = 0; i < dem_grid.size(); ++i)
 	{
-		out_tris[i]->id = i;
-		// destroy remaining DEM data
-		for (int d = 0; d < out_tris[i]->hosted_dem.size(); d++)
-		{
-			tmp_dem = out_tris[i]->hosted_dem[d];
-			delete tmp_dem; tmp_dem = nullptr;
-		}
+		tmp_dem = dem_grid[i];
+		delete tmp_dem; tmp_dem = nullptr;
 	}
+
 	return;
 }