From f136fdcde6d799d8328c3d2b5558437bc250cae9 Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Sun, 16 Mar 2025 00:52:57 +0100 Subject: [PATCH 01/17] WIP --- src/db/db/db.pro | 2 + src/db/db/dbQuadTree.cc | 24 ++ src/db/db/dbQuadTree.h | 494 ++++++++++++++++++++++++++ src/db/db/dbTriangle.cc | 20 +- src/db/db/dbTriangles.cc | 26 +- src/db/db/dbTriangles.h | 1 + src/db/unit_tests/dbQuadTreeTests.cc | 77 ++++ src/db/unit_tests/dbTriangleTests.cc | 36 ++ src/db/unit_tests/dbTrianglesTests.cc | 37 ++ src/db/unit_tests/unit_tests.pro | 1 + 10 files changed, 710 insertions(+), 8 deletions(-) create mode 100644 src/db/db/dbQuadTree.cc create mode 100644 src/db/db/dbQuadTree.h create mode 100644 src/db/unit_tests/dbQuadTreeTests.cc diff --git a/src/db/db/db.pro b/src/db/db/db.pro index 499046cc7..bbae9671a 100644 --- a/src/db/db/db.pro +++ b/src/db/db/db.pro @@ -81,6 +81,7 @@ SOURCES = \ dbPolygonGenerators.cc \ dbPropertiesFilter.cc \ dbPropertiesRepository.cc \ + dbQuadTree.cc \ dbReader.cc \ dbRecursiveInstanceIterator.cc \ dbRecursiveShapeIterator.cc \ @@ -319,6 +320,7 @@ HEADERS = \ dbPropertiesFilter.h \ dbPropertiesRepository.h \ dbPropertyConstraint.h \ + dbQuadTree.h \ dbReader.h \ dbRecursiveInstanceIterator.h \ dbRecursiveShapeIterator.h \ diff --git a/src/db/db/dbQuadTree.cc b/src/db/db/dbQuadTree.cc new file mode 100644 index 000000000..7e9992fbc --- /dev/null +++ b/src/db/db/dbQuadTree.cc @@ -0,0 +1,24 @@ + +/* + + 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 "dbQuadTree.h" + diff --git a/src/db/db/dbQuadTree.h b/src/db/db/dbQuadTree.h new file mode 100644 index 000000000..e3ab61ef2 --- /dev/null +++ b/src/db/db/dbQuadTree.h @@ -0,0 +1,494 @@ + +/* + + 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_dbQuadTree +#define HDR_dbQuadTree + +#include "dbBox.h" +#include + +namespace db +{ + +template +class quad_tree_node +{ +public: + typedef typename T::coord_type coord_type; + typedef db::box box_type; + typedef db::point point_type; + typedef db::vector vector_type; + typedef std::vector objects_vector; + + quad_tree_node (const point_type ¢er) + : m_split (false), m_center (center) + { + for (unsigned int i = 0; i < 4; ++i) { + m_q [i] = 0; + } + } + + ~quad_tree_node () + { + for (unsigned int i = 0; i < 4; ++i) { + delete m_q [i]; + m_q [i] = 0; + } + } + + const point_type ¢er () const + { + return m_center; + } + + void insert_top (const T &value, const box_type &total_box) + { + insert (value, propose_ucenter (total_box)); + } + + bool remove (const T &value) + { + box_type b = BC () (value); + + if (! m_split || b.contains (m_center)) { + for (auto i = m_objects.begin (); i != m_objects.end (); ++i) { + if (*i == value) { + m_objects.erase (i); + return true; + } + } + } + + for (unsigned int i = 0; i < 4; ++i) { + if (m_q[i] && b.inside (m_q[i]->box (m_center))) { + return b.remove (value); + } + } + + return false; + } + + const objects_vector &objects () const + { + return m_objects; + } + + box_type q_box (unsigned int n) const + { + if (m_q[n]) { + return m_q[n]->box (m_center); + } else { + return box_type (); + } + } + + quad_tree_node *node (unsigned int n) const + { + return m_q [n]; + } + + bool empty () const + { + if (m_objects.empty ()) { + for (unsigned int n = 0; n < 4; ++n) { + if (m_q[n] && ! m_q[n]->empty ()) { + return false; + } + } + return true; + } else { + return false; + } + } + + size_t size () const + { + size_t count = m_objects.size (); + for (unsigned int n = 0; n < 4; ++n) { + if (m_q[n]) { + count += m_q[n]->size (); + } + } + return count; + } + +private: + bool m_split; + point_type m_center; + quad_tree_node *m_q [4]; + objects_vector m_objects; + + box_type box (const point_type &ucenter) const + { + return box_type (ucenter, ucenter - (ucenter - m_center) * 2.0); + } + + box_type q (unsigned int n, const point_type &ucenter) const + { + // NOTE: with this definition the opposite quad index is 3 - n + + vector_type vx (std::abs (ucenter.x () - m_center.x ()), 0); + vector_type vy (0, std::abs (ucenter.y () - m_center.y ())); + switch (n) { + case 0: + return box_type (m_center - vx - vy, m_center); + case 1: + return box_type (m_center - vy, m_center + vx); + case 2: + return box_type (m_center - vx, m_center + vy); + default: + return box_type (m_center, m_center + vx + vy); + } + } + + void split (const point_type &ucenter) + { + m_split = true; + + objects_vector ov; + ov.swap (m_objects); + + for (auto o = ov.begin (); o != ov.end (); ++o) { + insert (*o, ucenter); + } + } + + void insert (const T &value, const point_type &ucenter) + { + if (! m_split && m_objects.size () + 1 < thr) { + + m_objects.push_back (value); + + } else { + + if (! m_split) { + split (ucenter); + } + + box_type b = BC () (value); + // @@@ should exclude m_center on box + if (b.contains (m_center)) { + m_objects.push_back (value); + return; + } + + for (unsigned int i = 0; i < 4; ++i) { + box_type bq = q (i, ucenter); + if (b.inside (bq)) { + if (! m_q[i]) { + m_q[i] = new quad_tree_node (bq.center ()); + } + m_q[i]->insert (value, m_center); + return; + } + } + + for (unsigned int i = 0; i < 4; ++i) { + if (m_q[i]) { + grow (m_center - (m_center - m_q[i]->center ()) * 2.0); + insert (value, ucenter); + return; + } + } + + tl_assert (false); + + } + } + + void grow (const point_type &ucenter) + { + for (unsigned int i = 0; i < 4; ++i) { + if (m_q[i]) { + quad_tree_node *n = m_q[i]; + m_q[i] = new quad_tree_node (q (i, ucenter).center ()); + m_q[i]->m_q[3 - i] = n; + } + } + } + + point_type propose_ucenter (const box_type &total_box) const + { + for (unsigned int i = 0; i < 4; ++i) { + if (m_q[i]) { + return m_center - (m_center - m_q[i]->center ()) * 2.0; + } + } + + coord_type dx = std::max (std::abs (total_box.left () - m_center.x ()), std::abs (total_box.right () - m_center.y ())); + coord_type dy = std::max (std::abs (total_box.bottom () - m_center.y ()), std::abs (total_box.top () - m_center.y ())); + return m_center - vector_type (dx, dy); + } +}; + +template +class quad_tree_iterator +{ +public: + typedef quad_tree_node quad_tree_node_type; + typedef typename T::coord_type coord_type; + typedef db::box box_type; + + quad_tree_iterator () + : m_s (), m_i (0) + { + // .. nothing yet .. + } + + quad_tree_iterator (const quad_tree_node_type *root, const S &s) + : m_s (s), m_i (0) + { + m_stack.push_back (std::make_pair (root, -1)); + validate (); + } + + bool at_end () const + { + return m_stack.empty (); + } + + quad_tree_iterator &operator++ () + { + ++m_i; + validate (); + return *this; + } + + const T &operator* () const + { + return m_stack.back ().first->objects () [m_i]; + } + + const T *operator-> () const + { + return (m_stack.back ().first->objects ().begin () + m_i).operator-> (); + } + +public: + S m_s; + std::vector > m_stack; + size_t m_i; + + void validate () + { + auto s = m_stack.end (); + tl_assert (s != m_stack.begin ()); + --s; + + const quad_tree_node_type *qn = s->first; + while (m_i < s->first->objects ().size () && ! m_s.select (s->first->objects () [m_i])) { + ++m_i; + } + + if (m_i < qn->objects ().size ()) { + return; + } + + m_i = 0; + + for (unsigned int n = 0; n < 4; ++n) { + box_type bq = qn->q_box (n); + if (! bq.empty () && m_s.select_quad (bq)) { + m_stack.back ().second = n; + m_stack.push_back (std::make_pair (qn->node (n), -1)); + validate (); + return; + } + } + + while (! m_stack.empty ()) { + + m_stack.pop_back (); + + int &n = m_stack.back ().second; + while (++n < 4) { + box_type bq = qn->q_box (n); + if (! bq.empty () && m_s.select_quad (bq)) { + m_stack.push_back (std::make_pair (qn->node (n), -1)); + validate (); + return; + } + } + + } + } +}; + +template +class quad_tree_always_sel +{ +public: + typedef typename T::coord_type coord_type; + typedef db::box box_type; + + bool select (const T &) const + { + return true; + } + + bool select_quad (const box_type &) const + { + return true; + } +}; + +template +class quad_tree_touching_sel +{ +public: + typedef typename T::coord_type coord_type; + typedef db::box box_type; + + quad_tree_touching_sel () + { + // .. nothing yet .. + } + + quad_tree_touching_sel (const box_type &box) + : m_box (box) + { + // .. nothing yet .. + } + + bool select (const T &value) const + { + return select_quad (BC () (value)); + } + + bool select_quad (const box_type &box) const + { + return m_box.touches (box); + } + +private: + box_type m_box; +}; + +template +class quad_tree_overlapping_sel +{ +public: + typedef typename T::coord_type coord_type; + typedef db::box box_type; + + quad_tree_overlapping_sel () + { + // .. nothing yet .. + } + + quad_tree_overlapping_sel (const box_type &box) + : m_box (box) + { + // .. nothing yet .. + } + + bool select (const T &value) const + { + return select_quad (BC () (value)); + } + + bool select_quad (const box_type &box) const + { + return m_box.overlaps (box); + } + +private: + box_type m_box; +}; + +template +class quad_tree +{ +public: + typedef quad_tree_node quad_tree_node_type; + typedef quad_tree_iterator > quad_tree_flat_iterator; + typedef quad_tree_iterator > quad_tree_touching_iterator; + typedef quad_tree_iterator > quad_tree_overlapping_iterator; + typedef typename T::coord_type coord_type; + typedef db::box box_type; + typedef db::point point_type; + typedef db::vector vector_type; + typedef std::vector objects_vector; + + quad_tree () + : m_root (point_type ()) + { + // .. nothing yet .. + } + + bool empty () const + { + return m_root.empty (); + } + + size_t size () const + { + return m_root.size (); + } + + void insert (const T &value) + { + box_type b = BC () (value); + if (b.empty ()) { + return; + } + + m_total_box += b; + m_root.insert_top (value, m_total_box); + } + + void erase (const T &value) + { + m_root.remove (value); + } + + quad_tree_flat_iterator begin () const + { + return quad_tree_flat_iterator (&m_root, quad_tree_always_sel ()); + } + + quad_tree_overlapping_iterator begin_overlapping (const box_type &box) const + { + if (m_total_box.overlaps (box)) { + return quad_tree_overlapping_iterator (&m_root, quad_tree_overlapping_sel (box)); + } else { + return quad_tree_overlapping_iterator (); + } + } + + quad_tree_touching_iterator begin_touching (const box_type &box) const + { + if (m_total_box.touches (box)) { + return quad_tree_touching_iterator (&m_root, quad_tree_touching_sel (box)); + } else { + return quad_tree_touching_iterator (); + } + } + +private: + box_type m_total_box; + quad_tree_node_type m_root; +}; + +} + +#endif diff --git a/src/db/db/dbTriangle.cc b/src/db/db/dbTriangle.cc index 1016ad3ec..0fa3018c9 100644 --- a/src/db/db/dbTriangle.cc +++ b/src/db/db/dbTriangle.cc @@ -507,15 +507,23 @@ Triangle::common_edge (const Triangle *other) const int Triangle::contains (const db::DPoint &point) const { + // the relative distance of the point from the edge + double alpha = db::epsilon; + + double d = db::vprod (*mp_v[2] - *mp_v[0], *mp_v[1] - *mp_v[0]); + if (fabs (d) < db::epsilon * db::epsilon) { + return -1; + } + int res = 1; - const Vertex *vl = mp_v[2];; + const Vertex *vl = mp_v[2]; for (int i = 0; i < 3; ++i) { - const Vertex *v = mp_v[i];; - int s = db::DEdge (*vl, *v).side_of (point); - if (s == 0) { - res = 0; - } else if (s > 0) { + const Vertex *v = mp_v[i]; + double n = db::vprod (point - *vl, *v - *vl) / d; + if (n < -alpha) { return -1; + } else if (n < alpha) { + res = 0; } vl = v; } diff --git a/src/db/db/dbTriangles.cc b/src/db/db/dbTriangles.cc index e6aff5ee5..5a0475a5d 100644 --- a/src/db/db/dbTriangles.cc +++ b/src/db/db/dbTriangles.cc @@ -355,6 +355,7 @@ Triangles::insert_point (db::DCoord x, db::DCoord y, std::list > *new_triangles) { + // @@@ printf("@@@ insert %d\n", (int)mp_triangles.size ()); fflush(stdout); std::vector tris = find_triangle_for_point (*vertex); // the new vertex is outside the domain @@ -412,7 +413,24 @@ Triangles::find_closest_edge (const db::DPoint &p, db::Vertex *vstart, bool insi { if (!vstart) { - if (! mp_triangles.empty ()) { + if (m_is_constrained && ! m_initial_segments.empty ()) { + + // Use the closest initial segment for the starting point + // TODO: m_initial_segments should be some search tree + + double d = -1; + for (auto i = m_initial_segments.begin (); i != m_initial_segments.end (); ++i) { + if (db::sprod_sign (p - i->first.p1 (), i->first.d ()) >= 0 && + db::sprod_sign (p - i->first.p2 (), i->first.d ()) < 0) { + double ds = std::abs (i->first.distance (p)); + if (d < 0.0 || ds < d) { + vstart = i->second; + d = ds; + } + } + } + + } else if (! mp_triangles.empty ()) { unsigned int ls = 0; size_t n = m_vertex_heap.size (); @@ -1319,6 +1337,7 @@ void Triangles::constrain (const std::vector > &contours) { tl_assert (! m_is_constrained); + m_initial_segments.clear (); std::vector > > resolved_edges; @@ -1329,8 +1348,10 @@ Triangles::constrain (const std::vector > &contours) if (vv == c->end ()) { vv = c->begin (); } - resolved_edges.push_back (std::make_pair (db::DEdge (**v, **vv), std::vector ())); + db::DEdge e (**v, **vv); + resolved_edges.push_back (std::make_pair (e, std::vector ())); resolved_edges.back ().second = ensure_edge (*v, *vv); + m_initial_segments.push_back (std::make_pair (e, *v)); } } @@ -1606,6 +1627,7 @@ Triangles::refine (const TriangulateParameters ¶meters) auto cr = (*t)->circumcircle(); auto center = cr.first; + // printf("@@@ %s -- %s\n", (*t)->to_string().c_str(), center.to_string ().c_str()); fflush(stdout); if ((*t)->contains (center) >= 0) { if (tl::verbosity () >= parameters.base_verbosity + 20) { diff --git a/src/db/db/dbTriangles.h b/src/db/db/dbTriangles.h index f86f77c1c..f15b2dc7a 100644 --- a/src/db/db/dbTriangles.h +++ b/src/db/db/dbTriangles.h @@ -311,6 +311,7 @@ private: std::vector m_returned_edges; tl::stable_vector m_vertex_heap; bool m_is_constrained; + std::vector > m_initial_segments; size_t m_level; size_t m_id; size_t m_flips, m_hops; diff --git a/src/db/unit_tests/dbQuadTreeTests.cc b/src/db/unit_tests/dbQuadTreeTests.cc new file mode 100644 index 000000000..fd74e164b --- /dev/null +++ b/src/db/unit_tests/dbQuadTreeTests.cc @@ -0,0 +1,77 @@ + +/* + + 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 "dbQuadTree.h" +#include "dbBoxConvert.h" +#include "tlUnitTest.h" +#include "tlString.h" + +typedef db::quad_tree, size_t (1)> my_quad_tree; + +std::string find_touching (const my_quad_tree &qt, const db::DBox &box) +{ + std::vector v; + auto i = qt.begin_touching (box); + while (! i.at_end ()) { + v.push_back (i->to_string ()); + ++i; + } + std::sort (v.begin (), v.end ()); + return tl::join (v, "/"); +} + +std::string find_overlapping (const my_quad_tree &qt, const db::DBox &box) +{ + std::vector v; + auto i = qt.begin_overlapping (box); + while (! i.at_end ()) { + v.push_back (i->to_string ()); + ++i; + } + std::sort (v.begin (), v.end ()); + return tl::join (v, "/"); +} + +TEST(basic) +{ + my_quad_tree tree; + EXPECT_EQ (tree.empty (), true); + EXPECT_EQ (tree.size (), size_t (0)); + + tree.insert (db::DBox ()); + EXPECT_EQ (tree.empty (), true); + EXPECT_EQ (tree.size (), size_t (0)); + + tree.insert (db::DBox (-1, -2, 3, 4)); + EXPECT_EQ (tree.empty (), false); + EXPECT_EQ (tree.size (), size_t (1)); + + EXPECT_EQ (find_touching (tree, db::DBox (-2, 0, -1, 0)), "..."); + + tree.insert (db::DBox (-1, -3, 3, 0)); + EXPECT_EQ (tree.empty (), false); + EXPECT_EQ (tree.size (), size_t (1)); + + EXPECT_EQ (find_touching (tree, db::DBox (-2, 0, -1, 0)), "..."); +} + diff --git a/src/db/unit_tests/dbTriangleTests.cc b/src/db/unit_tests/dbTriangleTests.cc index c93c93b8d..439eb1c65 100644 --- a/src/db/unit_tests/dbTriangleTests.cc +++ b/src/db/unit_tests/dbTriangleTests.cc @@ -302,6 +302,42 @@ TEST(Triangle_contains) } } +TEST(Triangle_contains_small) +{ + db::Vertex v1; + db::Vertex v2 (0.001, 0.002); + db::Vertex v3 (0.002, 0.001); + + TestableTriangleEdge s1 (&v1, &v2); + TestableTriangleEdge s2 (&v2, &v3); + TestableTriangleEdge s3 (&v3, &v1); + + { + db::Triangle tri (&s1, &s2, &s3); + EXPECT_EQ (tri.contains (db::DPoint (0, 0)), 0); + EXPECT_EQ (tri.contains (db::DPoint (-0.001, -0.002)), -1); + EXPECT_EQ (tri.contains (db::DPoint (0.0005, 0.001)), 0); + EXPECT_EQ (tri.contains (db::DPoint (0.0005, 0.002)), -1); + EXPECT_EQ (tri.contains (db::DPoint (0.0025, 0.001)), -1); + EXPECT_EQ (tri.contains (db::DPoint (0.001, -0.001)), -1); + EXPECT_EQ (tri.contains (db::DPoint (0.001, 0.001)), 1); + } + + s1.reverse (); + s2.reverse (); + s3.reverse (); + + { + db::Triangle tri2 (&s3, &s2, &s1); + EXPECT_EQ (tri2.contains(db::DPoint(0, 0)), 0); + EXPECT_EQ (tri2.contains(db::DPoint(0.0005, 0.001)), 0); + EXPECT_EQ (tri2.contains(db::DPoint(0.0005, 0.002)), -1); + EXPECT_EQ (tri2.contains(db::DPoint(0.0025, 0.001)), -1); + EXPECT_EQ (tri2.contains(db::DPoint(0.001, -0.001)), -1); + EXPECT_EQ (tri2.contains(db::DPoint(0.001, 0.001)), 1); + } +} + TEST(Triangle_circumcircle) { db::Vertex v1; diff --git a/src/db/unit_tests/dbTrianglesTests.cc b/src/db/unit_tests/dbTrianglesTests.cc index 0bc00893e..9b7dd021c 100644 --- a/src/db/unit_tests/dbTrianglesTests.cc +++ b/src/db/unit_tests/dbTrianglesTests.cc @@ -23,6 +23,7 @@ #include "dbTriangles.h" #include "dbWriter.h" +#include "dbRegionProcessors.h" #include "tlUnitTest.h" #include "tlStream.h" #include "tlFileUtils.h" @@ -917,3 +918,39 @@ TEST(triangulate_problematic) EXPECT_LT (tri.num_triangles (), size_t (490)); } +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::Triangles::TriangulateParameters param; + param.min_b = 0.5; + param.max_area = 0.0; + param.min_length = 2 * dbu; + + TestableTriangles tri; + db::DCplxTrans trans = db::DCplxTrans (dbu) * db::DCplxTrans (db::DTrans (db::DPoint () - poly.box ().center ())); + tri.triangulate (trans * poly, param); + + EXPECT_EQ (tri.check (false), true); + + // for debugging: + tri.dump ("debug.gds"); +} diff --git a/src/db/unit_tests/unit_tests.pro b/src/db/unit_tests/unit_tests.pro index a24390237..ecb3c7ccb 100644 --- a/src/db/unit_tests/unit_tests.pro +++ b/src/db/unit_tests/unit_tests.pro @@ -14,6 +14,7 @@ SOURCES = \ dbObjectWithPropertiesTests.cc \ dbPolygonNeighborhoodTests.cc \ dbPropertiesFilterTests.cc \ + dbQuadTreeTests.cc \ dbRecursiveInstanceIteratorTests.cc \ dbRegionCheckUtilsTests.cc \ dbTriangleTests.cc \ From 54b5d9f5d682c0831473e2687fe056c38a2ab820 Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Sun, 16 Mar 2025 13:35:51 +0100 Subject: [PATCH 02/17] WIP (quad tree) --- src/db/db/dbQuadTree.h | 156 +++++++++++++++++----- src/db/unit_tests/dbQuadTreeTests.cc | 187 ++++++++++++++++++++++++++- 2 files changed, 310 insertions(+), 33 deletions(-) diff --git a/src/db/db/dbQuadTree.h b/src/db/db/dbQuadTree.h index e3ab61ef2..353f7ddcd 100644 --- a/src/db/db/dbQuadTree.h +++ b/src/db/db/dbQuadTree.h @@ -24,6 +24,7 @@ #define HDR_dbQuadTree #include "dbBox.h" +#include "tlLog.h" #include namespace db @@ -38,6 +39,7 @@ public: typedef db::point point_type; typedef db::vector vector_type; typedef std::vector objects_vector; + typedef db::coord_traits coord_traits; quad_tree_node (const point_type ¢er) : m_split (false), m_center (center) @@ -65,23 +67,31 @@ public: insert (value, propose_ucenter (total_box)); } - bool remove (const T &value) + bool erase (const T &value) { box_type b = BC () (value); - if (! m_split || b.contains (m_center)) { + int n = quad_for (b); + + if (! m_split || n < 0) { + for (auto i = m_objects.begin (); i != m_objects.end (); ++i) { if (*i == value) { m_objects.erase (i); return true; } } - } - for (unsigned int i = 0; i < 4; ++i) { - if (m_q[i] && b.inside (m_q[i]->box (m_center))) { - return b.remove (value); + } else if (m_q[n]) { + + if (m_q[n]->erase (value)) { + if (m_q[n]->empty ()) { + delete m_q[n]; + m_q[n] = 0; + } + return true; } + } return false; @@ -131,12 +141,39 @@ public: return count; } + size_t levels () const + { + size_t l = 1; + for (unsigned int n = 0; n < 4; ++n) { + if (m_q[n]) { + l = std::max (l, m_q[n]->levels () + 1); + } + } + return l; + } + + bool check_top (const box_type &total_box) const + { + return check (propose_ucenter (total_box)); + } + private: bool m_split; point_type m_center; quad_tree_node *m_q [4]; objects_vector m_objects; + int quad_for (const box_type &box) const + { + int sx = coord_traits::less (box.right (), m_center.x ()) ? 0 : (coord_traits::less (m_center.x (), box.left ()) ? 1 : -1); + int sy = coord_traits::less (box.top (), m_center.y ()) ? 0 : (coord_traits::less (m_center.y (), box.bottom ()) ? 2 : -1); + if (sx < 0 || sy < 0) { + return -1; + } else { + return sx + sy; + } + } + box_type box (const point_type &ucenter) const { return box_type (ucenter, ucenter - (ucenter - m_center) * 2.0); @@ -185,33 +222,24 @@ private: } box_type b = BC () (value); - // @@@ should exclude m_center on box - if (b.contains (m_center)) { + int n = quad_for (b); + + if (n < 0) { m_objects.push_back (value); return; } - for (unsigned int i = 0; i < 4; ++i) { - box_type bq = q (i, ucenter); - if (b.inside (bq)) { - if (! m_q[i]) { - m_q[i] = new quad_tree_node (bq.center ()); - } - m_q[i]->insert (value, m_center); - return; + if (b.inside (box (ucenter))) { + if (! m_q[n]) { + box_type bq = q (n, ucenter); + m_q[n] = new quad_tree_node (bq.center ()); } + m_q[n]->insert (value, m_center); + } else { + grow (m_center - (m_center - ucenter) * 2.0); + insert (value, ucenter); } - for (unsigned int i = 0; i < 4; ++i) { - if (m_q[i]) { - grow (m_center - (m_center - m_q[i]->center ()) * 2.0); - insert (value, ucenter); - return; - } - } - - tl_assert (false); - } } @@ -238,6 +266,61 @@ private: coord_type dy = std::max (std::abs (total_box.bottom () - m_center.y ()), std::abs (total_box.top () - m_center.y ())); return m_center - vector_type (dx, dy); } + + bool check (const point_type &ucenter) const + { + bool result = true; + + box_type bq = box (ucenter); + + for (auto i = m_objects.begin (); i != m_objects.end (); ++i) { + box_type b = BC () (*i); + if (! b.inside (bq)) { + tl::error << "Box " << b.to_string () << " not inside quad box " << bq.to_string (); + result = false; + } + } + + if (m_split) { + + for (auto i = m_objects.begin (); i != m_objects.end (); ++i) { + box_type b = BC () (*i); + int n = quad_for (b); + if (n >= 0) { + tl::error << "Box " << b.to_string () << " on quad level not overlapping multiple quads"; + result = false; + } + } + + for (unsigned int n = 0; n < 4; ++n) { + if (m_q[n]) { + m_q[n]->check (m_center); + box_type bbq = m_q[n]->box (m_center); + if (bbq != q (n, ucenter)) { + tl::error << "Quad not centered (quad box is " << bbq.to_string () << ", should be " << q (n, ucenter).to_string (); + result = false; + } + } + } + + } else { + + if (m_objects.size () > thr) { + tl::error << "Non-split object count exceeds threshold " << m_objects.size () << " > " << thr; + result = false; + } + + for (unsigned int n = 0; n < 4; ++n) { + if (m_q[n]) { + tl::error << "Non-split node has child nodes"; + result = false; + } + } + + } + + return result; + } }; template @@ -315,9 +398,9 @@ public: } } - while (! m_stack.empty ()) { + m_stack.pop_back (); - m_stack.pop_back (); + while (! m_stack.empty ()) { int &n = m_stack.back ().second; while (++n < 4) { @@ -329,6 +412,8 @@ public: } } + m_stack.pop_back (); + } } }; @@ -415,6 +500,7 @@ private: box_type m_box; }; +// @@@ TODO: copy, assignment, move, swap template class quad_tree { @@ -445,6 +531,16 @@ public: return m_root.size (); } + size_t levels () const + { + return m_root.levels (); + } + + bool check () const + { + return m_root.check_top (m_total_box); + } + void insert (const T &value) { box_type b = BC () (value); @@ -456,9 +552,9 @@ public: m_root.insert_top (value, m_total_box); } - void erase (const T &value) + bool erase (const T &value) { - m_root.remove (value); + return m_root.erase (value); } quad_tree_flat_iterator begin () const diff --git a/src/db/unit_tests/dbQuadTreeTests.cc b/src/db/unit_tests/dbQuadTreeTests.cc index fd74e164b..62eec169e 100644 --- a/src/db/unit_tests/dbQuadTreeTests.cc +++ b/src/db/unit_tests/dbQuadTreeTests.cc @@ -28,6 +28,18 @@ typedef db::quad_tree, size_t (1)> my_quad_tree; +std::string find_all (const my_quad_tree &qt) +{ + std::vector v; + auto i = qt.begin (); + while (! i.at_end ()) { + v.push_back (i->to_string ()); + ++i; + } + std::sort (v.begin (), v.end ()); + return tl::join (v, "/"); +} + std::string find_touching (const my_quad_tree &qt, const db::DBox &box) { std::vector v; @@ -40,6 +52,20 @@ std::string find_touching (const my_quad_tree &qt, const db::DBox &box) return tl::join (v, "/"); } +std::string find_touching_from_all (const my_quad_tree &qt, const db::DBox &box) +{ + std::vector v; + auto i = qt.begin (); + while (! i.at_end ()) { + if (i->touches (box)) { + v.push_back (i->to_string ()); + } + ++i; + } + std::sort (v.begin (), v.end ()); + return tl::join (v, "/"); +} + std::string find_overlapping (const my_quad_tree &qt, const db::DBox &box) { std::vector v; @@ -52,26 +78,181 @@ std::string find_overlapping (const my_quad_tree &qt, const db::DBox &box) return tl::join (v, "/"); } +std::string find_overlapping_from_all (const my_quad_tree &qt, const db::DBox &box) +{ + std::vector v; + auto i = qt.begin (); + while (! i.at_end ()) { + if (i->overlaps (box)) { + v.push_back (i->to_string ()); + } + ++i; + } + std::sort (v.begin (), v.end ()); + return tl::join (v, "/"); +} + TEST(basic) { my_quad_tree tree; EXPECT_EQ (tree.empty (), true); EXPECT_EQ (tree.size (), size_t (0)); + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (tree.levels (), size_t (1)); tree.insert (db::DBox ()); EXPECT_EQ (tree.empty (), true); EXPECT_EQ (tree.size (), size_t (0)); + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (tree.levels (), size_t (1)); tree.insert (db::DBox (-1, -2, 3, 4)); EXPECT_EQ (tree.empty (), false); EXPECT_EQ (tree.size (), size_t (1)); + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (tree.levels (), size_t (1)); - EXPECT_EQ (find_touching (tree, db::DBox (-2, 0, -1, 0)), "..."); + EXPECT_EQ (find_all (tree), "(-1,-2;3,4)"); + + db::DBox bx; + + bx = db::DBox (-2, 0, -1, 0); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, -3, -1, -2); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, -3, -1, -2.5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 4, -1, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 4.5, -1, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 3, -1, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 3, -1.5, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); tree.insert (db::DBox (-1, -3, 3, 0)); EXPECT_EQ (tree.empty (), false); - EXPECT_EQ (tree.size (), size_t (1)); + EXPECT_EQ (tree.size (), size_t (2)); + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (tree.levels (), size_t (1)); - EXPECT_EQ (find_touching (tree, db::DBox (-2, 0, -1, 0)), "..."); + EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;3,0)"); + + bx = db::DBox (-2, 0, -1, 0); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, -3, -1, -2); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, -3, -1, -2.5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 4, -1, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 4.5, -1, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 3, -1, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 3, -1.5, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + + tree.insert (db::DBox (-1, -3, -0.5, -2)); + EXPECT_EQ (tree.empty (), false); + EXPECT_EQ (tree.size (), size_t (3)); + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (tree.levels (), size_t (2)); + + EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;3,0)"); + + bx = db::DBox (-2, 0, -1, 0); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, -3, -1, -2); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, -3, -1, -2.5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 4, -1, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 4.5, -1, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 3, -1, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 3, -1.5, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + + tree.insert (db::DBox (-1, -3, -0.5, 2)); + EXPECT_EQ (tree.empty (), false); + EXPECT_EQ (tree.size (), size_t (4)); + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (tree.levels (), size_t (2)); + + EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); + + bx = db::DBox (-2, 0, -1, 0); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, -3, -1, -2); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, -3, -1, -2.5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 4, -1, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 4.5, -1, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 3, -1, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); + bx = db::DBox (-2, 3, -1.5, 5); + EXPECT_EQ (find_touching (tree, bx), find_touching_from_all (tree, bx)); + EXPECT_EQ (find_overlapping (tree, bx), find_overlapping_from_all (tree, bx)); +} + +TEST(remove) +{ + my_quad_tree tree; + tree.insert (db::DBox (-1, -2, 3, 4)); + tree.insert (db::DBox (-1, -3, 3, 0)); + tree.insert (db::DBox (-1, -3, -0.5, -2)); + tree.insert (db::DBox (-1, -3, -0.5, 2)); + + EXPECT_EQ (tree.check (), true); + + EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); + + EXPECT_EQ (tree.erase (db::DBox (-1, -3, -0.5, -1)), false); + EXPECT_EQ (tree.erase (db::DBox (-1, -3, -0.5, -2)), true); + EXPECT_EQ (tree.check (), true); + + EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); + + while (! tree.empty ()) { + EXPECT_EQ (tree.erase (*tree.begin ()), true); + EXPECT_EQ (tree.check (), true); + } + + EXPECT_EQ (tree.size (), size_t (0)); + EXPECT_EQ (tree.levels (), size_t (1)); } From 5c8e0539eeecad7987a2678e7db89a389ca205e9 Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Sun, 16 Mar 2025 13:56:16 +0100 Subject: [PATCH 03/17] WIP (quad tree) --- src/db/db/dbQuadTree.h | 119 ++++++++++++++++++++++- src/db/unit_tests/dbQuadTreeTests.cc | 135 +++++++++++++++++++++++++++ 2 files changed, 249 insertions(+), 5 deletions(-) diff --git a/src/db/db/dbQuadTree.h b/src/db/db/dbQuadTree.h index 353f7ddcd..93674ea69 100644 --- a/src/db/db/dbQuadTree.h +++ b/src/db/db/dbQuadTree.h @@ -44,16 +44,77 @@ public: quad_tree_node (const point_type ¢er) : m_split (false), m_center (center) { - for (unsigned int i = 0; i < 4; ++i) { - m_q [i] = 0; - } + init (); } ~quad_tree_node () { + clear (); + } + + quad_tree_node (const quad_tree_node &other) + : m_split (false), m_center (center) + { + init (); + operator= (other); + } + + quad_tree_node &operator= (const quad_tree_node &other) + { + if (this != &other) { + clear (); + m_split = other.m_split; + m_center = other.m_center; + m_objects = other.m_objects; + for (unsigned int i = 0; i < 4; ++i) { + if (other.m_q[i]) { + m_q[i] = other.m_q[i]->clone (); + } + } + } + return *this; + } + + quad_tree_node (quad_tree_node &&other) + { + init (); + swap (other); + } + + quad_tree_node &operator= (quad_tree_node &&other) + { + swap (other); + return *this; + } + + void swap (quad_tree_node &other) + { + if (this != &other) { + std::swap (m_center, other.m_center); + std::swap (m_split, other.m_split); + m_objects.swap (other.m_objects); + for (unsigned int i = 0; i < 4; ++i) { + std::swap (m_q[i], other.m_q[i]); + } + } + } + + quad_tree_node *clone () const + { + quad_tree_node *node = new quad_tree_node (m_center); + *node = *this; + return node; + } + + void clear () + { + m_objects.clear (); + m_split = false; for (unsigned int i = 0; i < 4; ++i) { - delete m_q [i]; - m_q [i] = 0; + if (m_q[i]) { + delete m_q[i]; + } + m_q[i] = 0; } } @@ -163,6 +224,13 @@ private: quad_tree_node *m_q [4]; objects_vector m_objects; + void init () + { + for (unsigned int i = 0; i < 4; ++i) { + m_q[i] = 0; + } + } + int quad_for (const box_type &box) const { int sx = coord_traits::less (box.right (), m_center.x ()) ? 0 : (coord_traits::less (m_center.x (), box.left ()) ? 1 : -1); @@ -521,6 +589,47 @@ public: // .. nothing yet .. } + quad_tree (const quad_tree &other) + : m_root (point_type ()) + { + operator= (other); + } + + quad_tree &operator= (const quad_tree &other) + { + if (this != &other) { + m_root = other.m_root; + m_total_box = other.m_total_box; + } + return *this; + } + + quad_tree (quad_tree &&other) + : m_root (point_type ()) + { + swap (other); + } + + quad_tree &operator= (quad_tree &&other) + { + swap (other); + return *this; + } + + void clear () + { + m_root.clear (); + m_total_box = box_type (); + } + + void swap (quad_tree &other) + { + if (this != &other) { + m_root.swap (other.m_root); + std::swap (m_total_box, m_total_box); + } + } + bool empty () const { return m_root.empty (); diff --git a/src/db/unit_tests/dbQuadTreeTests.cc b/src/db/unit_tests/dbQuadTreeTests.cc index 62eec169e..6d39ea49b 100644 --- a/src/db/unit_tests/dbQuadTreeTests.cc +++ b/src/db/unit_tests/dbQuadTreeTests.cc @@ -256,3 +256,138 @@ TEST(remove) EXPECT_EQ (tree.levels (), size_t (1)); } +TEST(clear) +{ + my_quad_tree tree; + tree.insert (db::DBox (-1, -2, 3, 4)); + tree.insert (db::DBox (-1, -3, 3, 0)); + tree.insert (db::DBox (-1, -3, -0.5, -2)); + tree.insert (db::DBox (-1, -3, -0.5, 2)); + + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); + + tree.clear (); + + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (tree.empty (), true); + EXPECT_EQ (tree.size (), size_t (0)); + EXPECT_EQ (tree.levels (), size_t (1)); + EXPECT_EQ (find_all (tree), ""); +} + +TEST(copy) +{ + my_quad_tree tree; + tree.insert (db::DBox (-1, -2, 3, 4)); + tree.insert (db::DBox (-1, -3, 3, 0)); + tree.insert (db::DBox (-1, -3, -0.5, -2)); + tree.insert (db::DBox (-1, -3, -0.5, 2)); + + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); + EXPECT_EQ (tree.levels (), size_t (2)); + + my_quad_tree tree2 (tree); + tree.clear (); + + EXPECT_EQ (tree2.check (), true); + EXPECT_EQ (find_all (tree2), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); + EXPECT_EQ (tree2.levels (), size_t (2)); +} + +TEST(assign) +{ + my_quad_tree tree; + tree.insert (db::DBox (-1, -2, 3, 4)); + tree.insert (db::DBox (-1, -3, 3, 0)); + tree.insert (db::DBox (-1, -3, -0.5, -2)); + tree.insert (db::DBox (-1, -3, -0.5, 2)); + + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); + EXPECT_EQ (tree.levels (), size_t (2)); + + my_quad_tree tree2; + tree2 = tree; + tree.clear (); + + EXPECT_EQ (tree2.check (), true); + EXPECT_EQ (find_all (tree2), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); + EXPECT_EQ (tree2.levels (), size_t (2)); +} + +TEST(swap) +{ + my_quad_tree tree; + tree.insert (db::DBox (-1, -2, 3, 4)); + tree.insert (db::DBox (-1, -3, 3, 0)); + tree.insert (db::DBox (-1, -3, -0.5, -2)); + tree.insert (db::DBox (-1, -3, -0.5, 2)); + + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); + EXPECT_EQ (tree.levels (), size_t (2)); + + my_quad_tree tree2; + tree2.swap (tree); + + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (tree.empty (), true); + EXPECT_EQ (find_all (tree), ""); + EXPECT_EQ (tree.levels (), size_t (1)); + + EXPECT_EQ (tree2.check (), true); + EXPECT_EQ (find_all (tree2), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); + EXPECT_EQ (tree2.levels (), size_t (2)); +} + +TEST(move) +{ + my_quad_tree tree; + tree.insert (db::DBox (-1, -2, 3, 4)); + tree.insert (db::DBox (-1, -3, 3, 0)); + tree.insert (db::DBox (-1, -3, -0.5, -2)); + tree.insert (db::DBox (-1, -3, -0.5, 2)); + + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); + EXPECT_EQ (tree.levels (), size_t (2)); + + my_quad_tree tree2; + tree2 = std::move (tree); + + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (tree.empty (), true); + EXPECT_EQ (find_all (tree), ""); + EXPECT_EQ (tree.levels (), size_t (1)); + + EXPECT_EQ (tree2.check (), true); + EXPECT_EQ (find_all (tree2), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); + EXPECT_EQ (tree2.levels (), size_t (2)); +} + +TEST(move_ctor) +{ + my_quad_tree tree; + tree.insert (db::DBox (-1, -2, 3, 4)); + tree.insert (db::DBox (-1, -3, 3, 0)); + tree.insert (db::DBox (-1, -3, -0.5, -2)); + tree.insert (db::DBox (-1, -3, -0.5, 2)); + + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); + EXPECT_EQ (tree.levels (), size_t (2)); + + my_quad_tree tree2 (std::move (tree)); + + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (tree.empty (), true); + EXPECT_EQ (find_all (tree), ""); + EXPECT_EQ (tree.levels (), size_t (1)); + + EXPECT_EQ (tree2.check (), true); + EXPECT_EQ (find_all (tree2), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); + EXPECT_EQ (tree2.levels (), size_t (2)); +} + From f1f35ae2a4afc5dabb41bd40e06f14bbf5116717 Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Sun, 16 Mar 2025 14:46:03 +0100 Subject: [PATCH 04/17] WIP (quad tree) --- src/db/db/dbQuadTree.h | 33 ++++++----- src/db/unit_tests/dbQuadTreeTests.cc | 84 ++++++++++++++++++++++++++++ 2 files changed, 104 insertions(+), 13 deletions(-) diff --git a/src/db/db/dbQuadTree.h b/src/db/db/dbQuadTree.h index 93674ea69..546c007ff 100644 --- a/src/db/db/dbQuadTree.h +++ b/src/db/db/dbQuadTree.h @@ -290,22 +290,26 @@ private: } box_type b = BC () (value); - int n = quad_for (b); - - if (n < 0) { - m_objects.push_back (value); - return; - } - if (b.inside (box (ucenter))) { - if (! m_q[n]) { - box_type bq = q (n, ucenter); - m_q[n] = new quad_tree_node (bq.center ()); + + int n = quad_for (b); + if (n < 0) { + m_objects.push_back (value); + } else { + if (! m_q[n]) { + box_type bq = q (n, ucenter); + m_q[n] = new quad_tree_node (bq.center ()); + } + m_q[n]->insert (value, m_center); } - m_q[n]->insert (value, m_center); + } else { - grow (m_center - (m_center - ucenter) * 2.0); - insert (value, ucenter); + + tl_assert (m_q[0] || m_q[1] || m_q[2] || m_q[3]); + point_type new_ucenter = m_center - (m_center - ucenter) * 2.0; + grow (new_ucenter); + insert (value, new_ucenter); + } } @@ -317,6 +321,7 @@ private: if (m_q[i]) { quad_tree_node *n = m_q[i]; m_q[i] = new quad_tree_node (q (i, ucenter).center ()); + m_q[i]->m_split = true; m_q[i]->m_q[3 - i] = n; } } @@ -470,6 +475,8 @@ public: while (! m_stack.empty ()) { + qn = m_stack.back ().first; + int &n = m_stack.back ().second; while (++n < 4) { box_type bq = qn->q_box (n); diff --git a/src/db/unit_tests/dbQuadTreeTests.cc b/src/db/unit_tests/dbQuadTreeTests.cc index 6d39ea49b..1866907d1 100644 --- a/src/db/unit_tests/dbQuadTreeTests.cc +++ b/src/db/unit_tests/dbQuadTreeTests.cc @@ -26,6 +26,8 @@ #include "tlUnitTest.h" #include "tlString.h" +#include + typedef db::quad_tree, size_t (1)> my_quad_tree; std::string find_all (const my_quad_tree &qt) @@ -256,6 +258,58 @@ TEST(remove) EXPECT_EQ (tree.levels (), size_t (1)); } +TEST(grow) +{ + my_quad_tree tree; + tree.insert (db::DBox (-1, -2, 3, 4)); + tree.insert (db::DBox (-1, -3, 3, 0)); + tree.insert (db::DBox (-1, -3, -0.5, -2)); + tree.insert (db::DBox (-1, -3, -0.5, 2)); + EXPECT_EQ (tree.levels (), size_t (2)); + tree.insert (db::DBox (-100, -3, -99, 2)); + EXPECT_EQ (tree.levels (), size_t (8)); + + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)/(-100,-3;-99,2)"); + EXPECT_EQ (find_overlapping (tree, db::DBox (-100, -100, -90, 100)), "(-100,-3;-99,2)"); + + bool r = true; + while (r && ! tree.empty ()) { + r = tree.erase (*tree.begin ()); + EXPECT_EQ (r, true); + EXPECT_EQ (tree.check (), true); + } + + EXPECT_EQ (tree.size (), size_t (0)); + EXPECT_EQ (tree.levels (), size_t (1)); +} + +TEST(grow2) +{ + my_quad_tree tree; + tree.insert (db::DBox (-1, -2, 3, 4)); + tree.insert (db::DBox (-1, -3, 3, 0)); + tree.insert (db::DBox (-1, -3, -0.5, -2)); + tree.insert (db::DBox (-1, -3, -0.5, 2)); + EXPECT_EQ (tree.levels (), size_t (2)); + tree.insert (db::DBox (-100, -3, -99, -1)); + EXPECT_EQ (tree.levels (), size_t (8)); + + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)/(-100,-3;-99,-1)"); + EXPECT_EQ (find_overlapping (tree, db::DBox (-100, -100, -90, 100)), "(-100,-3;-99,-1)"); + + bool r = true; + while (r && ! tree.empty ()) { + r = tree.erase (*tree.begin ()); + EXPECT_EQ (r, true); + EXPECT_EQ (tree.check (), true); + } + + EXPECT_EQ (tree.size (), size_t (0)); + EXPECT_EQ (tree.levels (), size_t (1)); +} + TEST(clear) { my_quad_tree tree; @@ -391,3 +445,33 @@ TEST(move_ctor) EXPECT_EQ (tree2.levels (), size_t (2)); } +static double rvalue () +{ + return ((rand () % 1000000) - 5000) * 0.001; +} + +static db::DBox rbox () +{ + return db::DBox (db::DPoint (rvalue (), rvalue ()), db::DPoint (rvalue (), rvalue ())); +} + +TEST(many) +{ + return; // @@@ + my_quad_tree tree; + + unsigned int n = 1000; + + for (unsigned int i = 0; i < n; ++i) { + tree.insert (rbox ()); + } + + EXPECT_EQ (tree.check (), true); + + bool r = true; + while (r && ! tree.empty ()) { + r = tree.erase (*tree.begin ()); + EXPECT_EQ (r, true); + EXPECT_EQ (tree.check (), true); + } +} From 4369835e8b641f82d426fda584efc667b157a643 Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Sun, 16 Mar 2025 16:36:49 +0100 Subject: [PATCH 05/17] WIP (quad tree) --- src/db/db/dbQuadTree.h | 55 ++++++++++++----- src/db/unit_tests/dbQuadTreeTests.cc | 88 ++++++++++++++++------------ 2 files changed, 90 insertions(+), 53 deletions(-) diff --git a/src/db/db/dbQuadTree.h b/src/db/db/dbQuadTree.h index 546c007ff..31034827f 100644 --- a/src/db/db/dbQuadTree.h +++ b/src/db/db/dbQuadTree.h @@ -30,7 +30,7 @@ namespace db { -template +template class quad_tree_node { public: @@ -137,7 +137,7 @@ public: if (! m_split || n < 0) { for (auto i = m_objects.begin (); i != m_objects.end (); ++i) { - if (*i == value) { + if (CMP () (*i, value)) { m_objects.erase (i); return true; } @@ -290,7 +290,7 @@ private: } box_type b = BC () (value); - if (b.inside (box (ucenter))) { + if (inside (b, box (ucenter))) { int n = quad_for (b); if (n < 0) { @@ -337,7 +337,7 @@ private: coord_type dx = std::max (std::abs (total_box.left () - m_center.x ()), std::abs (total_box.right () - m_center.y ())); coord_type dy = std::max (std::abs (total_box.bottom () - m_center.y ()), std::abs (total_box.top () - m_center.y ())); - return m_center - vector_type (dx, dy); + return m_center - vector_type (dx * 2, dy * 2); } bool check (const point_type &ucenter) const @@ -348,10 +348,15 @@ private: for (auto i = m_objects.begin (); i != m_objects.end (); ++i) { box_type b = BC () (*i); - if (! b.inside (bq)) { + if (! inside (b, bq)) { tl::error << "Box " << b.to_string () << " not inside quad box " << bq.to_string (); result = false; } + if (coord_traits::equal (b.left (), ucenter.x ()) || coord_traits::equal (b.right (), ucenter.x ()) || + coord_traits::equal (b.bottom (), ucenter.y ()) || coord_traits::equal (b.top (), ucenter.y ())) { + tl::error << "Box " << b.to_string () << " touches center of upper-level quad " << ucenter.to_string (); + result = false; + } } if (m_split) { @@ -367,9 +372,11 @@ private: for (unsigned int n = 0; n < 4; ++n) { if (m_q[n]) { - m_q[n]->check (m_center); + if (! m_q[n]->check (m_center)) { + result = false; + } box_type bbq = m_q[n]->box (m_center); - if (bbq != q (n, ucenter)) { + if (! bbq.equal (q (n, ucenter))) { tl::error << "Quad not centered (quad box is " << bbq.to_string () << ", should be " << q (n, ucenter).to_string (); result = false; } @@ -394,13 +401,23 @@ private: return result; } + + static bool inside (const box_type &box, const box_type &in) + { + if (box.empty () || in.empty ()) { + return false; + } else { + return ! coord_traits::less (box.left (), in.left ()) && ! coord_traits::less (in.right (), box.right ()) && + ! coord_traits::less (box.bottom (), in.bottom ()) && ! coord_traits::less (in.top (), box.top ()); + } + } }; -template +template class quad_tree_iterator { public: - typedef quad_tree_node quad_tree_node_type; + typedef quad_tree_node quad_tree_node_type; typedef typename T::coord_type coord_type; typedef db::box box_type; @@ -575,15 +592,23 @@ private: box_type m_box; }; -// @@@ TODO: copy, assignment, move, swap -template +template +struct quad_tree_default_cmp +{ + bool operator() (const T &a, const T &b) const + { + return a == b; + } +}; + +template > class quad_tree { public: - typedef quad_tree_node quad_tree_node_type; - typedef quad_tree_iterator > quad_tree_flat_iterator; - typedef quad_tree_iterator > quad_tree_touching_iterator; - typedef quad_tree_iterator > quad_tree_overlapping_iterator; + typedef quad_tree_node quad_tree_node_type; + typedef quad_tree_iterator > quad_tree_flat_iterator; + typedef quad_tree_iterator > quad_tree_touching_iterator; + typedef quad_tree_iterator > quad_tree_overlapping_iterator; typedef typename T::coord_type coord_type; typedef db::box box_type; typedef db::point point_type; diff --git a/src/db/unit_tests/dbQuadTreeTests.cc b/src/db/unit_tests/dbQuadTreeTests.cc index 1866907d1..c0df73a44 100644 --- a/src/db/unit_tests/dbQuadTreeTests.cc +++ b/src/db/unit_tests/dbQuadTreeTests.cc @@ -28,9 +28,17 @@ #include -typedef db::quad_tree, size_t (1)> my_quad_tree; +struct MyQuadTreeCMP +{ + bool operator() (const db::DBox &a, const db::DBox &b) const + { + return a.equal (b); + } +}; -std::string find_all (const my_quad_tree &qt) +typedef db::quad_tree, size_t (1), MyQuadTreeCMP> MyQuadTree; + +std::string find_all (const MyQuadTree &qt) { std::vector v; auto i = qt.begin (); @@ -42,7 +50,7 @@ std::string find_all (const my_quad_tree &qt) return tl::join (v, "/"); } -std::string find_touching (const my_quad_tree &qt, const db::DBox &box) +std::string find_touching (const MyQuadTree &qt, const db::DBox &box) { std::vector v; auto i = qt.begin_touching (box); @@ -54,7 +62,7 @@ std::string find_touching (const my_quad_tree &qt, const db::DBox &box) return tl::join (v, "/"); } -std::string find_touching_from_all (const my_quad_tree &qt, const db::DBox &box) +std::string find_touching_from_all (const MyQuadTree &qt, const db::DBox &box) { std::vector v; auto i = qt.begin (); @@ -68,7 +76,7 @@ std::string find_touching_from_all (const my_quad_tree &qt, const db::DBox &box) return tl::join (v, "/"); } -std::string find_overlapping (const my_quad_tree &qt, const db::DBox &box) +std::string find_overlapping (const MyQuadTree &qt, const db::DBox &box) { std::vector v; auto i = qt.begin_overlapping (box); @@ -80,7 +88,7 @@ std::string find_overlapping (const my_quad_tree &qt, const db::DBox &box) return tl::join (v, "/"); } -std::string find_overlapping_from_all (const my_quad_tree &qt, const db::DBox &box) +std::string find_overlapping_from_all (const MyQuadTree &qt, const db::DBox &box) { std::vector v; auto i = qt.begin (); @@ -96,7 +104,7 @@ std::string find_overlapping_from_all (const my_quad_tree &qt, const db::DBox &b TEST(basic) { - my_quad_tree tree; + MyQuadTree tree; EXPECT_EQ (tree.empty (), true); EXPECT_EQ (tree.size (), size_t (0)); EXPECT_EQ (tree.check (), true); @@ -174,7 +182,7 @@ TEST(basic) EXPECT_EQ (tree.empty (), false); EXPECT_EQ (tree.size (), size_t (3)); EXPECT_EQ (tree.check (), true); - EXPECT_EQ (tree.levels (), size_t (2)); + EXPECT_EQ (tree.levels (), size_t (3)); EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;3,0)"); @@ -204,7 +212,7 @@ TEST(basic) EXPECT_EQ (tree.empty (), false); EXPECT_EQ (tree.size (), size_t (4)); EXPECT_EQ (tree.check (), true); - EXPECT_EQ (tree.levels (), size_t (2)); + EXPECT_EQ (tree.levels (), size_t (3)); EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); @@ -233,7 +241,7 @@ TEST(basic) TEST(remove) { - my_quad_tree tree; + MyQuadTree tree; tree.insert (db::DBox (-1, -2, 3, 4)); tree.insert (db::DBox (-1, -3, 3, 0)); tree.insert (db::DBox (-1, -3, -0.5, -2)); @@ -260,12 +268,12 @@ TEST(remove) TEST(grow) { - my_quad_tree tree; + MyQuadTree tree; tree.insert (db::DBox (-1, -2, 3, 4)); tree.insert (db::DBox (-1, -3, 3, 0)); tree.insert (db::DBox (-1, -3, -0.5, -2)); tree.insert (db::DBox (-1, -3, -0.5, 2)); - EXPECT_EQ (tree.levels (), size_t (2)); + EXPECT_EQ (tree.levels (), size_t (3)); tree.insert (db::DBox (-100, -3, -99, 2)); EXPECT_EQ (tree.levels (), size_t (8)); @@ -286,12 +294,12 @@ TEST(grow) TEST(grow2) { - my_quad_tree tree; + MyQuadTree tree; tree.insert (db::DBox (-1, -2, 3, 4)); tree.insert (db::DBox (-1, -3, 3, 0)); tree.insert (db::DBox (-1, -3, -0.5, -2)); tree.insert (db::DBox (-1, -3, -0.5, 2)); - EXPECT_EQ (tree.levels (), size_t (2)); + EXPECT_EQ (tree.levels (), size_t (3)); tree.insert (db::DBox (-100, -3, -99, -1)); EXPECT_EQ (tree.levels (), size_t (8)); @@ -312,7 +320,7 @@ TEST(grow2) TEST(clear) { - my_quad_tree tree; + MyQuadTree tree; tree.insert (db::DBox (-1, -2, 3, 4)); tree.insert (db::DBox (-1, -3, 3, 0)); tree.insert (db::DBox (-1, -3, -0.5, -2)); @@ -332,7 +340,7 @@ TEST(clear) TEST(copy) { - my_quad_tree tree; + MyQuadTree tree; tree.insert (db::DBox (-1, -2, 3, 4)); tree.insert (db::DBox (-1, -3, 3, 0)); tree.insert (db::DBox (-1, -3, -0.5, -2)); @@ -340,19 +348,19 @@ TEST(copy) EXPECT_EQ (tree.check (), true); EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); - EXPECT_EQ (tree.levels (), size_t (2)); + EXPECT_EQ (tree.levels (), size_t (3)); - my_quad_tree tree2 (tree); + MyQuadTree tree2 (tree); tree.clear (); EXPECT_EQ (tree2.check (), true); EXPECT_EQ (find_all (tree2), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); - EXPECT_EQ (tree2.levels (), size_t (2)); + EXPECT_EQ (tree2.levels (), size_t (3)); } TEST(assign) { - my_quad_tree tree; + MyQuadTree tree; tree.insert (db::DBox (-1, -2, 3, 4)); tree.insert (db::DBox (-1, -3, 3, 0)); tree.insert (db::DBox (-1, -3, -0.5, -2)); @@ -360,20 +368,20 @@ TEST(assign) EXPECT_EQ (tree.check (), true); EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); - EXPECT_EQ (tree.levels (), size_t (2)); + EXPECT_EQ (tree.levels (), size_t (3)); - my_quad_tree tree2; + MyQuadTree tree2; tree2 = tree; tree.clear (); EXPECT_EQ (tree2.check (), true); EXPECT_EQ (find_all (tree2), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); - EXPECT_EQ (tree2.levels (), size_t (2)); + EXPECT_EQ (tree2.levels (), size_t (3)); } TEST(swap) { - my_quad_tree tree; + MyQuadTree tree; tree.insert (db::DBox (-1, -2, 3, 4)); tree.insert (db::DBox (-1, -3, 3, 0)); tree.insert (db::DBox (-1, -3, -0.5, -2)); @@ -381,9 +389,9 @@ TEST(swap) EXPECT_EQ (tree.check (), true); EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); - EXPECT_EQ (tree.levels (), size_t (2)); + EXPECT_EQ (tree.levels (), size_t (3)); - my_quad_tree tree2; + MyQuadTree tree2; tree2.swap (tree); EXPECT_EQ (tree.check (), true); @@ -393,12 +401,12 @@ TEST(swap) EXPECT_EQ (tree2.check (), true); EXPECT_EQ (find_all (tree2), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); - EXPECT_EQ (tree2.levels (), size_t (2)); + EXPECT_EQ (tree2.levels (), size_t (3)); } TEST(move) { - my_quad_tree tree; + MyQuadTree tree; tree.insert (db::DBox (-1, -2, 3, 4)); tree.insert (db::DBox (-1, -3, 3, 0)); tree.insert (db::DBox (-1, -3, -0.5, -2)); @@ -406,9 +414,9 @@ TEST(move) EXPECT_EQ (tree.check (), true); EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); - EXPECT_EQ (tree.levels (), size_t (2)); + EXPECT_EQ (tree.levels (), size_t (3)); - my_quad_tree tree2; + MyQuadTree tree2; tree2 = std::move (tree); EXPECT_EQ (tree.check (), true); @@ -418,12 +426,12 @@ TEST(move) EXPECT_EQ (tree2.check (), true); EXPECT_EQ (find_all (tree2), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); - EXPECT_EQ (tree2.levels (), size_t (2)); + EXPECT_EQ (tree2.levels (), size_t (3)); } TEST(move_ctor) { - my_quad_tree tree; + MyQuadTree tree; tree.insert (db::DBox (-1, -2, 3, 4)); tree.insert (db::DBox (-1, -3, 3, 0)); tree.insert (db::DBox (-1, -3, -0.5, -2)); @@ -431,9 +439,9 @@ TEST(move_ctor) EXPECT_EQ (tree.check (), true); EXPECT_EQ (find_all (tree), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); - EXPECT_EQ (tree.levels (), size_t (2)); + EXPECT_EQ (tree.levels (), size_t (3)); - my_quad_tree tree2 (std::move (tree)); + MyQuadTree tree2 (std::move (tree)); EXPECT_EQ (tree.check (), true); EXPECT_EQ (tree.empty (), true); @@ -442,7 +450,7 @@ TEST(move_ctor) EXPECT_EQ (tree2.check (), true); EXPECT_EQ (find_all (tree2), "(-1,-2;3,4)/(-1,-3;-0.5,-2)/(-1,-3;-0.5,2)/(-1,-3;3,0)"); - EXPECT_EQ (tree2.levels (), size_t (2)); + EXPECT_EQ (tree2.levels (), size_t (3)); } static double rvalue () @@ -452,13 +460,16 @@ static double rvalue () static db::DBox rbox () { - return db::DBox (db::DPoint (rvalue (), rvalue ()), db::DPoint (rvalue (), rvalue ())); + db::DBox box; + while ((box = db::DBox (db::DPoint (rvalue (), rvalue ()), db::DPoint (rvalue (), rvalue ()))).empty ()) { + ; + } + return box; } TEST(many) { - return; // @@@ - my_quad_tree tree; + MyQuadTree tree; unsigned int n = 1000; @@ -467,6 +478,7 @@ TEST(many) } EXPECT_EQ (tree.check (), true); + EXPECT_EQ (tree.size (), size_t (n)); bool r = true; while (r && ! tree.empty ()) { From 477e2b5a31fd93b62426d1eca01de6576b5a1ae5 Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Sun, 16 Mar 2025 16:43:56 +0100 Subject: [PATCH 06/17] WIP (quad tree) --- src/db/unit_tests/dbQuadTreeTests.cc | 28 ++++++++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/src/db/unit_tests/dbQuadTreeTests.cc b/src/db/unit_tests/dbQuadTreeTests.cc index c0df73a44..26c06c196 100644 --- a/src/db/unit_tests/dbQuadTreeTests.cc +++ b/src/db/unit_tests/dbQuadTreeTests.cc @@ -50,7 +50,7 @@ std::string find_all (const MyQuadTree &qt) return tl::join (v, "/"); } -std::string find_touching (const MyQuadTree &qt, const db::DBox &box) +std::string find_touching (const MyQuadTree &qt, const db::DBox &box, bool report = false) { std::vector v; auto i = qt.begin_touching (box); @@ -58,6 +58,9 @@ std::string find_touching (const MyQuadTree &qt, const db::DBox &box) v.push_back (i->to_string ()); ++i; } + if (report) { + tl::info << v.size () << " items found."; + } std::sort (v.begin (), v.end ()); return tl::join (v, "/"); } @@ -76,7 +79,7 @@ std::string find_touching_from_all (const MyQuadTree &qt, const db::DBox &box) return tl::join (v, "/"); } -std::string find_overlapping (const MyQuadTree &qt, const db::DBox &box) +std::string find_overlapping (const MyQuadTree &qt, const db::DBox &box, bool report = false) { std::vector v; auto i = qt.begin_overlapping (box); @@ -84,6 +87,9 @@ std::string find_overlapping (const MyQuadTree &qt, const db::DBox &box) v.push_back (i->to_string ()); ++i; } + if (report) { + tl::info << v.size () << " items found."; + } std::sort (v.begin (), v.end ()); return tl::join (v, "/"); } @@ -472,6 +478,7 @@ TEST(many) MyQuadTree tree; unsigned int n = 1000; + unsigned int ntests = 100; for (unsigned int i = 0; i < n; ++i) { tree.insert (rbox ()); @@ -480,10 +487,27 @@ TEST(many) EXPECT_EQ (tree.check (), true); EXPECT_EQ (tree.size (), size_t (n)); + bool report = false; + for (unsigned int i = 0; i < ntests; ++i) { + if (report) { + tl::info << "Test iteration " << i << " ..."; + } + auto bx = rbox (); + EXPECT_EQ (find_overlapping (tree, bx, report), find_overlapping_from_all (tree, bx)); + EXPECT_EQ (find_touching (tree, bx, report), find_touching_from_all (tree, bx)); + bx = db::DBox (bx.center (), bx.center ()); + EXPECT_EQ (find_touching (tree, bx, report), find_touching_from_all (tree, bx)); + } + bool r = true; while (r && ! tree.empty ()) { r = tree.erase (*tree.begin ()); EXPECT_EQ (r, true); EXPECT_EQ (tree.check (), true); } + + EXPECT_EQ (tree.empty (), true); + EXPECT_EQ (tree.check (), true); + EXPECT_EQ (tree.levels (), size_t (1)); + EXPECT_EQ (tree.size (), size_t (0)); } From 4aeb94d42e70717bc7269da592ebf03e12ce798a Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Sun, 16 Mar 2025 17:13:06 +0100 Subject: [PATCH 07/17] WIP (quad tree) --- src/db/unit_tests/dbQuadTreeTests.cc | 78 ++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) diff --git a/src/db/unit_tests/dbQuadTreeTests.cc b/src/db/unit_tests/dbQuadTreeTests.cc index 26c06c196..219f8dc20 100644 --- a/src/db/unit_tests/dbQuadTreeTests.cc +++ b/src/db/unit_tests/dbQuadTreeTests.cc @@ -25,6 +25,7 @@ #include "dbBoxConvert.h" #include "tlUnitTest.h" #include "tlString.h" +#include "tlTimer.h" #include @@ -473,6 +474,13 @@ static db::DBox rbox () return box; } +static db::DBox rbox (double dim) +{ + db::DBox box; + db::DPoint c (rvalue (), rvalue ()); + return box = db::DBox (c, c).enlarged (db::DVector (dim * 0.5, dim * 0.5)); +} + TEST(many) { MyQuadTree tree; @@ -511,3 +519,73 @@ TEST(many) EXPECT_EQ (tree.levels (), size_t (1)); EXPECT_EQ (tree.size (), size_t (0)); } + +TEST(timing_insert) +{ + MyQuadTree tree; + + { + unsigned int n = 1000000; + tl::SelfTimer timer (tl::sprintf ("%d inserts ..", int (n))); + for (unsigned int i = 0; i < n; ++i) { + tree.insert (rbox ()); + } + } + + tree.clear (); + + { + unsigned int n = 2000000; + tl::SelfTimer timer (tl::sprintf ("%d inserts ..", int (n))); + for (unsigned int i = 0; i < n; ++i) { + tree.insert (rbox ()); + } + } +} + +TEST(timing_lookup) +{ + MyQuadTree tree; + + unsigned int n = 250000; + for (unsigned int i = 0; i < n; ++i) { + tree.insert (rbox (10.0)); + } + + unsigned int ntests = 100; + + std::vector > > tests; + for (unsigned int i = 0; i < ntests; ++i) { + db::DBox bx = rbox (10.0); + tests.push_back (std::make_pair (bx, std::make_pair (size_t (0), size_t (0)))); + } + + { + tl::SelfTimer timer (tl::sprintf ("%d tests (lookup) ..", int (ntests))); + for (auto t = tests.begin (); t != tests.end (); ++t) { + size_t n = 0; + for (auto i = tree.begin_touching (t->first); ! i.at_end (); ++i) { + ++n; + } + t->second.first = n; + } + } + + { + tl::SelfTimer timer (tl::sprintf ("%d tests (brute force) ..", int (ntests))); + for (auto t = tests.begin (); t != tests.end (); ++t) { + size_t n = 0; + for (auto i = tree.begin (); ! i.at_end (); ++i) { + if (i->touches (t->first)) { + ++n; + } + } + t->second.second = n; + } + } + + for (auto t = tests.begin (); t != tests.end (); ++t) { + EXPECT_EQ (t->second.first, t->second.second); + } +} + From f3037d11f3a771e16f4e65bfd45630194b7ec424 Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Sun, 16 Mar 2025 18:16:15 +0100 Subject: [PATCH 08/17] WIP (quad tree) --- src/db/db/dbQuadTree.h | 190 ++++++++++++++++++++------- src/db/unit_tests/dbQuadTreeTests.cc | 10 +- 2 files changed, 151 insertions(+), 49 deletions(-) diff --git a/src/db/db/dbQuadTree.h b/src/db/db/dbQuadTree.h index 31034827f..453ea78b8 100644 --- a/src/db/db/dbQuadTree.h +++ b/src/db/db/dbQuadTree.h @@ -30,6 +30,9 @@ namespace db { +/** + * @brief The quad tree node implementation class + */ template class quad_tree_node { @@ -42,9 +45,9 @@ public: typedef db::coord_traits coord_traits; quad_tree_node (const point_type ¢er) - : m_split (false), m_center (center) + : m_center (center) { - init (); + init (true); } ~quad_tree_node () @@ -53,9 +56,9 @@ public: } quad_tree_node (const quad_tree_node &other) - : m_split (false), m_center (center) + : m_center (center) { - init (); + init (true); operator= (other); } @@ -63,12 +66,14 @@ public: { if (this != &other) { clear (); - m_split = other.m_split; m_center = other.m_center; m_objects = other.m_objects; - for (unsigned int i = 0; i < 4; ++i) { - if (other.m_q[i]) { - m_q[i] = other.m_q[i]->clone (); + if (! other.is_leaf ()) { + init (false); + for (unsigned int i = 0; i < 4; ++i) { + if (other.m_q[i]) { + m_q[i] = other.m_q[i]->clone (); + } } } } @@ -77,7 +82,7 @@ public: quad_tree_node (quad_tree_node &&other) { - init (); + init (true); swap (other); } @@ -91,7 +96,6 @@ public: { if (this != &other) { std::swap (m_center, other.m_center); - std::swap (m_split, other.m_split); m_objects.swap (other.m_objects); for (unsigned int i = 0; i < 4; ++i) { std::swap (m_q[i], other.m_q[i]); @@ -109,12 +113,14 @@ public: void clear () { m_objects.clear (); - m_split = false; - for (unsigned int i = 0; i < 4; ++i) { - if (m_q[i]) { - delete m_q[i]; + if (! is_leaf ()) { + for (unsigned int i = 0; i < 4; ++i) { + if (m_q[i]) { + delete m_q[i]; + } + m_q[i] = 0; } - m_q[i] = 0; + init (true); } } @@ -134,7 +140,7 @@ public: int n = quad_for (b); - if (! m_split || n < 0) { + if (is_leaf () || n < 0) { for (auto i = m_objects.begin (); i != m_objects.end (); ++i) { if (CMP () (*i, value)) { @@ -165,7 +171,7 @@ public: box_type q_box (unsigned int n) const { - if (m_q[n]) { + if (! is_leaf () && m_q[n]) { return m_q[n]->box (m_center); } else { return box_type (); @@ -174,15 +180,18 @@ public: quad_tree_node *node (unsigned int n) const { + tl_assert (! is_leaf ()); return m_q [n]; } bool empty () const { if (m_objects.empty ()) { - for (unsigned int n = 0; n < 4; ++n) { - if (m_q[n] && ! m_q[n]->empty ()) { - return false; + if (! is_leaf ()) { + for (unsigned int n = 0; n < 4; ++n) { + if (m_q[n] && ! m_q[n]->empty ()) { + return false; + } } } return true; @@ -194,9 +203,11 @@ public: size_t size () const { size_t count = m_objects.size (); - for (unsigned int n = 0; n < 4; ++n) { - if (m_q[n]) { - count += m_q[n]->size (); + if (! is_leaf ()) { + for (unsigned int n = 0; n < 4; ++n) { + if (m_q[n]) { + count += m_q[n]->size (); + } } } return count; @@ -205,9 +216,11 @@ public: size_t levels () const { size_t l = 1; - for (unsigned int n = 0; n < 4; ++n) { - if (m_q[n]) { - l = std::max (l, m_q[n]->levels () + 1); + if (! is_leaf ()) { + for (unsigned int n = 0; n < 4; ++n) { + if (m_q[n]) { + l = std::max (l, m_q[n]->levels () + 1); + } } } return l; @@ -219,16 +232,23 @@ public: } private: - bool m_split; point_type m_center; quad_tree_node *m_q [4]; objects_vector m_objects; - void init () + void init (bool is_leaf) { for (unsigned int i = 0; i < 4; ++i) { m_q[i] = 0; } + if (is_leaf) { + m_q[0] = reinterpret_cast (1); + } + } + + bool is_leaf () const + { + return m_q[0] == reinterpret_cast (1); } int quad_for (const box_type &box) const @@ -267,7 +287,7 @@ private: void split (const point_type &ucenter) { - m_split = true; + init (false); objects_vector ov; ov.swap (m_objects); @@ -279,13 +299,13 @@ private: void insert (const T &value, const point_type &ucenter) { - if (! m_split && m_objects.size () + 1 < thr) { + if (is_leaf () && m_objects.size () + 1 < thr) { m_objects.push_back (value); } else { - if (! m_split) { + if (is_leaf ()) { split (ucenter); } @@ -321,7 +341,7 @@ private: if (m_q[i]) { quad_tree_node *n = m_q[i]; m_q[i] = new quad_tree_node (q (i, ucenter).center ()); - m_q[i]->m_split = true; + m_q[i]->init (false); m_q[i]->m_q[3 - i] = n; } } @@ -329,9 +349,11 @@ private: point_type propose_ucenter (const box_type &total_box) const { - for (unsigned int i = 0; i < 4; ++i) { - if (m_q[i]) { - return m_center - (m_center - m_q[i]->center ()) * 2.0; + if (! is_leaf ()) { + for (unsigned int i = 0; i < 4; ++i) { + if (m_q[i]) { + return m_center - (m_center - m_q[i]->center ()) * 2.0; + } } } @@ -359,7 +381,7 @@ private: } } - if (m_split) { + if (! is_leaf ()) { for (auto i = m_objects.begin (); i != m_objects.end (); ++i) { box_type b = BC () (*i); @@ -390,13 +412,6 @@ private: result = false; } - for (unsigned int n = 0; n < 4; ++n) { - if (m_q[n]) { - tl::error << "Non-split node has child nodes"; - result = false; - } - } - } return result; @@ -413,6 +428,9 @@ private: } }; +/** + * @brief The iterator implementation class + */ template class quad_tree_iterator { @@ -510,6 +528,9 @@ public: } }; +/** + * @brief The selector for implementing the all-iterator + */ template class quad_tree_always_sel { @@ -528,6 +549,9 @@ public: } }; +/** + * @brief The selector for implementing the touching iterator + */ template class quad_tree_touching_sel { @@ -560,6 +584,9 @@ private: box_type m_box; }; +/** + * @brief The selector for implementing the overlapping iterator + */ template class quad_tree_overlapping_sel { @@ -592,6 +619,9 @@ private: box_type m_box; }; +/** + * @brief The default compare function + */ template struct quad_tree_default_cmp { @@ -601,7 +631,22 @@ struct quad_tree_default_cmp } }; -template > +/** + * @brief A generic quad tree implementation + * + * In contrast to the box_tree implementation, this is a self-sorting implementation + * which is more generic. + * + * @param T The value to be stored + * @param BC The box converter + * @param thr The number of items per leaf node before splitting + * @param CMP The compare function (equality) + * + * T needs to have a type member named "coord_type". + * BC is a function of T and delivers a db::box for T + * CMP is a function that delivers a bool value for a pair of T (equality) + */ +template > class quad_tree { public: @@ -615,18 +660,29 @@ public: typedef db::vector vector_type; typedef std::vector objects_vector; + /** + * @brief Default constructor + * + * This creates an empty tree. + */ quad_tree () : m_root (point_type ()) { // .. nothing yet .. } + /** + * @brief Copy constructor + */ quad_tree (const quad_tree &other) : m_root (point_type ()) { operator= (other); } + /** + * @brief Assignment + */ quad_tree &operator= (const quad_tree &other) { if (this != &other) { @@ -636,52 +692,79 @@ public: return *this; } + /** + * @brief Move constructor + */ quad_tree (quad_tree &&other) : m_root (point_type ()) { swap (other); } + /** + * @brief Move assignment + */ quad_tree &operator= (quad_tree &&other) { swap (other); return *this; } + /** + * @brief Empties the tree + */ void clear () { m_root.clear (); m_total_box = box_type (); } + /** + * @brief Swaps the tree with another + */ void swap (quad_tree &other) { if (this != &other) { m_root.swap (other.m_root); - std::swap (m_total_box, m_total_box); + std::swap (m_total_box, other.m_total_box); } } + /** + * @brief Returns a value indicating whether the tree is empty + */ bool empty () const { return m_root.empty (); } + /** + * @brief Returns the number of items stored in the tree + */ size_t size () const { return m_root.size (); } + /** + * @brief Returns the number of quad levels (for testing) + */ size_t levels () const { return m_root.levels (); } + /** + * @brief Checks the tree for consistency (for testing) + */ bool check () const { return m_root.check_top (m_total_box); } + /** + * @brief Inserts an object into the tree + */ void insert (const T &value) { box_type b = BC () (value); @@ -693,16 +776,30 @@ public: m_root.insert_top (value, m_total_box); } + /** + * @brief Erases the given element from the tree + * + * @return true, if the element was found and erased. + * + * If multiple elements of the same kind are stored, the + * first one is erased. + */ bool erase (const T &value) { return m_root.erase (value); } + /** + * @brief begin iterator for all elements + */ quad_tree_flat_iterator begin () const { return quad_tree_flat_iterator (&m_root, quad_tree_always_sel ()); } + /** + * @brief begin iterator for all elements overlapping the given box + */ quad_tree_overlapping_iterator begin_overlapping (const box_type &box) const { if (m_total_box.overlaps (box)) { @@ -712,6 +809,9 @@ public: } } + /** + * @brief begin iterator for all elements touching the given box + */ quad_tree_touching_iterator begin_touching (const box_type &box) const { if (m_total_box.touches (box)) { diff --git a/src/db/unit_tests/dbQuadTreeTests.cc b/src/db/unit_tests/dbQuadTreeTests.cc index 219f8dc20..8ffe12327 100644 --- a/src/db/unit_tests/dbQuadTreeTests.cc +++ b/src/db/unit_tests/dbQuadTreeTests.cc @@ -530,6 +530,7 @@ TEST(timing_insert) for (unsigned int i = 0; i < n; ++i) { tree.insert (rbox ()); } + tl::info << "Quad levels: " << tree.levels (); } tree.clear (); @@ -540,6 +541,7 @@ TEST(timing_insert) for (unsigned int i = 0; i < n; ++i) { tree.insert (rbox ()); } + tl::info << "Quad levels: " << tree.levels (); } } @@ -547,16 +549,16 @@ TEST(timing_lookup) { MyQuadTree tree; - unsigned int n = 250000; + unsigned int n = 1000000; for (unsigned int i = 0; i < n; ++i) { - tree.insert (rbox (10.0)); + tree.insert (rbox (5.0)); } - unsigned int ntests = 100; + unsigned int ntests = 1000; std::vector > > tests; for (unsigned int i = 0; i < ntests; ++i) { - db::DBox bx = rbox (10.0); + db::DBox bx = rbox (5.0); tests.push_back (std::make_pair (bx, std::make_pair (size_t (0), size_t (0)))); } From 4e65b96cb769156e81dd12f90a893a029cd38eb9 Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Sun, 16 Mar 2025 21:10:15 +0100 Subject: [PATCH 09/17] WIP (quad tree) --- src/db/db/dbQuadTree.h | 32 +++++++++++++------------------- 1 file changed, 13 insertions(+), 19 deletions(-) diff --git a/src/db/db/dbQuadTree.h b/src/db/db/dbQuadTree.h index 453ea78b8..d17ee1af7 100644 --- a/src/db/db/dbQuadTree.h +++ b/src/db/db/dbQuadTree.h @@ -37,12 +37,11 @@ template class quad_tree_node { public: - typedef typename T::coord_type coord_type; - typedef db::box box_type; - typedef db::point point_type; - typedef db::vector vector_type; + typedef typename BC::box_type box_type; + typedef typename box_type::point_type point_type; + typedef typename box_type::vector_type vector_type; typedef std::vector objects_vector; - typedef db::coord_traits coord_traits; + typedef db::coord_traits coord_traits; quad_tree_node (const point_type ¢er) : m_center (center) @@ -357,8 +356,8 @@ private: } } - coord_type dx = std::max (std::abs (total_box.left () - m_center.x ()), std::abs (total_box.right () - m_center.y ())); - coord_type dy = std::max (std::abs (total_box.bottom () - m_center.y ()), std::abs (total_box.top () - m_center.y ())); + typename box_type::coord_type dx = std::max (std::abs (total_box.left () - m_center.x ()), std::abs (total_box.right () - m_center.y ())); + typename box_type::coord_type dy = std::max (std::abs (total_box.bottom () - m_center.y ()), std::abs (total_box.top () - m_center.y ())); return m_center - vector_type (dx * 2, dy * 2); } @@ -436,8 +435,7 @@ class quad_tree_iterator { public: typedef quad_tree_node quad_tree_node_type; - typedef typename T::coord_type coord_type; - typedef db::box box_type; + typedef typename BC::box_type box_type; quad_tree_iterator () : m_s (), m_i (0) @@ -535,8 +533,7 @@ template class quad_tree_always_sel { public: - typedef typename T::coord_type coord_type; - typedef db::box box_type; + typedef typename BC::box_type box_type; bool select (const T &) const { @@ -556,8 +553,7 @@ template class quad_tree_touching_sel { public: - typedef typename T::coord_type coord_type; - typedef db::box box_type; + typedef typename BC::box_type box_type; quad_tree_touching_sel () { @@ -591,8 +587,7 @@ template class quad_tree_overlapping_sel { public: - typedef typename T::coord_type coord_type; - typedef db::box box_type; + typedef typename BC::box_type box_type; quad_tree_overlapping_sel () { @@ -654,10 +649,9 @@ public: typedef quad_tree_iterator > quad_tree_flat_iterator; typedef quad_tree_iterator > quad_tree_touching_iterator; typedef quad_tree_iterator > quad_tree_overlapping_iterator; - typedef typename T::coord_type coord_type; - typedef db::box box_type; - typedef db::point point_type; - typedef db::vector vector_type; + typedef typename BC::box_type box_type; + typedef typename box_type::point_type point_type; + typedef typename box_type::vector_type vector_type; typedef std::vector objects_vector; /** From cd62f62140c1dfeeb0b7c4bcac52163164302a6d Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Sun, 16 Mar 2025 23:24:06 +0100 Subject: [PATCH 10/17] WIP (quad tree) --- src/db/db/dbQuadTree.h | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/db/db/dbQuadTree.h b/src/db/db/dbQuadTree.h index d17ee1af7..7cc9a13dc 100644 --- a/src/db/db/dbQuadTree.h +++ b/src/db/db/dbQuadTree.h @@ -373,9 +373,12 @@ private: tl::error << "Box " << b.to_string () << " not inside quad box " << bq.to_string (); result = false; } + point_type ucenter_opp = m_center + (m_center - ucenter); if (coord_traits::equal (b.left (), ucenter.x ()) || coord_traits::equal (b.right (), ucenter.x ()) || - coord_traits::equal (b.bottom (), ucenter.y ()) || coord_traits::equal (b.top (), ucenter.y ())) { - tl::error << "Box " << b.to_string () << " touches center of upper-level quad " << ucenter.to_string (); + coord_traits::equal (b.bottom (), ucenter.y ()) || coord_traits::equal (b.top (), ucenter.y ()) || + coord_traits::equal (b.left (), ucenter_opp.x ()) || coord_traits::equal (b.right (), ucenter_opp.x ()) || + coord_traits::equal (b.bottom (), ucenter_opp.y ()) || coord_traits::equal (b.top (), ucenter_opp.y ())) { + tl::error << "Box " << b.to_string () << " touches quad boundary " << ucenter.to_string () << " .. " << ucenter_opp.to_string (); result = false; } } @@ -421,8 +424,8 @@ private: if (box.empty () || in.empty ()) { return false; } else { - return ! coord_traits::less (box.left (), in.left ()) && ! coord_traits::less (in.right (), box.right ()) && - ! coord_traits::less (box.bottom (), in.bottom ()) && ! coord_traits::less (in.top (), box.top ()); + return coord_traits::less (in.left (), box.left ()) && coord_traits::less (box.right (), in.right ()) && + coord_traits::less (in.bottom (), box.bottom ()) && coord_traits::less (box.top (), in.top ()); } } }; From 6596008826580c8013adb9cacd2146a4f3393b36 Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Tue, 18 Mar 2025 00:10:00 +0100 Subject: [PATCH 11/17] instrumenting triangles implementation with Quad Tree, but without effect --- src/db/db/dbTriangles.cc | 120 ++++++++++++++++++++------ src/db/db/dbTriangles.h | 14 ++- src/db/unit_tests/dbTrianglesTests.cc | 1 + 3 files changed, 106 insertions(+), 29 deletions(-) diff --git a/src/db/db/dbTriangles.cc b/src/db/db/dbTriangles.cc index 5a0475a5d..b325579c9 100644 --- a/src/db/db/dbTriangles.cc +++ b/src/db/db/dbTriangles.cc @@ -44,9 +44,7 @@ Triangles::Triangles () Triangles::~Triangles () { - while (! mp_triangles.empty ()) { - remove_triangle (mp_triangles.begin ().operator-> ()); - } + clear (); } db::Vertex * @@ -88,6 +86,19 @@ Triangles::create_triangle (TriangleEdge *e1, TriangleEdge *e2, TriangleEdge *e3 db::Triangle *res = new db::Triangle (e1, e2, e3); res->set_id (++m_id); mp_triangles.push_back (res); + + m_triangle_qt.insert (res); + + // @@@ +#if 0 + if (mp_triangles.size () != m_triangle_qt.size ()) { + size_t a = mp_triangles.size (); + size_t b = m_triangle_qt.size (); + printf("@@@ %ld -- %ld\n", a, b); fflush(stdout); + } + tl_assert (mp_triangles.size () == m_triangle_qt.size ()); // @@@ +#endif + // @@@ return res; } @@ -99,7 +110,20 @@ Triangles::remove_triangle (db::Triangle *tri) edges [i] = tri->edge (i); } + bool removed = m_triangle_qt.erase (tri); + tl_assert (removed); + delete tri; + // @@@ +#if 0 + if (mp_triangles.size () != m_triangle_qt.size ()) { + size_t a = mp_triangles.size (); + size_t b = m_triangle_qt.size (); + printf("@@@ %ld -- %ld\n", a, b); fflush(stdout); + } + tl_assert (mp_triangles.size () == m_triangle_qt.size ()); // @@@ +#endif + // @@@ // clean up edges we do no longer need for (int i = 0; i < 3; ++i) { @@ -395,6 +419,8 @@ Triangles::insert (db::Vertex *vertex, std::list > *n std::vector Triangles::find_triangle_for_point (const db::DPoint &point) { +// @@@ +#if 1 // minimize distance search db::TriangleEdge *edge = find_closest_edge (point); std::vector res; @@ -405,7 +431,29 @@ Triangles::find_triangle_for_point (const db::DPoint &point) } } } + + // @@@ + std::set setb; + for (auto i = m_triangle_qt.begin_touching (db::DBox (point, point)); ! i.at_end (); ++i) { + if ((*i)->contains (point) >= 0) { + setb.insert (*i); + } + } + std::set seta (res.begin (), res.end ()); + if (seta != setb) { + tl_assert (false); + } + // @@@ return res; +#else + std::vector res; + for (auto i = m_triangle_qt.begin_touching (db::DBox (point, point)); ! i.at_end (); ++i) { + if ((*i)->contains (point) >= 0) { + res.push_back (*i); + } + } + return res; +#endif } db::TriangleEdge * @@ -413,24 +461,7 @@ Triangles::find_closest_edge (const db::DPoint &p, db::Vertex *vstart, bool insi { if (!vstart) { - if (m_is_constrained && ! m_initial_segments.empty ()) { - - // Use the closest initial segment for the starting point - // TODO: m_initial_segments should be some search tree - - double d = -1; - for (auto i = m_initial_segments.begin (); i != m_initial_segments.end (); ++i) { - if (db::sprod_sign (p - i->first.p1 (), i->first.d ()) >= 0 && - db::sprod_sign (p - i->first.p2 (), i->first.d ()) < 0) { - double ds = std::abs (i->first.distance (p)); - if (d < 0.0 || ds < d) { - vstart = i->second; - d = ds; - } - } - } - - } else if (! mp_triangles.empty ()) { + if (! mp_triangles.empty ()) { unsigned int ls = 0; size_t n = m_vertex_heap.size (); @@ -660,7 +691,7 @@ Triangles::split_triangle (db::Triangle *t, db::Vertex *vertex, std::list &tris, db::Vertex *vertex, db::TriangleEdge *split_edge, std::list > *new_triangles_out) +Triangles::split_triangles_on_edge (const std::vector &_tris /*@@@*/, db::Vertex *vertex, db::TriangleEdge *split_edge, std::list > *new_triangles_out) { TriangleEdge *s1 = create_edge (split_edge->v1 (), vertex); TriangleEdge *s2 = create_edge (split_edge->v2 (), vertex); @@ -669,6 +700,12 @@ Triangles::split_triangles_on_edge (const std::vector &tris, db: std::vector new_triangles; + std::vector tris; + tris.reserve (2); + for (auto t = split_edge->begin_triangles (); t != split_edge->end_triangles (); ++t) { + tris.push_back (t.operator-> ()); + } + for (auto t = tris.begin (); t != tris.end (); ++t) { (*t)->unlink (); @@ -1181,19 +1218,47 @@ Triangles::search_edges_crossing (Vertex *from, Vertex *to) } db::Vertex * -Triangles::find_vertex_for_point (const db::DPoint &pt) +Triangles::find_vertex_for_point (const db::DPoint &point) { - db::TriangleEdge *edge = find_closest_edge (pt); +// @@@ +#if 1 // minimize distance search + db::TriangleEdge *edge = find_closest_edge (point); if (!edge) { return 0; } db::Vertex *v = 0; - if (edge->v1 ()->equal (pt)) { + if (edge->v1 ()->equal (point)) { v = edge->v1 (); - } else if (edge->v2 ()->equal (pt)) { + } else if (edge->v2 ()->equal (point)) { v = edge->v2 (); } + // @@@ + db::Vertex *vv = 0; + for (auto i = m_triangle_qt.begin_touching (db::DBox (point, point)); ! i.at_end () && ! vv; ++i) { + db::Triangle *t = *i; + for (unsigned int i = 0; i < 3 && ! vv; ++i) { + if (t->vertex (i)->equal (point)) { + vv = t->vertex (i); + } + } + } + if (vv != v) { + tl_assert (false); // @@@ + } + // @@@ return v; +#else + for (auto i = m_triangle_qt.begin_touching (db::DBox (point, point)); ! i.at_end (); ++i) { + db::Triangle *t = *i; + for (unsigned int i = 0; i < 3; ++i) { + if (t->vertex (i)->equal (point)) { + return t->vertex (i); + } + } + } + return 0; +#endif +// @@@ } db::TriangleEdge * @@ -1337,7 +1402,6 @@ void Triangles::constrain (const std::vector > &contours) { tl_assert (! m_is_constrained); - m_initial_segments.clear (); std::vector > > resolved_edges; @@ -1351,7 +1415,6 @@ Triangles::constrain (const std::vector > &contours) db::DEdge e (**v, **vv); resolved_edges.push_back (std::make_pair (e, std::vector ())); resolved_edges.back ().second = ensure_edge (*v, *vv); - m_initial_segments.push_back (std::make_pair (e, *v)); } } @@ -1436,6 +1499,7 @@ void Triangles::clear () { mp_triangles.clear (); + m_triangle_qt.clear (); m_edges_heap.clear (); m_vertex_heap.clear (); m_returned_edges.clear (); diff --git a/src/db/db/dbTriangles.h b/src/db/db/dbTriangles.h index f15b2dc7a..96f342a1e 100644 --- a/src/db/db/dbTriangles.h +++ b/src/db/db/dbTriangles.h @@ -29,6 +29,7 @@ #include "dbTriangle.h" #include "dbBox.h" #include "dbRegion.h" +#include "dbQuadTree.h" #include "tlObjectCollection.h" #include "tlStableVector.h" @@ -306,12 +307,23 @@ protected: std::vector find_inside_circle (const db::DPoint ¢er, double radius) const; private: + struct TriangleBoxConvert + { + typedef db::DBox box_type; + box_type operator() (db::Triangle *t) const + { + return t ? t->bbox () : box_type (); + } + }; + + typedef db::quad_tree triangle_qt_type; + tl::list mp_triangles; + triangle_qt_type m_triangle_qt; tl::stable_vector m_edges_heap; std::vector m_returned_edges; tl::stable_vector m_vertex_heap; bool m_is_constrained; - std::vector > m_initial_segments; size_t m_level; size_t m_id; size_t m_flips, m_hops; diff --git a/src/db/unit_tests/dbTrianglesTests.cc b/src/db/unit_tests/dbTrianglesTests.cc index 9b7dd021c..c2209c3da 100644 --- a/src/db/unit_tests/dbTrianglesTests.cc +++ b/src/db/unit_tests/dbTrianglesTests.cc @@ -277,6 +277,7 @@ namespace { TEST(insert_many) { + return; // @@@ srand (0); TestableTriangles tris; From df631aa97013fe698287e2aaedbf2a18916fe60c Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Tue, 18 Mar 2025 00:19:15 +0100 Subject: [PATCH 12/17] Some minor refactoring --- src/db/db/dbQuadTree.h | 23 ++++---- src/db/db/dbTriangles.cc | 83 +-------------------------- src/db/db/dbTriangles.h | 8 +-- src/db/unit_tests/dbTrianglesTests.cc | 1 - 4 files changed, 13 insertions(+), 102 deletions(-) diff --git a/src/db/db/dbQuadTree.h b/src/db/db/dbQuadTree.h index 7cc9a13dc..45d1c0a53 100644 --- a/src/db/db/dbQuadTree.h +++ b/src/db/db/dbQuadTree.h @@ -128,15 +128,13 @@ public: return m_center; } - void insert_top (const T &value, const box_type &total_box) + void insert_top (const T &value, const box_type &total_box, const box_type &b) { - insert (value, propose_ucenter (total_box)); + insert (value, propose_ucenter (total_box), b); } - bool erase (const T &value) + bool erase (const T &value, const box_type &b) { - box_type b = BC () (value); - int n = quad_for (b); if (is_leaf () || n < 0) { @@ -150,7 +148,7 @@ public: } else if (m_q[n]) { - if (m_q[n]->erase (value)) { + if (m_q[n]->erase (value, b)) { if (m_q[n]->empty ()) { delete m_q[n]; m_q[n] = 0; @@ -292,11 +290,11 @@ private: ov.swap (m_objects); for (auto o = ov.begin (); o != ov.end (); ++o) { - insert (*o, ucenter); + insert (*o, ucenter, BC () (*o)); } } - void insert (const T &value, const point_type &ucenter) + void insert (const T &value, const point_type &ucenter, const box_type &b) { if (is_leaf () && m_objects.size () + 1 < thr) { @@ -308,7 +306,6 @@ private: split (ucenter); } - box_type b = BC () (value); if (inside (b, box (ucenter))) { int n = quad_for (b); @@ -319,7 +316,7 @@ private: box_type bq = q (n, ucenter); m_q[n] = new quad_tree_node (bq.center ()); } - m_q[n]->insert (value, m_center); + m_q[n]->insert (value, m_center, b); } } else { @@ -327,7 +324,7 @@ private: tl_assert (m_q[0] || m_q[1] || m_q[2] || m_q[3]); point_type new_ucenter = m_center - (m_center - ucenter) * 2.0; grow (new_ucenter); - insert (value, new_ucenter); + insert (value, new_ucenter, b); } @@ -770,7 +767,7 @@ public: } m_total_box += b; - m_root.insert_top (value, m_total_box); + m_root.insert_top (value, m_total_box, b); } /** @@ -783,7 +780,7 @@ public: */ bool erase (const T &value) { - return m_root.erase (value); + return m_root.erase (value, BC () (value)); } /** diff --git a/src/db/db/dbTriangles.cc b/src/db/db/dbTriangles.cc index b325579c9..246afa2e6 100644 --- a/src/db/db/dbTriangles.cc +++ b/src/db/db/dbTriangles.cc @@ -87,18 +87,6 @@ Triangles::create_triangle (TriangleEdge *e1, TriangleEdge *e2, TriangleEdge *e3 res->set_id (++m_id); mp_triangles.push_back (res); - m_triangle_qt.insert (res); - - // @@@ -#if 0 - if (mp_triangles.size () != m_triangle_qt.size ()) { - size_t a = mp_triangles.size (); - size_t b = m_triangle_qt.size (); - printf("@@@ %ld -- %ld\n", a, b); fflush(stdout); - } - tl_assert (mp_triangles.size () == m_triangle_qt.size ()); // @@@ -#endif - // @@@ return res; } @@ -110,20 +98,7 @@ Triangles::remove_triangle (db::Triangle *tri) edges [i] = tri->edge (i); } - bool removed = m_triangle_qt.erase (tri); - tl_assert (removed); - delete tri; - // @@@ -#if 0 - if (mp_triangles.size () != m_triangle_qt.size ()) { - size_t a = mp_triangles.size (); - size_t b = m_triangle_qt.size (); - printf("@@@ %ld -- %ld\n", a, b); fflush(stdout); - } - tl_assert (mp_triangles.size () == m_triangle_qt.size ()); // @@@ -#endif - // @@@ // clean up edges we do no longer need for (int i = 0; i < 3; ++i) { @@ -379,7 +354,6 @@ Triangles::insert_point (db::DCoord x, db::DCoord y, std::list > *new_triangles) { - // @@@ printf("@@@ insert %d\n", (int)mp_triangles.size ()); fflush(stdout); std::vector tris = find_triangle_for_point (*vertex); // the new vertex is outside the domain @@ -400,7 +374,7 @@ Triangles::insert (db::Vertex *vertex, std::list > *n if (! on_edges.empty ()) { if (on_edges.size () == size_t (1)) { - split_triangles_on_edge (tris, vertex, on_edges.front (), new_triangles); + split_triangles_on_edge (vertex, on_edges.front (), new_triangles); return vertex; } else { // the vertex is already present @@ -419,8 +393,6 @@ Triangles::insert (db::Vertex *vertex, std::list > *n std::vector Triangles::find_triangle_for_point (const db::DPoint &point) { -// @@@ -#if 1 // minimize distance search db::TriangleEdge *edge = find_closest_edge (point); std::vector res; @@ -432,28 +404,7 @@ Triangles::find_triangle_for_point (const db::DPoint &point) } } - // @@@ - std::set setb; - for (auto i = m_triangle_qt.begin_touching (db::DBox (point, point)); ! i.at_end (); ++i) { - if ((*i)->contains (point) >= 0) { - setb.insert (*i); - } - } - std::set seta (res.begin (), res.end ()); - if (seta != setb) { - tl_assert (false); - } - // @@@ return res; -#else - std::vector res; - for (auto i = m_triangle_qt.begin_touching (db::DBox (point, point)); ! i.at_end (); ++i) { - if ((*i)->contains (point) >= 0) { - res.push_back (*i); - } - } - return res; -#endif } db::TriangleEdge * @@ -691,7 +642,7 @@ Triangles::split_triangle (db::Triangle *t, db::Vertex *vertex, std::list &_tris /*@@@*/, db::Vertex *vertex, db::TriangleEdge *split_edge, std::list > *new_triangles_out) +Triangles::split_triangles_on_edge (db::Vertex *vertex, db::TriangleEdge *split_edge, std::list > *new_triangles_out) { TriangleEdge *s1 = create_edge (split_edge->v1 (), vertex); TriangleEdge *s2 = create_edge (split_edge->v2 (), vertex); @@ -1220,8 +1171,6 @@ Triangles::search_edges_crossing (Vertex *from, Vertex *to) db::Vertex * Triangles::find_vertex_for_point (const db::DPoint &point) { -// @@@ -#if 1 // minimize distance search db::TriangleEdge *edge = find_closest_edge (point); if (!edge) { return 0; @@ -1232,33 +1181,7 @@ Triangles::find_vertex_for_point (const db::DPoint &point) } else if (edge->v2 ()->equal (point)) { v = edge->v2 (); } - // @@@ - db::Vertex *vv = 0; - for (auto i = m_triangle_qt.begin_touching (db::DBox (point, point)); ! i.at_end () && ! vv; ++i) { - db::Triangle *t = *i; - for (unsigned int i = 0; i < 3 && ! vv; ++i) { - if (t->vertex (i)->equal (point)) { - vv = t->vertex (i); - } - } - } - if (vv != v) { - tl_assert (false); // @@@ - } - // @@@ return v; -#else - for (auto i = m_triangle_qt.begin_touching (db::DBox (point, point)); ! i.at_end (); ++i) { - db::Triangle *t = *i; - for (unsigned int i = 0; i < 3; ++i) { - if (t->vertex (i)->equal (point)) { - return t->vertex (i); - } - } - } - return 0; -#endif -// @@@ } db::TriangleEdge * @@ -1499,7 +1422,6 @@ void Triangles::clear () { mp_triangles.clear (); - m_triangle_qt.clear (); m_edges_heap.clear (); m_vertex_heap.clear (); m_returned_edges.clear (); @@ -1691,7 +1613,6 @@ Triangles::refine (const TriangulateParameters ¶meters) auto cr = (*t)->circumcircle(); auto center = cr.first; - // printf("@@@ %s -- %s\n", (*t)->to_string().c_str(), center.to_string ().c_str()); fflush(stdout); if ((*t)->contains (center) >= 0) { if (tl::verbosity () >= parameters.base_verbosity + 20) { diff --git a/src/db/db/dbTriangles.h b/src/db/db/dbTriangles.h index 96f342a1e..67c48b4c6 100644 --- a/src/db/db/dbTriangles.h +++ b/src/db/db/dbTriangles.h @@ -20,8 +20,6 @@ */ - - #ifndef HDR_dbTriangles #define HDR_dbTriangles @@ -29,7 +27,6 @@ #include "dbTriangle.h" #include "dbBox.h" #include "dbRegion.h" -#include "dbQuadTree.h" #include "tlObjectCollection.h" #include "tlStableVector.h" @@ -316,10 +313,7 @@ private: } }; - typedef db::quad_tree triangle_qt_type; - tl::list mp_triangles; - triangle_qt_type m_triangle_qt; tl::stable_vector m_edges_heap; std::vector m_returned_edges; tl::stable_vector m_vertex_heap; @@ -342,7 +336,7 @@ private: db::TriangleEdge *find_closest_edge (const db::DPoint &p, db::Vertex *vstart = 0, bool inside_only = false); db::Vertex *insert (db::Vertex *vertex, std::list > *new_triangles = 0); void split_triangle (db::Triangle *t, db::Vertex *vertex, std::list > *new_triangles_out); - void split_triangles_on_edge (const std::vector &tris, db::Vertex *vertex, db::TriangleEdge *split_edge, std::list > *new_triangles_out); + void split_triangles_on_edge (db::Vertex *vertex, db::TriangleEdge *split_edge, std::list > *new_triangles_out); void add_more_triangles (std::vector &new_triangles, db::TriangleEdge *incoming_edge, db::Vertex *from_vertex, db::Vertex *to_vertex, diff --git a/src/db/unit_tests/dbTrianglesTests.cc b/src/db/unit_tests/dbTrianglesTests.cc index c2209c3da..9b7dd021c 100644 --- a/src/db/unit_tests/dbTrianglesTests.cc +++ b/src/db/unit_tests/dbTrianglesTests.cc @@ -277,7 +277,6 @@ namespace { TEST(insert_many) { - return; // @@@ srand (0); TestableTriangles tris; From 1f5c2b5132ea14ffaff7ff712d1ccfde67062d46 Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Tue, 18 Mar 2025 14:44:27 +0100 Subject: [PATCH 13/17] Marked one test as long runner --- src/db/unit_tests/dbQuadTreeTests.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/db/unit_tests/dbQuadTreeTests.cc b/src/db/unit_tests/dbQuadTreeTests.cc index 8ffe12327..d92e07237 100644 --- a/src/db/unit_tests/dbQuadTreeTests.cc +++ b/src/db/unit_tests/dbQuadTreeTests.cc @@ -547,6 +547,8 @@ TEST(timing_insert) TEST(timing_lookup) { + test_is_long_runner (); + MyQuadTree tree; unsigned int n = 1000000; From d9343ee530ee228c978b6d974918a8cd190ac3f2 Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Tue, 18 Mar 2025 17:28:12 +0100 Subject: [PATCH 14/17] WIP --- src/db/db/dbTriangle.cc | 63 +++++++++++++++++---------- src/db/db/dbTriangle.h | 4 +- src/db/db/dbTriangles.cc | 10 +++-- src/db/unit_tests/dbTrianglesTests.cc | 46 ++++++++++++++++++- 4 files changed, 92 insertions(+), 31 deletions(-) diff --git a/src/db/db/dbTriangle.cc b/src/db/db/dbTriangle.cc index 0fa3018c9..afe396181 100644 --- a/src/db/db/dbTriangle.cc +++ b/src/db/db/dbTriangle.cc @@ -436,22 +436,40 @@ Triangle::bbox () const std::pair -Triangle::circumcircle () const +Triangle::circumcircle (bool *ok) const { - db::DVector v1 = *mp_v[0] - *mp_v[1]; - db::DVector v2 = *mp_v[0] - *mp_v[2]; - db::DVector n1 = db::DVector (v1.y (), -v1.x ()); - db::DVector n2 = db::DVector (v2.y (), -v2.x ()); + // see https://en.wikipedia.org/wiki/Circumcircle + // we set A=(0,0), so the formulas simplify - double p1s = v1.sq_length (); - double p2s = v2.sq_length (); + if (ok) { + *ok = true; + } - double s = db::vprod (v1, v2); - tl_assert (fabs (s) > db::epsilon); + db::DVector b = *mp_v[1] - *mp_v[0]; + db::DVector c = *mp_v[2] - *mp_v[0]; - db::DVector r = (n1 * p2s - n2 * p1s) * (0.5 / s); - db::DPoint center = *mp_v[0] + r; - double radius = r.length (); + double b2 = b.sq_length (); + double c2 = c.sq_length (); + + double sx = 0.5 * (b2 * c.y () - c2 * b.y ()); + double sy = 0.5 * (b.x () * c2 - c.x() * b2); + + double a1 = b.x() * c.y(); + double a2 = c.x() * b.y(); + double a = a1 - a2; + double a_abs = std::abs (a); + + if (a_abs < (std::abs (a1) + std::abs (a2)) * db::epsilon) { + if (ok) { + *ok = false; + return std::make_pair (db::DPoint (), 0.0); + } else { + tl_assert (false); + } + } + + double radius = sqrt (sx * sx + sy * sy) / a_abs; + db::DPoint center = *mp_v[0] + db::DVector (sx / a, sy / a); return std::make_pair (center, radius); } @@ -507,26 +525,22 @@ Triangle::common_edge (const Triangle *other) const int Triangle::contains (const db::DPoint &point) const { - // the relative distance of the point from the edge - double alpha = db::epsilon; - - double d = db::vprod (*mp_v[2] - *mp_v[0], *mp_v[1] - *mp_v[0]); - if (fabs (d) < db::epsilon * db::epsilon) { - return -1; - } + int vps = db::vprod_sign (*mp_v[2] - *mp_v[0], *mp_v[1] - *mp_v[0]); int res = 1; + const Vertex *vl = mp_v[2]; for (int i = 0; i < 3; ++i) { const Vertex *v = mp_v[i]; - double n = db::vprod (point - *vl, *v - *vl) / d; - if (n < -alpha) { + int n = db::vprod_sign (point - *vl, *v - *vl) * vps; + if (n < 0) { return -1; - } else if (n < alpha) { + } else if (n == 0) { res = 0; } vl = v; } + return res; } @@ -544,8 +558,9 @@ double Triangle::b () const { double lmin = min_edge_length (); - auto cr = circumcircle (); - return lmin / cr.second; + bool ok = false; + auto cr = circumcircle (&ok); + return ok ? lmin / cr.second : 0.0; } bool diff --git a/src/db/db/dbTriangle.h b/src/db/db/dbTriangle.h index dadcabc4a..29ff09abb 100644 --- a/src/db/db/dbTriangle.h +++ b/src/db/db/dbTriangle.h @@ -483,8 +483,10 @@ public: /** * @brief Gets the center point and radius of the circumcircle + * If ok is non-null, it will receive a boolean value indicating whether the circumcircle is valid. + * An invalid circumcircle is an indicator for a degenerated triangle with area 0 (or close to). */ - std::pair circumcircle () const; + std::pair circumcircle (bool *ok = 0) const; /** * @brief Gets the vertex opposite of the given edge diff --git a/src/db/db/dbTriangles.cc b/src/db/db/dbTriangles.cc index 246afa2e6..ae25fc574 100644 --- a/src/db/db/dbTriangles.cc +++ b/src/db/db/dbTriangles.cc @@ -937,13 +937,15 @@ Triangles::is_illegal_edge (db::TriangleEdge *edge) return false; } - auto lr = left->circumcircle (); - if (right->opposite (edge)->in_circle (lr.first, lr.second) > 0) { + bool ok = false; + + auto lr = left->circumcircle (&ok); + if (! ok || right->opposite (edge)->in_circle (lr.first, lr.second) > 0) { return true; } - auto rr = right->circumcircle(); - if (left->opposite (edge)->in_circle (rr.first, rr.second) > 0) { + auto rr = right->circumcircle(&ok); + if (! ok || left->opposite (edge)->in_circle (rr.first, rr.second) > 0) { return true; } diff --git a/src/db/unit_tests/dbTrianglesTests.cc b/src/db/unit_tests/dbTrianglesTests.cc index 9b7dd021c..a76f1e292 100644 --- a/src/db/unit_tests/dbTrianglesTests.cc +++ b/src/db/unit_tests/dbTrianglesTests.cc @@ -951,6 +951,48 @@ TEST(triangulate_thin) EXPECT_EQ (tri.check (false), true); - // for debugging: - tri.dump ("debug.gds"); + // for debugging: + // tri.dump ("debug.gds"); + + for (auto t = tri.begin (); t != tri.end (); ++t) { + EXPECT_GE (t->b (), param.min_b); + } + + EXPECT_GT (tri.num_triangles (), size_t (13000)); + EXPECT_LT (tri.num_triangles (), size_t (13200)); +} + +TEST(triangulate_issue1996) +{ + db::DPoint contour[] = { + db::DPoint (-8000, -8075), + db::DPoint (-8000, 8075), + db::DPoint (18000, 8075), + db::DPoint (18000, -8075) + }; + + db::DPolygon poly; + poly.assign_hull (contour + 0, contour + sizeof (contour) / sizeof (contour[0])); + + double dbu = 0.001; + + db::Triangles::TriangulateParameters param; + param.min_b = 0.5; + param.max_area = 500.0 * dbu * dbu; + + TestableTriangles tri; + db::DCplxTrans trans = db::DCplxTrans (dbu) * db::DCplxTrans (db::DTrans (db::DPoint () - poly.box ().center ())); + tri.triangulate (trans * poly, param); + + EXPECT_EQ (tri.check (false), true); + + // for debugging: + // tri.dump ("debug.gds"); + + for (auto t = tri.begin (); t != tri.end (); ++t) { + EXPECT_GE (t->b (), param.min_b); + } + + EXPECT_GT (tri.num_triangles (), size_t (13000)); + EXPECT_LT (tri.num_triangles (), size_t (13200)); } From 6f69efd4279488c8c0c829df92a664e2da9c85ec Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Thu, 20 Mar 2025 23:30:30 +0100 Subject: [PATCH 15/17] WIP --- src/db/db/dbTriangle.cc | 8 ++++- src/db/db/dbTriangles.cc | 78 +++++++++++++++++++++++++++++++++------- 2 files changed, 72 insertions(+), 14 deletions(-) diff --git a/src/db/db/dbTriangle.cc b/src/db/db/dbTriangle.cc index afe396181..13201bffb 100644 --- a/src/db/db/dbTriangle.cc +++ b/src/db/db/dbTriangle.cc @@ -525,7 +525,13 @@ Triangle::common_edge (const Triangle *other) const int Triangle::contains (const db::DPoint &point) const { - int vps = db::vprod_sign (*mp_v[2] - *mp_v[0], *mp_v[1] - *mp_v[0]); + auto c = *mp_v[2] - *mp_v[0]; + auto b = *mp_v[1] - *mp_v[0]; + + int vps = db::vprod_sign (c, b); + if (vps == 0) { + return db::vprod_sign (point - *mp_v[0], b) == 0 && db::vprod_sign (point - *mp_v[0], c) == 0 ? 0 : -1; + } int res = 1; diff --git a/src/db/db/dbTriangles.cc b/src/db/db/dbTriangles.cc index ae25fc574..0d49c53a8 100644 --- a/src/db/db/dbTriangles.cc +++ b/src/db/db/dbTriangles.cc @@ -36,6 +36,12 @@ namespace db { +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; +} + Triangles::Triangles () : m_is_constrained (false), m_level (0), m_id (0), m_flips (0), m_hops (0) { @@ -358,6 +364,12 @@ Triangles::insert (db::Vertex *vertex, std::list > *n // the new vertex is outside the domain if (tris.empty ()) { + // @@@ + if (m_is_constrained) { + dump("debug.gds"); // @@@ + find_triangle_for_point (*vertex); + } + // @@@ tl_assert (! m_is_constrained); insert_new_vertex (vertex, new_triangles); return vertex; @@ -365,26 +377,65 @@ Triangles::insert (db::Vertex *vertex, std::list > *n // 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) { db::TriangleEdge *e = tris.front ()->edge (i); if (e->side_of (*vertex) == 0) { - on_edges.push_back (e); + if (is_equal (*vertex, *e->v1 ()) || is_equal (*vertex, *e->v2 ())) { + on_vertex.push_back (e); + } else { + on_edges.push_back (e); + } } } - if (! on_edges.empty ()) { - if (on_edges.size () == size_t (1)) { - split_triangles_on_edge (vertex, on_edges.front (), new_triangles); - return vertex; - } else { - // the vertex is already present - tl_assert (on_edges.size () == size_t (2)); - return on_edges.front ()->common_vertex (on_edges [1]); - } + 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)); +// @@@ +auto v1 = on_edges.front()->v1(); +auto v2 = on_edges.front()->v2(); +unsigned int ns1 = 0, ns2 = 0; +for (auto e = v1->begin_edges (); e != v1->end_edges (); ++e) { + if ((*e)->is_segment()) { ++ns1; } +} +for (auto e = v2->begin_edges (); e != v2->end_edges (); ++e) { + if ((*e)->is_segment()) { ++ns2; } +} +std::string vs = vertex->to_string(); +std::string es = on_edges.front()->to_string(); +if (vs == "(-12.9999999999, 3.50126953125)" && es == "((-13, 3.5328125), (-12.9999999998, 3.4697265625))") { + printf("@@@ BANG!\n"); fflush(stdout); +} +// @@@ + split_triangles_on_edge (vertex, on_edges.front (), new_triangles); +// @@@ +unsigned int ns1p = 0, ns2p = 0; +for (auto e = v1->begin_edges (); e != v1->end_edges (); ++e) { + if ((*e)->is_segment()) { ++ns1p; } +} +for (auto e = v2->begin_edges (); e != v2->end_edges (); ++e) { + if ((*e)->is_segment()) { ++ns2p; } +} +if (ns1 != ns1p || ns2 != ns2p) { + printf("@@@ on '%s' '%s' - BANG!\n", vs.c_str(), es.c_str()); fflush(stdout); +} +tl_assert(ns1 == ns1p); +tl_assert(ns2 == ns2p); +// @@@ + 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); @@ -1178,9 +1229,9 @@ Triangles::find_vertex_for_point (const db::DPoint &point) return 0; } db::Vertex *v = 0; - if (edge->v1 ()->equal (point)) { + if (is_equal (*edge->v1 (), point)) { v = edge->v1 (); - } else if (edge->v2 ()->equal (point)) { + } else if (is_equal (*edge->v2 (), point)) { v = edge->v2 (); } return v; @@ -1194,7 +1245,7 @@ Triangles::find_edge_for_points (const db::DPoint &p1, const db::DPoint &p2) return 0; } for (auto e = v->begin_edges (); e != v->end_edges (); ++e) { - if ((*e)->other (v)->equal (p2)) { + if (is_equal (*(*e)->other (v), p2)) { return *e; } } @@ -1587,6 +1638,7 @@ Triangles::refine (const TriangulateParameters ¶meters) if (tl::verbosity () >= parameters.base_verbosity + 10) { tl::info << "Iteration " << nloop << " .."; } + printf("@@@ iteration %d ..\n", (int)nloop); fflush(stdout); std::list > to_consider; for (auto t = new_triangles.begin (); t != new_triangles.end (); ++t) { From a727ed0b1d284b39bbe832acef403969809a70b6 Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Sat, 22 Mar 2025 20:47:28 +0100 Subject: [PATCH 16/17] Fixed a compiler warning --- src/db/db/dbDeepTexts.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/db/db/dbDeepTexts.cc b/src/db/db/dbDeepTexts.cc index cabf79bfe..107a3f6ed 100644 --- a/src/db/db/dbDeepTexts.cc +++ b/src/db/db/dbDeepTexts.cc @@ -462,7 +462,7 @@ DeepTexts::apply_filter (const TextFilterBase &filter, bool with_true, bool with for (db::Shapes::shape_iterator si = s.begin (db::ShapeIterator::Texts); ! si.at_end (); ++si) { db::Text text; si->text (text); - if (filter.selected (text.transformed (*v), si->prop_id ())) { + if (filter.selected (text.transformed (tr), si->prop_id ())) { if (st_true) { st_true->insert (*si); } From 6228668fa12b40f734e11383ecbea3debf5746e4 Mon Sep 17 00:00:00 2001 From: Matthias Koefferlein Date: Sat, 22 Mar 2025 21:47:22 +0100 Subject: [PATCH 17/17] Fixing issue #1996: Providing a more robust triangulation --- src/db/db/dbTriangle.cc | 38 ++++++++++------ src/db/db/dbTriangles.cc | 62 ++++++++++----------------- src/db/unit_tests/dbTrianglesTests.cc | 14 +++--- 3 files changed, 57 insertions(+), 57 deletions(-) diff --git a/src/db/db/dbTriangle.cc b/src/db/db/dbTriangle.cc index 13201bffb..867187295 100644 --- a/src/db/db/dbTriangle.cc +++ b/src/db/db/dbTriangle.cc @@ -363,21 +363,33 @@ Triangle::Triangle (TriangleEdge *e1, TriangleEdge *e2, TriangleEdge *e3) } mp_v[2] = mp_e[1]->other (mp_v[1]); - // establish link to edges - for (int i = 0; i < 3; ++i) { - TriangleEdge *e = mp_e[i]; - int side_of = e->side_of (*mp_v[i == 0 ? 2 : i - 1]); - // NOTE: in the degenerated case, the triangle is not attached to an edge! - if (side_of < 0) { - e->set_left (this); - } else if (side_of > 0) { - e->set_right (this); - } + // 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); } - // enforce clockwise orientation - if (db::vprod_sign (*mp_v[2] - *mp_v[0], *mp_v[1] - *mp_v[0]) < 0) { - std::swap (mp_v[2], mp_v[1]); + // establish link to edges + for (int i = 0; i < 3; ++i) { + + TriangleEdge *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); + } + } } diff --git a/src/db/db/dbTriangles.cc b/src/db/db/dbTriangles.cc index 0d49c53a8..6ebbe6940 100644 --- a/src/db/db/dbTriangles.cc +++ b/src/db/db/dbTriangles.cc @@ -364,12 +364,6 @@ Triangles::insert (db::Vertex *vertex, std::list > *n // the new vertex is outside the domain if (tris.empty ()) { - // @@@ - if (m_is_constrained) { - dump("debug.gds"); // @@@ - find_triangle_for_point (*vertex); - } - // @@@ tl_assert (! m_is_constrained); insert_new_vertex (vertex, new_triangles); return vertex; @@ -397,37 +391,7 @@ Triangles::insert (db::Vertex *vertex, std::list > *n } else if (! on_edges.empty ()) { tl_assert (on_edges.size () == size_t (1)); -// @@@ -auto v1 = on_edges.front()->v1(); -auto v2 = on_edges.front()->v2(); -unsigned int ns1 = 0, ns2 = 0; -for (auto e = v1->begin_edges (); e != v1->end_edges (); ++e) { - if ((*e)->is_segment()) { ++ns1; } -} -for (auto e = v2->begin_edges (); e != v2->end_edges (); ++e) { - if ((*e)->is_segment()) { ++ns2; } -} -std::string vs = vertex->to_string(); -std::string es = on_edges.front()->to_string(); -if (vs == "(-12.9999999999, 3.50126953125)" && es == "((-13, 3.5328125), (-12.9999999998, 3.4697265625))") { - printf("@@@ BANG!\n"); fflush(stdout); -} -// @@@ split_triangles_on_edge (vertex, on_edges.front (), new_triangles); -// @@@ -unsigned int ns1p = 0, ns2p = 0; -for (auto e = v1->begin_edges (); e != v1->end_edges (); ++e) { - if ((*e)->is_segment()) { ++ns1p; } -} -for (auto e = v2->begin_edges (); e != v2->end_edges (); ++e) { - if ((*e)->is_segment()) { ++ns2p; } -} -if (ns1 != ns1p || ns2 != ns2p) { - printf("@@@ on '%s' '%s' - BANG!\n", vs.c_str(), es.c_str()); fflush(stdout); -} -tl_assert(ns1 == ns1p); -tl_assert(ns2 == ns2p); -// @@@ return vertex; } else if (tris.size () == size_t (1)) { @@ -1638,7 +1602,6 @@ Triangles::refine (const TriangulateParameters ¶meters) if (tl::verbosity () >= parameters.base_verbosity + 10) { tl::info << "Iteration " << nloop << " .."; } - printf("@@@ iteration %d ..\n", (int)nloop); fflush(stdout); std::list > to_consider; for (auto t = new_triangles.begin (); t != new_triangles.end (); ++t) { @@ -1667,7 +1630,28 @@ Triangles::refine (const TriangulateParameters ¶meters) auto cr = (*t)->circumcircle(); auto center = cr.first; - if ((*t)->contains (center) >= 0) { + 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); @@ -1677,7 +1661,7 @@ Triangles::refine (const TriangulateParameters ¶meters) } else { db::Vertex *vstart = 0; - for (int i = 0; i < 3; ++i) { + for (unsigned int i = 0; i < 3; ++i) { db::TriangleEdge *edge = (*t)->edge (i); vstart = (*t)->opposite (edge); if (edge->side_of (*vstart) * edge->side_of (center) < 0) { diff --git a/src/db/unit_tests/dbTrianglesTests.cc b/src/db/unit_tests/dbTrianglesTests.cc index a76f1e292..c3855ce29 100644 --- a/src/db/unit_tests/dbTrianglesTests.cc +++ b/src/db/unit_tests/dbTrianglesTests.cc @@ -700,6 +700,9 @@ TEST(triangulate_basic) EXPECT_GT (tri.num_triangles (), size_t (100)); EXPECT_LT (tri.num_triangles (), size_t (150)); + // for debugging: + // tri.dump ("debug.gds"); + param.min_b = 1.0; param.max_area = 0.1; @@ -914,8 +917,8 @@ TEST(triangulate_problematic) EXPECT_GE (t->b (), param.min_b); } - EXPECT_GT (tri.num_triangles (), size_t (470)); - EXPECT_LT (tri.num_triangles (), size_t (490)); + EXPECT_GT (tri.num_triangles (), size_t (540)); + EXPECT_LT (tri.num_triangles (), size_t (560)); } TEST(triangulate_thin) @@ -978,7 +981,7 @@ TEST(triangulate_issue1996) db::Triangles::TriangulateParameters param; param.min_b = 0.5; - param.max_area = 500.0 * dbu * dbu; + param.max_area = 5000.0 * dbu * dbu; TestableTriangles tri; db::DCplxTrans trans = db::DCplxTrans (dbu) * db::DCplxTrans (db::DTrans (db::DPoint () - poly.box ().center ())); @@ -990,9 +993,10 @@ TEST(triangulate_issue1996) // tri.dump ("debug.gds"); for (auto t = tri.begin (); t != tri.end (); ++t) { + EXPECT_LE (t->area (), param.max_area); EXPECT_GE (t->b (), param.min_b); } - EXPECT_GT (tri.num_triangles (), size_t (13000)); - EXPECT_LT (tri.num_triangles (), size_t (13200)); + EXPECT_GT (tri.num_triangles (), size_t (128000)); + EXPECT_LT (tri.num_triangles (), size_t (132000)); }