From d7193e972cbf33d384610626e8bfe535c6898c7b Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Sat, 12 Apr 2025 23:59:43 +0200 Subject: [PATCH] WIP --- src/db/db/dbTriangles.h | 19 +- src/db/unit_tests/dbTrianglesTests.cc | 296 +++++++++++++++----------- 2 files changed, 179 insertions(+), 136 deletions(-) diff --git a/src/db/db/dbTriangles.h b/src/db/db/dbTriangles.h index 31da3f48b..3bb34c87d 100644 --- a/src/db/db/dbTriangles.h +++ b/src/db/db/dbTriangles.h @@ -170,6 +170,17 @@ public: void triangulate (const db::DPolygon &poly, const TriangulateParameters ¶meters, const db::DCplxTrans &trans = db::DCplxTrans ()); void triangulate (const db::DPolygon &poly, const std::vector &vertexes, const TriangulateParameters ¶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. + */ + db::Vertex *insert_point (const db::DPoint &point, std::list > *new_triangles = 0); + /** * @brief Statistics: number of flips (fixing) */ @@ -214,14 +225,6 @@ protected: */ 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. - */ - db::Vertex *insert_point (const db::DPoint &point, std::list > *new_triangles = 0); - /** * @brief Inserts a new vertex as the given point * diff --git a/src/db/unit_tests/dbTrianglesTests.cc b/src/db/unit_tests/dbTrianglesTests.cc index fe09c3817..138248694 100644 --- a/src/db/unit_tests/dbTrianglesTests.cc +++ b/src/db/unit_tests/dbTrianglesTests.cc @@ -1081,7 +1081,7 @@ TEST(triangulate_with_vertexes) struct SortAngleAndEdgesByEdgeLength { - typedef std::list > angle_and_edges_list; + typedef std::list > angle_and_edges_list; bool operator() (const angle_and_edges_list::iterator &a, const angle_and_edges_list::iterator &b) const { @@ -1095,6 +1095,26 @@ struct SortAngleAndEdgesByEdgeLength } }; +// TODO: move to some generic header +template +struct less_compare_func +{ + bool operator() (const T &a, const T &b) const + { + return a.less (b); + } +}; + +// TODO: move to some generic header +template +struct equal_compare_func +{ + bool operator() (const T &a, const T &b) const + { + return a.equal (b); + } +}; + struct ConcaveCorner { ConcaveCorner () @@ -1113,13 +1133,13 @@ struct ConcaveCorner db::TriangleEdge *incoming, *outgoing; }; -db::TriangleEdge *find_outgoing_segment (db::Vertex *vertex, db::TriangleEdge *incoming, bool &is_concave) +db::TriangleEdge *find_outgoing_segment (db::Vertex *vertex, db::TriangleEdge *incoming, int &vp_max_sign) { db::Vertex *vfrom = incoming->other (vertex); db::DEdge e1 (*vfrom, *vertex); double vp_max = 0.0; - int vp_max_sign = 0; + vp_max_sign = 0; db::TriangleEdge *outgoing = 0; // Look for the outgoing edge. We pick the one which bends "most", favoring @@ -1147,8 +1167,6 @@ db::TriangleEdge *find_outgoing_segment (db::Vertex *vertex, db::TriangleEdge *i } - is_concave = (vp_max_sign > 0); - tl_assert (outgoing != 0); return outgoing; } @@ -1188,10 +1206,10 @@ void collect_concave_vertexes (db::Triangles &tris, std::vector & db::TriangleEdge *prev_segment = segment; - bool is_concave = false; - segment = find_outgoing_segment (vto, prev_segment, is_concave); + int vp_sign = 0; + segment = find_outgoing_segment (vto, prev_segment, vp_sign); - if (is_concave) { + if (vp_sign > 0) { concave_vertexes.push_back (ConcaveCorner (vto, prev_segment, segment)); } @@ -1202,122 +1220,94 @@ void collect_concave_vertexes (db::Triangles &tris, std::vector & } } -// Hertel-Mehlhorn :) -TEST(JoinTriangles) +std::pair +search_crossing_with_next_segment (const db::Vertex *v0, const db::DVector &direction) { -#if 1 - db::Point contour[] = { - db::Point (0, 0), - db::Point (0, 100), - db::Point (1000, 100), - db::Point (1000, 500), - db::Point (1100, 500), - db::Point (1100, 100), - db::Point (2100, 100), - db::Point (2100, 0) - }; -#else + auto vtri = v0->triangles (); // TODO: slow? + std::vector nvv, nvv_next; - db::Point contour[] = { - db::Point (0, 0), - db::Point (0, 100), - db::Point (1000, 100), - db::Point (1000, 500), - db::Point (1100, 500), - db::Point (1100, 100), - db::Point (2100, 100), - db::Point (2100, -1000), - db::Point (1050, -1000), - db::Point (1050, 0) - }; -#endif + for (auto it = vtri.begin (); it != vtri.end (); ++it) { - db::Polygon poly; - poly.assign_hull (contour + 0, contour + sizeof (contour) / sizeof (contour[0])); + // Search for a segment in the direction perpendicular to the edge + nvv.clear (); + nvv.push_back (v0); + const db::Triangle *t = *it; - // @@@ don't to anything if already convex + while (! nvv.empty ()) { - double dbu = 0.001; + nvv_next.clear (); - db::Triangles::TriangulateParameters param; - param.min_b = 0.0; + for (auto iv = nvv.begin (); iv != nvv.end (); ++iv) { - TestableTriangles tri; - db::CplxTrans trans = db::DCplxTrans (dbu) * db::CplxTrans (db::Trans (db::Point () - poly.box ().center ())); - trans = db::CplxTrans (dbu); // @@@ - tri.triangulate (poly, param, trans); + const db::Vertex *v = *iv; + const db::TriangleEdge *oe = t->opposite (v); + const db::Triangle *tt = oe->other (t); + const db::Vertex *v1 = oe->v1 (); + const db::Vertex *v2 = oe->v2 (); + if (db::sprod_sign (*v2 - *v, direction) >= 0 && db::sprod_sign (*v1 - *v, direction) >= 0 && + db::vprod_sign (*v2 - *v, direction) * db::vprod_sign (*v1 - *v, direction) < 0) { + + // this triangle covers the normal vector of e1 -> stop here or continue searching in that direction + if (oe->is_segment ()) { + auto cp = oe->edge ().cut_point (db::DEdge (*v0, *v0 + direction)); + if (cp.first) { + return std::make_pair (true, cp.second); + } + } else { + // continue searching in that direction + nvv_next.push_back (v1); + nvv_next.push_back (v2); + t = tt; + } + + break; + + } + + } + + nvv.swap (nvv_next); + + } + + } + + return std::make_pair (false, db::DPoint ()); +} + +void +hertel_mehlhorn_decomposition (db::Triangles &tri, bool diagonals_only, bool no_collinear_edges, std::list &polygons) +{ std::vector concave_vertexes; collect_concave_vertexes (tri, concave_vertexes); // @@@ return if no concave corners - // @@@ sort convex vertexes + // @@@ sort concave vertexes std::vector new_points; - // Cut off pieces from convex corners by creating connections to points perpendicular - // to the incoming and outgoing edges + if (! diagonals_only) { - for (auto cc = concave_vertexes.begin (); cc != concave_vertexes.end (); ++cc) { + // Create internal segments cutting off pieces orthogonal to the edges + // connecting the concave vertex. - auto vtri = cc->corner->triangles (); // @@@ slow? + for (auto cc = concave_vertexes.begin (); cc != concave_vertexes.end (); ++cc) { - std::vector nvv, nvv_next; + for (unsigned int ei = 0; ei < 2; ++ei) { - for (unsigned int ei = 0; ei < 2; ++ei) { - - db::DEdge ee; - db::Vertex *v0 = cc->corner; - if (ei == 0) { - ee = db::DEdge (*cc->incoming->other (v0), *v0); - } else { - ee = db::DEdge (*v0, *cc->outgoing->other (v0)); - } - - for (auto it = vtri.begin (); it != vtri.end (); ++it) { - - // Search for a segment in the direction perpendicular to the edge - nvv.clear (); - nvv.push_back (v0); - db::Triangle *t = *it; - - while (! nvv.empty ()) { - - nvv_next.clear (); - - for (auto iv = nvv.begin (); iv != nvv.end (); ++iv) { - - db::Vertex *v = *iv; - db::TriangleEdge *oe = t->opposite (v); - db::Triangle *tt = oe->other (t); - db::Vertex *v1 = oe->v1 (); - db::Vertex *v2 = oe->v2 (); - - if (db::vprod_sign (*v2 - *v, ee.d ()) >= 0 && db::vprod_sign (*v1 - *v, ee.d ()) >= 0 && - db::sprod_sign (*v2 - *v, ee.d ()) * db::sprod_sign (*v1 - *v, ee.d ()) < 0) { - - // this triangle covers the normal vector of e1 -> stop here or continue searching in that direction - if (oe->is_segment ()) { - auto cp = oe->edge ().cut_point (db::DEdge (*v0, *v0 + db::DVector (ee.dy (), -ee.dx ()))); - if (cp.first) { - new_points.push_back (cp.second); - } - } else { - // continue searching in that direction - nvv_next.push_back (v1); - nvv_next.push_back (v2); - t = tt; - } - - break; - - } - - } - - nvv.swap (nvv_next); + db::DEdge ee; + const db::Vertex *v0 = cc->corner; + if (ei == 0) { + ee = db::DEdge (*cc->incoming->other (v0), *v0); + } else { + ee = db::DEdge (*v0, *cc->outgoing->other (v0)); + } + auto cp = search_crossing_with_next_segment (v0, db::DVector (ee.dy (), -ee.dx ())); + if (cp.first) { + new_points.push_back (cp.second); } } @@ -1326,16 +1316,21 @@ TEST(JoinTriangles) } - // Insert the new points and make connections - - // @@@ TODO: what to do in case of equal new_points? - // -> sort, remove duplicates - for (auto p = new_points.begin (); p != new_points.end (); ++p) { - tri.insert_point (*p); - } + // eliminate duplicates and put the new points in some order if (! new_points.empty ()) { + + std::sort (new_points.begin (), new_points.end (), less_compare_func ()); + new_points.erase (std::unique (new_points.begin (), new_points.end (), equal_compare_func ()), new_points.end ()); + + // Insert the new points and make connections + for (auto p = new_points.begin (); p != new_points.end (); ++p) { + tri.insert_point (*p); + } + + // As the insertion invalidates the edges, we need to collect the concave vertexes again collect_concave_vertexes (tri, concave_vertexes); + } // Collect essential edges @@ -1345,27 +1340,25 @@ TEST(JoinTriangles) std::unordered_set essential_edges; - typedef std::list > angles_and_edges_list; + typedef std::list > angles_and_edges_list; angles_and_edges_list angles_and_edges; std::vector sorted_edges; for (auto cc = concave_vertexes.begin (); cc != concave_vertexes.end (); ++cc) { - std::cout << "@@@ cc=" << cc->corner->to_string () << ", in: " << cc->incoming->to_string () << ", out: " << cc->outgoing->to_string () << std::endl; // @@@ - angles_and_edges.clear (); - db::Vertex *v0 = cc->corner; + const db::Vertex *v0 = cc->corner; - db::TriangleEdge *e = cc->incoming; + const db::TriangleEdge *e = cc->incoming; while (e) { - db::Triangle *t = e->v2 () == v0 ? e->right () : e->left (); + const db::Triangle *t = e->v2 () == v0 ? e->right () : e->left (); tl_assert (t != 0); // @@@ make a method of triangle - db::TriangleEdge *en = 0; + const db::TriangleEdge *en = 0; for (unsigned int i = 0; i < 3; ++i) { - db::TriangleEdge *ee = t->edge (i); + const db::TriangleEdge *ee = t->edge (i); if (ee != e && (ee->v1 () == v0 || ee->v2 () == v0)) { en = ee; break; @@ -1379,7 +1372,6 @@ TEST(JoinTriangles) e = (en == cc->outgoing) ? 0 : en; angles_and_edges.push_back (std::make_pair (angle, e)); - std::cout << "@@@ [a,e] =" << angle << "," << (e ? e->to_string () : std::string ("ENDL")) << std::endl; // @@@ } @@ -1398,13 +1390,10 @@ TEST(JoinTriangles) angles_and_edges_list::iterator ii = *i; angles_and_edges_list::iterator iin = ii; ++iin; - std::cout << "@@@ checking [a,e] =" << ii->first << "," << iin->first << "," << ii->second->to_string () << std::endl; // @@@ - if (ii->first + iin->first < M_PI - db::epsilon) { - printf("@@@ %.12g, %.12g\n", ii->first + iin->first, M_PI); // @@@ + if (ii->first + iin->first < (no_collinear_edges ? M_PI - db::epsilon : M_PI + db::epsilon)) { // not an essential edge -> remove iin->first += ii->first; angles_and_edges.erase (ii); - std::cout << "@@@ -> not essential" << std::endl; // @@@ } } @@ -1416,8 +1405,6 @@ TEST(JoinTriangles) // Combine triangles, but don't cross essential edges - db::Region result; - std::unordered_set left_triangles; for (auto it = tri.begin (); it != tri.end (); ++it) { left_triangles.insert (it.operator-> ()); @@ -1483,11 +1470,64 @@ TEST(JoinTriangles) } while (v != v0); - db::DPolygon poly; - poly.assign_hull (polygon_points.begin (), polygon_points.end ()); - result.insert (trans.inverted () * poly); + polygons.push_back (db::DPolygon ()); + polygons.back ().assign_hull (polygon_points.begin (), polygon_points.end ()); } +} + + +// Hertel-Mehlhorn :) +TEST(JoinTriangles) +{ +#if 0 + db::Point contour[] = { + db::Point (0, 0), + db::Point (0, 100), + db::Point (1000, 100), + db::Point (1000, 500), + db::Point (1100, 500), + db::Point (1100, 100), + db::Point (2100, 100), + db::Point (2100, 0) + }; +#else + + db::Point contour[] = { + db::Point (0, 0), + db::Point (0, 100), + db::Point (1000, 100), + db::Point (1000, 500), + db::Point (1100, 500), + db::Point (1100, 100), + db::Point (2100, 100), + db::Point (2100, -1000), + db::Point (150, -1000), + db::Point (150, 0) + }; +#endif + + db::Polygon poly; + poly.assign_hull (contour + 0, contour + sizeof (contour) / sizeof (contour[0])); + + // @@@ don't to anything if already convex + + double dbu = 0.001; + + db::Triangles::TriangulateParameters param; + param.min_b = 0.0; + + TestableTriangles tri; + db::CplxTrans trans = db::CplxTrans (dbu); + tri.triangulate (poly, param, trans); + + std::list polygons; + hertel_mehlhorn_decomposition (tri, false, true, polygons); + + db::Region result; + for (auto p = polygons.begin (); p != polygons.end (); ++p) { + result.insert (trans.inverted () * *p); + } // @@@ tri.dump ("debugt.gds");