update readme

README.md

@@ -0,0 +1,30 @@
+### 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%;" />
+
+#### Compile
+
+```shell
+g++ demo.cpp
+```
+
+#### Run demo
+
+```shell
+./a.out
+```
+
+#### Visualization
+
+##### macOS
+
+```shell
+open -a Gmsh topo_TIN.msh
+```
+
+##### Linux/Unix
+
+```shell
+gmsh topo_TIN.msh
+```
+

demo.cpp

@@ -17,7 +17,7 @@ int main(int argc, char const *argv[])
 
 	std::vector<vertex2dc*> tin_vert;
 	std::vector<triangle*> tin_ele;
-	dem2tin(topo, 0, 1000, 0, 1000, 10, 10, tin_vert, tin_ele, 3.0);
+	dem2tin(topo, 0, 1000, 0, 1000, 10, 10, tin_vert, tin_ele, 1.0);
 
 	// Write a Gmsh's .msh file
 	std::ofstream outfile("topo_TIN.msh");

tin.h

@@ -145,8 +145,46 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
 
 	if (dem.size() != xnum*ynum) return;
 
-	vertex2dc *tmp_vert = nullptr;
+	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);
 
@@ -159,25 +197,65 @@ 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);
 
-	triangle *tmp_tri = nullptr;
-	std::vector<triangle*> cnst_tri;
-	std::vector<triangle*>::iterator t_iter;
+	for (int i = 0; i < 4; ++i)
+	{
+		// 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++;
+		}
 
-	tmp_tri = new triangle(out_verts[0], out_verts[1], out_verts[2]); // order the vertex anti-clock wise
+		// 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);
+		}
 
-	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;
+		// 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;
+		}
+	}
 
 	do // quit til the threshold is meet
 	{
@@ -266,6 +344,26 @@ 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;
 }