### Optimized Mesh::locate_triangle

`Added Mesh::delete_me to cleanup memory. Code needs to be refactored to remove naked new and delete during mesh refinement. I.e. move to pointer free management of nodes and edges. This will also improve cache coherency.`
parent 6a2f26e4
 ... ... @@ -9,7 +9,7 @@ Edge::Edge(Curve *c, bool direction) { Orientation = direction; } Edge::Edge(Curve *c, bool direction, const Point *v) { Edge::Edge(Curve *c, bool direction, Point const *v) { ConstraintCurve = c; Node = v; //Should have C->point(0.0) == *(v) || c->point(1.0) == *(v) Twin = this; ... ... @@ -22,9 +22,9 @@ bool Edge::is_protruding() const { /* See Chapter 1.4 of "Triangulations and Applications" by Øyvind Hjelle and Morten Dæhlen */ const Point *v0 = Prev->Node; const Point *v1 = Node; const Point *v2 = Next->Node; Point const *v0 = Prev->Node; Point const *v1 = Node; Point const *v2 = Next->Node; double v1x = v2->X - v1->X; double v1y = v2->Y - v1->Y; ... ... @@ -40,7 +40,7 @@ bool Edge::is_protruding() const { Edge * e = Next->Next; while (e != Prev) { const Point *p = e->Node; Point const *p = e->Node; double v2x = v2->X - p->X; double v2y = v2->Y - p->Y; ... ... @@ -75,10 +75,10 @@ bool Edge::is_locally_optimal() const { if (ConstraintCurve != nullptr) { return true; } else { const Point *p3 = Node; const Point *p2 = Prev->Node; const Point *p1 = Twin->Node; const Point *p4 = Twin->Prev->Node; Point const *p3 = Node; Point const *p2 = Prev->Node; Point const *p1 = Twin->Node; Point const *p4 = Twin->Prev->Node; double v1x = p3->X - p2->X; double v1y = p3->Y - p2->Y; ... ... @@ -115,7 +115,7 @@ bool Edge::is_valid() const { return value; } bool Edge::is_encroached(const Point *v2) const { bool Edge::is_encroached(Point const *v2) const { /* A constrained edge is encroached if a triangle and it's circumcenter lie on opposite sides of the edge. This is equivalent to a node being in the diameteral ball of the edge? ... ... @@ -127,8 +127,8 @@ bool Edge::is_encroached(const Point *v2) const { if (ConstraintCurve == nullptr) { return false; } else { const Point *v0 = base(); const Point *v1 = tip(); Point const *v0 = base(); Point const *v1 = tip(); double dx0 = v0->X - v2->X; double dy0 = v0->Y - v2->Y; ... ... @@ -143,7 +143,7 @@ bool Edge::is_encroached(const Point *v2) const { } } bool Edge::is_attached(const Point *p, Edge *&e) const { bool Edge::is_attached(Point const *p, Edge *&e) const { if (*this->tip() == *p) { return true; } ... ... @@ -238,12 +238,14 @@ bool Edge::recursive_swap() { } } void Edge::split_edge(std::vector &points, std::vector &edges) { void Edge::split_edge(std::vector &points, std::vector &edges) { /* Splits edge into two edges at the midpoint without creating any new triangles. Used for initial polygon refinement. */ // TODO: Refactor into split_constrainted_edge Point * p = new Point; points.push_back(p); ... ... @@ -284,56 +286,8 @@ void Edge::split_edge(std::vector &points, std::vector &e Next = e; Mark = false; } else { // Interior Edge Edge * e0 = new Edge; edges.push_back(e0); Edge * e1 = new Edge; edges.push_back(e1); // Constraint Curve e0->Orientation = Orientation; e1->Orientation = !Orientation; Twin->Orientation = !Orientation; if (Orientation) { e0->ConstraintCurve = c; } else { e0->ConstraintCurve = ConstraintCurve; ConstraintCurve = c; } Twin->ConstraintCurve = e0->ConstraintCurve; e1->ConstraintCurve = ConstraintCurve; // Connectivity e0->Node = p; e0->Next = Next; e0->Prev = this; e0->Twin = Twin; e0->Mark = false; e1->Node = p; e1->Next = Twin->Next; e1->Prev = Twin; e1->Twin = this; e1->Mark = false; if (Twin->Next != nullptr) { Twin->Next->Prev = e1; Twin->Next->Mark = false; } Twin->Next = e1; Twin->Twin = e0; Twin->Mark = false; if (Next != nullptr) { Next->Prev = e0; Next->Mark = false; } Next = e0; Twin = e1; Mark = false; } else { throw std::exception(); // Function should only be called on } } ... ...
 ... ... @@ -7,7 +7,7 @@ class Edge { public: friend class Mesh; friend bool are_intersecting(const Edge *e0, const Edge *e1); friend bool are_intersecting(Edge const *e0, Edge const *e1); // Constructors Edge() : Node(nullptr), Next(nullptr), Prev(nullptr), Twin(this), ConstraintCurve(nullptr), Orientation(true), Mark(false) {}; ... ... @@ -16,28 +16,28 @@ public: Edge(Curve *c, bool Orientation); Edge(Curve *c, bool Orientation, const Point *v); Edge(Curve *c, bool Orientation, Point const *v); // Accessors const Point *node() const { return Node; }; Point const *node() const { return Node; }; const Point *base() const { return Node; }; Point const *base() const { return Node; }; const Point *tip() const { return (next() == nullptr ? twin()->base() : next()->base()); }; Point const *tip() const { return (next() == nullptr ? twin()->base() : next()->base()); }; const Edge *next() const { return Next; }; Edge const *next() const { return Next; }; const Edge *twin() const { return Twin; }; Edge const *twin() const { return Twin; }; const Edge *prev() const { return Prev; }; Edge const *prev() const { return Prev; }; const Curve *constraint_curve() const { return ConstraintCurve; }; Curve const *constraint_curve() const { return ConstraintCurve; }; bool orientation() const { return Orientation; }; bool mark() const { return Mark; }; bool operator==(const Edge &e) const { bool operator==(Edge const &e) const { return (node() == e.node()) && (constraint_curve() == e.constraint_curve()) && (twin()->node() == e.twin()->node()); }; ... ... @@ -57,20 +57,20 @@ public: bool is_valid() const; bool is_encroached(const Point *p) const; bool is_encroached(Point const *p) const; bool is_attached(const Point *p, Edge *&e) const; bool is_attached(Point const *p, Edge *&e) const; bool swap(); bool recursive_swap(); void split_edge(std::vector &verts, std::vector &edges); void split_edge(std::vector &verts, std::vector &edges); void recursive_mark(); protected: const Point *Node; //Start of edge Point const *Node; //Start of edge Edge *Next; //In triangle Edge *Twin; //Adjacent triangle Edge *Prev; //In triangle ... ...
 ... ... @@ -25,9 +25,21 @@ size_t Mesh::num_edges() const { return count; } void Mesh::delete_me() { for (auto &ed : Edges) { delete ed; } Edges.clear(); for (auto &p : Points) { delete p; } Points.clear(); } // Topological Queries LocateTriangleResult Mesh::locate_triangle(const Point *p, Edge *&e) const { const Edge *ec = e; LocateTriangleResult Mesh::locate_triangle(Point const *p, Edge *&e) const { Edge const *ec = e; LocateTriangleResult result = locate_triangle(p, ec); ... ... @@ -36,121 +48,206 @@ LocateTriangleResult Mesh::locate_triangle(const Point *p, Edge *&e) const { return result; } LocateTriangleResult Mesh::locate_triangle(const Point *p, const Edge *&e) const { LocateTriangleResult Mesh::locate_triangle(Point const *p, Edge const *&e) const { double xp = p->X; double yp = p->Y; while (true) { const Point *p0 = e->node(); const Point *p1 = e->next()->node(); const Point *p2 = e->prev()->node(); const Point *p0 = e->node(); const Point *p1 = e->next()->node(); const Point *p2 = e->prev()->node(); if (*p == *p0) { return LocateTriangleResult::Point; } else if (*p == *p1) { e = e->next(); return LocateTriangleResult::Point; } else if (*p == *p2) { e = e->prev(); return LocateTriangleResult::Point; } if (*p == *p0) { return LocateTriangleResult::Point; } double dx0p = p0->X - xp; double dy0p = p0->Y - yp; double dx1p = p1->X - xp; double dy1p = p1->Y - yp; double dx2p = p2->X - xp; double dy2p = p2->Y - yp; double dx01 = p0->X - p1->X; double dy01 = p0->Y - p1->Y; double dx12 = p1->X - p2->X; double dy12 = p1->Y - p2->Y; double dx20 = p2->X - p0->X; double dy20 = p2->Y - p0->Y; double tol = FLT_EPSILON * (dx20 * dy01 - dy20 * dx01); double area01 = dx0p * dy1p - dx1p * dy0p; double area12 = dx1p * dy2p - dx2p * dy1p; double area20 = dx2p * dy0p - dx0p * dy2p; if (area01 > tol && area12 > tol && area20 > tol) { return LocateTriangleResult::Interior; } else if (area01 < -tol && e != e->twin()) { e = e->twin(); p2 = p1; p1 = p0; p0 = p2; dx2p = dx1p; dx1p = dx0p; dx0p = dx2p; dy2p = dy1p; dy1p = dy0p; dy0p = dy2p; dx01 = -dx01; dy01 = -dy01; } else if (area12 < -tol && e->next() != e->next()->twin()) { e = e->next()->twin(); p0 = p2; dx0p = dx2p; dy0p = dy2p; dx01 = -dx12; dy01 = -dy12; } else if (area20 < -tol && e->prev() != e->prev()->twin()) { e = e->prev()->twin(); p1 = p2; dx1p = dx2p; dy1p = dy2p; dx01 = -dx20; dy01 = -dy20; } else if (area01 > -tol && area12 > tol && area20 > tol) { e = e->twin(); return LocateTriangleResult::Interior; } else if (area01 > tol && area12 > -tol && area20 > tol) { e = e->next()->twin(); return LocateTriangleResult::Interior; } else if (area01 > tol && area12 > tol && area20 > -tol) { e = e->prev()->twin(); return LocateTriangleResult::Interior; } else if (area01 < -tol) { return LocateTriangleResult::Exterior; } else if (area12 < -tol) { e = e->next(); return LocateTriangleResult::Exterior; } else if (area20 < -tol) { e = e->prev(); return LocateTriangleResult::Exterior; } else { throw std::exception(); } if (*p == *p1) { e = e->next(); return LocateTriangleResult::Point; } while (true) { // After first iteration, area01 > 0 p2 = e->prev()->node(); if (*p == *p2) { e = e->prev(); return LocateTriangleResult::Point; } double dx0p = p0->X - xp; double dy0p = p0->Y - yp; double dp0 = sqrt(dx0p * dx0p + dy0p * dy0p); dx2p = p2->X - xp; dy2p = p2->Y - yp; double dx1p = p1->X - xp; double dy1p = p1->Y - yp; double dp1 = sqrt(dx1p * dx1p + dy1p * dy1p); dx12 = p1->X - p2->X; dy12 = p1->Y - p2->Y; double dx2p = p2->X - xp; double dy2p = p2->Y - yp; double dp2 = sqrt(dx2p * dx2p + dy2p * dy2p); dx20 = p2->X - p0->X; dy20 = p2->Y - p0->Y; double dx01 = p0->X - p1->X; double dy01 = p0->Y - p1->Y; double d01 = sqrt(dx01 * dx01 + dy01 * dy01); tol = FLT_EPSILON * (dx20 * dy01 - dy20 * dx01); double dx12 = p1->X - p2->X; double dy12 = p1->Y - p2->Y; double d12 = sqrt(dx12 * dx12 + dy12 * dy12); area12 = dx1p * dy2p - dx2p * dy1p; area20 = dx2p * dy0p - dx0p * dy2p; double dx20 = p2->X - p0->X; double dy20 = p2->Y - p0->Y; double d20 = sqrt(dx20 * dx20 + dy20 * dy20); if (area12 > tol && area20 > tol) { return LocateTriangleResult::Interior; } else if (area12 < -tol && e->next() != e->next()->twin()) { e = e->next()->twin(); p0 = p2; dx0p = dx2p; dy0p = dy2p; dx01 = -dx12; dy01 = -dy12; continue; } else if (area20 < -tol && e->prev() != e->prev()->twin()) { e = e->prev()->twin(); p1 = p2; dx1p = dx2p; dy1p = dy2p; dx01 = -dx20; dy01 = -dy20; continue; } else if (area12 > -tol && area20 > tol) { e = e->next()->twin(); return LocateTriangleResult::Interior; } else if (area12 > tol && area20 > -tol) { e = e->prev()->twin(); return LocateTriangleResult::Interior; } else if (area12 < -tol) { e = e->next(); return LocateTriangleResult::Exterior; } else if (area20 < -tol) { e = e->prev(); return LocateTriangleResult::Exterior; } else { throw std::exception(); } } } double tola = abs(dx20 * dy01 - dy20 * dx01); double tole = FLT_EPSILON * sqrt(tola); tola *= FLT_EPSILON; bool Mesh::in_triangle(Point const *p, Edge const *&e) const { double xp = p->X; double yp = p->Y; double area01 = dx0p * dy1p - dx1p * dy0p; double area12 = dx1p * dy2p - dx2p * dy1p; double area20 = dx2p * dy0p - dx0p * dy2p; Point const *p0 = e->node(); Point const *p1 = e->next()->node(); Point const *p2 = e->prev()->node(); if (area01 > tola && area12 > tola && area20 > tola) { return LocateTriangleResult::Interior; } double dx0p = p0->X - xp; double dy0p = p0->Y - yp; if (area01 < -tola) { if (e != e->twin()) { e = e->twin()->next(); continue; } } else if (area01 < tola) { if (abs(dp0 + dp1 - d01) < tole) { return LocateTriangleResult::Edge; } } double dx1p = p1->X - xp; double dy1p = p1->Y - yp; if (area12 < -tola) { if (e->next() != e->next()->twin()) { e = e->next()->twin()->next(); continue; } } else if (area12 < tola) { if (abs(dp1 + dp2 - d12) < tole) { e = e->next(); return LocateTriangleResult::Edge; } } double dx2p = p2->X - xp; double dy2p = p2->Y - yp; if (area20 < -tola) { if (e->prev() != e->prev()->twin()) { e = e->prev()->twin()->next(); continue; } } else if (area20 < tola) { if (abs(dp2 + dp0 - d20) < tole) { e = e->prev(); return LocateTriangleResult::Edge; } } double dx01 = p0->X - p1->X; double dy01 = p0->Y - p1->Y; // If not Interior, not on Edge, and cannot continue, point is exterior if (area01 < -tola) { return LocateTriangleResult::Exterior; } double dx12 = p1->X - p2->X; double dy12 = p1->Y - p2->Y; if (area12 < -tola) { e = e->next(); return LocateTriangleResult::Exterior; } double dx20 = p2->X - p0->X; double dy20 = p2->Y - p0->Y; if (area20 < -tola) { e = e->prev(); return LocateTriangleResult::Exterior; } double area012 = dx01 * dy12 - dy01 * dx12; throw std::exception(); } double tol = FLT_EPSILON * area012; double area01p = dx0p * dy1p - dx1p * dy0p; double area12p = dx1p * dy2p - dx2p * dy1p; double area20p = dx2p * dy0p - dx0p * dy2p; return (area01p > -tol && area12p > -tol && area20p > -tol); } // Point Insertion InsertPointResult Mesh::insert_point(const Point *vc, Edge *tri) { InsertPointResult Mesh::insert_point(Point const *vc, Edge *tri) { // Find triangle containing point LocateTriangleResult result = locate_triangle(vc, tri); ... ... @@ -161,6 +258,7 @@ InsertPointResult Mesh::insert_point(const Point *vc, Edge *tri) { // Test edges in current and adjacent triangles for encroachment // These are the only possible edges that are encroached due to empty circumcircle property // TODO: write is_encroached method std::vector test_edges; test_edges.reserve(9); ... ... @@ -191,16 +289,16 @@ InsertPointResult Mesh::insert_point(const Point *vc, Edge *tri) { Edge * next = tri->Next; Edge * prev = tri->Prev; const Point *vt = tri->Node; const Point *vn = next->Node; const Point *vp = prev->Node; Point const *vt = tri->Node; Point const *vn = next->Node; Point const *vp = prev->Node; Edge * e0 = new Edge; Edge * e1 = new Edge; Edge * e2 = new Edge; Edge * e3 = new Edge; Edge * e4 = new Edge; Edge * e5 = new Edge; Edge *e0 = new Edge; Edge *e1 = new Edge; Edge *e2 = new Edge; Edge *e3 = new Edge; Edge *e4 = new Edge; Edge *e5 = new Edge; Points.push_back(vc); Edges.push_back(e0); ... ... @@ -262,111 +360,14 @@ InsertPointResult Mesh::insert_point(const Point *vc, Edge *tri) { next->recursive_swap(); prev->recursive_swap(); return InsertPointResult::Success; } else if (result == LocateTriangleResult::Edge) { // Save edges for swapping Edge * next = tri->Next; Edge * prev = tri->Prev; Edge * twin = tri->Twin; Edge * twin_next = tri->Twin->Next; Edge * twin_prev = tri->Twin->Prev; // Get vertex pointers const Point *vt = tri->Node; const Point *vn = next->Node; const Point *vp = prev->Node; const Point *vtp = twin_prev->Node; // Allocate new edges Edge * e0 = new Edge; Edge * e1 = new Edge; Edge * e2 = new Edge; Edge * e3 = new Edge; Edge * e4 = new Edge; Edge * e5 = new Edge; // Save Points.push_back(vc); Edges.push_back(e0); Edges.push_back(e1); Edges.push_back(e2); Edges.push_back(e3); Edges.push_back(e4);