diff --git a/.gitignore b/.gitignore
index 259148f..6c3feba 100644
--- a/.gitignore
+++ b/.gitignore
@@ -30,3 +30,5 @@
 *.exe
 *.out
 *.app
+
+.DS_Store
\ No newline at end of file
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..e69de29
diff --git a/demo.cpp b/demo.cpp
new file mode 100644
index 0000000..9f3688f
--- /dev/null
+++ b/demo.cpp
@@ -0,0 +1,8 @@
+#include "tin.h"
+#include "iostream"
+
+int main(int argc, char const *argv[])
+{
+	
+	return 0;
+}
\ No newline at end of file
diff --git a/tin.h b/tin.h
new file mode 100644
index 0000000..59bf523
--- /dev/null
+++ b/tin.h
@@ -0,0 +1,242 @@
+/**
+ * @defgroup   TIN
+ *
+ * @brief      Generation of a Triangular Irregular Network (TIN) from a dense DEM grid
+ *
+ * @author     Yi Zhang
+ * @date       2021-09-15
+ */
+
+#ifndef _TIN_DELAUNAY_H
+#define _TIN_DELAUNAY_H
+#include "cmath"
+#include "vector"
+
+#define ZERO 1e-5
+
+// Start vertex definition
+struct vertex2dc
+{
+	unsigned int id; // index of the vertex
+	double x, y; // position of the vertex
+
+	vertex2dc() : x(NAN), y(NAN), id(0) {}
+	vertex2dc(double inx, double iny, unsigned int inid = 0) {set(inx, iny, inid);}
+	void set(double inx, double iny, unsigned int inid = 0)
+	{
+		x = inx; y = iny; 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;
+}
+// End vertex definition
+
+// Start edge definition
+struct edge
+{
+	vertex2dc *vert[2]; // vertex of the edge
+
+	edge() {vert[0] = vert[1] = nullptr;}
+	edge(vertex2dc *v0ptr, vertex2dc *v1ptr) {set(v0ptr, v1ptr);}
+	void set(vertex2dc *v0ptr, vertex2dc *v1ptr)
+	{
+		vert[0] = v0ptr; vert[1] = v1ptr;
+		return;
+	}
+};
+
+bool operator ==(const edge &a, const edge &b) // overload the == operator for edge type
+{
+	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 edge definition
+
+// Start triangle definition
+struct triangle
+{
+	vertex2dc *vert[3]; // vertex of the triangle
+	double cx, cy; // center of the triangle's circumcircle
+	double cr; // radius of the circumcircle
+
+	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;
+
+		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;
+	}
+};
+// End triangle definition
+
+/**
+ * @brief      2D Delaunay triangulation of some given points using the Bowyer-Watson algorithm.
+ *
+ * @param      in_verts  Input vertexes. Defined by the user.
+ * @param      out_tris  Output triangles. Compute by the function.
+ */
+void triangulation(std::vector<vertex2dc> &in_verts, std::vector<triangle> &out_tris)
+{
+	if (!out_tris.empty()) out_tris.clear();
+	if (in_verts.size() < 3) return;
+
+	// locate the surrounding box and initiate the staring two triangles
+	double xmin = in_verts[0].x, xmax = in_verts[0].x;
+	double ymin = in_verts[0].y, ymax = in_verts[0].y;
+	for (int i = 0; i < in_verts.size(); ++i)
+	{
+		xmin = std::min(xmin, in_verts[i].x);
+		xmax = std::max(xmax, in_verts[i].x);
+		ymin = std::min(ymin, in_verts[i].y);
+		ymax = std::max(ymax, in_verts[i].y);
+	}
+
+	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;
+
+	tmp_vert = new vertex2dc(midx - maxi_s, midy - maxi_s); // lower left corner
+	assit_vert.push_back(tmp_vert);
+
+	tmp_vert = new vertex2dc(midx + maxi_s, midy - maxi_s); // lower right corner
+	assit_vert.push_back(tmp_vert);
+
+	tmp_vert = new vertex2dc(midx + maxi_s, midy + maxi_s); // upper right corner
+	assit_vert.push_back(tmp_vert);
+
+	tmp_vert = new vertex2dc(midx - maxi_s, midy + maxi_s); // upper left corner
+	assit_vert.push_back(tmp_vert);
+
+	triangle *tmp_tri = nullptr;
+	std::vector<triangle*>  exist_tri, 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
+	exist_tri.push_back(tmp_tri);
+
+	tmp_tri = new triangle(assit_vert[0], assit_vert[2], assit_vert[3]); // order the vertex anti-clock wise
+	exist_tri.push_back(tmp_tri);
+
+	// loop all input vertice
+	bool removed;
+	double dist;
+	edge tmp_edge;
+	std::vector<edge> cnst_edge;
+	std::vector<edge>::iterator e_iter;
+	for (int i = 0; i < in_verts.size(); ++i)
+	{
+		// determine triangles that include the point and add the triangle to the cnst_tri and remove the triangle from exist_tri
+		// this is the 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 = exist_tri.begin(); t_iter != exist_tri.end(); )
+		{
+			tmp_tri = *t_iter;
+
+			dist = (tmp_tri->cx - in_verts[i].x) * (tmp_tri->cx - in_verts[i].x) + 
+				(tmp_tri->cy - in_verts[i].y) * (tmp_tri->cy - in_verts[i].y);
+
+			if ((dist - tmp_tri->cr) <= ZERO) // Points on the circumcircle are also included
+			{
+				t_iter = exist_tri.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 exist_tri
+		for (int c = 0; c < cnst_edge.size(); ++c)
+		{
+			tmp_tri = new triangle(cnst_edge[c].vert[0], cnst_edge[c].vert[1], &in_verts[i]); // order the vertex anti-clock wise
+			exist_tri.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;
+		}
+	}
+
+	// remove any triangles has an assistant vertex from exist_tri
+	for (t_iter = exist_tri.begin(); t_iter != exist_tri.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 = exist_tri.erase(t_iter);
+			delete tmp_tri; tmp_tri = nullptr;
+		}
+		else t_iter++;
+	}
+
+	// copy exist_tri to out_tris and destroy memories located
+	out_tris.resize(exist_tri.size());
+	for (int i = 0; i < exist_tri.size(); ++i)
+	{
+		out_tris[i].set(exist_tri[i]->vert[0], exist_tri[i]->vert[1], exist_tri[i]->vert[2]);
+		delete exist_tri[i]; exist_tri[i] = nullptr;
+	}
+
+	// destroy memories located for assit_vert
+	for (int i = 0; i < 4; ++i)
+	{
+		delete assit_vert[i]; assit_vert[i] = nullptr;
+	}
+	return;
+}
+
+#endif // _TIN_DELAUNAY_H
\ No newline at end of file