initial upload

.gitignore

@@ -30,3 +30,5 @@
 *.exe
 *.out
 *.app
+
+.DS_store

README.md

@@ -0,0 +1,18 @@
+## A C++ implementation of the 2D Delaunay triangulation using the Bowyer-Watson algorithm.
+
+### compile
+
+```shell
+g++ demo.cpp
+```
+
+### run
+
+```shell
+./a.out > demo.off
+```
+
+
+
+
+

delaunay.h

@@ -0,0 +1,236 @@
+/**
+ * @defgroup   DELAUNAY
+ *
+ * @brief      An implementation of the 2D Delaunay triangulation using the Bowyer-Watson algorithm.
+ *
+ * @author     Yi Zhang
+ * @date       2021-09-12
+ */
+
+#ifndef _BW_2D_DELAUNAY_H
+#define _BW_2D_DELAUNAY_H
+#include "cmath"
+#include "vector"
+
+#include "iostream"
+
+#define ZERO 1e-5
+
+// Start vertex definition
+struct vertex2dc
+{
+	unsigned int id;
+	double x, y;
+
+	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;
+	}
+};
+// End vertex definition
+
+// Start edge definition
+struct edge
+{
+	vertex2dc *vert[2];
+
+	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)
+{
+	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];
+	double cx, cy; // triangle's circular center
+	double cr; // triangle's circular radius
+
+	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 triangulation of given points using the Bowyer-Watson algorithm.
+ *
+ * @param      in_verts  Input vertexes
+ * @param      out_tris  Output triangles
+ */
+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 extra wide 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, remain_tri;
+
+	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
+		// Otherwise, add the triangle to the remain_tri
+		cnst_tri.clear();
+		remain_tri.clear();
+		for (int t = 0; t < exist_tri.size(); ++t)
+		{
+			dist = (exist_tri[t]->cx - in_verts[i].x) * (exist_tri[t]->cx - in_verts[i].x) + 
+				(exist_tri[t]->cy - in_verts[i].y) * (exist_tri[t]->cy - in_verts[i].y);
+
+			if ((dist - exist_tri[t]->cr) <= ZERO) // inside // think more later
+			{
+				cnst_tri.push_back(exist_tri[t]);
+			}
+			else
+			{
+				remain_tri.push_back(exist_tri[t]);
+			}
+		}
+
+		// 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(); ++e_iter)
+				{
+					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
+					}
+				}
+
+				if (!removed) // not a duplicate edge, add to the cnst_edge
+				{
+					cnst_edge.push_back(tmp_edge);
+				}
+			}
+		}
+
+		exist_tri.clear();
+		exist_tri.reserve(remain_tri.size() + cnst_edge.size());
+		// copy remain_tri to exist_tri
+		for (int r = 0; r < remain_tri.size(); ++r)
+		{
+			exist_tri.push_back(remain_tri[r]);
+		}
+
+		// 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);
+		}
+	}
+
+	// remove any triangles has an assistant vertex from exist_tri
+	std::vector<triangle*>::iterator t_iter;
+	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 memory 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 memory 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 memory located for assit_vert
+	for (int i = 0; i < 4; ++i)
+	{
+		delete assit_vert[i]; assit_vert[i] = nullptr;
+	}
+	return;
+}
+
+#endif // _BW_2D_DELAUNAY_H

demo.cpp

@@ -0,0 +1,36 @@
+#include "delaunay.h"
+#include "iostream"
+
+int main(int argc, char const *argv[])
+{
+	std::vector<vertex2dc> points(10);
+	points[0].set(-0.7, -0.3, 0);
+	points[1].set(0.1, -0.4, 1);
+	points[2].set(-0.1, -0.1, 2);
+	points[3].set(-0.4, 0.0, 3);
+	points[4].set(-0.4, -0.3, 4);
+	points[5].set(-0.2, -0.1, 5);
+	points[6].set(-0.2, -0.6, 6);
+	points[7].set(-0.2, -0.4, 7);
+	points[8].set(-0.5, -0.5, 8);
+	points[9].set(-0.6, -0.2, 9);
+
+	std::vector<triangle> elements;
+	triangulation(points, elements);
+
+	std::cout << "OFF\n";
+	std::cout << points.size() << " " << elements.size() << " 0\n";
+
+	for (int i = 0; i < points.size(); ++i)
+	{
+		std::cout << points[i].x << " " << points[i].y << " 0\n";
+	}
+
+	for (int i = 0; i < elements.size(); ++i)
+	{
+		std::cout << "3 " << elements[i].vert[0]->id << " " << elements[i].vert[1]->id << " " << elements[i].vert[2]->id << std::endl;
+	}
+	return 0;
+}
+
+

demo.off

@@ -0,0 +1,23 @@
+OFF
+10 11 0
+-0.7 -0.3 0
+0.1 -0.4 0
+-0.1 -0.1 0
+-0.4 0 0
+-0.4 -0.3 0
+-0.2 -0.1 0
+-0.2 -0.6 0
+-0.2 -0.4 0
+-0.5 -0.5 0
+-0.6 -0.2 0
+3 2 3 5
+3 3 4 5
+3 1 2 7
+3 6 1 7
+3 2 5 7
+3 5 4 7
+3 6 7 8
+3 7 4 8
+3 4 3 9
+3 0 8 9
+3 8 4 9