diff --git a/.gitignore b/.gitignore index 259148f..32192a8 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..abced24 --- /dev/null +++ b/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 +``` + + + + + diff --git a/delaunay.h b/delaunay.h new file mode 100644 index 0000000..3db46f2 --- /dev/null +++ b/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 &in_verts, std::vector &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 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 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 cnst_edge; + std::vector::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::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 \ No newline at end of file diff --git a/demo.cpp b/demo.cpp new file mode 100644 index 0000000..ff837b8 --- /dev/null +++ b/demo.cpp @@ -0,0 +1,36 @@ +#include "delaunay.h" +#include "iostream" + +int main(int argc, char const *argv[]) +{ + std::vector 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 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; +} + + diff --git a/demo.off b/demo.off new file mode 100644 index 0000000..11b7d72 --- /dev/null +++ b/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