diff --git a/README.md b/README.md index abced24..86db9a9 100644 --- a/README.md +++ b/README.md @@ -9,10 +9,5 @@ g++ demo.cpp ### run ```shell -./a.out > demo.off -``` - - - - - +./a.out +``` \ No newline at end of file diff --git a/delaunay.h b/delaunay.h index dcb6e70..160ea4a 100644 --- a/delaunay.h +++ b/delaunay.h @@ -12,8 +12,6 @@ #include "cmath" #include "vector" -#include "iostream" - #define ZERO 1e-5 // Start vertex definition @@ -39,63 +37,55 @@ bool operator ==(const vertex2dc &a, const vertex2dc &b) // overload the == oper } return false; } - -void circumcircle(vertex2dc *v0, vertex2dc *v1, vertex2dc *v2, double &cx, double &cy, double &cr) // calculate the circumcircle from three points -{ - 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; - - 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 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 { - 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; neigh[0] = neigh[1] = neigh[2] = nullptr;} + 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; - neigh[0] = neigh[1] = neigh[2] = nullptr; - circumcircle(vert[0], vert[1], vert[2], cx, cy, cr); + 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; } - - 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 @@ -105,7 +95,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; @@ -128,261 +118,119 @@ void triangulation(std::vector &in_verts, std::vector &out vertex2dc *tmp_vert = nullptr; std::vector assit_vert; - tmp_vert = new vertex2dc(midx - maxi_s, midy - maxi_s, 99990); // lower left corner + 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, 99991); // lower right corner + 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, 99992); // upper right corner + 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, 99993); // upper left corner + tmp_vert = new vertex2dc(midx - maxi_s, midy + maxi_s); // upper left corner assit_vert.push_back(tmp_vert); - triangle *old_tri = nullptr, *tmp_tri = nullptr; - triangle *cnst_tri[3], *old_neigh[6]; + 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 - out_tris.push_back(tmp_tri); + 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 - 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); + 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 the triangle that includes the new vertex and remove it from out_tris - for (t_iter = out_tris.begin(); t_iter != out_tris.end(); ) + // 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(); ) { - old_tri = *t_iter; - if (old_tri->bound_location(in_verts[i].x, in_verts[i].y)) + 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 = out_tris.erase(t_iter); - break; + t_iter = exist_tri.erase(t_iter); + cnst_tri.push_back(tmp_tri); } else t_iter++; } - // build three new triangles - for (int n = 0; n < 3; ++n) + // loop to remove duplicate edges + cnst_edge.clear(); + for (int c = 0; c < cnst_tri.size(); ++c) { - 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); - } + for (int e = 0; e < 3; ++e) + { + tmp_edge.set(cnst_tri[c]->vert[e], cnst_tri[c]->vert[(e+1)%3]); - // sort neighbors - for (int n = 0; n < 3; ++n) - { - if (old_tri->neigh[n] == nullptr) - { - 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 + removed = false; + for (e_iter = cnst_edge.begin(); e_iter != cnst_edge.end(); ) { - if (old_tri->neigh[n]->neigh[k] == old_tri) + if (tmp_edge == *e_iter) // duplicate edge, remove from cnst_edge { - old_tri->neigh[n]->neigh[k] = cnst_tri[n]; - break; + 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); } } } - // delete the old triangle - delete old_tri; old_tri = nullptr; - - // test if the cnst_tri need to be flipped - for (int n = 0; n < 3; ++n) + // construct new triangles and add to exist_tri + for (int c = 0; c < cnst_edge.size(); ++c) { - 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; - } - } - } + 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); } - std::cout << "Insert time: " << i + 1 << std::endl; - for (int e = 0; e < out_tris.size(); ++e) + // destroy memories used by cnst_edge + for (int c = 0; c < cnst_tri.size(); ++c) { - std::cout << out_tris[e]->vert[0]->id << " " << out_tris[e]->vert[1]->id << " " << out_tris[e]->vert[2]->id << std::endl; + tmp_tri = cnst_tri[c]; + delete tmp_tri; tmp_tri = nullptr; } - std::cout << "===================\n"; } - // remove any triangles has an assistant vertex from out_tris - for (t_iter = out_tris.begin(); t_iter != out_tris.end(); ) + // 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]) + 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]) { - 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 = 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); + t_iter = exist_tri.erase(t_iter); delete tmp_tri; tmp_tri = nullptr; } else t_iter++; } - // assign triangles index - for (int i = 0; i < out_tris.size(); i++) + // 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]->id = 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) { @@ -406,7 +254,7 @@ bool duplicated_vertex(const std::vector &in_verts) { for (int j = i+1; j < in_verts.size(); ++j) { - if (in_verts[i] == in_verts[j]) + if (in_verts[i] == in_verts[j] && in_verts[i].id != in_verts[j].id) { return true; } @@ -423,9 +271,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 false; + if (in_tris.empty()) return true; int count; double dist; @@ -434,10 +282,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); + 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) + if ((dist - in_tris[i].cr) <= ZERO) { count++; } diff --git a/delaunay_backup.h b/delaunay_backup.h deleted file mode 100644 index 9f169d3..0000000 --- a/delaunay_backup.h +++ /dev/null @@ -1,312 +0,0 @@ -/** - * @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; // 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, 99990); // lower left corner - assit_vert.push_back(tmp_vert); - - tmp_vert = new vertex2dc(midx + maxi_s, midy - maxi_s, 99991); // lower right corner - assit_vert.push_back(tmp_vert); - - tmp_vert = new vertex2dc(midx + maxi_s, midy + maxi_s, 99992); // upper right corner - assit_vert.push_back(tmp_vert); - - tmp_vert = new vertex2dc(midx - maxi_s, midy + maxi_s, 99993); // 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; - } - - std::cout << "Insert time: " << i + 1 << std::endl; - for (int e = 0; e < exist_tri.size(); ++e) - { - std::cout << exist_tri[e]->vert[0]->id << " " << exist_tri[e]->vert[1]->id << " " << exist_tri[e]->vert[2]->id << std::endl; - } - std::cout << "===================\n"; - } - - // 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 dde631d..4d68961 100644 --- a/demo.cpp +++ b/demo.cpp @@ -5,7 +5,7 @@ int main(int argc, char const *argv[]) { - std::vector points(7); + std::vector points(21); points[0].set(-0.8, -0.8, 0); points[1].set(0.4, -1.2, 1); points[2].set(1.2, -0.9, 2); @@ -13,7 +13,6 @@ int main(int argc, char const *argv[]) points[4].set(2.5, 0.5, 4); points[5].set(4.1, 0.7, 5); points[6].set(5.7, 1.8, 6); - /* points[7].set(5.1, 3.4, 7); points[8].set(2.5, 4.4, 8); points[9].set(1.2, 3.7, 9); @@ -28,26 +27,14 @@ 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)) { - std::cerr << "Duplicated vertice detected.\n"; + std::cerr << "Duplicated vertexes detected.\n"; return 0; } - std::vector elements; + std::vector elements; triangulation(points, elements); if (fully_delaunay(elements, points)) @@ -72,36 +59,12 @@ int main(int argc, char const *argv[]) outfile << i + 1 << " 2 0"; for (int j = 0; j < 3; j++) { - outfile << " " << elements[i]->vert[j]->id + 1; + 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 index 88efeda..fb98d1d 100644 --- a/demo.msh +++ b/demo.msh @@ -2,7 +2,7 @@ $MeshFormat 2.2 0 8 $EndMeshFormat $Nodes -7 +21 1 -0.8 -0.8 0.0 2 0.4 -1.2 0.0 3 1.2 -0.9 0.0 @@ -10,11 +10,51 @@ $Nodes 5 2.5 0.5 0.0 6 4.1 0.7 0.0 7 5.7 1.8 0.0 +8 5.1 3.4 0.0 +9 2.5 4.4 0.0 +10 1.2 3.7 0.0 +11 -1.2 3.9 0.0 +12 -3.2 5.1 0.0 +13 -4.3 2.9 0.0 +14 -3.1 0.7 0.0 +15 -1.3 0.6 0.0 +16 -2.1 2.9 0.0 +17 0.6 1.2 0.0 +18 0.1 2.4 0.0 +19 2.4 2.8 0.0 +20 3.5 1.8 0.0 +21 3.6 3.1 0.0 $EndNodes $Elements -4 -1 2 0 3 4 2 +30 +1 2 0 2 3 4 2 2 0 4 3 5 -3 2 0 7 5 6 -4 2 0 5 3 6 +3 2 0 5 3 6 +4 2 0 10 9 11 +5 2 0 11 9 12 +6 2 0 14 1 15 +7 2 0 11 12 16 +8 2 0 12 13 16 +9 2 0 13 14 16 +10 2 0 14 15 16 +11 2 0 1 2 17 +12 2 0 2 4 17 +13 2 0 4 5 17 +14 2 0 15 1 17 +15 2 0 10 11 18 +16 2 0 11 16 18 +17 2 0 16 15 18 +18 2 0 15 17 18 +19 2 0 9 10 19 +20 2 0 17 5 19 +21 2 0 10 18 19 +22 2 0 18 17 19 +23 2 0 6 7 20 +24 2 0 7 8 20 +25 2 0 5 6 20 +26 2 0 19 5 20 +27 2 0 8 9 21 +28 2 0 9 19 21 +29 2 0 19 20 21 +30 2 0 20 8 21 $EndElements diff --git a/demo.neigh b/demo.neigh deleted file mode 100644 index 5e310ec..0000000 --- a/demo.neigh +++ /dev/null @@ -1,5 +0,0 @@ -4 -1 2 -1 -1 -2 1 4 -1 -3 -1 4 -1 -4 2 -1 3 diff --git a/demo_back.off b/demo_back.off deleted file mode 100644 index 0d5e006..0000000 --- a/demo_back.off +++ /dev/null @@ -1,14 +0,0 @@ -OFF -7 5 0 --0.8 -0.8 0.0 -0.4 -1.2 0.0 -1.2 -0.9 0.0 -1.6 0.1 0.0 -2.5 0.5 0.0 -4.1 0.7 0.0 -5.7 1.8 0.0 -3 0 1 3 -3 1 2 3 -3 3 2 4 -3 4 2 5 -3 4 5 6 \ No newline at end of file diff --git a/demo_backup.cpp b/demo_backup.cpp deleted file mode 100644 index d3c62fa..0000000 --- a/demo_backup.cpp +++ /dev/null @@ -1,76 +0,0 @@ -#include "delaunay_backup.h" -#include "iostream" -#include "fstream" -#include "iomanip" - -int main(int argc, char const *argv[]) -{ - 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(1.6, 0.1, 3); - points[4].set(2.5, 0.5, 4); - points[5].set(4.1, 0.7, 5); - points[6].set(5.7, 1.8, 6); - /* - points[7].set(5.1, 3.4, 7); - points[8].set(2.5, 4.4, 8); - points[9].set(1.2, 3.7, 9); - points[10].set(-1.2, 3.9, 10); - points[11].set(-3.2, 5.1, 11); - points[12].set(-4.3, 2.9, 12); - points[13].set(-3.1, 0.7, 13); - points[14].set(-1.3, 0.6, 14); - points[15].set(-2.1, 2.9, 15); - points[16].set(0.6, 1.2, 16); - points[17].set(0.1, 2.4, 17); - 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)) - { - std::cerr << "Duplicated vertice detected.\n"; - return 0; - } - - std::vector elements; - triangulation(points, elements); - - if (fully_delaunay(elements, points)) - { - std::clog << "The triangulation is fully delaunay.\n"; - } - else std::clog << "The triangulation is not fully delaunay\n"; - - // Write a OFF file - std::cout << "OFF\n"; - std::cout << points.size() << " " << elements.size() << " 0\n"; - for (int i = 0; i < points.size(); i++) - { - std::cout << std::setprecision(16) << points[i].x << " " << points[i].y << " 0.0" << std::endl; - } - - 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/log.txt b/log.txt deleted file mode 100644 index fd45a0b..0000000 --- a/log.txt +++ /dev/null @@ -1,84 +0,0 @@ -Insert time: 1 -0 99991 99992 -99990 99991 0 -99992 99993 0 -99993 99990 0 -=================== -Insert time: 2 -99990 99991 1 -1 99993 0 -99993 99990 0 -0 99990 1 -99991 99992 1 -99992 99993 1 -=================== -Insert time: 3 -99990 99991 1 -1 99993 0 -99993 99990 0 -0 99990 1 -2 99993 1 -99991 99992 2 -99992 99993 2 -1 99991 2 -=================== -Insert time: 4 -99990 99991 1 -1 99993 0 -99993 99990 0 -0 99990 1 -2 3 1 -99991 99992 3 -1 99991 2 -99992 99993 3 -99993 1 3 -2 99991 3 -=================== -Insert time: 5 -99990 99991 1 -1 99993 0 -99993 99990 0 -0 99990 1 -2 3 1 -1 99991 2 -4 99993 3 -99993 1 3 -2 99991 4 -99991 99992 4 -99992 99993 4 -3 2 4 -=================== -Insert time: 6 -99990 99991 1 -1 99993 0 -99993 99990 0 -0 99990 1 -2 3 1 -1 99991 2 -4 99993 3 -99993 1 3 -2 99991 5 -99992 99993 4 -3 2 4 -99991 99992 5 -99992 4 5 -4 2 5 -=================== -Insert time: 7 -99990 99991 1 -1 99993 0 -99993 99990 0 -0 99990 1 -2 3 1 -1 99991 2 -4 99993 3 -99993 1 3 -2 99991 5 -99992 99993 4 -3 2 4 -6 4 5 -4 2 5 -99991 99992 6 -99992 4 6 -5 99991 6 -=================== diff --git a/log_backup.txt b/log_backup.txt deleted file mode 100644 index a97cf8a..0000000 --- a/log_backup.txt +++ /dev/null @@ -1,98 +0,0 @@ -Insert time: 1 -99990 99991 0 -99991 99992 0 -99992 99993 0 -99993 99990 0 -=================== -Insert time: 2 -99993 99990 0 -99990 99991 1 -0 99990 1 -99991 99992 1 -99992 99993 1 -99993 0 1 -=================== -Insert time: 3 -99993 99990 0 -99990 99991 1 -0 99990 1 -99991 99992 2 -1 99991 2 -99992 99993 2 -99993 0 2 -0 1 2 -=================== -Insert time: 4 -99993 99990 0 -99990 99991 1 -0 99990 1 -1 99991 2 -99991 99992 3 -2 99991 3 -99992 99993 3 -99993 0 3 -0 1 3 -1 2 3 -=================== -Insert time: 5 -99993 99990 0 -99990 99991 1 -0 99990 1 -1 99991 2 -99993 0 3 -0 1 3 -1 2 3 -99991 99992 4 -2 99991 4 -3 2 4 -99992 99993 4 -99993 3 4 -=================== -Insert time: 6 -99993 99990 0 -99990 99991 1 -0 99990 1 -1 99991 2 -99993 0 3 -0 1 3 -1 2 3 -3 2 4 -99992 99993 4 -99993 3 4 -99991 99992 5 -99992 4 5 -2 99991 5 -4 2 5 -=================== -Insert time: 7 -99993 99990 0 -99990 99991 1 -0 99990 1 -1 99991 2 -99993 0 3 -0 1 3 -1 2 3 -3 2 4 -99993 3 4 -2 99991 5 -4 2 5 -99992 99993 6 -99993 4 6 -99991 99992 6 -5 99991 6 -4 5 6 -=================== -OFF -7 5 0 --0.8 -0.8 0.0 -0.4 -1.2 0.0 -1.2 -0.9 0.0 -1.6 0.1 0.0 -2.5 0.5 0.0 -4.1 0.7 0.0 -5.7 1.8 0.0 -3 0 1 3 -3 1 2 3 -3 3 2 4 -3 4 2 5 -3 4 5 6