diff --git a/src/db/db/dbTriangle.cc b/src/db/db/dbTriangle.cc
index cb329ccc9..490f27aa1 100644
--- a/src/db/db/dbTriangle.cc
+++ b/src/db/db/dbTriangle.cc
@@ -119,7 +119,7 @@ Vertex::in_circle (const DPoint &point, const DPoint &center, double radius)
   double dy = point.y () - center.y ();
   double d2 = dx * dx + dy * dy;
   double r2 = radius * radius;
-  double delta = std::max (1.0, fabs (d2 + r2)) * db::epsilon;
+  double delta = fabs (d2 + r2) * db::epsilon;
   if (d2 < r2 - delta) {
     return 1;
   } else if (d2 < r2 + delta) {
@@ -427,15 +427,15 @@ Triangle::circumcircle () const
   db::DVector n1 = db::DVector (v1.y (), -v1.x ());
   db::DVector n2 = db::DVector (v2.y (), -v2.x ());
 
-  double p1s = vertex(0)->sq_distance (db::DPoint ());
-  double p2s = vertex(1)->sq_distance (db::DPoint ());
-  double p3s = vertex(2)->sq_distance (db::DPoint ());
+  double p1s = v1.sq_length ();
+  double p2s = v2.sq_length ();
 
   double s = db::vprod (v1, v2);
   tl_assert (fabs (s) > db::epsilon);
 
-  db::DPoint center = db::DPoint () + (n2 * (p1s - p2s) - n1 * (p1s - p3s)) * (0.5 / s);
-  double radius = (*vertex (0) - center).length ();
+  db::DVector r = (n1 * p2s - n2 * p1s) * (0.5 / s);
+  db::DPoint center = *vertex (0) + r;
+  double radius = r.length ();
 
   return std::make_pair (center, radius);
 }
diff --git a/src/db/db/dbTriangle.h b/src/db/db/dbTriangle.h
index 17e7aee63..89a4f7bc6 100644
--- a/src/db/db/dbTriangle.h
+++ b/src/db/db/dbTriangle.h
@@ -395,8 +395,8 @@ private:
   bool m_is_segment;
 
   //  no copying
-  TriangleEdge &operator= (const TriangleEdge &);
-  TriangleEdge (const TriangleEdge &);
+  // @@@ TriangleEdge &operator= (const TriangleEdge &);
+  // @@@ TriangleEdge (const TriangleEdge &);
 
   void set_left  (Triangle *t);
   void set_right (Triangle *t);
diff --git a/src/db/db/dbTriangles.cc b/src/db/db/dbTriangles.cc
index bdad2cd00..0d51148f9 100644
--- a/src/db/db/dbTriangles.cc
+++ b/src/db/db/dbTriangles.cc
@@ -314,19 +314,19 @@ Triangles::find_points_around (db::Vertex *vertex, double radius)
 }
 
 db::Vertex *
-Triangles::insert_point (const db::DPoint &point, std::vector<db::Triangle *> *new_triangles)
+Triangles::insert_point (const db::DPoint &point, std::list<tl::weak_ptr<db::Triangle> > *new_triangles)
 {
   return insert (create_vertex (point), new_triangles);
 }
 
 db::Vertex *
-Triangles::insert_point (db::DCoord x, db::DCoord y, std::vector<db::Triangle *> *new_triangles)
+Triangles::insert_point (db::DCoord x, db::DCoord y, std::list<tl::weak_ptr<db::Triangle> > *new_triangles)
 {
   return insert (create_vertex (x, y), new_triangles);
 }
 
 db::Vertex *
-Triangles::insert (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles)
+Triangles::insert (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles)
 {
   std::vector<db::Triangle *> tris = find_triangle_for_point (*vertex);
 
@@ -452,7 +452,7 @@ Triangles::find_closest_edge (const db::DPoint &p, db::Vertex *vstart, bool insi
 }
 
 void
-Triangles::insert_new_vertex (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles_out)
+Triangles::insert_new_vertex (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out)
 {
   if (mp_triangles.empty ()) {
 
@@ -547,7 +547,7 @@ Triangles::add_more_triangles (std::vector<db::Triangle *> &new_triangles,
 }
 
 void
-Triangles::split_triangle (db::Triangle *t, db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles_out)
+Triangles::split_triangle (db::Triangle *t, db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out)
 {
   t->unlink ();
 
@@ -577,7 +577,7 @@ Triangles::split_triangle (db::Triangle *t, db::Vertex *vertex, std::vector<db::
 }
 
 void
-Triangles::split_triangles_on_edge (const std::vector<db::Triangle *> &tris, db::Vertex *vertex, db::TriangleEdge *split_edge, std::vector<db::Triangle *> *new_triangles_out)
+Triangles::split_triangles_on_edge (const std::vector<db::Triangle *> &tris, db::Vertex *vertex, db::TriangleEdge *split_edge, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out)
 {
   TriangleEdge *s1 = create_edge (split_edge->v1 (), vertex);
   TriangleEdge *s2 = create_edge (split_edge->v2 (), vertex);
@@ -654,7 +654,7 @@ Triangles::find_inside_circle (const db::DPoint &center, double radius) const
 }
 
 void
-Triangles::remove (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles)
+Triangles::remove (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles)
 {
   if (vertex->begin_edges () == vertex->end_edges ()) {
     //  removing an orphan vertex -> ignore
@@ -666,7 +666,7 @@ Triangles::remove (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangle
 }
 
 void
-Triangles::remove_outside_vertex (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles_out)
+Triangles::remove_outside_vertex (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out)
 {
   auto to_remove = vertex->triangles ();
 
@@ -689,7 +689,7 @@ Triangles::remove_outside_vertex (db::Vertex *vertex, std::vector<db::Triangle *
 }
 
 void
-Triangles::remove_inside_vertex (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles_out)
+Triangles::remove_inside_vertex (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out)
 {
   std::set<db::Triangle *, TriangleLessFunc> triangles_to_fix;
 
@@ -789,7 +789,7 @@ Triangles::remove_inside_vertex (db::Vertex *vertex, std::vector<db::Triangle *>
 }
 
 int
-Triangles::fix_triangles (const std::vector<db::Triangle *> &tris, const std::vector<db::TriangleEdge *> &fixed_edges, std::vector<db::Triangle *> *new_triangles)
+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;
 
@@ -1387,4 +1387,172 @@ Triangles::create_constrained_delaunay (const db::Region &region, double dbu)
   constrain (edge_contours);
 }
 
+static bool is_skinny (db::Triangle *tri, const Triangles::TriangulateParameters &param)
+{
+  if (param.b < db::epsilon) {
+    return false;
+  } else {
+    auto cr = tri->circumcircle ();
+    double lmin = tri->edge (0)->d ().sq_length ();
+    for (int i = 1; i < 3; ++i) {
+      lmin = std::min (lmin, tri->edge (i)->d ().sq_length ());
+    }
+    double delta = (fabs (lmin / cr.second) + fabs (param.b)) * db::epsilon;
+    return lmin / cr.second < param.b - delta;
+  }
+}
+
+static bool is_invalid (db::Triangle *tri, const Triangles::TriangulateParameters &param)
+{
+  if (is_skinny (tri, param)) {
+    return true;
+  }
+
+  double amax = param.max_area;
+  if (param.max_area_border > db::epsilon) {
+    bool at_border = false;
+    for (int i = 0; i < 3; ++i) {
+      if (tri->edge (i)->is_segment ()) {
+        at_border = true;
+      }
+    }
+    if (at_border) {
+      amax = param.max_area_border;
+    }
+  }
+
+  if (amax > db::epsilon) {
+    double a = tri->area ();
+    double delta = (fabs (a) + fabs (amax)) * db::epsilon;
+    return tri->area () > amax + delta;
+  } else {
+    return false;
+  }
+}
+
+static unsigned int num_segments (db::Triangle *tri)
+{
+  unsigned int n = 0;
+  for (int i = 0; i < 3; ++i) {
+    if (tri->edge (i)->is_segment ()) {
+      ++n;
+    }
+  }
+  return n;
+}
+
+void
+Triangles::triangulate (const db::Region &region, const TriangulateParameters &parameters, double dbu)
+{
+  create_constrained_delaunay (region, dbu);
+
+  if (parameters.b < db::epsilon && parameters.max_area < db::epsilon && parameters.max_area_border < db::epsilon) {
+
+    //  no refinement requested - we're done.
+    remove_outside_triangles ();
+    return;
+
+  }
+
+  unsigned int nloop = 0;
+  std::list<tl::weak_ptr<db::Triangle> > new_triangles;
+  for (auto t = mp_triangles.begin (); t != mp_triangles.end (); ++t) {
+    new_triangles.push_back (t.operator-> ());
+  }
+
+  // @@@ TODO: break if iteration gets stuck
+  while (nloop < 20) { // @@@
+
+    ++nloop;
+    tl::info << "Iteration " << nloop << " ..";
+
+    std::list<tl::weak_ptr<db::Triangle> > to_consider;
+    for (auto t = new_triangles.begin (); t != new_triangles.end (); ++t) {
+      if (t->get () && ! (*t)->is_outside () && is_invalid (t->get (), parameters)) {
+        to_consider.push_back (*t);
+      }
+    }
+    tl::info << "@@@ to_consider " << to_consider.size();
+
+    if (to_consider.empty ()) {
+      break;
+    }
+
+    new_triangles.clear ();
+
+    for (auto t = to_consider.begin (); t != to_consider.end (); ++t) {
+
+      if (! t->get ()) {
+        //  triangle got removed during loop
+        continue;
+      }
+
+      auto cr = (*t)->circumcircle();
+      auto center = cr.first;
+
+      if ((*t)->contains (center) >= 0) {
+
+        //  heuristics #1: never insert a point into a triangle with more than two segments
+        if (num_segments (t->get ()) <= 1) {
+          //  @@@ print(f"Inserting in-triangle center {repr(center)} of {repr(tri)}")
+          insert_point (center, &new_triangles);
+        }
+
+      } else {
+
+        db::Vertex *vstart = 0;
+        for (int i = 0; i < 3; ++i) {
+          db::TriangleEdge *edge = (*t)->edge (i);
+          vstart = (*t)->opposite (edge);
+          if (edge->side_of (*vstart) * edge->side_of (center) < 0) {
+            break;
+          }
+        }
+
+        db::TriangleEdge *edge = find_closest_edge (center, vstart, true /*inside only*/);
+        tl_assert (edge != 0);
+
+        if (! edge->is_segment () || edge->side_of (*vstart) * edge->side_of (center) >= 0) {
+
+          // @@@ print(f"Inserting out-of-triangle center {repr(center)} of {repr(tri)}")
+          insert_point (center, &new_triangles);
+
+        } else {
+
+          double sr = edge->d ().length () * 0.5;
+          db::DPoint pnew = *edge->v1 () + edge->d () * 0.5;
+
+          // @@@ print(f"split edge {repr(edge)} at {repr(pnew)}")
+          db::Vertex *vnew = insert_point (pnew, &new_triangles);
+          auto vertexes_in_diametral_circle = find_points_around (vnew, sr);
+
+          std::vector<db::Vertex *> to_delete;
+          for (auto v = vertexes_in_diametral_circle.begin (); v != vertexes_in_diametral_circle.end (); ++v) {
+            bool has_segment = false;
+            for (auto e = (*v)->begin_edges (); e != (*v)->end_edges () && ! has_segment; ++e) {
+              has_segment = e->is_segment ();
+            }
+            if (! has_segment) {
+              to_delete.push_back (*v);
+            }
+          }
+
+          // @@@ print(f"  -> found {len(to_delete)} vertexes to remove")
+          for (auto v = to_delete.begin (); v != to_delete.end (); ++v) {
+            remove (*v, &new_triangles);
+          }
+
+        }
+
+      }
+
+    }
+
+    // @@@ tris.dump_as_gdstxt("debug2.txt")
+
+  }
+
+  remove_outside_triangles ();
+}
+
 }
diff --git a/src/db/db/dbTriangles.h b/src/db/db/dbTriangles.h
index ee1b8f17a..4e4743482 100644
--- a/src/db/db/dbTriangles.h
+++ b/src/db/db/dbTriangles.h
@@ -40,6 +40,30 @@ class Layout;
 class DB_PUBLIC Triangles
 {
 public:
+  struct TriangulateParameters
+  {
+    TriangulateParameters ()
+      : b (1.0),
+        max_area (0.0),
+        max_area_border (0.0)
+    { }
+
+    /**
+     *  @brief Min. readius-to-shortest edge ratio
+     */
+    double b;
+
+    /**
+     *  @brief Max area or zero for "no constraint"
+     */
+    double max_area;
+
+    /**
+     *  @brief Max area for border triangles or zero for "use max_area"
+     */
+    double max_area_border;
+  };
+
   typedef tl::shared_collection<db::Triangle> triangles_type;
   typedef triangles_type::const_iterator triangle_iterator;
 
@@ -82,7 +106,13 @@ public:
    */
   void clear ();
 
-  //  exposed for testing purposes:
+  /**
+   *  @brief Creates a refined Delaunay triangulation for the given region
+   */
+  void triangulate (const db::Region &region, const TriangulateParameters &parameters, double dbu = 1.0);
+
+
+  //  -- exposed for testing purposes --
 
   /**
    *  @brief Checks the triangle graph for consistency
@@ -113,7 +143,7 @@ public:
    *  If "new_triangles" is not null, it will receive the list of new triangles created during
    *  the remove step.
    */
-  db::Vertex *insert_point (const db::DPoint &point, std::vector<db::Triangle *> *new_triangles = 0);
+  db::Vertex *insert_point (const db::DPoint &point, std::list<tl::weak_ptr<db::Triangle> > *new_triangles = 0);
 
   /**
    *  @brief Inserts a new vertex as the given point
@@ -121,7 +151,7 @@ public:
    *  If "new_triangles" is not null, it will receive the list of new triangles created during
    *  the remove step.
    */
-  db::Vertex *insert_point (db::DCoord x, db::DCoord y, std::vector<db::Triangle *> *new_triangles = 0);
+  db::Vertex *insert_point (db::DCoord x, db::DCoord y, std::list<tl::weak_ptr<db::Triangle> > *new_triangles = 0);
 
   /**
    *  @brief Removes the given vertex
@@ -129,7 +159,7 @@ public:
    *  If "new_triangles" is not null, it will receive the list of new triangles created during
    *  the remove step.
    */
-  void remove (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles = 0);
+  void remove (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles = 0);
 
   /**
    *  @brief Flips the given edge
@@ -180,6 +210,11 @@ public:
    */
   void create_constrained_delaunay (const db::Region &region, double dbu = 1.0);
 
+  /**
+   *  @brief Returns a value indicating whether the edge is "illegal" (violates the Delaunay criterion)
+   */
+  static bool is_illegal_edge (db::TriangleEdge *edge);
+
   //  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;
@@ -198,21 +233,20 @@ private:
   db::Triangle *create_triangle (db::TriangleEdge *e1, db::TriangleEdge *e2, db::TriangleEdge *e3);
   void remove (db::Triangle *tri);
 
-  void remove_outside_vertex (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles = 0);
-  void remove_inside_vertex (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles_out = 0);
+  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::vector<db::Triangle *> *new_triangles);
-  static bool is_illegal_edge (db::TriangleEdge *edge);
+  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);
   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::vector<db::Triangle *> *new_triangles = 0);
-  void split_triangle (db::Triangle *t, db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles_out);
-  void split_triangles_on_edge (const std::vector<db::Triangle *> &tris, db::Vertex *vertex, db::TriangleEdge *split_edge, std::vector<db::Triangle *> *new_triangles_out);
+  db::Vertex *insert (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles = 0);
+  void split_triangle (db::Triangle *t, db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out);
+  void split_triangles_on_edge (const std::vector<db::Triangle *> &tris, db::Vertex *vertex, db::TriangleEdge *split_edge, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out);
   void add_more_triangles (std::vector<Triangle *> &new_triangles,
                                  db::TriangleEdge *incoming_edge,
                                  db::Vertex *from_vertex, db::Vertex *to_vertex,
                                  db::TriangleEdge *conn_edge);
-  void insert_new_vertex(db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles_out);
+  void insert_new_vertex(db::Vertex *vertex, std::list<tl::weak_ptr<Triangle> > *new_triangles_out);
   std::vector<db::TriangleEdge *> ensure_edge_inner (db::Vertex *from, db::Vertex *to);
   void join_edges (std::vector<TriangleEdge *> &edges);
 };
diff --git a/src/db/unit_tests/dbTrianglesTests.cc b/src/db/unit_tests/dbTrianglesTests.cc
index d915c2608..3807d6cb1 100644
--- a/src/db/unit_tests/dbTrianglesTests.cc
+++ b/src/db/unit_tests/dbTrianglesTests.cc
@@ -29,7 +29,7 @@
 #include <cstdlib>
 #include <cmath>
 
-TEST(Triangle_basic)
+TEST(Triangles_basic)
 {
   db::Triangles tris;
   tris.init_box (db::DBox (1, 0, 5, 4));
@@ -40,7 +40,7 @@ TEST(Triangle_basic)
   EXPECT_EQ (tris.check (), true);
 }
 
-TEST(Triangle_flip)
+TEST(Triangles_flip)
 {
   db::Triangles tris;
   tris.init_box (db::DBox (0, 0, 1, 1));
@@ -61,7 +61,7 @@ TEST(Triangle_flip)
   EXPECT_EQ (tris.check (), true);
 }
 
-TEST(Triangle_insert)
+TEST(Triangles_insert)
 {
   db::Triangles tris;
   tris.init_box (db::DBox (0, 0, 1, 1));
@@ -71,7 +71,7 @@ TEST(Triangle_insert)
   EXPECT_EQ (tris.check (), true);
 }
 
-TEST(Triangle_split_segment)
+TEST(Triangles_split_segment)
 {
   db::Triangles tris;
   tris.init_box (db::DBox (0, 0, 1, 1));
@@ -81,7 +81,7 @@ TEST(Triangle_split_segment)
   EXPECT_EQ (tris.check(), true);
 }
 
-TEST(Triangle_insert_vertex_twice)
+TEST(Triangles_insert_vertex_twice)
 {
   db::Triangles tris;
   tris.init_box (db::DBox (0, 0, 1, 1));
@@ -93,7 +93,7 @@ TEST(Triangle_insert_vertex_twice)
   EXPECT_EQ (tris.check(), true);
 }
 
-TEST(Triangle_insert_vertex_convex)
+TEST(Triangles_insert_vertex_convex)
 {
   db::Triangles tris;
   tris.insert_point (0.2, 0.2);
@@ -105,7 +105,7 @@ TEST(Triangle_insert_vertex_convex)
   EXPECT_EQ (tris.check(), true);
 }
 
-TEST(Triangle_insert_vertex_convex2)
+TEST(Triangles_insert_vertex_convex2)
 {
   db::Triangles tris;
   tris.insert_point (0.25, 0.1);
@@ -116,7 +116,7 @@ TEST(Triangle_insert_vertex_convex2)
   EXPECT_EQ (tris.check(), true);
 }
 
-TEST(Triangle_insert_vertex_convex3)
+TEST(Triangles_insert_vertex_convex3)
 {
   db::Triangles tris;
   tris.insert_point (0.25, 0.5);
@@ -127,7 +127,7 @@ TEST(Triangle_insert_vertex_convex3)
   EXPECT_EQ (tris.check(), true);
 }
 
-TEST(Triangle_search_edges_crossing)
+TEST(Triangles_search_edges_crossing)
 {
   db::Triangles tris;
   db::Vertex *v1 = tris.insert_point (0.2, 0.2);
@@ -147,6 +147,85 @@ TEST(Triangle_search_edges_crossing)
   EXPECT_EQ (std::find (xedges.begin (), xedges.end (), s2) != xedges.end (), true);
 }
 
+TEST(Triangles_illegal_edge1)
+{
+  db::Vertex v1 (0, 0);
+  db::Vertex v2 (1.6, 1.6);
+  db::Vertex v3 (1, 2);
+  db::Vertex v4 (2, 1);
+
+  {
+    db::TriangleEdge e1 (&v1, &v3);
+    db::TriangleEdge e2 (&v3, &v4);
+    db::TriangleEdge e3 (&v4, &v1);
+
+    db::Triangle t1 (&e1, &e2, &e3);
+
+    db::TriangleEdge ee1 (&v2, &v3);
+    db::TriangleEdge ee2 (&v4, &v2);
+
+    db::Triangle t2 (&ee1, &e2, &ee2);
+
+    EXPECT_EQ (db::Triangles::is_illegal_edge (&e2), true);
+  }
+
+  {
+    //  flipped
+    db::TriangleEdge e1 (&v1, &v2);
+    db::TriangleEdge e2 (&v2, &v3);
+    db::TriangleEdge e3 (&v3, &v1);
+
+    db::Triangle t1 (&e1, &e2, &e3);
+
+    db::TriangleEdge ee1 (&v1, &v4);
+    db::TriangleEdge ee2 (&v4, &v2);
+
+    db::Triangle t2 (&ee1, &ee2, &e1);
+
+    EXPECT_EQ (db::Triangles::is_illegal_edge (&e2), false);
+  }
+}
+
+TEST(Triangles_illegal_edge2)
+{
+  //  numerical border case
+  db::Vertex v1 (773.94756216690905, 114.45875269431208);
+  db::Vertex v2 (773.29574734131643, 113.47402096138073);
+  db::Vertex v3 (773.10652961562653, 114.25497975904504);
+  db::Vertex v4 (774.08856345337881, 113.60495072750861);
+
+  {
+    db::TriangleEdge e1 (&v1, &v2);
+    db::TriangleEdge e2 (&v2, &v4);
+    db::TriangleEdge e3 (&v4, &v1);
+
+    db::Triangle t1 (&e1, &e2, &e3);
+
+    db::TriangleEdge ee1 (&v2, &v3);
+    db::TriangleEdge ee2 (&v3, &v4);
+
+    db::Triangle t2 (&ee1, &ee2, &e2);
+
+    EXPECT_EQ (db::Triangles::is_illegal_edge (&e2), false);
+  }
+
+  {
+    //  flipped
+    db::TriangleEdge e1 (&v1, &v2);
+    db::TriangleEdge e2 (&v2, &v3);
+    db::TriangleEdge e3 (&v3, &v1);
+
+    db::Triangle t1 (&e1, &e2, &e3);
+
+    db::TriangleEdge ee1 (&v1, &v4);
+    db::TriangleEdge ee2 (&v4, &v2);
+
+    db::Triangle t2 (&ee1, &ee2, &e1);
+
+    EXPECT_EQ (db::Triangles::is_illegal_edge (&e1), true);
+  }
+}
+
 //  Returns a random float number between 0.0 and 1.0
 inline double flt_rand ()
 {
@@ -163,7 +242,7 @@ namespace {
   };
 }
 
-TEST(Triangle_test_heavy_insert)
+TEST(Triangles_test_heavy_insert)
 {
   tl::info << "Running test_heavy_insert " << tl::noendl;
 
@@ -219,7 +298,7 @@ TEST(Triangle_test_heavy_insert)
   tl::info << tl::endl << "done.";
 }
 
-TEST(Triangle_test_heavy_remove)
+TEST(Triangles_test_heavy_remove)
 {
   tl::info << "Running test_heavy_remove " << tl::noendl;
 
@@ -273,7 +352,7 @@ TEST(Triangle_test_heavy_remove)
   tl::info << tl::endl << "done.";
 }
 
-TEST(Triangle_test_ensure_edge)
+TEST(Triangles_test_ensure_edge)
 {
   srand (0);
 
@@ -339,7 +418,7 @@ bool safe_inside (const db::DBox &b1, const db::DBox &b2)
          (ct::less (b1.top (), b2.top ())       || ct::equal (b1.top (), b2.top ()));
 }
 
-TEST(Triangle_test_constrain)
+TEST(Triangles_test_constrain)
 {
   srand (0);
 
@@ -398,7 +477,7 @@ TEST(Triangle_test_constrain)
   EXPECT_EQ (tl::to_string (area_in), "0.25");
 }
 
-TEST(Triangle_test_heavy_constrain)
+TEST(Triangles_test_heavy_constrain)
 {
   tl::info << "Running test_heavy_constrain " << tl::noendl;
 
@@ -468,7 +547,7 @@ TEST(Triangle_test_heavy_constrain)
   tl::info << tl::endl << "done.";
 }
 
-TEST(Triangle_test_heavy_find_point_around)
+TEST(Triangles_test_heavy_find_point_around)
 {
   tl::info << "Running Triangle_test_heavy_find_point_around " << tl::noendl;
 
@@ -514,7 +593,7 @@ TEST(Triangle_test_heavy_find_point_around)
   tl::info << tl::endl << "done.";
 }
 
-TEST(Triangle_test_create_constrained_delaunay)
+TEST(Triangles_test_create_constrained_delaunay)
 {
   db::Region r;
   r.insert (db::Box (0, 0, 1000, 1000));
@@ -528,6 +607,8 @@ TEST(Triangle_test_create_constrained_delaunay)
   tri.create_constrained_delaunay (r);
   tri.remove_outside_triangles ();
 
+  EXPECT_EQ (tri.check (), true);
+
   EXPECT_EQ (tri.to_string (),
              "((1000, 0), (0, 0), (200, 200)), "
              "((0, 1000), (200, 200), (0, 0)), "
@@ -538,3 +619,45 @@ TEST(Triangle_test_create_constrained_delaunay)
              "((0, 1000), (800, 800), (200, 800)), "
              "((0, 1000), (200, 800), (200, 200))");
 }
+
+TEST(Triangles_test_triangulate)
+{
+  db::Region r;
+  r.insert (db::Box (0, 0, 1000, 1000));
+
+  db::Region r2;
+  r2.insert (db::Box (200, 200, 800, 800));
+
+  r -= r2;
+
+  db::Triangles::TriangulateParameters param;
+  param.max_area = 0.1;
+
+  db::Triangles tri;
+  tri.triangulate (r, param);
+
+  tri.dump ("debug.gds");
+
+  EXPECT_EQ (tri.check (), true);
+
+}
+
+#if 0
+
+assert tris.check(check_delaunay = False)
+
+amax = 0.0
+l2rmin = 2.0
+for tri in tris.triangles:
+if not tri.is_outside:
+  _, radius = tri.circumcircle()
+  lmin = min([math.sqrt(t.square(s.d())) for s in tri.edges()])
+  l2rmin = min(l2rmin, lmin / radius)
+  amax = max(amax, tri.area())
+print(f"max. area = {'%.5g'%amax}")
+print(f"l/R min = {'%.5g'%l2rmin}")
+
+tris.dump_as_gdstxt("out.txt")
+
+
+#endif