From 0ba513f2c7a24867d300722050c1f122cf00cf29 Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Sat, 18 Sep 2021 08:56:06 +0800 Subject: [PATCH] tmp update --- delaunay.h | 337 +++++++++++++++++++++++++++++++++++-------------- delaunay_tmp.h | 303 ++++++++++++++++++++++++++++++++++++++++++++ demo.cpp | 65 ++++++++-- demo.msh | 23 ++++ demo.neigh | 8 ++ demo.off | 53 -------- 6 files changed, 631 insertions(+), 158 deletions(-) create mode 100644 delaunay_tmp.h create mode 100644 demo.msh create mode 100644 demo.neigh delete mode 100644 demo.off diff --git a/delaunay.h b/delaunay.h index 2b2e0e5..b9f0c07 100644 --- a/delaunay.h +++ b/delaunay.h @@ -12,6 +12,8 @@ #include "cmath" #include "vector" +#include "iostream" + #define ZERO 1e-5 // Start vertex definition @@ -37,55 +39,63 @@ bool operator ==(const vertex2dc &a, const vertex2dc &b) // overload the == oper } return false; } -// End vertex definition -// Start edge definition -struct edge +void circumcircle(vertex2dc *v0, vertex2dc *v1, vertex2dc *v2, double &cx, double &cy, double &cr) // calculate the circumcircle from three points { - vertex2dc *vert[2]; // vertex of the edge + double s = 0.5 / ((v1->x - v0->x) * (v2->y - v0->y) - (v1->y - v0->y) * (v2->x - v0->x)); + double m = v1->x*v1->x - v0->x*v0->x + v1->y*v1->y - v0->y*v0->y; + double u = v2->x*v2->x - v0->x*v0->x + v2->y*v2->y - v0->y*v0->y; - 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; + cx = ((v2->y - v0->y)*m + (v0->y - v1->y)*u)*s; + cy = ((v0->x - v2->x)*m + (v1->x - v0->x)*u)*s; + cr = (v0->x - cx)*(v0->x - cx) + (v0->y - cy)*(v0->y - cy); // not need to calculate the squared root here + return; } -// End edge definition +// End vertex definition // Start triangle definition struct triangle { + int id; vertex2dc *vert[3]; // vertex of the triangle + triangle *neigh[3]; // neighbors 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() {vert[0] = vert[1] = vert[2] = nullptr; neigh[0] = neigh[1] = neigh[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; + neigh[0] = neigh[1] = neigh[2] = nullptr; - 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 + circumcircle(vert[0], vert[1], vert[2], cx, cy, cr); return; } + + void set_neighbor(triangle *n0ptr, triangle *n1ptr, triangle *n2ptr) + { + neigh[0] = n0ptr; neigh[1] = n1ptr; neigh[2] = n2ptr; + return; + } + + 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++) + { + l1x = vert[(i+1)%3]->x - vert[i]->x; + l1y = vert[(i+1)%3]->y - vert[i]->y; + l2x = inx - vert[i]->x; + l2y = iny - vert[i]->y; + + if ((l1x*l2y - l1y*l2x) < 0) // This condition includes points on the triangle's edge + { + return false; + } + } + return true; + } }; // End triangle definition @@ -95,7 +105,7 @@ struct triangle * @param in_verts Input vertexes. Defined by the user. * @param out_tris Output triangles. Compute by the function. */ -void triangulation(std::vector &in_verts, std::vector &out_tris) +void triangulation(std::vector &in_verts, std::vector &out_tris) { if (!out_tris.empty()) out_tris.clear(); if (in_verts.size() < 3) return; @@ -130,107 +140,242 @@ void triangulation(std::vector &in_verts, std::vector &out_ 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; + triangle *old_tri = nullptr, *tmp_tri = nullptr; + triangle *cnst_tri[3], *old_neigh[6]; std::vector::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); + 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 - exist_tri.push_back(tmp_tri); + out_tris.push_back(tmp_tri); + + out_tris[0]->set_neighbor(nullptr, nullptr, out_tris[1]); + out_tris[1]->set_neighbor(out_tris[0], nullptr, nullptr); // 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 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(); ) + // determine the triangle that includes the new vertex and remove it from out_tris + for (t_iter = out_tris.begin(); t_iter != out_tris.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 + old_tri = *t_iter; + if (old_tri->bound_location(in_verts[i].x, in_verts[i].y)) { - t_iter = exist_tri.erase(t_iter); - cnst_tri.push_back(tmp_tri); + t_iter = out_tris.erase(t_iter); + break; } else t_iter++; } - // loop to remove duplicate edges - cnst_edge.clear(); - for (int c = 0; c < cnst_tri.size(); ++c) + // build three new triangles + for (int n = 0; n < 3; ++n) { - for (int e = 0; e < 3; ++e) + tmp_tri = new triangle(old_tri->vert[n], old_tri->vert[(n+1)%3], &in_verts[i]); + cnst_tri[n] = tmp_tri; + out_tris.push_back(tmp_tri); + } + + // sort neighbors + for (int n = 0; n < 3; ++n) + { + if (old_tri->neigh[n] == nullptr) { - 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(); ) + cnst_tri[n]->set_neighbor(nullptr, cnst_tri[(n+1)%3], cnst_tri[(n+2)%3]); + } + else + { + cnst_tri[n]->set_neighbor(old_tri->neigh[n], cnst_tri[(n+1)%3], cnst_tri[(n+2)%3]); + for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle { - if (tmp_edge == *e_iter) // duplicate edge, remove from cnst_edge + if (old_tri->neigh[n]->neigh[k] == old_tri) { - e_iter = cnst_edge.erase(e_iter); - removed = true; - break; // no need to search more + old_tri->neigh[n]->neigh[k] = cnst_tri[n]; + break; } - 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); - } + // delete the old triangle + delete old_tri; old_tri = nullptr; - // destroy memories used by cnst_edge - for (int c = 0; c < cnst_tri.size(); ++c) + // test if the cnst_tri need to be flipped + for (int n = 0; n < 3; ++n) { - tmp_tri = cnst_tri[c]; - delete tmp_tri; tmp_tri = nullptr; + if (cnst_tri[n]->neigh[0] != nullptr) // must has neighbor on this side + { + old_tri = cnst_tri[n]->neigh[0]; + for (int v = 0; v < 3; ++v) + { + tmp_vert = old_tri->vert[v]; + if (tmp_vert != cnst_tri[n]->vert[0] && tmp_vert != cnst_tri[n]->vert[1]) // find the opposite vertex + { + //dist = (cnst_tri[n]->cx - tmp_vert->x) * (cnst_tri[n]->cx - tmp_vert->x) + + // (cnst_tri[n]->cy - tmp_vert->y) * (cnst_tri[n]->cy - tmp_vert->y); + + //if ((dist - cnst_tri[n]->cr) <= ZERO) // need to be flipped + + dist = (old_tri->cx - cnst_tri[n]->vert[2]->x) * (old_tri->cx - cnst_tri[n]->vert[2]->x) + + (old_tri->cy - cnst_tri[n]->vert[2]->y) * (old_tri->cy - cnst_tri[n]->vert[2]->y); + + if ((dist - old_tri->cr) <= ZERO) // need to be flipped + { + // record the original neighbors + old_neigh[0] = cnst_tri[n]->neigh[0]; + old_neigh[1] = cnst_tri[n]->neigh[1]; + old_neigh[2] = cnst_tri[n]->neigh[2]; + old_neigh[3] = old_tri->neigh[0]; + old_neigh[4] = old_tri->neigh[1]; + old_neigh[5] = old_tri->neigh[2]; + + cnst_tri[n]->set(cnst_tri[n]->vert[0], tmp_vert, cnst_tri[n]->vert[2]); // flip + + if (v == 0) + { + old_tri->set(old_tri->vert[0], old_tri->vert[1], cnst_tri[n]->vert[2]); //flip + + // Sort neighbors + cnst_tri[n]->set_neighbor(old_neigh[5], old_tri, old_neigh[2]); + old_tri->set_neighbor(old_neigh[3], old_neigh[1], cnst_tri[n]); + + for (int k = 0; k < 3; ++k) + { + if (cnst_tri[n]->neigh[0] != nullptr && cnst_tri[n]->neigh[0]->neigh[k] == old_tri) + { + cnst_tri[n]->neigh[0]->neigh[k] = cnst_tri[n]; + break; + } + } + + for (int k = 0; k < 3; ++k) + { + if (old_tri->neigh[1] != nullptr && old_tri->neigh[1]->neigh[k] == cnst_tri[n]) + { + old_tri->neigh[1]->neigh[k] = old_tri; + break; + } + } + } + else if (v == 1) + { + old_tri->set(cnst_tri[n]->vert[2], old_tri->vert[1], old_tri->vert[2]); //flip + + // Sort neighbors + cnst_tri[n]->set_neighbor(old_neigh[3], old_tri, old_neigh[2]); + old_tri->set_neighbor(cnst_tri[n], old_neigh[4], old_neigh[1]); + + for (int k = 0; k < 3; ++k) + { + if (cnst_tri[n]->neigh[0] != nullptr && cnst_tri[n]->neigh[0]->neigh[k] == old_tri) + { + cnst_tri[n]->neigh[0]->neigh[k] = cnst_tri[n]; + break; + } + } + + for (int k = 0; k < 3; ++k) + { + if (old_tri->neigh[2] != nullptr && old_tri->neigh[2]->neigh[k] == cnst_tri[n]) + { + old_tri->neigh[2]->neigh[k] = old_tri; + break; + } + } + } + else + { + old_tri->set(old_tri->vert[0], cnst_tri[n]->vert[2], old_tri->vert[2]); //flip + + // Sort neighbors + cnst_tri[n]->set_neighbor(old_neigh[4], old_tri, old_neigh[2]); + old_tri->set_neighbor(old_neigh[1], cnst_tri[n], old_neigh[5]); + + for (int k = 0; k < 3; ++k) + { + if (cnst_tri[n]->neigh[0] != nullptr && cnst_tri[n]->neigh[0]->neigh[k] == old_tri) + { + cnst_tri[n]->neigh[0]->neigh[k] = cnst_tri[n]; + break; + } + } + + for (int k = 0; k < 3; ++k) + { + if (old_tri->neigh[0] != nullptr && old_tri->neigh[0]->neigh[k] == cnst_tri[n]) + { + old_tri->neigh[0]->neigh[k] = old_tri; + break; + } + } + } + } + break; + } + } + } } } - // remove any triangles has an assistant vertex from exist_tri - for (t_iter = exist_tri.begin(); t_iter != exist_tri.end(); ) + // 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]) + 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]) { + for (int k = 0; k < 3; ++k) + { + if (tmp_tri->neigh[1] != nullptr && tmp_tri->neigh[1]->neigh[k] == tmp_tri) + { + tmp_tri->neigh[1]->neigh[k] = nullptr; + break; + } + } + // destroy the memories located and remove from the vector - t_iter = exist_tri.erase(t_iter); + t_iter = out_tris.erase(t_iter); + delete tmp_tri; tmp_tri = nullptr; + } + else if (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]) + { + for (int k = 0; k < 3; ++k) + { + if (tmp_tri->neigh[2] != nullptr && tmp_tri->neigh[2]->neigh[k] == tmp_tri) + { + tmp_tri->neigh[2]->neigh[k] = nullptr; + break; + } + } + + // destroy the memories located and remove from the vector + t_iter = out_tris.erase(t_iter); + delete tmp_tri; tmp_tri = nullptr; + } + else if (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]) + { + for (int k = 0; k < 3; ++k) + { + if (tmp_tri->neigh[0] != nullptr && tmp_tri->neigh[0]->neigh[k] == tmp_tri) + { + tmp_tri->neigh[0]->neigh[k] = nullptr; + break; + } + } + + // 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++; } - // 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) + // assign triangles index + for (int i = 0; i < out_tris.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; + out_tris[i]->id = i; } - + // destroy memories located for assit_vert for (int i = 0; i < 4; ++i) { @@ -271,9 +416,9 @@ bool duplicated_vertex(const std::vector &in_verts) * * @return If the triangulation is fully delaunay */ -bool fully_delaunay(const std::vector &in_tris, const std::vector &in_verts) +bool fully_delaunay(const std::vector &in_tris, const std::vector &in_verts) { - if (in_tris.empty()) return true; + if (in_tris.empty()) return false; int count; double dist; @@ -282,10 +427,10 @@ bool fully_delaunay(const std::vector &in_tris, const std::vectorcx - in_verts[j].x) * (in_tris[i]->cx - in_verts[j].x) + + (in_tris[i]->cy - in_verts[j].y) * (in_tris[i]->cy - in_verts[j].y); - if ((dist - in_tris[i].cr) <= ZERO) + if ((dist - in_tris[i]->cr) <= ZERO) { count++; } diff --git a/delaunay_tmp.h b/delaunay_tmp.h new file mode 100644 index 0000000..2b2e0e5 --- /dev/null +++ b/delaunay_tmp.h @@ -0,0 +1,303 @@ +/** + * @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" + +#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 &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 four times bigger 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; + std::vector::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 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 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; +} + +/** + * @brief Check for duplicated vertex + * + * @param[in] in_verts Input vertexes + * + * @return If there is duplicated vertex + */ +bool duplicated_vertex(const std::vector &in_verts) +{ + if (in_verts.empty()) return false; + + for (int i = 0; i < in_verts.size()-1; ++i) + { + for (int j = i+1; j < in_verts.size(); ++j) + { + if (in_verts[i] == in_verts[j]) + { + return true; + } + } + } + return false; +} + +/** + * @brief Check to see if the triangulation is fully delaunay + * + * @param[in] in_tris Input triangles + * @param[in] in_verts Input vertexes + * + * @return If the triangulation is fully delaunay + */ +bool fully_delaunay(const std::vector &in_tris, const std::vector &in_verts) +{ + if (in_tris.empty()) return true; + + int count; + double dist; + for (int i = 0; i < in_tris.size(); ++i) + { + count = 0; + for (int j = 0; j < in_verts.size(); ++j) + { + dist = (in_tris[i].cx - in_verts[j].x) * (in_tris[i].cx - in_verts[j].x) + + (in_tris[i].cy - in_verts[j].y) * (in_tris[i].cy - in_verts[j].y); + + if ((dist - in_tris[i].cr) <= ZERO) + { + count++; + } + } + + if (count > 3) + { + return false; + } + } + + return true; +} + +#endif // _BW_2D_DELAUNAY_H \ No newline at end of file diff --git a/demo.cpp b/demo.cpp index e142f77..a11bf29 100644 --- a/demo.cpp +++ b/demo.cpp @@ -1,8 +1,11 @@ #include "delaunay.h" #include "iostream" +#include "fstream" +#include "iomanip" int main(int argc, char const *argv[]) { + /* std::vector points(21); points[0].set(-0.8, -0.8, 0); points[1].set(0.4, -1.2, 1); @@ -25,6 +28,16 @@ int main(int argc, char const *argv[]) points[18].set(2.4, 2.8, 18); points[19].set(3.5, 1.8, 19); points[20].set(3.6, 3.1, 20); + */ + + std::vector points(7); + points[0].set(-0.8, -0.8, 0); + points[1].set(0.4, -1.2, 1); + points[2].set(1.2, 0.9, 2); + points[3].set(-0.4, 0.5, 3); + points[4].set(0.2, -0.15, 4); + points[5].set(0.5, 0.375, 5); + points[6].set(0.7, -0.15, 6); if (duplicated_vertex(points)) { @@ -32,7 +45,7 @@ int main(int argc, char const *argv[]) return 0; } - std::vector elements; + std::vector elements; triangulation(points, elements); if (fully_delaunay(elements, points)) @@ -41,17 +54,51 @@ int main(int argc, char const *argv[]) } else std::clog << "The triangulation is not fully delaunay\n"; - std::cout << "OFF\n"; - std::cout << points.size() << " " << elements.size() << " 0\n"; - - for (int i = 0; i < points.size(); ++i) + // Write a Gmsh's .msh file + std::ofstream outfile("demo.msh"); + outfile << "$MeshFormat" << std::endl << "2.2 0 8" << std::endl << "$EndMeshFormat "<id << " " << elements[i].vert[1]->id << " " << elements[i].vert[2]->id << std::endl; + outfile << i + 1 << " 2 0"; + for (int j = 0; j < 3; j++) + { + outfile << " " << elements[i]->vert[j]->id + 1; + } + outfile << std::endl; + } + outfile << "$EndElements"<< std::endl; + outfile.close(); + + // write a neighbor file + outfile.open("demo.neigh"); + outfile << elements.size() << std::endl; + for (int i = 0; i < elements.size(); i++) + { + outfile << i + 1; + for (int j = 0; j < 3; j++) + { + if (elements[i]->neigh[j] != nullptr) + { + outfile << " " << elements[i]->neigh[j]->id + 1; + } + else outfile << " -1"; + } + outfile << std::endl; + } + outfile.close(); + + // destroy allocated memories + for (int i = 0; i < elements.size(); i++) + { + delete elements[i]; } return 0; } diff --git a/demo.msh b/demo.msh new file mode 100644 index 0000000..fa79819 --- /dev/null +++ b/demo.msh @@ -0,0 +1,23 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +7 +1 -0.8 -0.8 0.0 +2 0.4 -1.2 0.0 +3 1.2 0.9 0.0 +4 -0.4 0.5 0.0 +5 0.2 -0.15 0.0 +6 0.5 0.375 0.0 +7 0.7 -0.15 0.0 +$EndNodes +$Elements +7 +1 2 0 1 2 5 +2 2 0 4 1 5 +3 2 0 3 4 6 +4 2 0 4 5 6 +5 2 0 5 2 7 +6 2 0 3 6 7 +7 2 0 6 5 7 +$EndElements diff --git a/demo.neigh b/demo.neigh new file mode 100644 index 0000000..5ce3eec --- /dev/null +++ b/demo.neigh @@ -0,0 +1,8 @@ +7 +1 -1 5 2 +2 -1 1 4 +3 -1 4 6 +4 2 7 3 +5 1 -1 7 +6 3 7 -1 +7 4 5 6 diff --git a/demo.off b/demo.off deleted file mode 100644 index 9f68ef3..0000000 --- a/demo.off +++ /dev/null @@ -1,53 +0,0 @@ -OFF -21 30 0 --0.8 -0.8 0 -0.4 -1.2 0 -1.2 -0.9 0 -1.6 0.1 0 -2.5 0.5 0 -4.1 0.7 0 -5.7 1.8 0 -5.1 3.4 0 -2.5 4.4 0 -1.2 3.7 0 --1.2 3.9 0 --3.2 5.1 0 --4.3 2.9 0 --3.1 0.7 0 --1.3 0.6 0 --2.1 2.9 0 -0.6 1.2 0 -0.1 2.4 0 -2.4 2.8 0 -3.5 1.8 0 -3.6 3.1 0 -3 1 2 3 -3 3 2 4 -3 4 2 5 -3 9 8 10 -3 10 8 11 -3 13 0 14 -3 10 11 15 -3 11 12 15 -3 12 13 15 -3 13 14 15 -3 0 1 16 -3 1 3 16 -3 3 4 16 -3 14 0 16 -3 9 10 17 -3 10 15 17 -3 15 14 17 -3 14 16 17 -3 8 9 18 -3 16 4 18 -3 9 17 18 -3 17 16 18 -3 5 6 19 -3 6 7 19 -3 4 5 19 -3 18 4 19 -3 7 8 20 -3 8 18 20 -3 18 19 20 -3 19 7 20