WIP

src/db/db/dbTriangle.cc

@@ -436,22 +436,40 @@ Triangle::bbox () const
 
 
 std::pair<db::DPoint, double>
-Triangle::circumcircle () const
+Triangle::circumcircle (bool *ok) const
 {
-  db::DVector v1 = *mp_v[0] - *mp_v[1];
-  db::DVector v2 = *mp_v[0] - *mp_v[2];
-  db::DVector n1 = db::DVector (v1.y (), -v1.x ());
-  db::DVector n2 = db::DVector (v2.y (), -v2.x ());
+  //  see https://en.wikipedia.org/wiki/Circumcircle
+  //  we set A=(0,0), so the formulas simplify
 
-  double p1s = v1.sq_length ();
-  double p2s = v2.sq_length ();
+  if (ok) {
+    *ok = true;
+  }
 
-  double s = db::vprod (v1, v2);
-  tl_assert (fabs (s) > db::epsilon);
+  db::DVector b = *mp_v[1] - *mp_v[0];
+  db::DVector c = *mp_v[2] - *mp_v[0];
 
-  db::DVector r = (n1 * p2s - n2 * p1s) * (0.5 / s);
-  db::DPoint center = *mp_v[0] + r;
-  double radius = r.length ();
+  double b2 = b.sq_length ();
+  double c2 = c.sq_length ();
+
+  double sx = 0.5 * (b2 * c.y () - c2 * b.y ());
+  double sy = 0.5 * (b.x () * c2 - c.x() * b2);
+
+  double a1 = b.x() * c.y();
+  double a2 = c.x() * b.y();
+  double a = a1 - a2;
+  double a_abs = std::abs (a);
+
+  if (a_abs < (std::abs (a1) + std::abs (a2)) * db::epsilon) {
+    if (ok) {
+      *ok = false;
+      return std::make_pair (db::DPoint (), 0.0);
+    } else {
+      tl_assert (false);
+    }
+  }
+
+  double radius = sqrt (sx * sx + sy * sy) / a_abs;
+  db::DPoint center = *mp_v[0] + db::DVector (sx / a, sy / a);
 
   return std::make_pair (center, radius);
 }
@@ -507,26 +525,22 @@ Triangle::common_edge (const Triangle *other) const
 int
 Triangle::contains (const db::DPoint &point) const
 {
-  //  the relative distance of the point from the edge
-  double alpha = db::epsilon;
-
-  double d = db::vprod (*mp_v[2] - *mp_v[0], *mp_v[1] - *mp_v[0]);
-  if (fabs (d) < db::epsilon * db::epsilon) {
-    return -1;
-  }
+  int vps = db::vprod_sign (*mp_v[2] - *mp_v[0], *mp_v[1] - *mp_v[0]);
 
   int res = 1;
+
   const Vertex *vl = mp_v[2];
   for (int i = 0; i < 3; ++i) {
     const Vertex *v = mp_v[i];
-    double n = db::vprod (point - *vl, *v - *vl) / d;
-    if (n < -alpha) {
+    int n = db::vprod_sign (point - *vl, *v - *vl) * vps;
+    if (n < 0) {
       return -1;
-    } else if (n < alpha) {
+    } else if (n == 0) {
       res = 0;
     }
     vl = v;
   }
+
   return res;
 }
 
@@ -544,8 +558,9 @@ double
 Triangle::b () const
 {
   double lmin = min_edge_length ();
-  auto cr = circumcircle ();
-  return lmin / cr.second;
+  bool ok = false;
+  auto cr = circumcircle (&ok);
+  return ok ? lmin / cr.second : 0.0;
 }
 
 bool

src/db/db/dbTriangle.h

@@ -483,8 +483,10 @@ public:
 
   /**
    *  @brief Gets the center point and radius of the circumcircle
+   *  If ok is non-null, it will receive a boolean value indicating whether the circumcircle is valid.
+   *  An invalid circumcircle is an indicator for a degenerated triangle with area 0 (or close to).
    */
-  std::pair<db::DPoint, double> circumcircle () const;
+  std::pair<db::DPoint, double> circumcircle (bool *ok = 0) const;
 
   /**
    *  @brief Gets the vertex opposite of the given edge

src/db/db/dbTriangles.cc

@@ -937,13 +937,15 @@ Triangles::is_illegal_edge (db::TriangleEdge *edge)
     return false;
   }
 
-  auto lr = left->circumcircle ();
-  if (right->opposite (edge)->in_circle (lr.first, lr.second) > 0) {
+  bool ok = false;
+
+  auto lr = left->circumcircle (&ok);
+  if (! ok || right->opposite (edge)->in_circle (lr.first, lr.second) > 0) {
     return true;
   }
 
-  auto rr = right->circumcircle();
-  if (left->opposite (edge)->in_circle (rr.first, rr.second) > 0) {
+  auto rr = right->circumcircle(&ok);
+  if (! ok || left->opposite (edge)->in_circle (rr.first, rr.second) > 0) {
     return true;
   }
 

src/db/unit_tests/dbTrianglesTests.cc

@@ -951,6 +951,48 @@ TEST(triangulate_thin)
 
   EXPECT_EQ (tri.check (false), true);
 
-  // for debugging:
-  tri.dump ("debug.gds");
+  //  for debugging:
+  //  tri.dump ("debug.gds");
+
+  for (auto t = tri.begin (); t != tri.end (); ++t) {
+    EXPECT_GE (t->b (), param.min_b);
+  }
+
+  EXPECT_GT (tri.num_triangles (), size_t (13000));
+  EXPECT_LT (tri.num_triangles (), size_t (13200));
+}
+
+TEST(triangulate_issue1996)
+{
+  db::DPoint contour[] = {
+    db::DPoint (-8000, -8075),
+    db::DPoint (-8000, 8075),
+    db::DPoint (18000, 8075),
+    db::DPoint (18000, -8075)
+  };
+
+  db::DPolygon poly;
+  poly.assign_hull (contour + 0, contour + sizeof (contour) / sizeof (contour[0]));
+
+  double dbu = 0.001;
+
+  db::Triangles::TriangulateParameters param;
+  param.min_b = 0.5;
+  param.max_area = 500.0 * dbu * dbu;
+
+  TestableTriangles tri;
+  db::DCplxTrans trans = db::DCplxTrans (dbu) * db::DCplxTrans (db::DTrans (db::DPoint () - poly.box ().center ()));
+  tri.triangulate (trans * poly, param);
+
+  EXPECT_EQ (tri.check (false), true);
+
+  //  for debugging:
+  //  tri.dump ("debug.gds");
+
+  for (auto t = tri.begin (); t != tri.end (); ++t) {
+    EXPECT_GE (t->b (), param.min_b);
+  }
+
+  EXPECT_GT (tri.num_triangles (), size_t (13000));
+  EXPECT_LT (tri.num_triangles (), size_t (13200));
 }