update src

README.md

@@ -1,6 +1,6 @@
 ### Generation of a Triangular Irregular Network (TIN) from a dense DEM grid
 
-<img src="topo_TIN_err1e0.png" alt="topo_TIN_err1e0" style="zoom:24%;" />
+<img src="topo_TIN.png" alt="topo_TIN" style="zoom:24%;" />
 
 #### Compile
 

tin.h

@@ -12,8 +12,6 @@
 #include "cmath"
 #include "vector"
 
-#include "iostream"
-
 #define ZERO 1e-5
 
 // Start vertex definition
@@ -40,6 +38,16 @@ bool operator ==(const vertex2dc &a, const vertex2dc &b) // overload the == oper
 	}
 	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;
+}
 // End vertex definition
 
 // Start edge definition
@@ -90,7 +98,7 @@ struct triangle
 		return;
 	}
 
-	bool bound_location(double inx, double iny)
+	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++)
@@ -108,7 +116,7 @@ struct triangle
 		return true;
 	}
 
-	double interpolate(double inx, double iny)
+	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));
@@ -138,6 +146,7 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
 	if (!out_verts.empty()) out_verts.clear();
 	if (!out_tris.empty())  out_tris.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;
@@ -145,46 +154,8 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
 
 	if (dem.size() != xnum*ynum) return;
 
-	double midx = 0.5*(xmin + xmax);
-	double midy = 0.5*(ymin + ymax);
-	double maxi_s = std::max(xmax - xmin, ymax - ymin); // use an four times bigger rectangle to include all points
-
 	vertex2dc *tmp_vert = nullptr;
-	std::vector<vertex2dc*> assit_vert, corner_vert;
 
-	tmp_vert = new vertex2dc(midx - maxi_s, midy - maxi_s, 0.0); // lower left corner
-	assit_vert.push_back(tmp_vert);
-
-	tmp_vert = new vertex2dc(midx + maxi_s, midy - maxi_s, 0.0); // lower right corner
-	assit_vert.push_back(tmp_vert);
-
-	tmp_vert = new vertex2dc(midx + maxi_s, midy + maxi_s, 0.0); // upper right corner
-	assit_vert.push_back(tmp_vert);
-
-	tmp_vert = new vertex2dc(midx - maxi_s, midy + maxi_s, 0.0); // upper left corner
-	assit_vert.push_back(tmp_vert);
-
-	triangle *tmp_tri = nullptr;
-	std::vector<triangle*> cnst_tri;
-	std::vector<triangle*>::iterator t_iter;
-
-	tmp_tri = new triangle(assit_vert[0], assit_vert[1], assit_vert[2]); // order the vertex anti-clock wise
-	out_tris.push_back(tmp_tri);
-
-	tmp_tri = new triangle(assit_vert[0], assit_vert[2], assit_vert[3]); // order the vertex anti-clock wise
-	out_tris.push_back(tmp_tri);
-
-	int now_maxi_id;
-	double now_x, now_y, now_err;
-	double now_maxi_err;
-
-	bool removed;
-	double dist;
-	edge tmp_edge;
-	std::vector<edge> cnst_edge;
-	std::vector<edge>::iterator e_iter;
-
-	// We first add the four corner points of the DEM to the TIN
 	tmp_vert = new vertex2dc(xmin, ymin, dem[0], out_verts.size()); // lower left corner
 	out_verts.push_back(tmp_vert);
 
@@ -197,66 +168,32 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
 	tmp_vert = new vertex2dc(xmin, ymax, dem[xnum*(ynum-1)], out_verts.size()); // upper left corner
 	out_verts.push_back(tmp_vert);
 
-	for (int i = 0; i < 4; ++i)
+	triangle *tmp_tri = nullptr;
+	std::vector<triangle*> cnst_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
 	{
-		// determine triangles that include the point and add the triangle to the cnst_tri and remove the triangle from out_tris
-		// this is also a part that could take a lot of time if we are working with a large amount of points. We will fix this later
-		cnst_tri.clear();
-		for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
-		{
-			tmp_tri = *t_iter;
-			dist = (tmp_tri->cx - out_verts[i]->x) * (tmp_tri->cx - out_verts[i]->x) 
-				+ (tmp_tri->cy - out_verts[i]->y) * (tmp_tri->cy - out_verts[i]->y);
-			if ((dist - tmp_tri->cr) <= ZERO) // Points on the circumcircle are also included
-			{
-				t_iter = out_tris.erase(t_iter);
-				cnst_tri.push_back(tmp_tri);
-			}
-			else t_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
-		for (int c = 0; c < cnst_edge.size(); ++c)
-		{
-			tmp_tri = new triangle(cnst_edge[c].vert[0], cnst_edge[c].vert[1], out_verts[i]); // order the vertex anti-clock wise
-			out_tris.push_back(tmp_tri);
-		}
-
-		// destroy memories used by cnst_edge
-		for (int c = 0; c < cnst_tri.size(); ++c)
-		{
-			tmp_tri = cnst_tri[c];
-			delete tmp_tri; tmp_tri = nullptr;
-		}
+		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);
 	}
 
+	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);
+	}
+
+	int now_maxi_id;
+	double now_x, now_y, now_err;
+	double now_maxi_err;
+
+	bool removed;
+	double dist;
+	edge tmp_edge;
+	std::vector<edge> cnst_edge;
+	std::vector<edge>::iterator e_iter;
+
 	do // quit til the threshold is meet
 	{
 		// loop all DEM data to find the location with maximal error
@@ -287,8 +224,8 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
 		tmp_vert = new vertex2dc(now_x, now_y, dem[now_maxi_id], out_verts.size());
 		out_verts.push_back(tmp_vert);
 
-		// determine triangles that include the point and add the triangle to the cnst_tri and remove the triangle from out_tris
-		// this is also a part that could take a lot of time if we are working with a large amount of points. We will fix this later
+		// determine triangles that include the point and add the triangle to the cnst_tri and remove it from out_tris
+		// this is also a part that could take a lot of time if we are working with a large amount of points. We will fix it later
 		cnst_tri.clear();
 		for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
 		{
@@ -332,8 +269,11 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
 		// construct new triangles and add to out_tris
 		for (int c = 0; c < cnst_edge.size(); ++c)
 		{
-			tmp_tri = new triangle(cnst_edge[c].vert[0], cnst_edge[c].vert[1], out_verts.back()); // order the vertex anti-clock wise
-			out_tris.push_back(tmp_tri);
+			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);
+			}
 		}
 
 		// destroy memories used by cnst_edge
@@ -344,26 +284,6 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
 		}
 	} while (now_maxi_err >= maxi_err);
 
-	// remove any triangles has an assistant vertex from out_tris
-	for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
-	{
-		tmp_tri = *t_iter;
-		if (tmp_tri->vert[0] == assit_vert[0] || tmp_tri->vert[0] == assit_vert[1] || tmp_tri->vert[0] == assit_vert[2] || tmp_tri->vert[0] == assit_vert[3] || 
-			tmp_tri->vert[1] == assit_vert[0] || tmp_tri->vert[1] == assit_vert[1] || tmp_tri->vert[1] == assit_vert[2] || tmp_tri->vert[1] == assit_vert[3] || 
-			tmp_tri->vert[2] == assit_vert[0] || tmp_tri->vert[2] == assit_vert[1] || tmp_tri->vert[2] == assit_vert[2] || tmp_tri->vert[2] == assit_vert[3])
-		{
-			// destroy the memories located and remove from the vector
-			t_iter = out_tris.erase(t_iter);
-			delete tmp_tri; tmp_tri = nullptr;
-		}
-		else t_iter++;
-	}
-
-	// destroy memories located for assit_vert
-	for (int i = 0; i < 4; ++i)
-	{
-		delete assit_vert[i]; assit_vert[i] = nullptr;
-	}
 	return;
 }