diff --git a/src/db/db/db.pro b/src/db/db/db.pro index 732d42aa1..888a22dcc 100644 --- a/src/db/db/db.pro +++ b/src/db/db/db.pro @@ -70,13 +70,14 @@ SOURCES = \ dbNetlistSpiceReaderExpressionParser.cc \ dbObject.cc \ dbObjectWithProperties.cc \ + dbPLC.cc \ + dbPLCTriangulation.cc \ dbPath.cc \ dbPCellDeclaration.cc \ dbPCellHeader.cc \ dbPCellVariant.cc \ dbPoint.cc \ dbPolygon.cc \ - dbPolygonGraph.cc \ dbPolygonNeighborhood.cc \ dbPolygonTools.cc \ dbPolygonGenerators.cc \ @@ -309,13 +310,14 @@ HEADERS = \ dbObject.h \ dbObjectTag.h \ dbObjectWithProperties.h \ + dbPLC.h \ + dbPLCTriangulation.h \ dbPath.h \ dbPCellDeclaration.h \ dbPCellHeader.h \ dbPCellVariant.h \ dbPoint.h \ dbPolygon.h \ - dbPolygonGraph.h \ dbPolygonNeighborhood.h \ dbPolygonTools.h \ dbPolygonGenerators.h \ diff --git a/src/db/db/dbPolygonGraph.cc b/src/db/db/dbPLC.cc similarity index 68% rename from src/db/db/dbPolygonGraph.cc rename to src/db/db/dbPLC.cc index f745b0da2..559bf4572 100644 --- a/src/db/db/dbPolygonGraph.cc +++ b/src/db/db/dbPLC.cc @@ -21,7 +21,7 @@ */ -#include "dbPolygonGraph.h" +#include "dbPLC.h" #include "dbLayout.h" #include "dbWriter.h" #include "tlStream.h" @@ -36,28 +36,43 @@ namespace db { +namespace plc +{ + // ------------------------------------------------------------------------------------- -// GVertex implementation +// Vertex implementation -GVertex::GVertex () - : DPoint (), m_is_precious (false) +Vertex::Vertex (Graph *graph) + : DPoint (), mp_graph (graph), m_is_precious (false) { // .. nothing yet .. } -GVertex::GVertex (const db::DPoint &p) - : DPoint (p), m_is_precious (false) +Vertex::Vertex (Graph *graph, const db::DPoint &p) + : DPoint (p), mp_graph (graph), m_is_precious (false) { // .. nothing yet .. } -GVertex::GVertex (const GVertex &v) - : DPoint (), m_is_precious (false) +Vertex::Vertex (Graph *graph, const Vertex &v) + : DPoint (), mp_graph (graph), m_is_precious (false) { operator= (v); } -GVertex &GVertex::operator= (const GVertex &v) +Vertex::Vertex (Graph *graph, db::DCoord x, db::DCoord y) + : DPoint (x, y), mp_graph (graph), m_is_precious (false) +{ + // .. nothing yet .. +} + +Vertex::Vertex (const Vertex &v) + : DPoint (v), mp_graph (v.mp_graph), m_is_precious (v.m_is_precious) +{ + // NOTE: edges are not copied! +} + +Vertex &Vertex::operator= (const Vertex &v) { if (this != &v) { // NOTE: edges are not copied! @@ -67,15 +82,8 @@ GVertex &GVertex::operator= (const GVertex &v) return *this; } -GVertex::GVertex (db::DCoord x, db::DCoord y) - : DPoint (x, y), m_is_precious (false) -{ - // .. nothing yet .. -} - -#if 0 // @@@ bool -GVertex::is_outside () const +Vertex::is_outside () const { for (auto e = mp_edges.begin (); e != mp_edges.end (); ++e) { if ((*e)->is_outside ()) { @@ -84,13 +92,12 @@ GVertex::is_outside () const } return false; } -#endif -std::vector -GVertex::polygons () const +std::vector +Vertex::polygons () const { - std::set seen; - std::vector res; + std::set seen; + std::vector res; for (auto e = mp_edges.begin (); e != mp_edges.end (); ++e) { for (auto t = (*e)->begin_polygons (); t != (*e)->end_polygons (); ++t) { if (seen.insert (t.operator-> ()).second) { @@ -102,7 +109,7 @@ GVertex::polygons () const } bool -GVertex::has_edge (const GPolygonEdge *edge) const +Vertex::has_edge (const Edge *edge) const { for (auto e = mp_edges.begin (); e != mp_edges.end (); ++e) { if (*e == edge) { @@ -113,7 +120,7 @@ GVertex::has_edge (const GPolygonEdge *edge) const } size_t -GVertex::num_edges (int max_count) const +Vertex::num_edges (int max_count) const { if (max_count < 0) { // NOTE: this can be slow for a std::list, so we have max_count to limit this effort @@ -128,7 +135,7 @@ GVertex::num_edges (int max_count) const } std::string -GVertex::to_string (bool with_id) const +Vertex::to_string (bool with_id) const { std::string res = tl::sprintf ("(%.12g, %.12g)", x (), y()); if (with_id) { @@ -138,7 +145,7 @@ GVertex::to_string (bool with_id) const } int -GVertex::in_circle (const DPoint &point, const DPoint ¢er, double radius) +Vertex::in_circle (const DPoint &point, const DPoint ¢er, double radius) { double dx = point.x () - center.x (); double dy = point.y () - center.y (); @@ -155,34 +162,34 @@ GVertex::in_circle (const DPoint &point, const DPoint ¢er, double radius) } // ------------------------------------------------------------------------------------- -// GPolygonEdge implementation +// Edge implementation -GPolygonEdge::GPolygonEdge () - : mp_v1 (0), mp_v2 (0), mp_left (), mp_right (), m_level (0), m_id (0), m_is_segment (false) +Edge::Edge (Graph *graph) + : mp_graph (graph), mp_v1 (0), mp_v2 (0), mp_left (), mp_right (), m_level (0), m_id (0), m_is_segment (false) { // .. nothing yet .. } -GPolygonEdge::GPolygonEdge (GVertex *v1, GVertex *v2) - : mp_v1 (v1), mp_v2 (v2), mp_left (), mp_right (), m_level (0), m_id (0), m_is_segment (false) +Edge::Edge (Graph *graph, Vertex *v1, Vertex *v2) + : mp_graph (graph), mp_v1 (v1), mp_v2 (v2), mp_left (), mp_right (), m_level (0), m_id (0), m_is_segment (false) { // .. nothing yet .. } void -GPolygonEdge::set_left (GPolygon *t) +Edge::set_left (Polygon *t) { mp_left = t; } void -GPolygonEdge::set_right (GPolygon *t) +Edge::set_right (Polygon *t) { mp_right = t; } void -GPolygonEdge::link () +Edge::link () { mp_v1->mp_edges.push_back (this); m_ec_v1 = --mp_v1->mp_edges.end (); @@ -192,7 +199,7 @@ GPolygonEdge::link () } void -GPolygonEdge::unlink () +Edge::unlink () { if (mp_v1) { mp_v1->remove_edge (m_ec_v1); @@ -203,8 +210,8 @@ GPolygonEdge::unlink () mp_v1 = mp_v2 = 0; } -GPolygon * -GPolygonEdge::other (const GPolygon *t) const +Polygon * +Edge::other (const Polygon *t) const { if (t == mp_left) { return mp_right; @@ -216,8 +223,8 @@ GPolygonEdge::other (const GPolygon *t) const return 0; } -GVertex * -GPolygonEdge::other (const GVertex *t) const +Vertex * +Edge::other (const Vertex *t) const { if (t == mp_v1) { return mp_v2; @@ -230,13 +237,13 @@ GPolygonEdge::other (const GVertex *t) const } bool -GPolygonEdge::has_vertex (const GVertex *v) const +Edge::has_vertex (const Vertex *v) const { return mp_v1 == v || mp_v2 == v; } -GVertex * -GPolygonEdge::common_vertex (const GPolygonEdge *other) const +Vertex * +Edge::common_vertex (const Edge *other) const { if (has_vertex (other->v1 ())) { return (other->v1 ()); @@ -248,7 +255,7 @@ GPolygonEdge::common_vertex (const GPolygonEdge *other) const } std::string -GPolygonEdge::to_string (bool with_id) const +Edge::to_string (bool with_id) const { std::string res = std::string ("(") + mp_v1->to_string (with_id) + ", " + mp_v2->to_string (with_id) + ")"; if (with_id) { @@ -258,7 +265,7 @@ GPolygonEdge::to_string (bool with_id) const } double -GPolygonEdge::distance (const db::DEdge &e, const db::DPoint &p) +Edge::distance (const db::DEdge &e, const db::DPoint &p) { double l = db::sprod (p - e.p1 (), e.d ()) / e.d ().sq_length (); db::DPoint pp; @@ -273,27 +280,27 @@ GPolygonEdge::distance (const db::DEdge &e, const db::DPoint &p) } bool -GPolygonEdge::crosses (const db::DEdge &e, const db::DEdge &other) +Edge::crosses (const db::DEdge &e, const db::DEdge &other) { return e.side_of (other.p1 ()) * e.side_of (other.p2 ()) < 0 && other.side_of (e.p1 ()) * other.side_of (e.p2 ()) < 0; } bool -GPolygonEdge::crosses_including (const db::DEdge &e, const db::DEdge &other) +Edge::crosses_including (const db::DEdge &e, const db::DEdge &other) { return e.side_of (other.p1 ()) * e.side_of (other.p2 ()) <= 0 && other.side_of (e.p1 ()) * other.side_of (e.p2 ()) <= 0; } db::DPoint -GPolygonEdge::intersection_point (const db::DEdge &e, const db::DEdge &other) +Edge::intersection_point (const db::DEdge &e, const db::DEdge &other) { return e.intersect_point (other).second; } bool -GPolygonEdge::point_on (const db::DEdge &edge, const db::DPoint &point) +Edge::point_on (const db::DEdge &edge, const db::DPoint &point) { if (edge.side_of (point) != 0) { return false; @@ -302,42 +309,72 @@ GPolygonEdge::point_on (const db::DEdge &edge, const db::DPoint &point) } } -#if 0 // @@@ bool -GPolygonEdge::is_for_outside_polygons () const +Edge::can_flip () const +{ + if (! left () || ! right ()) { + return false; + } + + const Vertex *v1 = left ()->opposite (this); + const Vertex *v2 = right ()->opposite (this); + return crosses (db::DEdge (*v1, *v2)); +} + +bool +Edge::can_join_via (const Vertex *vertex) const +{ + if (! left () || ! right ()) { + return false; + } + + tl_assert (has_vertex (vertex)); + const Vertex *v1 = left ()->opposite (this); + const Vertex *v2 = right ()->opposite (this); + return db::DEdge (*v1, *v2).side_of (*vertex) == 0; +} + +bool +Edge::is_outside () const +{ + return left () == 0 || right () == 0; +} + +bool +Edge::is_for_outside_triangles () const { return (left () && left ()->is_outside ()) || (right () && right ()->is_outside ()); } -#endif bool -GPolygonEdge::has_polygon (const GPolygon *t) const +Edge::has_polygon (const Polygon *t) const { return t != 0 && (left () == t || right () == t); } // ------------------------------------------------------------------------------------- -// GPolygon implementation +// Polygon implementation -GPolygon::GPolygon () - : m_id (0) +Polygon::Polygon (Graph *graph) + : mp_graph (graph), m_is_outside (false), m_id (0) { // .. nothing yet .. } void -GPolygon::init () +Polygon::init () { m_id = 0; + m_is_outside = false; if (mp_e.empty ()) { return; } - std::vector e; + std::vector e; e.swap (mp_e); - std::multimap v2e; + std::multimap v2e; for (auto i = e.begin (); i != e.end (); ++i) { if (i != e.begin ()) { @@ -379,7 +416,7 @@ GPolygon::init () // establish clockwise order of the vertexes double area = 0.0; - const db::GVertex *vm1 = vertex (-1), *v0; + const Vertex *vm1 = vertex (-1), *v0; for (auto i = mp_v.begin (); i != mp_v.end (); ++i) { v0 = *i; area += db::vprod (*vm1 - db::DPoint (), *v0 - *vm1); @@ -394,8 +431,8 @@ GPolygon::init () // link the polygon to the edges for (size_t i = 0; i < size (); ++i) { - db::GVertex *v = mp_v[i]; - db::GPolygonEdge *e = mp_e[i]; + Vertex *v = mp_v[i]; + Edge *e = mp_e[i]; if (e->v1 () == v) { e->set_right (this); } else { @@ -404,13 +441,62 @@ GPolygon::init () } } -GPolygon::~GPolygon () +Polygon::Polygon (Graph *graph, Edge *e1, Edge *e2, Edge *e3) + : mp_graph (graph), m_is_outside (false), m_id (0) +{ + mp_e.resize (3, 0); + mp_v.resize (3, 0); + + mp_e[0] = e1; + mp_v[0] = e1->v1 (); + mp_v[1] = e1->v2 (); + + if (e2->has_vertex (mp_v[1])) { + mp_e[1] = e2; + mp_e[2] = e3; + } else { + mp_e[1] = e3; + mp_e[2] = e2; + } + mp_v[2] = mp_e[1]->other (mp_v[1]); + + // enforce clockwise orientation + int s = db::vprod_sign (*mp_v[2] - *mp_v[0], *mp_v[1] - *mp_v[0]); + if (s < 0) { + std::swap (mp_v[2], mp_v[1]); + } else if (s == 0) { + // Triangle is not orientable + tl_assert (false); + } + + // establish link to edges + for (int i = 0; i < 3; ++i) { + + Edge *e = mp_e[i]; + + unsigned int i1 = 0; + for ( ; e->v1 () != mp_v[i1] && i1 < 3; ++i1) + ; + unsigned int i2 = 0; + for ( ; e->v2 () != mp_v[i2] && i2 < 3; ++i2) + ; + + if ((i1 + 1) % 3 == i2) { + e->set_right (this); + } else { + e->set_left (this); + } + + } +} + +Polygon::~Polygon () { unlink (); } void -GPolygon::unlink () +Polygon::unlink () { for (auto e = mp_e.begin (); e != mp_e.end (); ++e) { if ((*e)->left () == this) { @@ -423,7 +509,7 @@ GPolygon::unlink () } std::string -GPolygon::to_string (bool with_id) const +Polygon::to_string (bool with_id) const { std::string res = "("; for (int i = 0; i < int (size ()); ++i) { @@ -441,13 +527,13 @@ GPolygon::to_string (bool with_id) const } double -GPolygon::area () const +Polygon::area () const { return fabs (db::vprod (mp_e[0]->d (), mp_e[1]->d ())) * 0.5; } db::DBox -GPolygon::bbox () const +Polygon::bbox () const { db::DBox box; for (auto i = mp_v.begin (); i != mp_v.end (); ++i) { @@ -456,8 +542,77 @@ GPolygon::bbox () const return box; } -GPolygonEdge * -GPolygon::find_edge_with (const GVertex *v1, const GVertex *v2) const +std::pair +Polygon::circumcircle (bool *ok) const +{ + tl_assert (mp_v.size () == 3); + + // see https://en.wikipedia.org/wiki/Circumcircle + // we set A=(0,0), so the formulas simplify + + if (ok) { + *ok = true; + } + + db::DVector b = *mp_v[1] - *mp_v[0]; + db::DVector c = *mp_v[2] - *mp_v[0]; + + 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); +} + +Vertex * +Polygon::opposite (const Edge *edge) const +{ + tl_assert (mp_v.size () == 3); + + for (int i = 0; i < 3; ++i) { + Vertex *v = mp_v[i]; + if (! edge->has_vertex (v)) { + return v; + } + } + tl_assert (false); +} + +Edge * +Polygon::opposite (const Vertex *vertex) const +{ + tl_assert (mp_v.size () == 3); + + for (int i = 0; i < 3; ++i) { + Edge *e = mp_e[i]; + if (! e->has_vertex (vertex)) { + return e; + } + } + tl_assert (false); +} + +Edge * +Polygon::find_edge_with (const Vertex *v1, const Vertex *v2) const { for (auto e = mp_e.begin (); e != mp_e.end (); ++e) { if ((*e)->has_vertex (v1) && (*e)->has_vertex (v2)) { @@ -467,8 +622,8 @@ GPolygon::find_edge_with (const GVertex *v1, const GVertex *v2) const tl_assert (false); } -GPolygonEdge * -GPolygon::common_edge (const GPolygon *other) const +Edge * +Polygon::common_edge (const Polygon *other) const { for (auto e = mp_e.begin (); e != mp_e.end (); ++e) { if ((*e)->other (this) == other) { @@ -478,10 +633,11 @@ GPolygon::common_edge (const GPolygon *other) const return 0; } -#if 0 // @@@ int -GPolygon::contains (const db::DPoint &point) const +Polygon::contains (const db::DPoint &point) const { + tl_assert (mp_v.size () == 3); + auto c = *mp_v[2] - *mp_v[0]; auto b = *mp_v[1] - *mp_v[0]; @@ -492,9 +648,9 @@ GPolygon::contains (const db::DPoint &point) const int res = 1; - const GVertex *vl = mp_v[2]; + const Vertex *vl = mp_v[2]; for (int i = 0; i < 3; ++i) { - const GVertex *v = mp_v[i]; + const Vertex *v = mp_v[i]; int n = db::vprod_sign (point - *vl, *v - *vl) * vps; if (n < 0) { return -1; @@ -506,10 +662,9 @@ GPolygon::contains (const db::DPoint &point) const return res; } -#endif double -GPolygon::min_edge_length () const +Polygon::min_edge_length () const { double lmin = mp_e[0]->d ().length (); for (auto e = mp_e.begin (); e != mp_e.end (); ++e) { @@ -518,19 +673,17 @@ GPolygon::min_edge_length () const return lmin; } -#if 0 // @@@ double -GPolygon::b () const +Polygon::b () const { double lmin = min_edge_length (); bool ok = false; auto cr = circumcircle (&ok); return ok ? lmin / cr.second : 0.0; } -#endif bool -GPolygon::has_segment () const +Polygon::has_segment () const { for (auto e = mp_e.begin (); e != mp_e.end (); ++e) { if ((*e)->is_segment ()) { @@ -541,7 +694,7 @@ GPolygon::has_segment () const } unsigned int -GPolygon::num_segments () const +Polygon::num_segments () const { unsigned int n = 0; for (auto e = mp_e.begin (); e != mp_e.end (); ++e) { @@ -554,49 +707,42 @@ GPolygon::num_segments () const // ----------------------------------------------------------------------------------- -static inline bool is_equal (const db::DPoint &a, const db::DPoint &b) -{ - return std::abs (a.x () - b.x ()) < std::max (1.0, (std::abs (a.x ()) + std::abs (b.x ()))) * db::epsilon && - std::abs (a.y () - b.y ()) < std::max (1.0, (std::abs (a.y ()) + std::abs (b.y ()))) * db::epsilon; -} - -PolygonGraph::PolygonGraph () +Graph::Graph () : m_id (0) - // @@@: m_is_constrained (false), m_level (0), m_id (0), m_flips (0), m_hops (0) { // .. nothing yet .. } -PolygonGraph::~PolygonGraph () +Graph::~Graph () { clear (); } -db::GVertex * -PolygonGraph::create_vertex (double x, double y) +Vertex * +Graph::create_vertex (double x, double y) { - m_vertex_heap.push_back (db::GVertex (x, y)); + m_vertex_heap.push_back (Vertex (this, x, y)); return &m_vertex_heap.back (); } -db::GVertex * -PolygonGraph::create_vertex (const db::DPoint &pt) +Vertex * +Graph::create_vertex (const db::DPoint &pt) { - m_vertex_heap.push_back (pt); + m_vertex_heap.push_back (Vertex (this, pt)); return &m_vertex_heap.back (); } -db::GPolygonEdge * -PolygonGraph::create_edge (db::GVertex *v1, db::GVertex *v2) +Edge * +Graph::create_edge (Vertex *v1, Vertex *v2) { - db::GPolygonEdge *edge = 0; + Edge *edge = 0; if (! m_returned_edges.empty ()) { edge = m_returned_edges.back (); m_returned_edges.pop_back (); - *edge = db::GPolygonEdge (v1, v2); + *edge = Edge (this, v1, v2); } else { - m_edges_heap.push_back (db::GPolygonEdge (v1, v2)); + m_edges_heap.push_back (Edge (this, v1, v2)); edge = &m_edges_heap.back (); } @@ -605,11 +751,21 @@ PolygonGraph::create_edge (db::GVertex *v1, db::GVertex *v2) return edge; } -void -PolygonGraph::remove_polygon (db::GPolygon *poly) +Polygon * +Graph::create_triangle (Edge *e1, Edge *e2, Edge *e3) { - std::vector edges; - edges.reserve (poly->size ()); + Polygon *res = new Polygon (this, e1, e2, e3); + res->set_id (++m_id); + mp_polygons.push_back (res); + + return res; +} + +void +Graph::remove_polygon (Polygon *poly) +{ + std::vector edges; + edges.resize (poly->size (), 0); for (int i = 0; i < int (poly->size ()); ++i) { edges [i] = poly->edge (i); } @@ -625,63 +781,8 @@ PolygonGraph::remove_polygon (db::GPolygon *poly) } } -#if 0 // @@@ -void -PolygonGraph::convex_decompose (const db::DPolygon &polygon) -{ - clear (); - - if (polygon.begin_edge ().at_end ()) { - return; - } - - std::vector edges; - - for (unsigned int c = 0; c < polygon.holes () + 1; ++c) { - - const db::DSimplePolygon::contour_type &ctr = polygon.contour (c); - - db::GVertex *v0 = 0, *vv, *v; - size_t n = ctr.size (); - for (size_t i = 0; i < n; ++i) { - - db::DPoint pm1 = ctr [i > 0 ? i - 1 : n - 1]; - db::DPoint pp1 = ctr [i + 1 < n ? i + 1 : 0]; - db::DPoint p = ctr [i]; - - bool is_convex = db::vprod_sign (p - pm1, pp1 - p); - // @@@ - - v = create_vertex (p.x (), p.y ()); - if (! v0) { - v0 = v; - } else { - edges.push_back (create_edge (vv, v)); - } - - vv = v; - - } - - if (v0 && v0 != v) { - edges.push_back (create_edge (v, v0)); - } - - } -} -#endif - -void -PolygonGraph::convex_decompose (const db::DPolygon &poly) -{ - - // @@@ - - -} - std::string -PolygonGraph::to_string () +Graph::to_string () { std::string res; for (auto t = mp_polygons.begin (); t != mp_polygons.end (); ++t) { @@ -694,7 +795,7 @@ PolygonGraph::to_string () } db::DBox -PolygonGraph::bbox () const +Graph::bbox () const { db::DBox box; for (auto t = mp_polygons.begin (); t != mp_polygons.end (); ++t) { @@ -705,7 +806,7 @@ PolygonGraph::bbox () const #if 0 // @@@ bool -PolygonGraph::check (bool check_delaunay) const +Graph::check (bool check_delaunay) const { bool res = true; @@ -759,7 +860,7 @@ PolygonGraph::check (bool check_delaunay) const tl::error << "(check error) edges " << e->to_string (true) << " vertex 2 not found in adjacent polygon " << t->to_string (true); res = false; } - db::GVertex *vopp = t->opposite (e.operator-> ()); + Vertex *vopp = t->opposite (e.operator-> ()); double sgn = (e->left () == t.operator-> ()) ? 1.0 : -1.0; double vp = db::vprod (e->d(), *vopp - *e->v1 ()); // positive if on left side if (vp * sgn <= 0.0) { @@ -803,7 +904,7 @@ PolygonGraph::check (bool check_delaunay) const #endif db::Layout * -PolygonGraph::to_layout (bool decompose_by_id) const +Graph::to_layout (bool decompose_by_id) const { db::Layout *layout = new db::Layout (); layout->dbu (0.001); @@ -850,7 +951,7 @@ PolygonGraph::to_layout (bool decompose_by_id) const } void -PolygonGraph::dump (const std::string &path, bool decompose_by_id) const +Graph::dump (const std::string &path, bool decompose_by_id) const { std::unique_ptr ly (to_layout (decompose_by_id)); @@ -860,36 +961,19 @@ PolygonGraph::dump (const std::string &path, bool decompose_by_id) const db::Writer writer (opt); writer.write (*ly, stream); - tl::info << "PolygonGraph written to " << path; + tl::info << "Graph written to " << path; } void -PolygonGraph::clear () +Graph::clear () { mp_polygons.clear (); m_edges_heap.clear (); m_vertex_heap.clear (); m_returned_edges.clear (); - // @@@m_is_constrained = false; - // @@@m_level = 0; m_id = 0; } -template -void -PolygonGraph::make_contours (const Poly &poly, const Trans &trans, std::vector > &edge_contours) -{ - edge_contours.push_back (std::vector ()); - for (auto pt = poly.begin_hull (); pt != poly.end_hull (); ++pt) { - edge_contours.back ().push_back (insert_point (trans * *pt)); - } +} // namespace plc - for (unsigned int h = 0; h < poly.holes (); ++h) { - edge_contours.push_back (std::vector ()); - for (auto pt = poly.begin_hole (h); pt != poly.end_hole (h); ++pt) { - edge_contours.back ().push_back (insert_point (trans * *pt)); - } - } -} - -} +} // namespace db diff --git a/src/db/db/dbPolygonGraph.h b/src/db/db/dbPLC.h similarity index 64% rename from src/db/db/dbPolygonGraph.h rename to src/db/db/dbPLC.h index d296b886f..8d410d0b8 100644 --- a/src/db/db/dbPolygonGraph.h +++ b/src/db/db/dbPLC.h @@ -20,8 +20,8 @@ */ -#ifndef HDR_dbPolygonGraph -#define HDR_dbPolygonGraph +#ifndef HDR_dbPLC +#define HDR_dbPLC #include "dbCommon.h" #include "dbTriangle.h" @@ -41,8 +41,38 @@ namespace db class Layout; -class GPolygon; -class GPolygonEdge; +namespace plc +{ + +/** + * @brief A framework for piecewise linear curves + * + * This framework implements classes for dealing with piecewise linear + * curves. It is the basis for triangulation and polygon decomposition + * algorithms. + * + * The core class is the PLCGraph which is a collection of vertices, + * edges and edge loops (polygons). Vertices, edges and polygons form + * graphs. + * + * A "vertex" (db::plc::Vertex) is a point. A point connects two or + * more edges. + * A vertex has some attributes: + * * 'precious': if set, the vertex is not removed during triangulation + * for example. + * + * An "edge" (db::plc::Edge) is a line connecting two vertexes. The + * edge runs from vertex v1 to vertex v2. An edge separates two + * polygons (left and right, as seen in the run direction). + * + * A "segment" as an edge that is part of an original polygon outline. + * + * A "polygon" (db::plc::Polygon) is a loop of edges. + */ + +class Polygon; +class Edge; +class Graph; /** * @brief A class representing a vertex in a Delaunay triangulation graph @@ -51,40 +81,71 @@ class GPolygonEdge; * an integer value that can be used in traversal algorithms * ("level") */ -class DB_PUBLIC GVertex +class DB_PUBLIC Vertex : public db::DPoint { public: - typedef std::list edges_type; + typedef std::list edges_type; typedef edges_type::const_iterator edges_iterator; typedef edges_type::iterator edges_iterator_non_const; - GVertex (); - GVertex (const DPoint &p); - GVertex (const GVertex &v); - GVertex (db::DCoord x, db::DCoord y); + Vertex (const Vertex &v); + Vertex &operator= (const Vertex &v); - GVertex &operator= (const GVertex &v); - -#if 0 // @@@ + /** + * @brief Gets a value indicating whether any of the attached edges is "outside" + */ bool is_outside () const; -#endif - std::vector polygons () const; + /** + * @brief Gets a list of polygons that are attached to this vertex + */ + std::vector polygons() const; + + /** + * @brief Gets the graph object this vertex belongs to + */ + Graph *graph () const { return mp_graph; } + + /** + * @brief Iterates the edges on this vertex (begin) + */ edges_iterator begin_edges () const { return mp_edges.begin (); } + + /** + * @brief Iterates the edges on this vertex (end) + */ edges_iterator end_edges () const { return mp_edges.end (); } + + /** + * @brief Returns the number of edges attached to this vertex + */ size_t num_edges (int max_count = -1) const; - bool has_edge (const GPolygonEdge *edge) const; + /** + * @brief Returns a value indicating whether the given edge is attached to this vertex + */ + bool has_edge (const Edge *edge) const; + /** + * @brief Sets a valid indicating whether the vertex is precious + * + * "precious" vertexes are not removed during triangulation for example. + */ void set_is_precious (bool f) { m_is_precious = f; } + + /** + * @brief Gets a valid indicating whether the vertex is precious + */ bool is_precious () const { return m_is_precious; } + /** + * @brief Returns a string representation of the vertex + */ std::string to_string (bool with_id = false) const; /** * @brief Returns 1 is the point is inside the circle, 0 if on the circle and -1 if outside - * TODO: Move to db::DPoint */ static int in_circle (const db::DPoint &point, const db::DPoint ¢er, double radius); @@ -97,13 +158,20 @@ public: } private: - friend class GPolygonEdge; + friend class Edge; + friend class Graph; + + Vertex (Graph *graph); + Vertex (Graph *graph, const DPoint &p); + Vertex (Graph *graph, const Vertex &v); + Vertex (Graph *graph, db::DCoord x, db::DCoord y); void remove_edge (const edges_iterator_non_const &ec) { mp_edges.erase (ec); } + Graph *mp_graph; edges_type mp_edges; bool m_is_precious; }; @@ -111,15 +179,15 @@ private: /** * @brief A class representing an edge in the Delaunay triangulation graph */ -class DB_PUBLIC GPolygonEdge +class DB_PUBLIC Edge { public: - class GPolygonIterator + class PolygonIterator { public: - typedef GPolygon value_type; - typedef GPolygon &reference; - typedef GPolygon *pointer; + typedef Polygon value_type; + typedef Polygon &reference; + typedef Polygon *pointer; reference operator*() const { @@ -131,17 +199,17 @@ public: return m_index ? mp_edge->right () : mp_edge->left (); } - bool operator== (const GPolygonIterator &other) const + bool operator== (const PolygonIterator &other) const { return m_index == other.m_index; } - bool operator!= (const GPolygonIterator &other) const + bool operator!= (const PolygonIterator &other) const { return !operator== (other); } - GPolygonIterator &operator++ () + PolygonIterator &operator++ () { while (++m_index < 2 && operator-> () == 0) ; @@ -149,9 +217,9 @@ public: } private: - friend class GPolygonEdge; + friend class Edge; - GPolygonIterator (const GPolygonEdge *edge) + PolygonIterator (const Edge *edge) : mp_edge (edge), m_index (0) { if (! edge) { @@ -162,48 +230,72 @@ public: } } - const GPolygonEdge *mp_edge; + const Edge *mp_edge; unsigned int m_index; }; - GPolygonEdge (); - GPolygonEdge (GVertex *v1, GVertex *v2); + /** + * @brief Gets the first vertex ("from") + */ + Vertex *v1 () const { return mp_v1; } - GVertex *v1 () const { return mp_v1; } - GVertex *v2 () const { return mp_v2; } + /** + * @brief Gets the first vertex ("to") + */ + Vertex *v2 () const { return mp_v2; } + /** + * @brief Reverses the edge + */ void reverse () { std::swap (mp_v1, mp_v2); std::swap (mp_left, mp_right); } - GPolygon *left () const { return mp_left; } - GPolygon *right () const { return mp_right; } + /** + * @brief Gets the polygon on the left side (can be null) + */ + Polygon *left () const { return mp_left; } - GPolygonIterator begin_polygons () const + /** + * @brief Gets the polygon on the right side (can be null) + */ + Polygon *right () const { return mp_right; } + + /** + * @brief Iterates the polygons (one or two, begin iterator) + */ + PolygonIterator begin_polygons () const { - return GPolygonIterator (this); + return PolygonIterator (this); } - GPolygonIterator end_polygons () const + /** + * @brief Iterates the polygons (end iterator) + */ + PolygonIterator end_polygons () const { - return GPolygonIterator (0); + return PolygonIterator (0); } - void set_level (size_t l) { m_level = l; } - size_t level () const { return m_level; } - - void set_id (size_t id) { m_id = id; } - size_t id () const { return m_id; } - - void set_is_segment (bool is_seg) { m_is_segment = is_seg; } + /** + * @brief Gets a value indicating whether the edge is a segment + */ bool is_segment () const { return m_is_segment; } + /** + * @brief Gets the edge ID (a unique identifier) + */ + size_t id () const { return m_id; } + + /** + * @brief Gets a string representation of the edge + */ std::string to_string (bool with_id = false) const; /** - * @brief Converts to an db::DEdge + * @brief Converts to a db::DEdge */ db::DEdge edge () const { @@ -254,7 +346,7 @@ public: * "crosses" is true, if both edges share at least one point which is not an endpoint * of one of the edges. */ - bool crosses (const db::GPolygonEdge &other) const + bool crosses (const Edge &other) const { return crosses (edge (), other.edge ()); } @@ -279,7 +371,7 @@ public: * @brief Returns a value indicating whether this edge crosses the other one * "crosses" is true, if both edges share at least one point. */ - bool crosses_including (const db::GPolygonEdge &other) const + bool crosses_including (const Edge &other) const { return crosses_including (edge (), other.edge ()); } @@ -301,7 +393,7 @@ public: /** * @brief Gets the intersection point */ - db::DPoint intersection_point (const GPolygonEdge &other) const + db::DPoint intersection_point (const Edge &other) const { return intersection_point (edge (), other.edge ()); } @@ -353,24 +445,23 @@ public: /** * @brief Gets the other triangle for the given one */ - GPolygon *other (const GPolygon *) const; + Polygon *other (const Polygon *) const; /** * @brief Gets the other vertex for the given one */ - GVertex *other (const GVertex *) const; + Vertex *other (const Vertex *) const; /** * @brief Gets a value indicating whether the edge has the given vertex */ - bool has_vertex (const GVertex *) const; + bool has_vertex (const Vertex *) const; /** * @brief Gets the common vertex of the other edge and this edge or null if there is no common vertex */ - GVertex *common_vertex (const GPolygonEdge *other) const; + Vertex *common_vertex (const Edge *other) const; -#if 0 // @@@ /** * @brief Returns a value indicating whether this edge can be flipped */ @@ -379,81 +470,93 @@ public: /** * @brief Returns a value indicating whether the edge separates two triangles that can be joined into one (via the given vertex) */ - bool can_join_via (const GVertex *vertex) const; + bool can_join_via (const Vertex *vertex) const; + + /** + * @brief Returns a value indicating whether this edge belongs to outside triangles + */ + bool is_for_outside_triangles () const; /** * @brief Returns a value indicating whether this edge is an outside edge (no other triangles) */ bool is_outside () const; - /** - * @brief Returns a value indicating whether this edge belongs to outside triangles - */ - bool is_for_outside_triangles () const; -#endif // @@@ - /** * @brief Returns a value indicating whether t is attached to this edge */ - bool has_polygon (const GPolygon *t) const; + bool has_polygon (const Polygon *t) const; protected: void unlink (); void link (); private: - friend class GPolygon; - friend class PolygonGraph; + friend class Polygon; + friend class Graph; + friend class Triangulation; - GVertex *mp_v1, *mp_v2; - GPolygon *mp_left, *mp_right; - GVertex::edges_iterator_non_const m_ec_v1, m_ec_v2; + void set_level (size_t l) { m_level = l; } + size_t level () const { return m_level; } + + void set_id (size_t id) { m_id = id; } + + void set_is_segment (bool is_seg) { m_is_segment = is_seg; } + + Edge (Graph *graph); + Edge (Graph *graph, Vertex *v1, Vertex *v2); + + Graph *mp_graph; + Vertex *mp_v1, *mp_v2; + Polygon *mp_left, *mp_right; + Vertex::edges_iterator_non_const m_ec_v1, m_ec_v2; size_t m_level; size_t m_id; bool m_is_segment; - void set_left (GPolygon *t); - void set_right (GPolygon *t); + void set_left (Polygon *t); + void set_right (Polygon *t); }; /** - * @brief A compare function that compares triangles by ID + * @brief A compare function that compares edges by ID * * The ID acts as a more predicable unique ID for the object in sets and maps. */ -struct GPolygonEdgeLessFunc +struct EdgeLessFunc { - bool operator () (GPolygonEdge *a, GPolygonEdge *b) const + bool operator () (Edge *a, Edge *b) const { return a->id () < b->id (); } }; /** - * @brief A class representing a triangle + * @brief A class representing a polygon */ -class DB_PUBLIC GPolygon - : public tl::list_node, public tl::Object +class DB_PUBLIC Polygon + : public tl::list_node, public tl::Object { public: - GPolygon (); + Polygon (Graph *graph); + Polygon (Graph *graph, Edge *e1, Edge *e2, Edge *e3); template - GPolygon (Iter from, Iter to) - : mp_e (from, to) + Polygon (Graph *graph, Iter from, Iter to) + : mp_graph (graph), mp_e (from, to) { init (); } - ~GPolygon (); + ~Polygon (); void unlink (); void set_id (size_t id) { m_id = id; } size_t id () const { return m_id; } - // @@@bool is_outside () const { return m_is_outside; } - // @@@void set_outside (bool o) { m_is_outside = o; } + bool is_outside () const { return m_is_outside; } + void set_outside (bool o) { m_is_outside = o; } std::string to_string (bool with_id = false) const; @@ -469,7 +572,7 @@ public: * @brief Gets the nth vertex (n wraps around and can be negative) * The vertexes are oriented clockwise. */ - inline GVertex *vertex (int n) const + inline Vertex *vertex (int n) const { size_t sz = size (); tl_assert (sz > 0); @@ -483,7 +586,7 @@ public: /** * @brief Gets the nth edge (n wraps around and can be negative) */ - inline GPolygonEdge *edge (int n) const + inline Edge *edge (int n) const { size_t sz = size (); tl_assert (sz > 0); @@ -504,46 +607,50 @@ public: */ db::DBox bbox () const; -#if 0 // @@@ /** * @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). + * + * This method only applies to triangles. */ std::pair circumcircle (bool *ok = 0) const; /** * @brief Gets the vertex opposite of the given edge + * + * This method only applies to triangles. */ - GVertex *opposite (const GPolygonEdge *edge) const; + Vertex *opposite (const Edge *edge) const; /** * @brief Gets the edge opposite of the given vertex + * + * This method only applies to triangles. */ - GPolygonEdge *opposite (const GVertex *vertex) const; -#endif + Edge *opposite (const Vertex *vertex) const; /** * @brief Gets the edge with the given vertexes */ - GPolygonEdge *find_edge_with (const GVertex *v1, const GVertex *v2) const; + Edge *find_edge_with (const Vertex *v1, const Vertex *v2) const; /** * @brief Finds the common edge for both polygons */ - GPolygonEdge *common_edge (const GPolygon *other) const; + Edge *common_edge (const Polygon *other) const; -#if 0 // @@@ /** * @brief Returns a value indicating whether the point is inside (1), on the polygon (0) or outside (-1) + * + * This method only applies to triangles currently. */ int contains (const db::DPoint &point) const; -#endif /** * @brief Gets a value indicating whether the triangle has the given vertex */ - inline bool has_vertex (const db::GVertex *v) const + inline bool has_vertex (const Vertex *v) const { for (auto i = mp_v.begin (); i != mp_v.end (); ++i) { if (*i == v) { @@ -556,7 +663,7 @@ public: /** * @brief Gets a value indicating whether the triangle has the given edge */ - inline bool has_edge (const db::GPolygonEdge *e) const + inline bool has_edge (const Edge *e) const { for (auto i = mp_e.begin (); i != mp_e.end (); ++i) { if (*i == e) { @@ -571,12 +678,12 @@ public: */ double min_edge_length () const; -#if 0 // @@@ /** * @brief Returns the min edge length to circumcircle radius ratio + * + * This method only applies to triangles currently. */ double b () const; -#endif /** * @brief Returns a value indicating whether the polygon borders to a segment @@ -589,16 +696,17 @@ public: unsigned int num_segments () const; private: - // @@@ bool m_is_outside; - std::vector mp_e; - std::vector mp_v; + Graph *mp_graph; + bool m_is_outside; + std::vector mp_e; + std::vector mp_v; size_t m_id; void init (); // no copying - GPolygon &operator= (const GPolygon &); - GPolygon (const GPolygon &); + Polygon &operator= (const Polygon &); + Polygon (const Polygon &); }; /** @@ -606,87 +714,29 @@ private: * * The ID acts as a more predicable unique ID for the object in sets and maps. */ -struct GPolygonLessFunc +struct PolygonLessFunc { - bool operator () (GPolygon *a, GPolygon *b) const + bool operator () (Polygon *a, Polygon *b) const { return a->id () < b->id (); } }; -class DB_PUBLIC PolygonGraph +/** + * @brief A class representing the polygon graph + * + * A polygon graph is the main container, holding vertexes, edges and polygons. + * The graph can be of "triangles" type, in which case it is guaranteed to only + * hold triangles (polygons with 3 vertexes). + */ +class DB_PUBLIC Graph { public: -#if 0 // @@@ - struct TriangulateParameters - { - TriangulateParameters () - : min_b (1.0), - min_length (0.0), - max_area (0.0), - max_area_border (0.0), - max_iterations (std::numeric_limits::max ()), - base_verbosity (30), - mark_triangles (false) - { } - - /** - * @brief Min. readius-to-shortest edge ratio - */ - double min_b; - - /** - * @brief Min. edge length - * - * This parameter does not provide a guarantee about a minimume edge length, but - * helps avoiding ever-reducing triangle splits in acute corners of the input polygon. - * Splitting of edges stops when the edge is less than the min length. - */ - double min_length; - - /** - * @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; - - /** - * @brief Max number of iterations - */ - size_t max_iterations; - - /** - * @brief The verbosity level above which triangulation reports details - */ - int base_verbosity; - - /** - * @brief If true, final triangles are marked using the "id" integer as a bit field - * - * This provides information about the result quality. - * - * Bit 0: skinny triangle - * Bit 1: bad-quality (skinny or area too large) - * Bit 2: non-Delaunay (in the strict sense) - */ - bool mark_triangles; - }; -#endif - - typedef tl::list polygons_type; + typedef tl::list polygons_type; typedef polygons_type::const_iterator polygon_iterator; - PolygonGraph (); - ~PolygonGraph (); - - /** - * @brief Creates a convex decomposition for the given polygon - */ - void convex_decompose (const DPolygon &poly); + Graph (); + ~Graph (); /** * @brief Returns a string representation of the polygon graph. @@ -718,15 +768,6 @@ public: */ void clear (); -protected: -#if 0 // @@@ - /** - * @brief Checks the polygon graph for consistency - * This method is for testing purposes mainly. - */ - bool check (bool check_delaunay = true) const; -#endif - /** * @brief Dumps the polygon graph to a GDS file at the given path * This method is for testing purposes mainly. @@ -743,35 +784,43 @@ protected: */ db::Layout *to_layout (bool decompose_by_id = false) const; -private: - tl::list mp_polygons; - tl::stable_vector m_edges_heap; - std::vector m_returned_edges; - tl::stable_vector m_vertex_heap; -// @@@ bool m_is_constrained; -// @@@ size_t m_level; - size_t m_id; -// @@@ size_t m_flips, m_hops; - - db::GVertex *create_vertex (double x, double y); - db::GVertex *create_vertex (const db::DPoint &pt); - db::GPolygonEdge *create_edge (db::GVertex *v1, db::GVertex *v2); +protected: + Vertex *create_vertex (double x, double y); + Vertex *create_vertex (const db::DPoint &pt); + Edge *create_edge (Vertex *v1, Vertex *v2); template - db::GPolygon * + Polygon * create_polygon (Iter from, Iter to) { - db::GPolygon *res = new db::GPolygon (from ,to); + Polygon *res = new Polygon (this, from ,to); res->set_id (++m_id); mp_polygons.push_back (res); return res; } - void remove_polygon (db::GPolygon *tri); - template void make_contours (const Poly &poly, const Trans &trans, std::vector > &contours); + Polygon *create_triangle (Edge *e1, Edge *e2, Edge *e3); + + void remove_polygon (Polygon *p); + +private: + friend class Triangulation; + friend class ConvexDecomposition; + + tl::list mp_polygons; + tl::stable_vector m_edges_heap; + std::vector m_returned_edges; + tl::stable_vector m_vertex_heap; + size_t m_id; + + tl::list &polygons () { return mp_polygons; } + tl::stable_vector &edges () { return m_edges_heap; } + tl::stable_vector &vertexes () { return m_vertex_heap; } }; -} +} // namespace plc + +} // namespace db #endif diff --git a/src/db/db/dbPLCTriangulation.cc b/src/db/db/dbPLCTriangulation.cc new file mode 100644 index 000000000..879cd9e7d --- /dev/null +++ b/src/db/db/dbPLCTriangulation.cc @@ -0,0 +1,1649 @@ + +/* + + KLayout Layout Viewer + Copyright (C) 2006-2025 Matthias Koefferlein + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +*/ + + +#include "dbPLCTriangulation.h" +#include "dbLayout.h" +#include "dbWriter.h" +#include "tlStream.h" +#include "tlLog.h" +#include "tlTimer.h" + +#include +#include +#include +#include + +namespace db +{ + +namespace plc +{ + +static inline bool is_equal (const db::DPoint &a, const db::DPoint &b) +{ + return std::abs (a.x () - b.x ()) < std::max (1.0, (std::abs (a.x ()) + std::abs (b.x ()))) * db::epsilon && + std::abs (a.y () - b.y ()) < std::max (1.0, (std::abs (a.y ()) + std::abs (b.y ()))) * db::epsilon; +} + +Triangulation::Triangulation (Graph *graph) +{ + mp_graph = graph; + clear (); +} + +void +Triangulation::clear () +{ + mp_graph->clear (); + + m_is_constrained = false; + m_level = 0; + m_id = 0; + m_flips = m_hops = 0; +} + +void +Triangulation::init_box (const db::DBox &box) +{ + double xmin = box.left (), xmax = box.right (); + double ymin = box.bottom (), ymax = box.top (); + + Vertex *vbl = mp_graph->create_vertex (xmin, ymin); + Vertex *vtl = mp_graph->create_vertex (xmin, ymax); + Vertex *vbr = mp_graph->create_vertex (xmax, ymin); + Vertex *vtr = mp_graph->create_vertex (xmax, ymax); + + Edge *sl = mp_graph->create_edge (vbl, vtl); + Edge *sd = mp_graph->create_edge (vtl, vbr); + Edge *sb = mp_graph->create_edge (vbr, vbl); + + Edge *sr = mp_graph->create_edge (vbr, vtr); + Edge *st = mp_graph->create_edge (vtr, vtl); + + mp_graph->create_triangle (sl, sd, sb); + mp_graph->create_triangle (sd, sr, st); +} + +bool +Triangulation::check (bool check_delaunay) const +{ + bool res = true; + + for (auto t = mp_graph->polygons ().begin (); t != mp_graph->polygons ().end (); ++t) { + if (t->size () != 3) { + res = false; + tl::error << "(check error) not a triangle: " << t->to_string (); + } + } + + if (! res) { + return false; + } + + if (check_delaunay) { + for (auto t = mp_graph->polygons ().begin (); t != mp_graph->polygons ().end (); ++t) { + auto cp = t->circumcircle (); + auto vi = find_inside_circle (cp.first, cp.second); + if (! vi.empty ()) { + res = false; + tl::error << "(check error) triangle does not meet Delaunay criterion: " << t->to_string (); + for (auto v = vi.begin (); v != vi.end (); ++v) { + tl::error << " vertex inside circumcircle: " << (*v)->to_string (true); + } + } + } + } + + for (auto t = mp_graph->polygons ().begin (); t != mp_graph->polygons ().end (); ++t) { + for (int i = 0; i < 3; ++i) { + if (! t->edge (i)->has_polygon (t.operator-> ())) { + tl::error << "(check error) edges " << t->edge (i)->to_string (true) + << " attached to triangle " << t->to_string (true) << " does not refer to this triangle"; + res = false; + } + } + } + + for (auto e = mp_graph->edges ().begin (); e != mp_graph->edges ().end (); ++e) { + + if (!e->left () && !e->right ()) { + continue; + } + + if (e->left () && e->right ()) { + if (e->left ()->is_outside () != e->right ()->is_outside () && ! e->is_segment ()) { + tl::error << "(check error) edge " << e->to_string (true) << " splits an outside and inside triangle, but is not a segment"; + res = false; + } + } + + for (auto t = e->begin_polygons (); t != e->end_polygons (); ++t) { + if (! t->has_edge (e.operator-> ())) { + tl::error << "(check error) edge " << e->to_string (true) << " not found in adjacent triangle " << t->to_string (true); + res = false; + } + if (! t->has_vertex (e->v1 ())) { + tl::error << "(check error) edges " << e->to_string (true) << " vertex 1 not found in adjacent triangle " << t->to_string (true); + res = false; + } + if (! t->has_vertex (e->v2 ())) { + tl::error << "(check error) edges " << e->to_string (true) << " vertex 2 not found in adjacent triangle " << t->to_string (true); + res = false; + } + Vertex *vopp = t->opposite (e.operator-> ()); + double sgn = (e->left () == t.operator-> ()) ? 1.0 : -1.0; + double vp = db::vprod (e->d(), *vopp - *e->v1 ()); // positive if on left side + if (vp * sgn <= 0.0) { + const char * side_str = sgn > 0.0 ? "left" : "right"; + tl::error << "(check error) external point " << vopp->to_string (true) << " not on " << side_str << " side of edge " << e->to_string (true); + res = false; + } + } + + if (! e->v1 ()->has_edge (e.operator-> ())) { + tl::error << "(check error) edge " << e->to_string (true) << " vertex 1 does not list this edge"; + res = false; + } + if (! e->v2 ()->has_edge (e.operator-> ())) { + tl::error << "(check error) edge " << e->to_string (true) << " vertex 2 does not list this edge"; + res = false; + } + + } + + for (auto v = mp_graph->vertexes ().begin (); v != mp_graph->vertexes ().end (); ++v) { + unsigned int num_outside_edges = 0; + for (auto e = v->begin_edges (); e != v->end_edges (); ++e) { + if ((*e)->is_outside ()) { + ++num_outside_edges; + } + } + if (num_outside_edges > 0 && num_outside_edges != 2) { + tl::error << "(check error) vertex " << v->to_string (true) << " has " << num_outside_edges << " outside edges (can only be 2)"; + res = false; + for (auto e = v->begin_edges (); e != v->end_edges (); ++e) { + if ((*e)->is_outside ()) { + tl::error << " Outside edge is " << (*e)->to_string (true); + } + } + } + } + + return res; +} + +std::vector +Triangulation::find_points_around (Vertex *vertex, double radius) +{ + std::set seen; + seen.insert (vertex); + + std::vector res; + std::vector new_vertexes, next_vertexes; + new_vertexes.push_back (vertex); + + while (! new_vertexes.empty ()) { + next_vertexes.clear (); + for (auto v = new_vertexes.begin (); v != new_vertexes.end (); ++v) { + for (auto e = (*v)->begin_edges (); e != (*v)->end_edges (); ++e) { + Vertex *ov = (*e)->other (*v); + if (ov->in_circle (*vertex, radius) == 1 && seen.insert (ov).second) { + next_vertexes.push_back (ov); + res.push_back (ov); + } + } + } + new_vertexes.swap (next_vertexes); + } + + return res; +} + +Vertex * +Triangulation::insert_point (const db::DPoint &point, std::list > *new_triangles) +{ + return insert (mp_graph->create_vertex (point), new_triangles); +} + +Vertex * +Triangulation::insert_point (db::DCoord x, db::DCoord y, std::list > *new_triangles) +{ + return insert (mp_graph->create_vertex (x, y), new_triangles); +} + +Vertex * +Triangulation::insert (Vertex *vertex, std::list > *new_triangles) +{ + std::vector tris = find_triangle_for_point (*vertex); + + // the new vertex is outside the domain + if (tris.empty ()) { + tl_assert (! m_is_constrained); + insert_new_vertex (vertex, new_triangles); + return vertex; + } + + // check, if the new vertex is on an edge (may be edge between triangles or edge on outside) + std::vector on_edges; + std::vector on_vertex; + for (int i = 0; i < 3; ++i) { + Edge *e = tris.front ()->edge (i); + if (e->side_of (*vertex) == 0) { + if (is_equal (*vertex, *e->v1 ()) || is_equal (*vertex, *e->v2 ())) { + on_vertex.push_back (e); + } else { + on_edges.push_back (e); + } + } + } + + if (! on_vertex.empty ()) { + + tl_assert (on_vertex.size () == size_t (2)); + return on_vertex.front ()->common_vertex (on_vertex [1]); + + } else if (! on_edges.empty ()) { + + tl_assert (on_edges.size () == size_t (1)); + split_triangles_on_edge (vertex, on_edges.front (), new_triangles); + return vertex; + + } else if (tris.size () == size_t (1)) { + + // the new vertex is inside one triangle + split_triangle (tris.front (), vertex, new_triangles); + return vertex; + + } + + tl_assert (false); +} + +std::vector +Triangulation::find_triangle_for_point (const db::DPoint &point) +{ + Edge *edge = find_closest_edge (point); + + std::vector res; + if (edge) { + for (auto t = edge->begin_polygons (); t != edge->end_polygons (); ++t) { + if (t->contains (point) >= 0) { + res.push_back (t.operator-> ()); + } + } + } + + return res; +} + +Edge * +Triangulation::find_closest_edge (const db::DPoint &p, Vertex *vstart, bool inside_only) +{ + if (!vstart) { + + if (! mp_graph->polygons ().empty ()) { + + unsigned int ls = 0; + size_t n = mp_graph->vertexes ().size (); + size_t m = n; + + // A simple heuristics that takes a sqrt(N) sample from the + // vertexes to find a good starting point + + vstart = mp_graph->polygons ().begin ()->vertex (0); + double dmin = vstart->distance (p); + + while (ls * ls < m) { + m /= 2; + for (size_t i = m / 2; i < n; i += m) { + ++ls; + // NOTE: this assumes the heap is not too loaded with orphan vertexes + Vertex *v = (mp_graph->vertexes ().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; + + } + + } + + db::DEdge line (*vstart, p); + + double d = -1.0; + Edge *edge = 0; + Vertex *v = vstart; + + while (v) { + + Vertex *vnext = 0; + + for (auto e = v->begin_edges (); e != v->end_edges (); ++e) { + + if (inside_only) { + // NOTE: in inside mode we stay on the line of sight as we don't + // want to walk around outside pockets. + if (! (*e)->is_segment () && (*e)->is_for_outside_triangles ()) { + continue; + } + if (! (*e)->crosses_including (line)) { + continue; + } + } + + double ds = (*e)->distance (p); + + if (d < 0.0) { + + d = ds; + edge = *e; + vnext = edge->other (v); + + } else if (fabs (ds - d) < std::max (1.0, fabs (ds) + fabs (d)) * db::epsilon) { + + // this differentiation selects the edge which bends further towards + // the target point if both edges share a common point and that + // is the one the determines the distance. + Vertex *cv = edge->common_vertex (*e); + if (cv) { + db::DVector edge_d = *edge->other (cv) - *cv; + db::DVector e_d = *(*e)->other(cv) - *cv; + db::DVector r = p - *cv; + double edge_sp = db::sprod (r, edge_d) / edge_d.length (); + double s_sp = db::sprod (r, e_d) / e_d.length (); + if (s_sp > edge_sp + db::epsilon) { + edge = *e; + vnext = edge->other (v); + } + } + + } else if (ds < d) { + + d = ds; + edge = *e; + vnext = edge->other (v); + + } + + } + + ++m_hops; + + v = vnext; + + } + + return edge; +} + +void +Triangulation::insert_new_vertex (Vertex *vertex, std::list > *new_triangles_out) +{ + if (mp_graph->polygons ().empty ()) { + + tl_assert (mp_graph->vertexes ().size () <= size_t (3)); // fails if vertexes were created but not inserted. + + if (mp_graph->vertexes ().size () == 3) { + + std::vector vv; + for (auto v = mp_graph->vertexes ().begin (); v != mp_graph->vertexes ().end (); ++v) { + vv.push_back (v.operator-> ()); + } + + // form the first triangle + Edge *s1 = mp_graph->create_edge (vv[0], vv[1]); + Edge *s2 = mp_graph->create_edge (vv[1], vv[2]); + Edge *s3 = mp_graph->create_edge (vv[2], vv[0]); + + if (db::vprod_sign (s1->d (), s2->d ()) == 0) { + // avoid degenerate Triangles to happen here + tl_assert (false); + } else { + Polygon *t = mp_graph->create_triangle (s1, s2, s3); + if (new_triangles_out) { + new_triangles_out->push_back (t); + } + } + + } + + return; + + } + + std::vector new_triangles; + + // Find closest edge + Edge *closest_edge = find_closest_edge (*vertex); + tl_assert (closest_edge != 0); + + Edge *s1 = mp_graph->create_edge (vertex, closest_edge->v1 ()); + Edge *s2 = mp_graph->create_edge (vertex, closest_edge->v2 ()); + + Polygon *t = mp_graph->create_triangle (s1, closest_edge, s2); + new_triangles.push_back (t); + + add_more_triangles (new_triangles, closest_edge, closest_edge->v1 (), vertex, s1); + add_more_triangles (new_triangles, closest_edge, closest_edge->v2 (), vertex, s2); + + if (new_triangles_out) { + new_triangles_out->insert (new_triangles_out->end (), new_triangles.begin (), new_triangles.end ()); + } + + fix_triangles (new_triangles, std::vector (), new_triangles_out); +} + +void +Triangulation::add_more_triangles (std::vector &new_triangles, + Edge *incoming_edge, + Vertex *from_vertex, Vertex *to_vertex, + Edge *conn_edge) +{ + while (true) { + + Edge *next_edge = 0; + + for (auto e = from_vertex->begin_edges (); e != from_vertex->end_edges (); ++e) { + if (! (*e)->has_vertex (to_vertex) && (*e)->is_outside ()) { + // TODO: remove and break + tl_assert (next_edge == 0); + next_edge = *e; + } + } + + tl_assert (next_edge != 0); + Vertex *next_vertex = next_edge->other (from_vertex); + + db::DVector d_from_to = *to_vertex - *from_vertex; + Vertex *incoming_vertex = incoming_edge->other (from_vertex); + if (db::vprod_sign(*from_vertex - *incoming_vertex, d_from_to) * db::vprod_sign(*from_vertex - *next_vertex, d_from_to) >= 0) { + return; + } + + Edge *next_conn_edge = mp_graph->create_edge (next_vertex, to_vertex); + Polygon *t = mp_graph->create_triangle (next_conn_edge, next_edge, conn_edge); + new_triangles.push_back (t); + + incoming_edge = next_edge; + conn_edge = next_conn_edge; + from_vertex = next_vertex; + + } +} + +void +Triangulation::split_triangle (Polygon *t, Vertex *vertex, std::list > *new_triangles_out) +{ + t->unlink (); + + std::map v2new_edges; + std::vector new_edges; + for (int i = 0; i < 3; ++i) { + Vertex *v = t->vertex (i); + Edge *e = mp_graph->create_edge (v, vertex); + v2new_edges[v] = e; + new_edges.push_back (e); + } + + std::vector new_triangles; + for (int i = 0; i < 3; ++i) { + Edge *e = t->edge (i); + Polygon *new_triangle = mp_graph->create_triangle (e, v2new_edges[e->v1 ()], v2new_edges[e->v2 ()]); + if (new_triangles_out) { + new_triangles_out->push_back (new_triangle); + } + new_triangle->set_outside (t->is_outside ()); + new_triangles.push_back (new_triangle); + } + + mp_graph->remove_polygon (t); + + fix_triangles (new_triangles, new_edges, new_triangles_out); +} + +void +Triangulation::split_triangles_on_edge (Vertex *vertex, Edge *split_edge, std::list > *new_triangles_out) +{ + Edge *s1 = mp_graph->create_edge (split_edge->v1 (), vertex); + Edge *s2 = mp_graph->create_edge (split_edge->v2 (), vertex); + s1->set_is_segment (split_edge->is_segment ()); + s2->set_is_segment (split_edge->is_segment ()); + + std::vector new_triangles; + + std::vector tris; + tris.reserve (2); + for (auto t = split_edge->begin_polygons (); t != split_edge->end_polygons (); ++t) { + tris.push_back (t.operator-> ()); + } + + for (auto t = tris.begin (); t != tris.end (); ++t) { + + (*t)->unlink (); + + Vertex *ext_vertex = (*t)->opposite (split_edge); + Edge *new_edge = mp_graph->create_edge (ext_vertex, vertex); + + for (int i = 0; i < 3; ++i) { + + Edge *e = (*t)->edge (i); + if (e->has_vertex (ext_vertex)) { + + Edge *partial = e->has_vertex (split_edge->v1 ()) ? s1 : s2; + Polygon *new_triangle = mp_graph->create_triangle (new_edge, partial, e); + + if (new_triangles_out) { + new_triangles_out->push_back (new_triangle); + } + new_triangle->set_outside ((*t)->is_outside ()); + new_triangles.push_back (new_triangle); + + } + + } + + } + + for (auto t = tris.begin (); t != tris.end (); ++t) { + mp_graph->remove_polygon (*t); + } + + std::vector fixed_edges; + fixed_edges.push_back (s1); + fixed_edges.push_back (s2); + fix_triangles (new_triangles, fixed_edges, new_triangles_out); +} + +std::vector +Triangulation::find_touching (const db::DBox &box) const +{ + // NOTE: this is a naive, slow implementation for test purposes + std::vector res; + for (auto v = mp_graph->vertexes ().begin (); v != mp_graph->vertexes ().end (); ++v) { + if (v->begin_edges () != v->end_edges ()) { + if (box.contains (*v)) { + res.push_back (const_cast (v.operator-> ())); + } + } + } + return res; +} + +std::vector +Triangulation::find_inside_circle (const db::DPoint ¢er, double radius) const +{ + // NOTE: this is a naive, slow implementation for test purposes + std::vector res; + for (auto v = mp_graph->vertexes ().begin (); v != mp_graph->vertexes ().end (); ++v) { + if (v->begin_edges () != v->end_edges ()) { + if (v->in_circle (center, radius) == 1) { + res.push_back (const_cast (v.operator-> ())); + } + } + } + return res; +} + +void +Triangulation::remove (Vertex *vertex, std::list > *new_triangles) +{ + if (vertex->begin_edges () == vertex->end_edges ()) { + // removing an orphan vertex -> ignore + } else if (vertex->is_outside ()) { + remove_outside_vertex (vertex, new_triangles); + } else { + remove_inside_vertex (vertex, new_triangles); + } +} + +void +Triangulation::remove_outside_vertex (Vertex *vertex, std::list > *new_triangles_out) +{ + auto to_remove = vertex->polygons (); + + std::vector outer_edges; + for (auto t = to_remove.begin (); t != to_remove.end (); ++t) { + outer_edges.push_back ((*t)->opposite (vertex)); + } + + for (auto t = to_remove.begin (); t != to_remove.end (); ++t) { + (*t)->unlink (); + } + + auto new_triangles = fill_concave_corners (outer_edges); + + for (auto t = to_remove.begin (); t != to_remove.end (); ++t) { + mp_graph->remove_polygon (*t); + } + + fix_triangles (new_triangles, std::vector (), new_triangles_out); +} + +void +Triangulation::remove_inside_vertex (Vertex *vertex, std::list > *new_triangles_out) +{ + std::set triangles_to_fix; + + bool make_new_triangle = true; + + while (vertex->num_edges (4) > 3) { + + Edge *to_flip = 0; + for (auto e = vertex->begin_edges (); e != vertex->end_edges () && to_flip == 0; ++e) { + if ((*e)->can_flip ()) { + to_flip = *e; + } + } + if (! to_flip) { + break; + } + + // NOTE: in the "can_join" case zero-area triangles are created which we will sort out later + triangles_to_fix.erase (to_flip->left ()); + triangles_to_fix.erase (to_flip->right ()); + + auto pp = flip (to_flip); + triangles_to_fix.insert (pp.first.first); + triangles_to_fix.insert (pp.first.second); + + } + + if (vertex->num_edges (4) > 3) { + + tl_assert (vertex->num_edges (5) == 4); + + // This case can happen if two edges attached to the vertex are collinear + // in this case choose the "join" strategy + Edge *jseg = 0; + for (auto e = vertex->begin_edges (); e != vertex->end_edges () && !jseg; ++e) { + if ((*e)->can_join_via (vertex)) { + jseg = *e; + } + } + tl_assert (jseg != 0); + + Vertex *v1 = jseg->left ()->opposite (jseg); + Edge *s1 = jseg->left ()->opposite (vertex); + Vertex *v2 = jseg->right ()->opposite (jseg); + Edge *s2 = jseg->right ()->opposite (vertex); + + Edge *jseg_opp = 0; + for (auto e = vertex->begin_edges (); e != vertex->end_edges () && !jseg_opp; ++e) { + if (!(*e)->has_polygon (jseg->left ()) && !(*e)->has_polygon (jseg->right ())) { + jseg_opp = *e; + } + } + + Edge *s1opp = jseg_opp->left ()->opposite (vertex); + Edge *s2opp = jseg_opp->right ()->opposite (vertex); + + Edge *new_edge = mp_graph->create_edge (v1, v2); + Polygon *t1 = mp_graph->create_triangle (s1, s2, new_edge); + Polygon *t2 = mp_graph->create_triangle (s1opp, s2opp, new_edge); + + triangles_to_fix.insert (t1); + triangles_to_fix.insert (t2); + + make_new_triangle = false; + + } + + auto to_remove = vertex->polygons (); + + std::vector outer_edges; + for (auto t = to_remove.begin (); t != to_remove.end (); ++t) { + outer_edges.push_back ((*t)->opposite (vertex)); + } + + if (make_new_triangle) { + + tl_assert (outer_edges.size () == size_t (3)); + + Polygon *nt = mp_graph->create_triangle (outer_edges[0], outer_edges[1], outer_edges[2]); + triangles_to_fix.insert (nt); + + } + + for (auto t = to_remove.begin (); t != to_remove.end (); ++t) { + triangles_to_fix.erase (*t); + mp_graph->remove_polygon (*t); + } + + if (new_triangles_out) { + for (auto t = triangles_to_fix.begin (); t != triangles_to_fix.end (); ++t) { + new_triangles_out->push_back (*t); + } + } + + std::vector to_fix_a (triangles_to_fix.begin (), triangles_to_fix.end ()); + fix_triangles (to_fix_a, std::vector (), new_triangles_out); +} + +void +Triangulation::fix_triangles (const std::vector &tris, const std::vector &fixed_edges, std::list > *new_triangles) +{ + m_level += 1; + for (auto e = fixed_edges.begin (); e != fixed_edges.end (); ++e) { + (*e)->set_level (m_level); + } + + std::set queue, todo; + + for (auto t = tris.begin (); t != tris.end (); ++t) { + for (int i = 0; i < 3; ++i) { + Edge *e = (*t)->edge (i); + if (e->level () < m_level && ! e->is_segment ()) { + queue.insert (e); + } + } + } + + while (! queue.empty ()) { + + todo.clear (); + todo.swap (queue); + + // NOTE: we cannot be sure that already treated edges will not become + // illegal by neighbor edges flipping .. + // for s in todo: + // s.level = self.level + + for (auto e = todo.begin (); e != todo.end (); ++e) { + + if (is_illegal_edge (*e)) { + + queue.erase (*e); + + auto pp = flip (*e); + auto t1 = pp.first.first; + auto t2 = pp.first.second; + auto s12 = pp.second; + + if (new_triangles) { + new_triangles->push_back (t1); + new_triangles->push_back (t2); + } + + ++m_flips; + tl_assert (! is_illegal_edge (s12)); // TODO: remove later! + + for (int i = 0; i < 3; ++i) { + Edge *s1 = t1->edge (i); + if (s1->level () < m_level && ! s1->is_segment ()) { + queue.insert (s1); + } + } + + for (int i = 0; i < 3; ++i) { + Edge *s2 = t2->edge (i); + if (s2->level () < m_level && ! s2->is_segment ()) { + queue.insert (s2); + } + } + + } + + } + + } +} + +bool +Triangulation::is_illegal_edge (Edge *edge) +{ + Polygon *left = edge->left (); + Polygon *right = edge->right (); + if (!left || !right) { + return false; + } + + 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(&ok); + if (! ok || left->opposite (edge)->in_circle (rr.first, rr.second) > 0) { + return true; + } + + return false; +} + +std::pair, Edge *> +Triangulation::flip (Edge *edge) +{ + Polygon *t1 = edge->left (); + Polygon *t2 = edge->right (); + + bool outside = t1->is_outside (); + tl_assert (t1->is_outside () == outside); + + // prepare for the new triangle to replace this one + t1->unlink (); + t2->unlink (); + + Vertex *t1_vext = t1->opposite (edge); + Edge *t1_sext1 = t1->find_edge_with (t1_vext, edge->v1 ()); + Edge *t1_sext2 = t1->find_edge_with (t1_vext, edge->v2 ()); + + Vertex *t2_vext = t2->opposite (edge); + Edge *t2_sext1 = t2->find_edge_with (t2_vext, edge->v1 ()); + Edge *t2_sext2 = t2->find_edge_with (t2_vext, edge->v2 ()); + + Edge *s_new = mp_graph->create_edge (t1_vext, t2_vext); + + Polygon *t1_new = mp_graph->create_triangle (s_new, t1_sext1, t2_sext1); + t1_new->set_outside (outside); + Polygon *t2_new = mp_graph->create_triangle (s_new, t1_sext2, t2_sext2); + t2_new->set_outside (outside); + + mp_graph->remove_polygon (t1); + mp_graph->remove_polygon (t2); + + return std::make_pair (std::make_pair (t1_new, t2_new), s_new); +} + +std::vector +Triangulation::fill_concave_corners (const std::vector &edges) +{ + std::vector res; + std::vector points, terminals; + + std::map > vertex2edge; + for (auto e = edges.begin (); e != edges.end (); ++e) { + + auto i = vertex2edge.insert (std::make_pair ((*e)->v1 (), std::vector ())); + if (i.second) { + points.push_back ((*e)->v1 ()); + } + i.first->second.push_back (*e); + + i = vertex2edge.insert (std::make_pair ((*e)->v2 (), std::vector ())); + if (i.second) { + points.push_back ((*e)->v2 ()); + } + i.first->second.push_back (*e); + + } + + while (points.size () > size_t (2)) { + + terminals.clear (); + for (auto p = points.begin (); p != points.end (); ++p) { + if (vertex2edge [*p].size () == 1) { + terminals.push_back (*p); + } + } + tl_assert (terminals.size () == size_t (2)); + Vertex *v = terminals[0]; + + bool any_connected = false; + Vertex *vp = 0; + + std::set to_remove; + + while (vertex2edge [v].size () >= size_t (2) || ! vp) { + + Edge *seg = 0; + std::vector &ee = vertex2edge [v]; + for (auto e = ee.begin (); e != ee.end (); ++e) { + if (! (*e)->has_vertex (vp)) { + seg = (*e); + break; + } + } + + tl_assert (seg != 0); + Polygon *tri = seg->left () ? seg->left () : seg->right (); + Vertex *vn = seg->other (v); + + std::vector &een = vertex2edge [vn]; + if (een.size () < size_t (2)) { + break; + } + tl_assert (een.size () == size_t (2)); + + Edge *segn = 0; + for (auto e = een.begin (); e != een.end (); ++e) { + if (! (*e)->has_vertex (v)) { + segn = (*e); + break; + } + } + + tl_assert (segn != 0); + Vertex *vnn = segn->other (vn); + std::vector &eenn = vertex2edge [vnn]; + + // NOTE: tri can be None in case a lonely edge stays after removing + // attached triangles + if (! tri || seg->side_of (*vnn) * seg->side_of (*tri->opposite (seg)) < 0) { + + // is concave + Edge *new_edge = mp_graph->create_edge (v, vnn); + for (auto e = ee.begin (); e != ee.end (); ++e) { + if (*e == seg) { + ee.erase (e); + break; + } + } + ee.push_back (new_edge); + + for (auto e = eenn.begin (); e != eenn.end (); ++e) { + if (*e == segn) { + eenn.erase (e); + break; + } + } + eenn.push_back (new_edge); + + vertex2edge.erase (vn); + to_remove.insert (vn); + + Polygon *new_triangle = mp_graph->create_triangle (seg, segn, new_edge); + res.push_back (new_triangle); + any_connected = true; + + } else { + + vp = v; + v = vn; + + } + + } + + if (! any_connected) { + break; + } + + std::vector::iterator wp = points.begin (); + for (auto p = points.begin (); p != points.end (); ++p) { + if (to_remove.find (*p) == to_remove.end ()) { + *wp++ = *p; + } + } + points.erase (wp, points.end ()); + + } + + return res; +} + +std::vector +Triangulation::search_edges_crossing (Vertex *from, Vertex *to) +{ + Vertex *v = from; + Vertex *vv = to; + db::DEdge edge (*from, *to); + + Polygon *current_triangle = 0; + Edge *next_edge = 0; + + std::vector result; + + for (auto e = v->begin_edges (); e != v->end_edges () && ! next_edge; ++e) { + for (auto t = (*e)->begin_polygons (); t != (*e)->end_polygons (); ++t) { + Edge *os = t->opposite (v); + if (os->has_vertex (vv)) { + return result; + } + if (os->crosses (edge)) { + result.push_back (os); + current_triangle = t.operator-> (); + next_edge = os; + break; + } + } + } + + tl_assert (current_triangle != 0); + + while (true) { + + current_triangle = next_edge->other (current_triangle); + + // Note that we're convex, so there has to be a path across triangles + tl_assert (current_triangle != 0); + + Edge *cs = next_edge; + next_edge = 0; + for (int i = 0; i < 3; ++i) { + Edge *e = current_triangle->edge (i); + if (e != cs) { + if (e->has_vertex (vv)) { + return result; + } + if (e->crosses (edge)) { + result.push_back (e); + next_edge = e; + break; + } + } + } + + tl_assert (next_edge != 0); + + } +} + +Vertex * +Triangulation::find_vertex_for_point (const db::DPoint &point) +{ + Edge *edge = find_closest_edge (point); + if (!edge) { + return 0; + } + Vertex *v = 0; + if (is_equal (*edge->v1 (), point)) { + v = edge->v1 (); + } else if (is_equal (*edge->v2 (), point)) { + v = edge->v2 (); + } + return v; +} + +Edge * +Triangulation::find_edge_for_points (const db::DPoint &p1, const db::DPoint &p2) +{ + Vertex *v = find_vertex_for_point (p1); + if (!v) { + return 0; + } + for (auto e = v->begin_edges (); e != v->end_edges (); ++e) { + if (is_equal (*(*e)->other (v), p2)) { + return *e; + } + } + return 0; +} + +std::vector +Triangulation::ensure_edge_inner (Vertex *from, Vertex *to) +{ + auto crossed_edges = search_edges_crossing (from, to); + std::vector result; + + if (crossed_edges.empty ()) { + + // no crossing edge - there should be a edge already + Edge *res = find_edge_for_points (*from, *to); + tl_assert (res != 0); + result.push_back (res); + + } else if (crossed_edges.size () == 1) { + + // can be solved by flipping + auto pp = flip (crossed_edges.front ()); + Edge *res = pp.second; + tl_assert (res->has_vertex (from) && res->has_vertex (to)); + result.push_back (res); + + } else { + + // split edge close to center + db::DPoint split_point; + double d = -1.0; + double l_half = 0.25 * (*to - *from).sq_length (); + for (auto e = crossed_edges.begin (); e != crossed_edges.end (); ++e) { + db::DPoint p = (*e)->intersection_point (db::DEdge (*from, *to)); + double dp = fabs ((p - *from).sq_length () - l_half); + if (d < 0.0 || dp < d) { + dp = d; + split_point = p; + } + } + + Vertex *split_vertex = insert_point (split_point); + + result = ensure_edge_inner (from, split_vertex); + + auto result2 = ensure_edge_inner (split_vertex, to); + result.insert (result.end (), result2.begin (), result2.end ()); + + } + + return result; +} + +std::vector +Triangulation::ensure_edge (Vertex *from, Vertex *to) +{ +#if 0 + // NOTE: this should not be required if the original segments are non-overlapping + // TODO: this is inefficient + for v in self.vertexes: + if edge.point_on(v): + return self.ensure_edge(Edge(edge.p1, v)) + self.ensure_edge(Edge(v, edge.p2)) +#endif + + auto edges = ensure_edge_inner (from, to); + for (auto e = edges.begin (); e != edges.end (); ++e) { + // mark the edges as fixed "forever" so we don't modify them when we ensure other edges + (*e)->set_level (std::numeric_limits::max ()); + } + return edges; +} + +void +Triangulation::join_edges (std::vector &edges) +{ + // edges are supposed to be ordered + for (size_t i = 1; i < edges.size (); ) { + + Edge *s1 = edges [i - 1]; + Edge *s2 = edges [i]; + tl_assert (s1->is_segment () == s2->is_segment ()); + Vertex *cp = s1->common_vertex (s2); + tl_assert (cp != 0); + + std::vector join_edges; + for (auto e = cp->begin_edges (); e != cp->end_edges (); ++e) { + if (*e != s1 && *e != s2) { + if ((*e)->can_join_via (cp)) { + join_edges.push_back (*e); + } else { + join_edges.clear (); + break; + } + } + } + + if (! join_edges.empty ()) { + + tl_assert (join_edges.size () <= 2); + + Edge *new_edge = mp_graph->create_edge (s1->other (cp), s2->other (cp)); + new_edge->set_is_segment (s1->is_segment ()); + + for (auto js = join_edges.begin (); js != join_edges.end (); ++js) { + + Polygon *t1 = (*js)->left (); + Polygon *t2 = (*js)->right (); + Edge *tedge1 = t1->opposite (cp); + Edge *tedge2 = t2->opposite (cp); + t1->unlink (); + t2->unlink (); + Polygon *tri = mp_graph->create_triangle (tedge1, tedge2, new_edge); + tri->set_outside (t1->is_outside ()); + mp_graph->remove_polygon (t1); + mp_graph->remove_polygon (t2); + } + + edges [i - 1] = new_edge; + edges.erase (edges.begin () + i); + + } else { + ++i; + } + + } +} + +void +Triangulation::constrain (const std::vector > &contours) +{ + tl_assert (! m_is_constrained); + + std::vector > > resolved_edges; + + for (auto c = contours.begin (); c != contours.end (); ++c) { + for (auto v = c->begin (); v != c->end (); ++v) { + auto vv = v; + ++vv; + if (vv == c->end ()) { + vv = c->begin (); + } + db::DEdge e (**v, **vv); + resolved_edges.push_back (std::make_pair (e, std::vector ())); + resolved_edges.back ().second = ensure_edge (*v, *vv); + } + } + + for (auto tri = mp_graph->polygons ().begin (); tri != mp_graph->polygons ().end (); ++tri) { + tri->set_outside (false); + for (int i = 0; i < 3; ++i) { + tri->edge (i)->set_is_segment (false); + } + } + + std::set new_tri; + + for (auto re = resolved_edges.begin (); re != resolved_edges.end (); ++re) { + auto edge = re->first; + auto edges = re->second; + for (auto e = edges.begin (); e != edges.end (); ++e) { + (*e)->set_is_segment (true); + Polygon *outer_tri = 0; + int d = db::sprod_sign (edge.d (), (*e)->d ()); + if (d > 0) { + outer_tri = (*e)->left (); + } + if (d < 0) { + outer_tri = (*e)->right (); + } + if (outer_tri) { + new_tri.insert (outer_tri); + outer_tri->set_outside (true); + } + } + } + + while (! new_tri.empty ()) { + + std::set next_tris; + + for (auto tri = new_tri.begin (); tri != new_tri.end (); ++tri) { + for (int i = 0; i < 3; ++i) { + auto e = (*tri)->edge (i); + if (! e->is_segment ()) { + auto ot = e->other (*tri); + if (ot && ! ot->is_outside ()) { + next_tris.insert (ot); + ot->set_outside (true); + } + } + } + } + + new_tri.swap (next_tris); + + } + + // join edges where possible + for (auto re = resolved_edges.begin (); re != resolved_edges.end (); ++re) { + auto edges = re->second; + join_edges (edges); + } + + m_is_constrained = true; +} + +void +Triangulation::remove_outside_triangles () +{ + tl_assert (m_is_constrained); + + // NOTE: don't remove while iterating + std::vector to_remove; + for (auto tri = mp_graph->begin (); tri != mp_graph->end (); ++tri) { + if (tri->is_outside ()) { + to_remove.push_back (const_cast (tri.operator-> ())); + } + } + + for (auto t = to_remove.begin (); t != to_remove.end (); ++t) { + mp_graph->remove_polygon (*t); + } +} + + +template +void +Triangulation::make_contours (const Poly &poly, const Trans &trans, std::vector > &edge_contours) +{ + edge_contours.push_back (std::vector ()); + for (auto pt = poly.begin_hull (); pt != poly.end_hull (); ++pt) { + edge_contours.back ().push_back (insert_point (trans * *pt)); + } + + for (unsigned int h = 0; h < poly.holes (); ++h) { + edge_contours.push_back (std::vector ()); + for (auto pt = poly.begin_hole (h); pt != poly.end_hole (h); ++pt) { + edge_contours.back ().push_back (insert_point (trans * *pt)); + } + } +} + +void +Triangulation::create_constrained_delaunay (const db::Region ®ion, const CplxTrans &trans) +{ + clear (); + + std::vector > edge_contours; + + for (auto p = region.begin_merged (); ! p.at_end (); ++p) { + make_contours (*p, trans, edge_contours); + } + + constrain (edge_contours); +} + +void +Triangulation::create_constrained_delaunay (const db::Polygon &p, const std::vector &vertexes, const CplxTrans &trans) +{ + clear (); + + for (auto v = vertexes.begin (); v != vertexes.end (); ++v) { + insert_point (trans * *v)->set_is_precious (true); + } + + std::vector > edge_contours; + make_contours (p, trans, edge_contours); + + constrain (edge_contours); +} + +void +Triangulation::create_constrained_delaunay (const db::DPolygon &p, const std::vector &vertexes, const DCplxTrans &trans) +{ + clear (); + + for (auto v = vertexes.begin (); v != vertexes.end (); ++v) { + insert_point (trans * *v)->set_is_precious (true); + } + + std::vector > edge_contours; + make_contours (p, trans, edge_contours); + + constrain (edge_contours); +} + +static bool is_skinny (const Polygon *tri, const TriangulationParameters ¶m) +{ + if (param.min_b < db::epsilon) { + return false; + } else { + double b = tri->b (); + double delta = (b + param.min_b) * db::epsilon; + return b < param.min_b - delta; + } +} + +static bool is_invalid (const Polygon *tri, const TriangulationParameters ¶m) +{ + if (is_skinny (tri, param)) { + return true; + } + + double amax = param.max_area; + if (param.max_area_border > db::epsilon) { + if (tri->has_segment ()) { + amax = param.max_area_border; + } + } + + if (amax > db::epsilon) { + double a = tri->area (); + double delta = (a + amax) * db::epsilon; + return tri->area () > amax + delta; + } else { + return false; + } +} + +void +Triangulation::triangulate (const db::Region ®ion, const TriangulationParameters ¶meters, double dbu) +{ + tl::SelfTimer timer (tl::verbosity () > parameters.base_verbosity, "Triangles::triangulate"); + + create_constrained_delaunay (region, db::CplxTrans (dbu)); + refine (parameters); +} + +void +Triangulation::triangulate (const db::Region ®ion, const TriangulationParameters ¶meters, const db::CplxTrans &trans) +{ + tl::SelfTimer timer (tl::verbosity () > parameters.base_verbosity, "Triangles::triangulate"); + + create_constrained_delaunay (region, trans); + refine (parameters); +} + +void +Triangulation::triangulate (const db::Polygon &poly, const TriangulationParameters ¶meters, double dbu) +{ + triangulate (poly, std::vector (), parameters, dbu); +} + +void +Triangulation::triangulate (const db::Polygon &poly, const std::vector &vertexes, const TriangulationParameters ¶meters, double dbu) +{ + tl::SelfTimer timer (tl::verbosity () > parameters.base_verbosity, "Triangles::triangulate"); + + create_constrained_delaunay (poly, vertexes, db::CplxTrans (dbu)); + refine (parameters); +} + +void +Triangulation::triangulate (const db::Polygon &poly, const TriangulationParameters ¶meters, const db::CplxTrans &trans) +{ + triangulate (poly, std::vector (), parameters, trans); +} + +void +Triangulation::triangulate (const db::Polygon &poly, const std::vector &vertexes, const TriangulationParameters ¶meters, const db::CplxTrans &trans) +{ + tl::SelfTimer timer (tl::verbosity () > parameters.base_verbosity, "Triangles::triangulate"); + + create_constrained_delaunay (poly, vertexes, trans); + refine (parameters); +} + +void +Triangulation::triangulate (const db::DPolygon &poly, const TriangulationParameters ¶meters, const DCplxTrans &trans) +{ + triangulate (poly, std::vector (), parameters, trans); +} + +void +Triangulation::triangulate (const db::DPolygon &poly, const std::vector &vertexes, const TriangulationParameters ¶meters, const DCplxTrans &trans) +{ + tl::SelfTimer timer (tl::verbosity () > parameters.base_verbosity, "Triangles::triangulate"); + + create_constrained_delaunay (poly, vertexes, trans); + refine (parameters); +} + +void +Triangulation::refine (const TriangulationParameters ¶meters) +{ + if (parameters.min_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 > new_triangles; + for (auto t = mp_graph->polygons ().begin (); t != mp_graph->polygons ().end (); ++t) { + new_triangles.push_back (t.operator-> ()); + } + + // TODO: break if iteration gets stuck + while (nloop < parameters.max_iterations) { + + ++nloop; + if (tl::verbosity () >= parameters.base_verbosity + 10) { + tl::info << "Iteration " << nloop << " .."; + } + + std::list > 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); + } + } + + if (to_consider.empty ()) { + break; + } + + if (tl::verbosity () >= parameters.base_verbosity + 10) { + tl::info << to_consider.size() << " triangles to consider"; + } + + 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; + + int s = (*t)->contains (center); + if (s >= 0) { + + if (s > 0) { + + double snap = 1e-3; + + // Snap the center to a segment center if "close" to it. + // This avoids generating very skinny triangles that can't be fixed as the + // segment cannot be flipped. This a part of the issue #1996 problem. + for (unsigned int i = 0; i < 3; ++i) { + if ((*t)->edge (i)->is_segment ()) { + auto e = (*t)->edge (i)->edge (); + auto c = e.p1 () + e.d () * 0.5; + if (c.distance (center) < e.length () * 0.5 * snap - db::epsilon) { + center = c; + break; + } + } + } + + } + + if (tl::verbosity () >= parameters.base_verbosity + 20) { + tl::info << "Inserting in-triangle center " << center.to_string () << " of " << (*t)->to_string (true); + } + insert_point (center, &new_triangles); + + } else { + + Vertex *vstart = 0; + for (unsigned int i = 0; i < 3; ++i) { + Edge *edge = (*t)->edge (i); + vstart = (*t)->opposite (edge); + if (edge->side_of (*vstart) * edge->side_of (center) < 0) { + break; + } + } + + Edge *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) { + + if (tl::verbosity () >= parameters.base_verbosity + 20) { + tl::info << "Inserting out-of-triangle center " << center << " of " << (*t)->to_string (true); + } + insert_point (center, &new_triangles); + + } else { + + double sr = edge->d ().length () * 0.5; + if (sr >= parameters.min_length) { + + db::DPoint pnew = *edge->v1 () + edge->d () * 0.5; + + if (tl::verbosity () >= parameters.base_verbosity + 20) { + tl::info << "Split edge " << edge->to_string (true) << " at " << pnew.to_string (); + } + Vertex *vnew = insert_point (pnew, &new_triangles); + auto vertexes_in_diametral_circle = find_points_around (vnew, sr); + + std::vector 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 && ! (*v)->is_precious ()) { + to_delete.push_back (*v); + } + } + + if (tl::verbosity () >= parameters.base_verbosity + 20) { + tl::info << " -> found " << to_delete.size () << " vertexes to remove"; + } + for (auto v = to_delete.begin (); v != to_delete.end (); ++v) { + remove (*v, &new_triangles); + } + + } + + } + + } + + } + + } + + if (tl::verbosity () >= parameters.base_verbosity + 20) { + tl::info << "Finishing .."; + } + + if (parameters.mark_triangles) { + + for (auto t = mp_graph->begin (); t != mp_graph->end (); ++t) { + + size_t id = 0; + + if (! t->is_outside ()) { + + if (is_skinny (t.operator-> (), parameters)) { + id |= 1; + } + if (is_invalid (t.operator-> (), parameters)) { + id |= 2; + } + auto cp = t->circumcircle (); + auto vi = find_inside_circle (cp.first, cp.second); + if (! vi.empty ()) { + id |= 4; + } + + } + + (const_cast (t.operator->()))->set_id (id); + + } + + } + + remove_outside_triangles (); +} + +} // namespace plc + +} // namespace db diff --git a/src/db/db/dbPLCTriangulation.h b/src/db/db/dbPLCTriangulation.h new file mode 100644 index 000000000..b3e941f89 --- /dev/null +++ b/src/db/db/dbPLCTriangulation.h @@ -0,0 +1,306 @@ + +/* + + KLayout Layout Viewer + Copyright (C) 2006-2025 Matthias Koefferlein + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +*/ + +#ifndef HDR_dbPLCTriangulation +#define HDR_dbPLCTriangulation + +#include "dbCommon.h" +#include "dbPLC.h" + +#include +#include +#include +#include + +namespace db +{ + +namespace plc +{ + +struct DB_PUBLIC TriangulationParameters +{ + TriangulationParameters () + : min_b (1.0), + min_length (0.0), + max_area (0.0), + max_area_border (0.0), + max_iterations (std::numeric_limits::max ()), + base_verbosity (30), + mark_triangles (false) + { } + + /** + * @brief Min. readius-to-shortest edge ratio + */ + double min_b; + + /** + * @brief Min. edge length + * + * This parameter does not provide a guarantee about a minimume edge length, but + * helps avoiding ever-reducing triangle splits in acute corners of the input polygon. + * Splitting of edges stops when the edge is less than the min length. + */ + double min_length; + + /** + * @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; + + /** + * @brief Max number of iterations + */ + size_t max_iterations; + + /** + * @brief The verbosity level above which triangulation reports details + */ + int base_verbosity; + + /** + * @brief If true, final triangles are marked using the "id" integer as a bit field + * + * This provides information about the result quality. + * + * Bit 0: skinny triangle + * Bit 1: bad-quality (skinny or area too large) + * Bit 2: non-Delaunay (in the strict sense) + */ + bool mark_triangles; +}; + +/** + * @brief A Triangulation algorithm + * + * This class implements a constrained refined Delaunay triangulation using Chew's algorithm. + */ +class DB_PUBLIC Triangulation +{ +public: + /** + * @brief The constructor + * + * The graph will be one filled by the triangulation. + */ + Triangulation (Graph *graph); + + /** + * @brief Clears the triangulation + */ + void clear (); + + /** + * @brief Initializes the triangle collection with a box + * Two triangles will be created. + */ + void init_box (const db::DBox &box); + + /** + * @brief Creates a refined Delaunay triangulation for the given region + * + * The database unit should be chosen in a way that target area values are "in the order of 1". + * For inputs featuring acute angles (angles < ~25 degree), the parameters should defined a min + * edge length ("min_length"). + * "min_length" should be at least 1e-4. If a min edge length is given, the max area constaints + * may not be satisfied. + * + * Edges in the input should not be shorter than 1e-4. + */ + void triangulate (const db::Region ®ion, const TriangulationParameters ¶meters, double dbu = 1.0); + + // more versions + void triangulate (const db::Region ®ion, const TriangulationParameters ¶meters, const db::CplxTrans &trans = db::CplxTrans ()); + void triangulate (const db::Polygon &poly, const TriangulationParameters ¶meters, double dbu = 1.0); + void triangulate (const db::Polygon &poly, const std::vector &vertexes, const TriangulationParameters ¶meters, double dbu = 1.0); + void triangulate (const db::Polygon &poly, const TriangulationParameters ¶meters, const db::CplxTrans &trans = db::CplxTrans ()); + void triangulate (const db::Polygon &poly, const std::vector &vertexes, const TriangulationParameters ¶meters, const db::CplxTrans &trans = db::CplxTrans ()); + + /** + * @brief Triangulates a floating-point polygon + */ + void triangulate (const db::DPolygon &poly, const TriangulationParameters ¶meters, const db::DCplxTrans &trans = db::DCplxTrans ()); + void triangulate (const db::DPolygon &poly, const std::vector &vertexes, const TriangulationParameters ¶meters, const db::DCplxTrans &trans = db::DCplxTrans ()); + + /** + * @brief Inserts a new vertex as the given point + * + * If "new_triangles" is not null, it will receive the list of new triangles created during + * the remove step. + * + * This method can be called after "triangulate" to add new points and adjust the triangulation. + * Inserting new points will maintain the (constrained) Delaunay condition. + */ + Vertex *insert_point (const db::DPoint &point, std::list > *new_triangles = 0); + + /** + * @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; + } + +protected: + /** + * @brief Checks the polygon graph for consistency + * This method is for testing purposes mainly. + */ + bool check (bool check_delaunay = true) const; + + /** + * @brief Finds the points within (not "on") a circle of radius "radius" around the given vertex. + */ + std::vector find_points_around (Vertex *vertex, double radius); + + /** + * @brief Inserts a new vertex as the given point + * + * If "new_triangles" is not null, it will receive the list of new triangles created during + * the remove step. + */ + Vertex *insert_point (db::DCoord x, db::DCoord y, std::list > *new_triangles = 0); + + /** + * @brief Removes the given vertex + * + * If "new_triangles" is not null, it will receive the list of new triangles created during + * the remove step. + */ + void remove (Vertex *vertex, std::list > *new_triangles = 0); + + /** + * @brief Flips the given edge + */ + std::pair, Edge *> flip (Edge *edge); + + /** + * @brief Finds all edges that cross the given one for a convex triangulation + * + * Requirements: + * * self must be a convex triangulation + * * edge must not contain another vertex from the triangulation except p1 and p2 + */ + std::vector search_edges_crossing (Vertex *from, Vertex *to); + + /** + * @brief Finds the edge for two given points + */ + Edge *find_edge_for_points (const db::DPoint &p1, const db::DPoint &p2); + + /** + * @brief Finds the vertex for a point + */ + Vertex *find_vertex_for_point (const db::DPoint &pt); + + /** + * @brief Ensures all points between from an to are connected by edges and makes these segments + */ + std::vector ensure_edge (Vertex *from, Vertex *to); + + /** + * @brief Given a set of contours with edges, mark outer triangles + * + * The edges must be made from existing vertexes. Edge orientation is + * clockwise. + * + * This will also mark triangles as outside ones. + */ + void constrain (const std::vector > &contours); + + /** + * @brief Removes the outside triangles. + */ + void remove_outside_triangles (); + + /** + * @brief Creates a constrained Delaunay triangulation from the given Region + */ + void create_constrained_delaunay (const db::Region ®ion, const db::CplxTrans &trans = db::CplxTrans ()); + + /** + * @brief Creates a constrained Delaunay triangulation from the given Polygon + */ + void create_constrained_delaunay (const db::Polygon &poly, const std::vector &vertexes, const db::CplxTrans &trans = db::CplxTrans ()); + + /** + * @brief Creates a constrained Delaunay triangulation from the given DPolygon + */ + void create_constrained_delaunay (const db::DPolygon &poly, const std::vector &vertexes, const DCplxTrans &trans); + + /** + * @brief Returns a value indicating whether the edge is "illegal" (violates the Delaunay criterion) + */ + static bool is_illegal_edge (Edge *edge); + + // NOTE: these functions are SLOW and intended to test purposes only + std::vector find_touching (const db::DBox &box) const; + std::vector find_inside_circle (const db::DPoint ¢er, double radius) const; + +private: + Graph *mp_graph; + bool m_is_constrained; + size_t m_level; + size_t m_id; + size_t m_flips, m_hops; + + template void make_contours (const Poly &poly, const Trans &trans, std::vector > &contours); + + void remove_outside_vertex (Vertex *vertex, std::list > *new_triangles = 0); + void remove_inside_vertex (Vertex *vertex, std::list > *new_triangles_out = 0); + std::vector fill_concave_corners (const std::vector &edges); + void fix_triangles (const std::vector &tris, const std::vector &fixed_edges, std::list > *new_triangles); + std::vector find_triangle_for_point (const db::DPoint &point); + Edge *find_closest_edge (const db::DPoint &p, Vertex *vstart = 0, bool inside_only = false); + Vertex *insert (Vertex *vertex, std::list > *new_triangles = 0); + void split_triangle (Polygon *t, Vertex *vertex, std::list > *new_triangles_out); + void split_triangles_on_edge (Vertex *vertex, Edge *split_edge, std::list > *new_triangles_out); + void add_more_triangles (std::vector &new_triangles, + Edge *incoming_edge, + Vertex *from_vertex, Vertex *to_vertex, + Edge *conn_edge); + void insert_new_vertex(Vertex *vertex, std::list > *new_triangles_out); + std::vector ensure_edge_inner (Vertex *from, Vertex *to); + void join_edges (std::vector &edges); + void refine (const TriangulationParameters ¶m); +}; + +} // namespace plc + +} // namespace db + +#endif + diff --git a/src/db/unit_tests/dbPolygonGraphTests.cc b/src/db/unit_tests/dbPLCGraphTest.cc similarity index 76% rename from src/db/unit_tests/dbPolygonGraphTests.cc rename to src/db/unit_tests/dbPLCGraphTest.cc index edf152942..e02449401 100644 --- a/src/db/unit_tests/dbPolygonGraphTests.cc +++ b/src/db/unit_tests/dbPLCGraphTest.cc @@ -21,31 +21,21 @@ */ -#include "dbPolygonGraph.h" +#include "dbPLC.h" #include "tlUnitTest.h" #include #include -class TestablePolygonGraph - : public db::PolygonGraph -{ -public: - using db::PolygonGraph::PolygonGraph; - // @@@ using db::PolygonGraph::check; - using db::PolygonGraph::dump; -}; - - TEST(basic) { db::DBox box (0, 0, 100.0, 200.0); - TestablePolygonGraph pg; + db::plc::Graph plc; // @@@ pg.insert_polygon (db::DSimplePolygon (box)); // @@@ - tl::info << pg.to_string (); - pg.dump ("debug.gds"); // @@@ + tl::info << plc.to_string (); + plc.dump ("debug.gds"); // @@@ } diff --git a/src/db/unit_tests/dbPLCTriangulationTests.cc b/src/db/unit_tests/dbPLCTriangulationTests.cc new file mode 100644 index 000000000..0d58ffc04 --- /dev/null +++ b/src/db/unit_tests/dbPLCTriangulationTests.cc @@ -0,0 +1,1117 @@ + +/* + + KLayout Layout Viewer + Copyright (C) 2006-2025 Matthias Koefferlein + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +*/ + + +#include "dbPLCTriangulation.h" +#include "dbWriter.h" +#include "dbRegionProcessors.h" +#include "tlUnitTest.h" +#include "tlStream.h" +#include "tlFileUtils.h" + +#include +#include +#include +#include + +class TestableTriangulation + : public db::plc::Triangulation +{ +public: + using db::plc::Triangulation::Triangulation; + using db::plc::Triangulation::check; + using db::plc::Triangulation::flip; + using db::plc::Triangulation::insert_point; + using db::plc::Triangulation::search_edges_crossing; + using db::plc::Triangulation::find_edge_for_points; + using db::plc::Triangulation::find_points_around; + using db::plc::Triangulation::find_inside_circle; + using db::plc::Triangulation::create_constrained_delaunay; + using db::plc::Triangulation::is_illegal_edge; + using db::plc::Triangulation::find_vertex_for_point; + using db::plc::Triangulation::remove; + using db::plc::Triangulation::ensure_edge; + using db::plc::Triangulation::constrain; + using db::plc::Triangulation::remove_outside_triangles; +}; + +class TestableGraph + : public db::plc::Graph +{ +public: + using db::plc::Graph::Graph; + using db::plc::Graph::create_vertex; + using db::plc::Graph::create_edge; + using db::plc::Graph::create_triangle; +}; + +TEST(basic) +{ + db::plc::Graph plc; + TestableTriangulation tris (&plc); + tris.init_box (db::DBox (1, 0, 5, 4)); + + EXPECT_EQ (plc.bbox ().to_string (), "(1,0;5,4)"); + EXPECT_EQ (plc.to_string (), "((1, 0), (1, 4), (5, 0)), ((1, 4), (5, 4), (5, 0))"); + + EXPECT_EQ (tris.check (), true); + + tris.clear (); + + EXPECT_EQ (plc.bbox ().to_string (), "()"); + EXPECT_EQ (plc.to_string (), ""); + + EXPECT_EQ (tris.check (), true); +} + +TEST(flip) +{ + db::plc::Graph plc; + TestableTriangulation tris (&plc); + tris.init_box (db::DBox (0, 0, 1, 1)); + EXPECT_EQ (plc.to_string (), "((0, 0), (0, 1), (1, 0)), ((0, 1), (1, 1), (1, 0))"); + + EXPECT_EQ (plc.num_polygons (), size_t (2)); + EXPECT_EQ (tris.check (), true); + + const db::plc::Polygon &t1 = *plc.begin (); + db::plc::Edge *diag_segment; + for (int i = 0; i < 3; ++i) { + diag_segment = t1.edge (i); + if (diag_segment->side_of (db::DPoint (0.5, 0.5)) == 0) { + break; + } + } + tris.flip (diag_segment); + EXPECT_EQ (plc.to_string (), "((1, 1), (0, 0), (0, 1)), ((1, 1), (1, 0), (0, 0))"); + EXPECT_EQ (tris.check (), true); +} + +TEST(insert) +{ + db::plc::Graph plc; + TestableTriangulation tris (&plc); + tris.init_box (db::DBox (0, 0, 1, 1)); + + tris.insert_point (0.2, 0.2); + EXPECT_EQ (plc.to_string (), "((0, 0), (0, 1), (0.2, 0.2)), ((1, 0), (0, 0), (0.2, 0.2)), ((1, 1), (0.2, 0.2), (0, 1)), ((1, 1), (1, 0), (0.2, 0.2))"); + EXPECT_EQ (tris.check (), true); +} + +TEST(split_segment) +{ + db::plc::Graph plc; + TestableTriangulation tris (&plc); + tris.init_box (db::DBox (0, 0, 1, 1)); + + tris.insert_point (0.5, 0.5); + EXPECT_EQ (plc.to_string (), "((1, 1), (1, 0), (0.5, 0.5)), ((1, 1), (0.5, 0.5), (0, 1)), ((0, 0), (0, 1), (0.5, 0.5)), ((0, 0), (0.5, 0.5), (1, 0))"); + EXPECT_EQ (tris.check(), true); +} + +TEST(insert_vertex_twice) +{ + db::plc::Graph plc; + TestableTriangulation tris (&plc); + tris.init_box (db::DBox (0, 0, 1, 1)); + + tris.insert_point (0.5, 0.5); + // inserted a vertex twice does not change anything + tris.insert_point (0.5, 0.5); + EXPECT_EQ (plc.to_string (), "((1, 1), (1, 0), (0.5, 0.5)), ((1, 1), (0.5, 0.5), (0, 1)), ((0, 0), (0, 1), (0.5, 0.5)), ((0, 0), (0.5, 0.5), (1, 0))"); + EXPECT_EQ (tris.check(), true); +} + +TEST(insert_vertex_convex) +{ + db::plc::Graph plc; + TestableTriangulation tris (&plc); + tris.insert_point (0.2, 0.2); + tris.insert_point (0.2, 0.8); + tris.insert_point (0.6, 0.5); + tris.insert_point (0.7, 0.5); + tris.insert_point (0.6, 0.4); + EXPECT_EQ (plc.to_string (), "((0.2, 0.2), (0.2, 0.8), (0.6, 0.5)), ((0.2, 0.8), (0.7, 0.5), (0.6, 0.5)), ((0.6, 0.4), (0.6, 0.5), (0.7, 0.5)), ((0.6, 0.4), (0.2, 0.2), (0.6, 0.5))"); + EXPECT_EQ (tris.check(), true); +} + +TEST(insert_vertex_convex2) +{ + db::plc::Graph plc; + TestableTriangulation tris (&plc); + tris.insert_point (0.25, 0.1); + tris.insert_point (0.1, 0.4); + tris.insert_point (0.4, 0.15); + tris.insert_point (1, 0.7); + EXPECT_EQ (plc.to_string (), "((0.25, 0.1), (0.1, 0.4), (0.4, 0.15)), ((1, 0.7), (0.4, 0.15), (0.1, 0.4))"); + EXPECT_EQ (tris.check(), true); +} + +TEST(insert_vertex_convex3) +{ + db::plc::Graph plc; + TestableTriangulation tris (&plc); + tris.insert_point (0.25, 0.5); + tris.insert_point (0.25, 0.55); + tris.insert_point (0.15, 0.8); + tris.insert_point (1, 0.4); + EXPECT_EQ (plc.to_string (), "((0.25, 0.5), (0.15, 0.8), (0.25, 0.55)), ((1, 0.4), (0.25, 0.5), (0.25, 0.55)), ((0.15, 0.8), (1, 0.4), (0.25, 0.55))"); + EXPECT_EQ (tris.check(), true); +} + +TEST(search_edges_crossing) +{ + db::plc::Graph plc; + TestableTriangulation tris (&plc); + db::plc::Vertex *v1 = tris.insert_point (0.2, 0.2); + db::plc::Vertex *v2 = tris.insert_point (0.2, 0.8); + db::plc::Vertex *v3 = tris.insert_point (0.6, 0.5); + /*db::plc::Vertex *v4 =*/ tris.insert_point (0.7, 0.5); + db::plc::Vertex *v5 = tris.insert_point (0.6, 0.4); + db::plc::Vertex *v6 = tris.insert_point (0.7, 0.2); + EXPECT_EQ (tris.check(), true); + + auto xedges = tris.search_edges_crossing (v2, v6); + + EXPECT_EQ (xedges.size (), size_t (2)); + auto s1 = tris.find_edge_for_points (*v1, *v3); + auto s2 = tris.find_edge_for_points (*v1, *v5); + EXPECT_EQ (std::find (xedges.begin (), xedges.end (), s1) != xedges.end (), true); + EXPECT_EQ (std::find (xedges.begin (), xedges.end (), s2) != xedges.end (), true); +} + +TEST(illegal_edge1) +{ + TestableGraph plc; + + db::plc::Vertex *v1 = plc.create_vertex (0, 0); + db::plc::Vertex *v2 = plc.create_vertex (1.6, 1.6); + db::plc::Vertex *v3 = plc.create_vertex (1, 2); + db::plc::Vertex *v4 = plc.create_vertex (2, 1); + + { + db::plc::Edge *e1 = plc.create_edge (v1, v3); + db::plc::Edge *e2 = plc.create_edge (v3, v4); + db::plc::Edge *e3 = plc.create_edge (v4, v1); + + plc.create_triangle (e1, e2, e3); + + db::plc::Edge *ee1 = plc.create_edge (v2, v3); + db::plc::Edge *ee2 = plc.create_edge (v4, v2); + + plc.create_triangle (ee1, e2, ee2); + + EXPECT_EQ (TestableTriangulation::is_illegal_edge (e2), true); + } + + { + // flipped + db::plc::Edge *e1 = plc.create_edge (v1, v2); + db::plc::Edge *e2 = plc.create_edge (v2, v3); + db::plc::Edge *e3 = plc.create_edge (v3, v1); + + plc.create_triangle (e1, e2, e3); + + db::plc::Edge *ee1 = plc.create_edge (v1, v4); + db::plc::Edge *ee2 = plc.create_edge (v4, v2); + + plc.create_triangle (ee1, ee2, e1); + + EXPECT_EQ (TestableTriangulation::is_illegal_edge (e2), false); + } +} + +TEST(illegal_edge2) +{ + TestableGraph plc; + + // numerical border case + db::plc::Vertex *v1 = plc.create_vertex (773.94756216690905, 114.45875269431208); + db::plc::Vertex *v2 = plc.create_vertex (773.29574734131643, 113.47402096138073); + db::plc::Vertex *v3 = plc.create_vertex (773.10652961562653, 114.25497975904504); + db::plc::Vertex *v4 = plc.create_vertex (774.08856345337881, 113.60495072750861); + + { + db::plc::Edge *e1 = plc.create_edge (v1, v2); + db::plc::Edge *e2 = plc.create_edge (v2, v4); + db::plc::Edge *e3 = plc.create_edge (v4, v1); + + plc.create_triangle (e1, e2, e3); + + db::plc::Edge *ee1 = plc.create_edge (v2, v3); + db::plc::Edge *ee2 = plc.create_edge (v3, v4); + + plc.create_triangle (ee1, ee2, e2); + + EXPECT_EQ (TestableTriangulation::is_illegal_edge (e2), false); + } + + { + // flipped + db::plc::Edge *e1 = plc.create_edge (v1, v2); + db::plc::Edge *e2 = plc.create_edge (v2, v3); + db::plc::Edge *e3 = plc.create_edge (v3, v1); + + plc.create_triangle (e1, e2, e3); + + db::plc::Edge *ee1 = plc.create_edge (v1, v4); + db::plc::Edge *ee2 = plc.create_edge (v4, v2); + + plc.create_triangle (ee1, ee2, e1); + + EXPECT_EQ (TestableTriangulation::is_illegal_edge (e1), false); + } +} + +// Returns a random float number between 0.0 and 1.0 +inline double flt_rand () +{ + return rand () * (1.0 / double (RAND_MAX)); +} + +namespace { + struct PointLessOp + { + bool operator() (const db::DPoint &a, const db::DPoint &b) const + { + return a.less (b); + } + }; +} + +TEST(insert_many) +{ + srand (0); + + db::plc::Graph plc; + TestableTriangulation tris (&plc); + double res = 65536.0; + + db::DBox bbox; + + unsigned int n = 200000; + for (unsigned int i = 0; i < n; ++i) { + double x = round (flt_rand () * res) * 0.0001; + double y = round (flt_rand () * res) * 0.0001; + tris.insert_point (x, y); + } + + // slow: EXPECT_EQ (tris.check (), true); + EXPECT_LT (double (tris.flips ()) / double (n), 3.1); + EXPECT_LT (double (tris.hops ()) / double (n), 23.0); +} + +TEST(heavy_insert) +{ + tl::info << "Running test_heavy_insert " << tl::noendl; + + for (unsigned int l = 0; l < 100; ++l) { + + srand (l); + tl::info << "." << tl::noendl; + + db::plc::Graph plc; + TestableTriangulation tris (&plc); + double res = 128.0; + + unsigned int n = rand () % 190 + 10; + + db::DBox bbox; + std::map vmap; + + for (unsigned int i = 0; i < n; ++i) { + double x = round (flt_rand () * res) * (1.0 / res); + double y = round (flt_rand () * res) * (1.0 / res); + db::plc::Vertex *v = tris.insert_point (x, y); + bbox += db::DPoint (x, y); + vmap.insert (std::make_pair (*v, false)); + } + + // not strictly true, but very likely with at least 10 vertexes: + EXPECT_GT (plc.num_polygons (), size_t (0)); + EXPECT_EQ (plc.bbox ().to_string (), bbox.to_string ()); + + bool ok = true; + for (auto t = plc.begin (); t != plc.end (); ++t) { + for (int i = 0; i < 3; ++i) { + auto f = vmap.find (*t->vertex (i)); + if (f == vmap.end ()) { + tl::error << "Could not identify triangle vertex " << t->vertex (i)->to_string () << " as inserted vertex"; + ok = false; + } else { + f->second = true; + } + } + } + for (auto m = vmap.begin (); m != vmap.end (); ++m) { + if (!m->second) { + tl::error << "Could not identify vertex " << m->first.to_string () << " with a triangle"; + ok = false; + } + } + EXPECT_EQ (ok, true); + + EXPECT_EQ (tris.check(), true); + + } + + tl::info << tl::endl << "done."; +} + +TEST(heavy_remove) +{ + tl::info << "Running test_heavy_remove " << tl::noendl; + + for (unsigned int l = 0; l < 100; ++l) { + + srand (l); + tl::info << "." << tl::noendl; + + db::plc::Graph plc; + TestableTriangulation tris (&plc); + double res = 128.0; + + unsigned int n = rand () % 190 + 10; + + for (unsigned int i = 0; i < n; ++i) { + double x = round (flt_rand () * res) * (1.0 / res); + double y = round (flt_rand () * res) * (1.0 / res); + tris.insert_point (x, y); + } + + EXPECT_EQ (tris.check(), true); + + std::set vset; + std::vector vertexes; + for (auto t = plc.begin (); t != plc.end (); ++t) { + for (int i = 0; i < 3; ++i) { + db::plc::Vertex *v = t->vertex (i); + if (vset.insert (v).second) { + vertexes.push_back (v); + } + } + } + + while (! vertexes.empty ()) { + + unsigned int n = rand () % (unsigned int) vertexes.size (); + db::plc::Vertex *v = vertexes [n]; + tris.remove (v); + vertexes.erase (vertexes.begin () + n); + + // just a few times as it wastes time otherwise + if (vertexes.size () % 10 == 0) { + EXPECT_EQ (tris.check (), true); + } + + } + + EXPECT_EQ (plc.num_polygons (), size_t (0)); + + } + + tl::info << tl::endl << "done."; +} + +TEST(ensure_edge) +{ + srand (0); + + db::plc::Graph plc; + TestableTriangulation tris (&plc); + double res = 128.0; + + db::DEdge ee[] = { + db::DEdge (0.25, 0.25, 0.25, 0.75), + db::DEdge (0.25, 0.75, 0.75, 0.75), + db::DEdge (0.75, 0.75, 0.75, 0.25), + db::DEdge (0.75, 0.25, 0.25, 0.25) + }; + + for (unsigned int i = 0; i < 200; ++i) { + double x = round (flt_rand () * res) * (1.0 / res); + double y = round (flt_rand () * res) * (1.0 / res); + bool ok = true; + for (unsigned int j = 0; j < sizeof (ee) / sizeof (ee[0]); ++j) { + if (ee[j].side_of (db::DPoint (x, y)) == 0) { + --i; + ok = false; + } + } + if (ok) { + tris.insert_point (x, y); + } + } + + for (unsigned int i = 0; i < sizeof (ee) / sizeof (ee[0]); ++i) { + tris.insert_point (ee[i].p1 ()); + } + + EXPECT_EQ (tris.check (), true); + + for (unsigned int i = 0; i < sizeof (ee) / sizeof (ee[0]); ++i) { + tris.ensure_edge (tris.find_vertex_for_point (ee[i].p1 ()), tris.find_vertex_for_point (ee[i].p2 ())); + } + + EXPECT_EQ (tris.check (false), true); + + double area_in = 0.0; + db::DBox clip_box; + for (unsigned int i = 0; i < sizeof (ee) / sizeof (ee[0]); ++i) { + clip_box += ee[i].p1 (); + } + for (auto t = plc.begin (); t != plc.end (); ++t) { + if (clip_box.overlaps (t->bbox ())) { + EXPECT_EQ (t->bbox ().inside (clip_box), true); + area_in += t->area (); + } + } + + EXPECT_EQ (tl::to_string (area_in), "0.25"); +} + +static bool safe_inside (const db::DBox &b1, const db::DBox &b2) +{ + typedef db::coord_traits ct; + + return (ct::less (b2.left (), b1.left ()) || ct::equal (b2.left (), b1.left ())) && + (ct::less (b1.right (), b2.right ()) || ct::equal (b1.right (), b2.right ())) && + (ct::less (b2.bottom (), b1.bottom ()) || ct::equal (b2.bottom (), b1.bottom ())) && + (ct::less (b1.top (), b2.top ()) || ct::equal (b1.top (), b2.top ())); +} + +TEST(constrain) +{ + srand (0); + + db::plc::Graph plc; + TestableTriangulation tris (&plc); + double res = 128.0; + + db::DEdge ee[] = { + db::DEdge (0.25, 0.25, 0.25, 0.75), + db::DEdge (0.25, 0.75, 0.75, 0.75), + db::DEdge (0.75, 0.75, 0.75, 0.25), + db::DEdge (0.75, 0.25, 0.25, 0.25) + }; + + for (unsigned int i = 0; i < 200; ++i) { + double x = round (flt_rand () * res) * (1.0 / res); + double y = round (flt_rand () * res) * (1.0 / res); + bool ok = true; + for (unsigned int j = 0; j < sizeof (ee) / sizeof (ee[0]); ++j) { + if (ee[j].side_of (db::DPoint (x, y)) == 0) { + --i; + ok = false; + } + } + if (ok) { + tris.insert_point (x, y); + } + } + + std::vector contour; + for (unsigned int i = 0; i < sizeof (ee) / sizeof (ee[0]); ++i) { + contour.push_back (tris.insert_point (ee[i].p1 ())); + } + std::vector > contours; + contours.push_back (contour); + + EXPECT_EQ (tris.check (), true); + + tris.constrain (contours); + EXPECT_EQ (tris.check (false), true); + + tris.remove_outside_triangles (); + + EXPECT_EQ (tris.check (), true); + + double area_in = 0.0; + db::DBox clip_box; + for (unsigned int i = 0; i < sizeof (ee) / sizeof (ee[0]); ++i) { + clip_box += ee[i].p1 (); + } + for (auto t = plc.begin (); t != plc.end (); ++t) { + EXPECT_EQ (clip_box.overlaps (t->bbox ()), true); + EXPECT_EQ (safe_inside (t->bbox (), clip_box), true); + area_in += t->area (); + } + + EXPECT_EQ (tl::to_string (area_in), "0.25"); +} + +TEST(heavy_constrain) +{ + tl::info << "Running test_heavy_constrain " << tl::noendl; + + for (unsigned int l = 0; l < 100; ++l) { + + srand (l); + tl::info << "." << tl::noendl; + + db::plc::Graph plc; + TestableTriangulation tris (&plc); + double res = 128.0; + + db::DEdge ee[] = { + db::DEdge (0.25, 0.25, 0.25, 0.75), + db::DEdge (0.25, 0.75, 0.75, 0.75), + db::DEdge (0.75, 0.75, 0.75, 0.25), + db::DEdge (0.75, 0.25, 0.25, 0.25) + }; + + unsigned int n = rand () % 150 + 50; + + for (unsigned int i = 0; i < n; ++i) { + double x = round (flt_rand () * res) * (1.0 / res); + double y = round (flt_rand () * res) * (1.0 / res); + bool ok = true; + for (unsigned int j = 0; j < sizeof (ee) / sizeof (ee[0]); ++j) { + if (ee[j].side_of (db::DPoint (x, y)) == 0) { + --i; + ok = false; + } + } + if (ok) { + tris.insert_point (x, y); + } + } + + std::vector contour; + for (unsigned int i = 0; i < sizeof (ee) / sizeof (ee[0]); ++i) { + contour.push_back (tris.insert_point (ee[i].p1 ())); + } + std::vector > contours; + contours.push_back (contour); + + EXPECT_EQ (tris.check (), true); + + tris.constrain (contours); + EXPECT_EQ (tris.check (false), true); + + tris.remove_outside_triangles (); + + EXPECT_EQ (tris.check (), true); + + double area_in = 0.0; + db::DBox clip_box; + for (unsigned int i = 0; i < sizeof (ee) / sizeof (ee[0]); ++i) { + clip_box += ee[i].p1 (); + } + for (auto t = plc.begin (); t != plc.end (); ++t) { + EXPECT_EQ (clip_box.overlaps (t->bbox ()), true); + EXPECT_EQ (safe_inside (t->bbox (), clip_box), true); + area_in += t->area (); + } + + EXPECT_EQ (tl::to_string (area_in), "0.25"); + + } + + tl::info << tl::endl << "done."; +} + +TEST(heavy_find_point_around) +{ + tl::info << "Running Triangle_test_heavy_find_point_around " << tl::noendl; + + for (unsigned int l = 0; l < 100; ++l) { + + srand (l); + tl::info << "." << tl::noendl; + + db::plc::Graph plc; + TestableTriangulation tris (&plc); + double res = 128.0; + + unsigned int n = rand () % 190 + 10; + + std::vector vertexes; + + for (unsigned int i = 0; i < n; ++i) { + double x = round (flt_rand () * res) * (1.0 / res); + double y = round (flt_rand () * res) * (1.0 / res); + vertexes.push_back (tris.insert_point (x, y)); + } + + EXPECT_EQ (tris.check(), true); + + for (int i = 0; i < 100; ++i) { + + unsigned int nv = rand () % (unsigned int) vertexes.size (); + auto vertex = vertexes [nv]; + + double r = round (flt_rand () * res) * (1.0 / res); + auto p1 = tris.find_points_around (vertex, r); + auto p2 = tris.find_inside_circle (*vertex, r); + + std::set sp1 (p1.begin (), p1.end ()); + std::set sp2 (p2.begin (), p2.end ()); + sp2.erase (vertex); + + EXPECT_EQ (sp1 == sp2, true); + + } + + } + + tl::info << tl::endl << "done."; +} + +TEST(create_constrained_delaunay) +{ + 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::plc::Graph plc; + TestableTriangulation tris (&plc); + tris.create_constrained_delaunay (r); + tris.remove_outside_triangles (); + + EXPECT_EQ (tris.check (), true); + + EXPECT_EQ (plc.to_string (), + "((1000, 0), (0, 0), (200, 200)), " + "((0, 1000), (200, 200), (0, 0)), " + "((1000, 0), (200, 200), (800, 200)), " + "((1000, 0), (800, 200), (1000, 1000)), " + "((800, 200), (800, 800), (1000, 1000)), " + "((0, 1000), (1000, 1000), (800, 800)), " + "((0, 1000), (800, 800), (200, 800)), " + "((0, 1000), (200, 800), (200, 200))"); +} + +TEST(triangulate_basic) +{ + db::Region r; + r.insert (db::Box (0, 0, 10000, 10000)); + + db::Region r2; + r2.insert (db::Box (2000, 2000, 8000, 8000)); + + r -= r2; + + db::plc::TriangulationParameters param; + param.min_b = 1.2; + param.max_area = 1.0; + + db::plc::Graph plc; + TestableTriangulation tri (&plc); + tri.triangulate (r, param, 0.001); + + EXPECT_EQ (tri.check (), true); + + for (auto t = plc.begin (); t != plc.end (); ++t) { + EXPECT_LE (t->area (), param.max_area); + EXPECT_GE (t->b (), param.min_b); + } + + EXPECT_GT (plc.num_polygons (), size_t (100)); + EXPECT_LT (plc.num_polygons (), size_t (150)); + + // for debugging: + // tri.dump ("debug.gds"); + + param.min_b = 1.0; + param.max_area = 0.1; + + tri.triangulate (r, param, 0.001); + + EXPECT_EQ (tri.check (), true); + + for (auto t = plc.begin (); t != plc.end (); ++t) { + EXPECT_LE (t->area (), param.max_area); + EXPECT_GE (t->b (), param.min_b); + } + + EXPECT_GT (plc.num_polygons (), size_t (900)); + EXPECT_LT (plc.num_polygons (), size_t (1000)); +} + +static void read_polygons (const std::string &path, db::Region ®ion, double dbu) +{ + tl::InputStream is (path); + tl::TextInputStream ti (is); + + unsigned int nvert = 0, nedges = 0; + + { + tl::Extractor ex (ti.get_line ().c_str ()); + ex.read (nvert); + ex.read (nedges); + } + + std::vector v; + auto dbu_trans = db::CplxTrans (dbu).inverted (); + for (unsigned int i = 0; i < nvert; ++i) { + double x = 0, y = 0; + tl::Extractor ex (ti.get_line ().c_str ()); + ex.read (x); + ex.read (y); + v.push_back (dbu_trans * db::DPoint (x, y)); + } + + unsigned int nstart = 0; + bool new_contour = true; + std::vector contour; + + for (unsigned int i = 0; i < nedges; ++i) { + + unsigned int n1 = 0, n2 = 0; + + tl::Extractor ex (ti.get_line ().c_str ()); + ex.read (n1); + ex.read (n2); + + if (new_contour) { + nstart = n1; + new_contour = false; + } + + contour.push_back (v[n1]); + + if (n2 == nstart) { + // finish contour + db::SimplePolygon sp; + sp.assign_hull (contour.begin (), contour.end ()); + region.insert (sp); + new_contour = true; + contour.clear (); + } else if (n2 <= n1) { + tl::error << "Invalid polygon wrap in line " << ti.line_number (); + tl_assert (false); + } + + } +} + +TEST(triangulate_geo) +{ + double dbu = 0.001; + + db::Region r; + read_polygons (tl::combine_path (tl::testsrc (), "testdata/algo/triangles1.txt"), r, dbu); + + // for debugging purposes dump the inputs + if (false) { + + db::Layout layout = db::Layout (); + layout.dbu (dbu); + db::Cell &top = layout.cell (layout.add_cell ("DUMP")); + unsigned int l1 = layout.insert_layer (db::LayerProperties (1, 0)); + r.insert_into (&layout, top.cell_index (), l1); + + { + tl::OutputStream stream ("input.gds"); + db::SaveLayoutOptions opt; + db::Writer writer (opt); + writer.write (layout, stream); + } + + } + + db::plc::TriangulationParameters param; + param.min_b = 1.0; + param.max_area = 0.1; + param.min_length = 0.001; + + db::plc::Graph plc; + TestableTriangulation tri (&plc); + tri.triangulate (r, param, dbu); + + EXPECT_EQ (tri.check (false), true); + + // for debugging: + // tri.dump ("debug.gds"); + + size_t n_skinny = 0; + for (auto t = plc.begin (); t != plc.end (); ++t) { + EXPECT_LE (t->area (), param.max_area); + if (t->b () < param.min_b) { + ++n_skinny; + } + } + + EXPECT_LT (n_skinny, size_t (20)); + EXPECT_GT (plc.num_polygons (), size_t (29000)); + EXPECT_LT (plc.num_polygons (), size_t (30000)); +} + +TEST(triangulate_analytic) +{ + double dbu = 0.0001; + + double star1 = 9.0, star2 = 5.0; + double r = 1.0; + int n = 100; + + auto dbu_trans = db::CplxTrans (dbu).inverted (); + + std::vector contour1, contour2; + for (int i = 0; i < n; ++i) { + double a = -M_PI * 2.0 * double (i) / double (n); // "-" for clockwise orientation + double rr, x, y; + rr = r * (1.0 + 0.4 * cos (star1 * a)); + x = rr * cos (a); + y = rr * sin (a); + contour1.push_back (dbu_trans * db::DPoint (x, y)); + rr = r * (0.1 + 0.03 * cos (star2 * a)); + x = rr * cos (a); + y = rr * sin (a); + contour2.push_back (dbu_trans * db::DPoint (x, y)); + } + + db::Region rg; + + db::SimplePolygon sp1; + sp1.assign_hull (contour1.begin (), contour1.end ()); + db::SimplePolygon sp2; + sp2.assign_hull (contour2.begin (), contour2.end ()); + + rg = db::Region (sp1) - db::Region (sp2); + + db::plc::TriangulationParameters param; + param.min_b = 1.0; + param.max_area = 0.01; + + db::plc::Graph plc; + TestableTriangulation tri (&plc); + tri.triangulate (rg, param, dbu); + + EXPECT_EQ (tri.check (false), true); + + // for debugging: + // tri.dump ("debug.gds"); + + for (auto t = plc.begin (); t != plc.end (); ++t) { + EXPECT_LE (t->area (), param.max_area); + EXPECT_GE (t->b (), param.min_b); + } + + EXPECT_GT (plc.num_polygons (), size_t (1250)); + EXPECT_LT (plc.num_polygons (), size_t (1300)); +} + +TEST(triangulate_problematic) +{ + db::DPoint contour[] = { + db::DPoint (129145.00000, -30060.80000), + db::DPoint (129145.00000, -28769.50000), + db::DPoint (129159.50000, -28754.90000), // this is a very short edge <-- from here. + db::DPoint (129159.60000, -28754.80000), // <-- to here. + db::DPoint (129159.50000, -28754.70000), + db::DPoint (129366.32200, -28547.90000), + db::DPoint (130958.54600, -26955.84600), + db::DPoint (131046.25000, -27043.55000), + db::DPoint (130152.15000, -27937.65000), + db::DPoint (130152.15000, -30060.80000) + }; + + db::DPolygon poly; + poly.assign_hull (contour + 0, contour + sizeof (contour) / sizeof (contour[0])); + + db::plc::TriangulationParameters param; + param.min_b = 1.0; + param.max_area = 100000.0; + param.min_length = 0.002; + + db::plc::Graph plc; + TestableTriangulation tri (&plc); + tri.triangulate (poly, param); + + EXPECT_EQ (tri.check (false), true); + + // for debugging: + // tri.dump ("debug.gds"); + + for (auto t = plc.begin (); t != plc.end (); ++t) { + EXPECT_LE (t->area (), param.max_area); + EXPECT_GE (t->b (), param.min_b); + } + + EXPECT_GT (plc.num_polygons (), size_t (540)); + EXPECT_LT (plc.num_polygons (), size_t (560)); +} + +TEST(triangulate_thin) +{ + db::DPoint contour[] = { + db::DPoint (18790, 58090), + db::DPoint (18790, 58940), + db::DPoint (29290, 58940), + db::DPoint (29290, 58090) + }; + + db::DPoint hole[] = { + db::DPoint (18791, 58091), + db::DPoint (29289, 58091), + db::DPoint (29289, 58939), + db::DPoint (18791, 58939) + }; + + db::DPolygon poly; + poly.assign_hull (contour + 0, contour + sizeof (contour) / sizeof (contour[0])); + poly.insert_hole (hole + 0, hole + sizeof (hole) / sizeof (hole[0])); + + double dbu = 0.001; + + db::plc::TriangulationParameters param; + param.min_b = 0.5; + param.max_area = 0.0; + param.min_length = 2 * dbu; + + db::plc::Graph plc; + TestableTriangulation tri (&plc); + 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 = plc.begin (); t != plc.end (); ++t) { + EXPECT_GE (t->b (), param.min_b); + } + + EXPECT_GT (plc.num_polygons (), size_t (13000)); + EXPECT_LT (plc.num_polygons (), 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::plc::TriangulationParameters param; + param.min_b = 0.5; + param.max_area = 5000.0 * dbu * dbu; + + db::plc::Graph plc; + TestableTriangulation tri (&plc); + 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 = plc.begin (); t != plc.end (); ++t) { + EXPECT_LE (t->area (), param.max_area); + EXPECT_GE (t->b (), param.min_b); + } + + EXPECT_GT (plc.num_polygons (), size_t (128000)); + EXPECT_LT (plc.num_polygons (), size_t (132000)); +} + +TEST(triangulate_with_vertexes) +{ + db::Point contour[] = { + db::Point (0, 0), + db::Point (0, 100), + db::Point (1000, 100), + db::Point (1000, 0) + }; + + db::Polygon poly; + poly.assign_hull (contour + 0, contour + sizeof (contour) / sizeof (contour[0])); + + double dbu = 0.001; + + db::plc::TriangulationParameters param; + param.min_b = 0.0; + param.max_area = 0.0; + + std::vector vertexes; + + db::plc::Graph plc; + TestableTriangulation tri (&plc); + db::CplxTrans trans = db::DCplxTrans (dbu) * db::CplxTrans (db::Trans (db::Point () - poly.box ().center ())); + tri.triangulate (poly, param, trans); + + EXPECT_EQ (plc.to_string (), "((-0.5, -0.05), (-0.5, 0.05), (0.5, 0.05)), ((0.5, -0.05), (-0.5, -0.05), (0.5, 0.05))"); + + vertexes.clear (); + + // outside vertexes are ignored, but lead to a different triangulation + vertexes.push_back (db::Point (50, 150)); + tri.triangulate (poly, vertexes, param, trans); + + EXPECT_EQ (plc.to_string (), "((-0.5, -0.05), (-0.133333333333, 0.05), (0.5, -0.05)), ((0.5, 0.05), (0.5, -0.05), (-0.133333333333, 0.05)), ((-0.133333333333, 0.05), (-0.5, -0.05), (-0.5, 0.05))"); + + for (auto v = vertexes.begin (); v != vertexes.end (); ++v) { + auto *vp = tri.find_vertex_for_point (trans * *v); + EXPECT_EQ (vp, 0); + } + + vertexes.clear (); + vertexes.push_back (db::Point (50, 50)); + tri.triangulate (poly, vertexes, param, trans); + + EXPECT_EQ (plc.to_string (), "((-0.45, 0), (-0.5, -0.05), (-0.5, 0.05)), ((0.5, 0.05), (-0.45, 0), (-0.5, 0.05)), ((-0.45, 0), (0.5, -0.05), (-0.5, -0.05)), ((-0.45, 0), (0.5, 0.05), (0.5, -0.05))"); + + for (auto v = vertexes.begin (); v != vertexes.end (); ++v) { + auto *vp = tri.find_vertex_for_point (trans * *v); + if (! vp) { + tl::warn << "Vertex not present in output: " << v->to_string (); + EXPECT_EQ (1, 0); + } + } + + // aggressive triangulation + param.min_b = 1.0; + param.max_area = 20 * 20 * dbu * dbu; + + tri.triangulate (poly, vertexes, param, trans); + + EXPECT_GT (plc.num_polygons (), size_t (380)); + EXPECT_LT (plc.num_polygons (), size_t (400)); + + for (auto t = plc.begin (); t != plc.end (); ++t) { + EXPECT_LE (t->area (), param.max_area); + EXPECT_GE (t->b (), param.min_b); + } + + for (auto v = vertexes.begin (); v != vertexes.end (); ++v) { + auto *vp = tri.find_vertex_for_point (trans * *v); + if (! vp) { + tl::warn << "Vertex not present in output: " << v->to_string (); + EXPECT_EQ (1, 0); + } + } +} diff --git a/src/db/unit_tests/dbTrianglesTests.cc b/src/db/unit_tests/dbTrianglesTests.cc index 138248694..cd5779f2d 100644 --- a/src/db/unit_tests/dbTrianglesTests.cc +++ b/src/db/unit_tests/dbTrianglesTests.cc @@ -461,7 +461,7 @@ TEST(ensure_edge) EXPECT_EQ (tl::to_string (area_in), "0.25"); } -bool safe_inside (const db::DBox &b1, const db::DBox &b2) +static bool safe_inside (const db::DBox &b1, const db::DBox &b2) { typedef db::coord_traits ct; diff --git a/src/db/unit_tests/unit_tests.pro b/src/db/unit_tests/unit_tests.pro index 047e46ea6..ad9f9e382 100644 --- a/src/db/unit_tests/unit_tests.pro +++ b/src/db/unit_tests/unit_tests.pro @@ -12,7 +12,8 @@ SOURCES = \ dbFillToolTests.cc \ dbLogTests.cc \ dbObjectWithPropertiesTests.cc \ - dbPolygonGraphTests.cc \ + dbPLCGraphTest.cc \ + dbPLCTriangulationTests.cc \ dbPolygonNeighborhoodTests.cc \ dbPropertiesFilterTests.cc \ dbQuadTreeTests.cc \