Performance optimization of Triangle insert

src/db/db/dbTriangle.cc

@@ -309,31 +309,24 @@ Triangle::Triangle ()
 }
 
 Triangle::Triangle (TriangleEdge *e1, TriangleEdge *e2, TriangleEdge *e3)
-  : m_is_outside (false), mp_e1 (e1), mp_e2 (e2), mp_e3 (e3), m_id (0)
+  : m_is_outside (false), mp_e1 (e1), m_id (0)
 {
-  mp_v1 = e1->v1 ();
-  mp_v2 = e1->other (mp_v1);
+  mp_v1 = mp_e1->v1 ();
+  mp_v2 = mp_e1->other (mp_v1);
 
   if (e2->has_vertex (mp_v2)) {
-    mp_v3 = e2->other (mp_v2);
-    tl_assert (e3->other (mp_v3) == mp_v1);
+    mp_e2 = e2;
+    mp_e3 = e3;
   } else {
-    mp_v3 = e3->other (mp_v2);
-    tl_assert (e2->other (mp_v3) == mp_v1);
-  }
-
-  //  enforce clockwise orientation
-  if (db::vprod_sign (*mp_v3 - *mp_v1, *mp_v2 - *mp_v1) < 0) {
-    std::swap (mp_v3, mp_v2);
+    mp_e2 = e3;
+    mp_e3 = e2;
   }
+  mp_v3 = mp_e2->other (mp_v2);
 
   //  establish link to edges
   for (int i = 0; i < 3; ++i) {
     TriangleEdge *e = edge (i);
-    int side_of = 0;
-    for (int j = 0; j < 3; ++j) {
-      side_of += e->side_of (*vertex (j));
-    }
+    int side_of = e->side_of (*vertex (i - 1));
     //  NOTE: in the degenerated case, the triangle is not attached to an edge!
     if (side_of < 0) {
       e->set_left (this);
@@ -341,6 +334,11 @@ Triangle::Triangle (TriangleEdge *e1, TriangleEdge *e2, TriangleEdge *e3)
       e->set_right (this);
     }
   }
+
+  //  enforce clockwise orientation
+  if (db::vprod_sign (*mp_v3 - *mp_v1, *mp_v2 - *mp_v1) < 0) {
+    std::swap (mp_v3, mp_v2);
+  }
 }
 
 void

src/db/db/dbTriangles.cc

@@ -34,7 +34,7 @@ namespace db
 {
 
 Triangles::Triangles ()
-  : m_is_constrained (false), m_level (0)
+  : m_is_constrained (false), m_level (0), m_flips (0), m_hops (0)
 {
   //  .. nothing yet ..
 }
@@ -385,10 +385,38 @@ db::TriangleEdge *
 Triangles::find_closest_edge (const db::DPoint &p, db::Vertex *vstart, bool inside_only)
 {
   if (!vstart) {
-    if (! mp_triangles.empty ()) {
-      vstart = mp_triangles.front ()->vertex (0);
+
+    if (! mp_edges.empty ()) {
+
+      unsigned int ls = 0;
+      size_t n = m_vertex_heap.size ();
+      size_t m = n;
+
+      //  A sample heuristics that takes a sqrt(N) sample from the
+      //  vertexes to find a good starting point
+
+      vstart = mp_edges.front ()->v1 ();
+      double dmin = vstart->distance (p);
+
+      while (ls * ls < m) {
+        m /= 2;
+        for (size_t i = m / 2; i < n; i += m) {
+          ++ls;
+          db::Vertex *v = (m_vertex_heap.begin () + i).operator-> ();
+          if (v->begin_edges () != v->end_edges ()) {
+            double d = v->distance (p);
+            if (d < dmin) {
+              vstart = v;
+              dmin = d;
+            }
+          }
+        }
+      }
+
     } else {
+
       return 0;
+
     }
   }
 
@@ -445,6 +473,8 @@ Triangles::find_closest_edge (const db::DPoint &p, db::Vertex *vstart, bool insi
 
     }
 
+    ++m_hops;
+
     v = vnext;
 
   }
@@ -789,11 +819,9 @@ Triangles::remove_inside_vertex (db::Vertex *vertex, std::list<tl::weak_ptr<db::
   fix_triangles (to_fix_a, std::vector<db::TriangleEdge *> (), new_triangles_out);
 }
 
-int
+void
 Triangles::fix_triangles (const std::vector<db::Triangle *> &tris, const std::vector<db::TriangleEdge *> &fixed_edges, std::list<tl::weak_ptr<db::Triangle> > *new_triangles)
 {
-  int flips = 0;
-
   m_level += 1;
   for (auto e = fixed_edges.begin (); e != fixed_edges.end (); ++e) {
     (*e)->set_level (m_level);
@@ -836,7 +864,7 @@ Triangles::fix_triangles (const std::vector<db::Triangle *> &tris, const std::ve
           new_triangles->push_back (t2);
         }
 
-        ++flips;
+        ++m_flips;
         tl_assert (! is_illegal_edge (s12)); // @@@ remove later!
 
         for (int i = 0; i < 3; ++i) {
@@ -858,8 +886,6 @@ Triangles::fix_triangles (const std::vector<db::Triangle *> &tris, const std::ve
     }
 
   }
-
-  return flips;
 }
 
 bool

src/db/db/dbTriangles.h

@@ -31,6 +31,7 @@
 #include "dbRegion.h"
 
 #include "tlObjectCollection.h"
+#include "tlStableVector.h"
 
 namespace db
 {
@@ -230,6 +231,22 @@ public:
    */
   static bool is_illegal_edge (db::TriangleEdge *edge);
 
+  /**
+   *  @brief Statistics: number of flips (fixing)
+   */
+  size_t flips () const
+  {
+    return m_flips;
+  }
+
+  /**
+   *  @brief Statistics: number of hops (searching)
+   */
+  size_t hops () const
+  {
+    return m_hops;
+  }
+
   //  NOTE: these functions are SLOW and intended to test purposes only
   std::vector<db::Vertex *> find_touching (const db::DBox &box) const;
   std::vector<db::Vertex *> find_inside_circle (const db::DPoint &center, double radius) const;
@@ -237,10 +254,11 @@ public:
 private:
   tl::shared_collection<db::Triangle> mp_triangles;
   tl::weak_collection<db::TriangleEdge> mp_edges;
-  std::list<db::Vertex> m_vertex_heap;
+  tl::stable_vector<db::Vertex> m_vertex_heap;
   bool m_is_constrained;
   size_t m_level;
   size_t m_id;
+  size_t m_flips, m_hops;
 
   db::Vertex *create_vertex (double x, double y);
   db::Vertex *create_vertex (const db::DPoint &pt);
@@ -251,7 +269,7 @@ private:
   void remove_outside_vertex (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles = 0);
   void remove_inside_vertex (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out = 0);
   std::vector<db::Triangle *> fill_concave_corners (const std::vector<TriangleEdge *> &edges);
-  int fix_triangles(const std::vector<db::Triangle *> &tris, const std::vector<db::TriangleEdge *> &fixed_edges, std::list<tl::weak_ptr<db::Triangle> > *new_triangles);
+  void fix_triangles (const std::vector<db::Triangle *> &tris, const std::vector<db::TriangleEdge *> &fixed_edges, std::list<tl::weak_ptr<db::Triangle> > *new_triangles);
   std::vector<db::Triangle *> find_triangle_for_point (const db::DPoint &point);
   db::TriangleEdge *find_closest_edge (const db::DPoint &p, db::Vertex *vstart = 0, bool inside_only = false);
   db::Vertex *insert (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles = 0);

src/db/unit_tests/dbTrianglesTests.cc

@@ -258,7 +258,9 @@ TEST(insert_many)
     tris.insert_point (x, y);
   }
 
-  tris.dump ("debug.gds");
+  tl::info << "avg. flips = " << double (tris.flips ()) / double (n);
+  tl::info << "avg. hops = " << double (tris.hops ()) / double (n);
+  // @@@ tris.dump ("debug.gds");
 }
 
 TEST(heavy_insert)