From 6a58b01e8b53e69ad3b602488a833ce7f265ac1a Mon Sep 17 00:00:00 2001 From: Thomas Christiansen Date: Fri, 27 Jun 2025 10:51:10 -0600 Subject: [PATCH 1/4] Add a method to GmPtSearch that lets you find points within a box. --- xmsgrid/geometry/GmPtSearch.cpp | 86 ++++++++++++++++++++++++++++----- xmsgrid/geometry/GmPtSearch.h | 1 + xmsgrid/geometry/GmPtSearch.t.h | 1 + 3 files changed, 76 insertions(+), 12 deletions(-) diff --git a/xmsgrid/geometry/GmPtSearch.cpp b/xmsgrid/geometry/GmPtSearch.cpp index 00b8d510..55584254 100644 --- a/xmsgrid/geometry/GmPtSearch.cpp +++ b/xmsgrid/geometry/GmPtSearch.cpp @@ -193,6 +193,8 @@ class GmPtSearchImpl : public GmPtSearch const Pt3d& a_pt, double a_dist, std::vector& a_nearest) const override; + virtual void PtsInBoxInRtree(const Pt3d& a_min, const Pt3d& a_max, std::vector& a_nearest) const override; + void PtsInBoxInRtree(int a_ptIdx, const Pt3d& a_min, const Pt3d& a_max, std::vector& a_nearest) const; virtual void NearestPtsToPt(const Pt3d& a_pt, int a_numPtsToFind, @@ -501,21 +503,46 @@ bool GmPtSearchImpl::PtInRTree(const Pt3d& a_pt, const double a_tol) return !n.empty(); } // GmPtSearchImpl::PtInRTree //------------------------------------------------------------------------------ -/// \brief Finds the nearest points to the input a_pt with an option to search -/// quadrants/octants. Similar to NearestPtsToPt but this method will always -/// pass the fSatisfies class to the RTree. This method is used to find the -/// nearest points to a point that is included in the RTree. -/// \param a_ptIdx The index of the point a_pt. a_pt is in the RTree and we -/// don't want it in the results. -/// \param a_pt The location we are interested in. -/// \param a_distance The distance from the point where we want to find -/// additional points -/// \param a_nearest Vector that is filled in by this method. +/// \brief Finds points in the RTree that are active and near a provided point. +/// This method defines distance in terms of a box rather than a circle. +/// It works by creating a square box centered on the point, with all edges +/// a_distance from the point. It then finds all the points in that box, +/// and filters out the point with the provided index from the result. This +/// means that all points within an a_distance radius of a_pt are found, +/// but there may be some extras near the corners that are further than +/// a_distance from a_pt. +/// \param a_ptIdx The index of a point that should be excluded from the +/// result. Typically the index in the RTree where a_pt appears. If -1, no +/// points will be filtered from the result, and a_pt will be included. +/// \param a_pt The location to find points near. +/// \param a_distance The distance from a_pt to look for other points. +/// \param a_nearest Will be filled in with the indexes of any found points. +/// a_ptIdx will never be included in this result. //------------------------------------------------------------------------------ void GmPtSearchImpl::PtsWithinDistanceToPtInRtree(int a_ptIdx, const Pt3d& a_pt, double a_distance, std::vector& a_nearest) const +{ + Pt3d bMin(a_pt - a_distance), bMax(a_pt + a_distance); + PtsInBoxInRtree(a_ptIdx, bMin, bMax, a_nearest); +} // GmPtSearchImpl::PtsWithinDistanceToPtInRtree + +//------------------------------------------------------------------------------ +/// \brief Finds points in the RTree that are active and in or on a given box. +/// +/// This overload is mainly to give the other overload and +/// PtsWithinDistanceToPtInRtree a common implementation and reduce +/// code duplication. It isn't expected to be useful externally. +/// +/// \param a_ptIdx The index of a point that should be excluded. This index +/// will not be inserted into a_nearest. +/// \param a_pt The location to find points near. +/// \param a_distance The distance from a_pt to look for other points. +/// \param a_nearest Will be filled in with the indexes of any found points. +/// a_ptIdx will never be included in this result. +//------------------------------------------------------------------------------ +void GmPtSearchImpl::PtsInBoxInRtree(int a_ptIdx, const Pt3d& a_min, const Pt3d& a_max, std::vector& a_nearest) const { fSatisfies fsat(m_rTree->size()); if (!m_activity.empty()) @@ -523,7 +550,7 @@ void GmPtSearchImpl::PtsWithinDistanceToPtInRtree(int a_ptIdx, if (a_ptIdx > -1) fsat.m_bits.set(a_ptIdx); - Pt3d bMin(a_pt - a_distance), bMax(a_pt + a_distance); + Pt3d bMin = a_min, bMax = a_max; if (m_2dSearch) { bMin.z = -1; @@ -541,7 +568,20 @@ void GmPtSearchImpl::PtsWithinDistanceToPtInRtree(int a_ptIdx, { a_nearest.push_back((int)nearest[i]); } -} // GmPtSearchImpl::PtsWithinDistanceToPtInRtree +} // GmPtSearchImpl::PtsInBoxInRtree + +//------------------------------------------------------------------------------ +/// \brief Finds points in the RTree that are active and in or on a given box. +/// \param a_pt The location to find points near. +/// \param a_distance The distance from a_pt to look for other points. +/// \param a_nearest Will be filled in with the indexes of any found points. +/// a_ptIdx will never be included in this result. +//------------------------------------------------------------------------------ +void GmPtSearchImpl::PtsInBoxInRtree(const Pt3d& a_min, const Pt3d& a_max, std::vector& a_nearest) const +{ + PtsInBoxInRtree(-1, a_min, a_max, a_nearest); +} // GmPtSearchImpl::PtsInBoxInRtree + //------------------------------------------------------------------------------ /// \brief Sets activity on the points in the Rtree so that points can be /// ignored when interpolating. @@ -1048,6 +1088,28 @@ void PtSearchUnitTests::testPtsWithinDist() p.PtsWithinDistanceToPtInRtree(2, (*pts)[2], 1.0, n); TS_ASSERT_EQUALS_VEC(base, n); } + +//------------------------------------------------------------------------------ +/// \brief testing Point in a box. +//------------------------------------------------------------------------------ +void PtSearchUnitTests::testPtsInBox() +{ + BSHP> pts(new std::vector()); + *pts = {{0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {1, 1, 0}, {1, -1, 0}}; + + GmPtSearchImpl p(true); + p.PtsToSearch(pts); + + std::vector n; + p.PtsInBoxInRtree(Pt3d(-50.0, -100.0), Pt3d(-10.0, -10.0), n); + std::vector base = {}; + TS_ASSERT_EQUALS_VEC(base, n); + + p.PtsInBoxInRtree(Pt3d(0.0, -1.0), Pt3d(1.0, 0.0), n); + base = {0, 1, 4}; + TS_ASSERT_EQUALS_VEC(base, n); +} // PtSearchUnitTests::testPtsInBox + //------------------------------------------------------------------------------ /// \brief testing VectorThatGrows functionality //------------------------------------------------------------------------------ diff --git a/xmsgrid/geometry/GmPtSearch.h b/xmsgrid/geometry/GmPtSearch.h index f594a8df..a5d6de94 100644 --- a/xmsgrid/geometry/GmPtSearch.h +++ b/xmsgrid/geometry/GmPtSearch.h @@ -52,6 +52,7 @@ class GmPtSearch const Pt3d& a_pt, double a_dist, VecInt& a_nearest) const = 0; + virtual void PtsInBoxInRtree(const Pt3d& a_min, const Pt3d& a_max, std::vector& a_nearest) const = 0; virtual void SetActivity(boost::dynamic_bitset& a_activity) = 0; virtual boost::dynamic_bitset GetActivity() = 0; diff --git a/xmsgrid/geometry/GmPtSearch.t.h b/xmsgrid/geometry/GmPtSearch.t.h index 38bd1989..3aa84616 100644 --- a/xmsgrid/geometry/GmPtSearch.t.h +++ b/xmsgrid/geometry/GmPtSearch.t.h @@ -35,6 +35,7 @@ class PtSearchUnitTests : public CxxTest::TestSuite void testActivity2d(); void testActivity3d(); void testPtsWithinDist(); + void testPtsInBox(); void testVectorThatGrows(); }; //----- Function prototypes ---------------------------------------------------- From a07f38fee7bbbccad2c762b10d9da8ee44a534d9 Mon Sep 17 00:00:00 2001 From: Thomas Christiansen Date: Thu, 17 Jul 2025 13:54:47 -0600 Subject: [PATCH 2/4] Add a point in box 3d test. --- xmsgrid/geometry/GmPtSearch.cpp | 31 ++++++++++++++++++++++++++++++- xmsgrid/geometry/GmPtSearch.t.h | 1 + 2 files changed, 31 insertions(+), 1 deletion(-) diff --git a/xmsgrid/geometry/GmPtSearch.cpp b/xmsgrid/geometry/GmPtSearch.cpp index 55584254..134e5148 100644 --- a/xmsgrid/geometry/GmPtSearch.cpp +++ b/xmsgrid/geometry/GmPtSearch.cpp @@ -1095,7 +1095,7 @@ void PtSearchUnitTests::testPtsWithinDist() void PtSearchUnitTests::testPtsInBox() { BSHP> pts(new std::vector()); - *pts = {{0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {1, 1, 0}, {1, -1, 0}}; + *pts = {{0, 0, -1}, {1, 0, 1}, {2, 0, 0}, {1, 1, 0}, {1, -1, 0}}; GmPtSearchImpl p(true); p.PtsToSearch(pts); @@ -1108,6 +1108,35 @@ void PtSearchUnitTests::testPtsInBox() p.PtsInBoxInRtree(Pt3d(0.0, -1.0), Pt3d(1.0, 0.0), n); base = {0, 1, 4}; TS_ASSERT_EQUALS_VEC(base, n); + + p.PtsInBoxInRtree(Pt3d(1.0, 0.0), Pt3d(0.0, -1.0), n); + base = {0, 1, 4}; + TS_ASSERT_EQUALS_VEC(base, n); +} // PtSearchUnitTests::testPtsInBox + +//------------------------------------------------------------------------------ +/// \brief testing Point in a box. +//------------------------------------------------------------------------------ +void PtSearchUnitTests::testPtsInBox3d() +{ + BSHP> pts(new std::vector()); + *pts = {{0, 0, -1}, {1, 0, 1}, {2, 0, 0}, {1, 1, 0}, {1, -1, 0}}; + + GmPtSearchImpl p(false); + p.PtsToSearch(pts); + + std::vector n; + p.PtsInBoxInRtree(Pt3d(-50.0, -100.0), Pt3d(-10.0, -10.0), n); + std::vector base = {}; + TS_ASSERT_EQUALS_VEC(base, n); + + p.PtsInBoxInRtree(Pt3d(0.0, -1.0, -1.0), Pt3d(1.0, 0.0, 1.0), n); + base = {0, 1, 4}; + TS_ASSERT_EQUALS_VEC(base, n); + + p.PtsInBoxInRtree(Pt3d(0.0, -1.0, 0.0), Pt3d(1.0, 0.0, 0.0), n); + base = {4}; + TS_ASSERT_EQUALS_VEC(base, n); } // PtSearchUnitTests::testPtsInBox //------------------------------------------------------------------------------ diff --git a/xmsgrid/geometry/GmPtSearch.t.h b/xmsgrid/geometry/GmPtSearch.t.h index 3aa84616..e522476e 100644 --- a/xmsgrid/geometry/GmPtSearch.t.h +++ b/xmsgrid/geometry/GmPtSearch.t.h @@ -36,6 +36,7 @@ class PtSearchUnitTests : public CxxTest::TestSuite void testActivity3d(); void testPtsWithinDist(); void testPtsInBox(); + void testPtsInBox3d(); void testVectorThatGrows(); }; //----- Function prototypes ---------------------------------------------------- From e54c868313ba2541dce6debaf974b35f64ded746 Mon Sep 17 00:00:00 2001 From: Thomas Christiansen Date: Tue, 22 Jul 2025 14:31:58 -0600 Subject: [PATCH 3/4] Fix issue 0015785: HydroAS arcs snap weirdly --- test_files/bug-0015785.xmugrid | 86 ++++++ .../GmMultiPolyIntersectionSorterTerse.cpp | 52 +++- .../GmMultiPolyIntersectionSorterTerse.h | 1 + xmsgrid/geometry/GmMultiPolyIntersector.cpp | 289 +++++++++++------- xmsgrid/geometry/GmMultiPolyIntersector.t.h | 1 + 5 files changed, 318 insertions(+), 111 deletions(-) create mode 100644 test_files/bug-0015785.xmugrid diff --git a/test_files/bug-0015785.xmugrid b/test_files/bug-0015785.xmugrid new file mode 100644 index 00000000..9b2d01d0 --- /dev/null +++ b/test_files/bug-0015785.xmugrid @@ -0,0 +1,86 @@ +ASCII XmUGrid Version 2 +LOCATIONS 48 + POINT 0 -38528.08 254321.0 769.744019 + POINT 1 -38528.0 254321.25 769.744019 + POINT 2 -38528.0 254321.0 769.744019 + POINT 3 -36899.68 262325.09 769.744019 + POINT 4 -36892.32 262327.386 769.744019 + POINT 5 -36897.87 262322.32 769.744019 + POINT 6 -40659.55 247159.57 675.431478 + POINT 7 -40665.6009 247157.7 675.431478 + POINT 8 -40663.43 247164.9 675.431478 + POINT 9 -38526.34 254328.197 723.4776 + POINT 10 -38526.568 254330.333 723.418762 + POINT 11 -38525.608 254328.434 723.4776 + POINT 12 -38526.757 254333.772 724.232239 + POINT 13 -38525.3923 254328.592 724.625953 + POINT 14 -38525.534 254330.438 723.418762 + POINT 15 -38525.709 254333.817 724.399628 + POINT 16 -38524.6954 254328.775 724.307859 + POINT 17 -38524.7269 254330.534 724.372841 + POINT 18 -38524.7873 254333.899 724.49722 + POINT 19 -38523.9986 254328.958 723.989764 + POINT 20 -38523.9341 254330.627 723.935229 + POINT 21 -38524.2333 254333.949 724.219912 + POINT 22 -38523.3775 254329.029 723.866525 + POINT 23 -38523.1413 254330.721 723.497616 + POINT 24 -38523.6793 254333.999 723.942604 + POINT 25 -38522.7565 254329.101 723.743286 + POINT 26 -38522.3485 254330.815 723.060004 + POINT 27 -38523.1339 254334.047 723.512122 + POINT 28 -38522.1354 254329.172 723.620047 + POINT 29 -38521.5557 254330.909 722.622391 + POINT 30 -38522.6326 254334.092 723.081258 + POINT 31 -38520.7629 254331.003 722.184779 + POINT 32 -38520.4547 254329.366 721.993129 + POINT 33 -38521.3661 254334.206 722.559801 + POINT 34 -38519.8416 254331.112 721.720529 + POINT 35 -38519.1657 254329.513 721.800792 + POINT 36 -38520.2036 254334.309 722.415745 + POINT 37 -38518.552 254331.353 721.256279 + POINT 38 -38517.8767 254329.659 721.608454 + POINT 39 -38519.0418 254334.413 722.27193 + POINT 40 -38517.0778 254331.439 720.32778 + POINT 41 -38516.83 254329.95 720.0238 + POINT 42 -38517.5943 254334.543 720.961488 + POINT 43 -38515.934 254331.574 720.144667 + POINT 44 -38515.064 254330.029 719.84848 + POINT 45 -38516.3064 254334.658 720.592489 + POINT 46 -38514.7903 254331.71 719.961554 + POINT 47 -38515.0185 254334.773 720.223491 +CELL_STREAM 191 + CELL 0 TRIANGLE 3 1 0 2 + CELL 1 TRIANGLE 3 4 3 5 + CELL 2 TRIANGLE 3 7 6 8 + CELL 3 QUAD 4 9 11 14 10 + CELL 4 QUAD 4 12 10 14 15 + CELL 5 TRIANGLE 3 11 13 14 + CELL 6 QUAD 4 16 17 14 13 + CELL 7 QUAD 4 17 18 15 14 + CELL 8 QUAD 4 19 20 17 16 + CELL 9 TRIANGLE 3 20 18 17 + CELL 10 TRIANGLE 3 20 21 18 + CELL 11 QUAD 4 20 19 22 23 + CELL 12 QUAD 4 24 20 23 27 + CELL 13 TRIANGLE 3 20 24 21 + CELL 14 TRIANGLE 3 25 23 22 + CELL 15 QUAD 4 26 23 25 28 + CELL 16 QUAD 4 27 23 26 30 + CELL 17 TRIANGLE 3 29 26 28 + CELL 18 TRIANGLE 3 30 26 29 + CELL 19 TRIANGLE 3 29 28 32 + CELL 20 TRIANGLE 3 31 29 32 + CELL 21 TRIANGLE 3 29 31 33 + CELL 22 TRIANGLE 3 29 33 30 + CELL 23 QUAD 4 32 35 34 31 + CELL 24 QUAD 4 36 33 31 34 + CELL 25 QUAD 4 35 38 37 34 + CELL 26 QUAD 4 36 34 37 39 + CELL 27 TRIANGLE 3 38 40 37 + CELL 28 QUAD 4 39 37 40 42 + CELL 29 TRIANGLE 3 40 38 41 + CELL 30 TRIANGLE 3 43 40 41 + CELL 31 QUAD 4 45 42 40 43 + CELL 32 TRIANGLE 3 43 41 44 + CELL 33 TRIANGLE 3 44 46 43 + CELL 34 QUAD 4 45 43 46 47 diff --git a/xmsgrid/geometry/GmMultiPolyIntersectionSorterTerse.cpp b/xmsgrid/geometry/GmMultiPolyIntersectionSorterTerse.cpp index 72a97624..b518e2bb 100644 --- a/xmsgrid/geometry/GmMultiPolyIntersectionSorterTerse.cpp +++ b/xmsgrid/geometry/GmMultiPolyIntersectionSorterTerse.cpp @@ -70,8 +70,9 @@ void GmMultiPolyIntersectionSorterTerse::Sort( RemoveCornerTouches(); RemoveDuplicateEdges(); + RemoveIntersectionsWithoutMatch(); SwapAdjacents(); - RemoveCornerTouches(); + IntersectionsToPolyIdsAndTValues(polyids, tvalues, a_pts); FixArrays(polyids, tvalues, a_pts); } // GmMultiPolyIntersectionSorterTerse::Sort @@ -264,6 +265,55 @@ void GmMultiPolyIntersectionSorterTerse::SwapAdjacents() { } } } // GmMultiPolyIntersectionSorterTerse::SwapAdjacents +//------------------------------------------------------------------------------ +/// \brief Remove an entry/exit without a matching exit/entry. +/// See test GmMultiPolyIntersector2IntermediateTests::testBug15785. +//------------------------------------------------------------------------------ +void GmMultiPolyIntersectionSorterTerse::RemoveIntersectionsWithoutMatch() +{ + VecBool hasMatch(m_d->m_ixs.size(), false); + if (hasMatch.size() < 2) + { + return; + } + + for (int i = 0; i < m_d->m_ixs.size(); i++) + { + if (hasMatch[i]) + { + continue; + } + + double t = m_d->m_ixs[i].m_t; + int cell = m_d->m_ixs[i].m_i; + + for (int j = i + 1; j < m_d->m_ixs.size(); j++) + { + // If it's already paired, or not part of a future block. + if (hasMatch[j] || m_d->m_ixs[j].m_t <= t) + { + continue; + } + + if (m_d->m_ixs[j].m_i == cell) + { + hasMatch[i] = hasMatch[j] = true; + break; + } + } + } + + std::vector oldIx = m_d->m_ixs; + m_d->m_ixs.clear(); + for (int i = 0; i < oldIx.size(); i++) + { + if (hasMatch[i]) + { + m_d->m_ixs.push_back(oldIx[i]); + } + } +} // GmMultiPolyIntersectionSorterTerse::RemoveIntersectionsWithoutMatch + //------------------------------------------------------------------------------ /// \brief Add endpoint polygon IDs that may be removed as duplicates. /// \param a_tChange: indexes of t-value changes. diff --git a/xmsgrid/geometry/GmMultiPolyIntersectionSorterTerse.h b/xmsgrid/geometry/GmMultiPolyIntersectionSorterTerse.h index bbb7e5ed..0cf169a8 100644 --- a/xmsgrid/geometry/GmMultiPolyIntersectionSorterTerse.h +++ b/xmsgrid/geometry/GmMultiPolyIntersectionSorterTerse.h @@ -40,6 +40,7 @@ class GmMultiPolyIntersectionSorterTerse void RemoveCornerTouches(); void RemoveDuplicateEdges(); void SwapAdjacents(); + void RemoveIntersectionsWithoutMatch(); void AddMissingEndpointIds(const std::vector& a_tChange); void IntersectionsToPolyIdsAndTValuesFor2(std::vector &polyids, std::vector &tvalues, diff --git a/xmsgrid/geometry/GmMultiPolyIntersector.cpp b/xmsgrid/geometry/GmMultiPolyIntersector.cpp index e38434f9..985ecbb1 100644 --- a/xmsgrid/geometry/GmMultiPolyIntersector.cpp +++ b/xmsgrid/geometry/GmMultiPolyIntersector.cpp @@ -2330,8 +2330,8 @@ void GmMultiPolyIntersector2IntermediateTests::testBug13273() VecPt3d locs = xmGrid->GetLocations(); for (auto& p : locs) p.z = 0.0; - iRunTest(8494183.82, 822079.23, 8494283.14, 822106.67, - locs, cellPolys, expectedIds, expectedTvals, expectedPoints); + iRunTest(8494183.82, 822079.23, 8494283.14, 822106.67, locs, cellPolys, expectedIds, + expectedTvals, expectedPoints); } // GmMultiPolyIntersector2IntermediateTests::testBug13273 //------------------------------------------------------------------------------ @@ -2358,74 +2358,58 @@ void GmMultiPolyIntersector2IntermediateTests::testPointOnPolygonVertex() } } - VecPt3d segmentPoints = { - {380.0, 0.0, 0.0}, - {361.3333333333333, 6.666666666666667, 0.0}, - {342.6666666666667, 13.333333333333334, 0.0}, - {324.0, 20.0, 0.0}, - {305.3333333333333, 26.666666666666668, 0.0}, - {286.6666666666667, 33.33333333333333, 0.0}, - {268.0, 40.0, 0.0}, - {249.33333333333334, 46.666666666666664, 0.0}, - {230.66666666666666, 53.333333333333336, 0.0}, - {212.0, 60.0, 0.0}, - {193.33333333333334, 66.66666666666666, 0.0}, - {174.66666666666669, 73.33333333333333, 0.0}, - {156.00000000000003, 80.0, 0.0}, - {137.3333333333334, 86.66666666666664, 0.0}, - {118.66666666666669, 93.33333333333333, 0.0}, - {100.0, 100.0, 0.0} - }; - - VecInt2d expectedPolyIds = { - {167, -1}, - {166, -1}, - {165, -1}, - {164, -1}, - {163, -1}, - {162, -1}, - {161, -1}, - {160, -1}, - {159, -1}, - {158, -1}, - {157, -1}, - {156, -1}, - {155, -1}, - {308, -1}, - {308, -1} - }; - - VecPt3d2d expectedPoints = { - {{380, 0.0, 0.0}, {361.333, 6.66667, 0.0}}, - {{361.333, 6.66667, 0.0}, {342.667, 13.3333, 0.0}}, - {{342.667, 13.3333, 0.0}, {324, 20, 0.0}}, - {{324, 20, 0.0}, {305.333, 26.6667, 0.0}}, - {{305.333, 26.6667, 0.0}, {286.667, 33.3333, 0.0}}, - {{286.667, 33.3333, 0.0}, {268, 40, 0.0}}, - {{268, 40, 0.0}, {249.333,46.6667, 0.0}}, - {{249.333, 46.6667, 0.0}, {230.667, 53.3333, 0.0}}, - {{230.667, 53.3333, 0.0}, {212, 60, 0.0}}, - {{212, 60, 0.0}, {193.333, 66.6667, 0.0}}, - {{193.333, 66.6667, 0.0}, {174.667, 73.3333, 0.0}}, - {{174.667, 73.3333, 0.0}, {156, 80, 0.0}}, - {{156, 80, 0.0}, {137.333, 86.6667, 0.0}}, - {{137.333,86.6667,0.0}, {118.667, 93.3333, 0.0}}, - {{118.667, 93.3333, 0.0}, {117.733, 93.6667, 0.0}} - }; + VecPt3d segmentPoints = {{380.0, 0.0, 0.0}, + {361.3333333333333, 6.666666666666667, 0.0}, + {342.6666666666667, 13.333333333333334, 0.0}, + {324.0, 20.0, 0.0}, + {305.3333333333333, 26.666666666666668, 0.0}, + {286.6666666666667, 33.33333333333333, 0.0}, + {268.0, 40.0, 0.0}, + {249.33333333333334, 46.666666666666664, 0.0}, + {230.66666666666666, 53.333333333333336, 0.0}, + {212.0, 60.0, 0.0}, + {193.33333333333334, 66.66666666666666, 0.0}, + {174.66666666666669, 73.33333333333333, 0.0}, + {156.00000000000003, 80.0, 0.0}, + {137.3333333333334, 86.66666666666664, 0.0}, + {118.66666666666669, 93.33333333333333, 0.0}, + {100.0, 100.0, 0.0}}; + + VecInt2d expectedPolyIds = {{167, -1}, {166, -1}, {165, -1}, {164, -1}, {163, -1}, + {162, -1}, {161, -1}, {160, -1}, {159, -1}, {158, -1}, + {157, -1}, {156, -1}, {155, -1}, {308, -1}, {308, -1}}; + + VecPt3d2d expectedPoints = {{{380, 0.0, 0.0}, {361.333, 6.66667, 0.0}}, + {{361.333, 6.66667, 0.0}, {342.667, 13.3333, 0.0}}, + {{342.667, 13.3333, 0.0}, {324, 20, 0.0}}, + {{324, 20, 0.0}, {305.333, 26.6667, 0.0}}, + {{305.333, 26.6667, 0.0}, {286.667, 33.3333, 0.0}}, + {{286.667, 33.3333, 0.0}, {268, 40, 0.0}}, + {{268, 40, 0.0}, {249.333, 46.6667, 0.0}}, + {{249.333, 46.6667, 0.0}, {230.667, 53.3333, 0.0}}, + {{230.667, 53.3333, 0.0}, {212, 60, 0.0}}, + {{212, 60, 0.0}, {193.333, 66.6667, 0.0}}, + {{193.333, 66.6667, 0.0}, {174.667, 73.3333, 0.0}}, + {{174.667, 73.3333, 0.0}, {156, 80, 0.0}}, + {{156, 80, 0.0}, {137.333, 86.6667, 0.0}}, + {{137.333, 86.6667, 0.0}, {118.667, 93.3333, 0.0}}, + {{118.667, 93.3333, 0.0}, {117.733, 93.6667, 0.0}}}; VecDbl2d expectedTValues = {}; expectedTValues.resize(15, {0.0, 1.0}); expectedTValues[14] = {0.0, 0.05}; - BSHP sorter = - BSHP(new GmMultiPolyIntersectionSorterTerse()); - BSHP mpi = GmMultiPolyIntersector::New(xmGrid->GetLocations(), cellPolys, sorter, 0); + BSHP sorter = + BSHP(new GmMultiPolyIntersectionSorterTerse()); + BSHP mpi = + GmMultiPolyIntersector::New(xmGrid->GetLocations(), cellPolys, sorter, 0); for (size_t i = 0; i < segmentPoints.size() - 1; ++i) { VecInt polyIds; VecDbl tValues; VecPt3d points; - mpi->TraverseLineSegment(segmentPoints[i].x, segmentPoints[i].y, segmentPoints[i + 1].x, segmentPoints[i + 1].y, polyIds, tValues, points); + mpi->TraverseLineSegment(segmentPoints[i].x, segmentPoints[i].y, segmentPoints[i + 1].x, + segmentPoints[i + 1].y, polyIds, tValues, points); TS_ASSERT_EQUALS(expectedPolyIds[i], polyIds); TS_ASSERT_DELTA_VECPT3D(expectedPoints[i], points, 1.0e-3); TS_ASSERT_DELTA_VEC(expectedTValues[i], tValues, 1.0e-5); @@ -2448,43 +2432,25 @@ void GmMultiPolyIntersector2IntermediateTests::testPointOnPolygonVertex() {176.00000000000003, 80.0, 0.0}, {157.3333333333334, 86.66666666666664, 0.0}, {138.66666666666669, 93.33333333333333, 0.0}, - {120.0, 100.0, 0.0} - }; - - expectedPolyIds = { - {}, - {13, -1}, - {12, -1}, - {11, -1}, - {10, -1}, - {9, -1}, - {8, -1}, - {7, -1}, - {6, -1}, - {5, -1}, - {4, -1}, - {3, -1}, - {2, -1}, - {1, -1}, - {0, -1} - }; - expectedPoints = { - {}, - {{381.333, 6.66667, 0.0}, {362.667, 13.3333, 0.0}}, - {{362.667, 13.3333, 0.0}, {344, 20, 0.0}}, - {{344, 20, 0.0}, {325.333, 26.6667, 0.0}}, - {{325.333, 26.6667, 0.0}, {306.667, 33.3333, 0.0}}, - {{306.667, 33.3333, 0.0}, {288, 40, 0.0}}, - {{288, 40, 0.0}, {269.333, 46.6667, 0.0}}, - {{269.333, 46.6667, 0.0}, {250.667, 53.3333, 0.0}}, - {{250.667, 53.3333, 0.0}, {232, 60, 0.0}}, - {{232, 60, 0.0}, {213.333, 66.6667, 0.0}}, - {{213.333, 66.6667, 0.0}, {194.667, 73.3333, 0.0}}, - {{194.667, 73.3333, 0.0}, {176, 80, 0.0}}, - {{176, 80, 0.0}, {157.333, 86.6667, 0.0}}, - {{157.333, 86.6667, 0.0}, {138.667, 93.3333, 0.0}}, - {{138.667, 93.3333, 0.0}, {120, 100, 0.0}} - }; + {120.0, 100.0, 0.0}}; + + expectedPolyIds = {{}, {13, -1}, {12, -1}, {11, -1}, {10, -1}, {9, -1}, {8, -1}, {7, -1}, + {6, -1}, {5, -1}, {4, -1}, {3, -1}, {2, -1}, {1, -1}, {0, -1}}; + expectedPoints = {{}, + {{381.333, 6.66667, 0.0}, {362.667, 13.3333, 0.0}}, + {{362.667, 13.3333, 0.0}, {344, 20, 0.0}}, + {{344, 20, 0.0}, {325.333, 26.6667, 0.0}}, + {{325.333, 26.6667, 0.0}, {306.667, 33.3333, 0.0}}, + {{306.667, 33.3333, 0.0}, {288, 40, 0.0}}, + {{288, 40, 0.0}, {269.333, 46.6667, 0.0}}, + {{269.333, 46.6667, 0.0}, {250.667, 53.3333, 0.0}}, + {{250.667, 53.3333, 0.0}, {232, 60, 0.0}}, + {{232, 60, 0.0}, {213.333, 66.6667, 0.0}}, + {{213.333, 66.6667, 0.0}, {194.667, 73.3333, 0.0}}, + {{194.667, 73.3333, 0.0}, {176, 80, 0.0}}, + {{176, 80, 0.0}, {157.333, 86.6667, 0.0}}, + {{157.333, 86.6667, 0.0}, {138.667, 93.3333, 0.0}}, + {{138.667, 93.3333, 0.0}, {120, 100, 0.0}}}; expectedTValues.assign(15, {0.0, 1.0}); expectedTValues[0] = {}; for (size_t i = 0; i < segmentPoints.size() - 1; ++i) @@ -2492,7 +2458,8 @@ void GmMultiPolyIntersector2IntermediateTests::testPointOnPolygonVertex() VecInt polyIds; VecDbl tValues; VecPt3d points; - mpi->TraverseLineSegment(segmentPoints[i].x, segmentPoints[i].y, segmentPoints[i + 1].x, segmentPoints[i + 1].y, polyIds, tValues, points); + mpi->TraverseLineSegment(segmentPoints[i].x, segmentPoints[i].y, segmentPoints[i + 1].x, + segmentPoints[i + 1].y, polyIds, tValues, points); TS_ASSERT_EQUALS(expectedPolyIds[i], polyIds); TS_ASSERT_DELTA_VECPT3D(expectedPoints[i], points, 1.0e-3); TS_ASSERT_DELTA_VEC(expectedTValues[i], tValues, 1.0e-5); @@ -2523,24 +2490,22 @@ void GmMultiPolyIntersector2IntermediateTests::testPointsNearEdgePoints() } } - VecPt3d segmentPoints = { - {958452.39285714, 498553.35714286, 7531.0}, - {958450.41714286, 498552.80285714, 7531.0} - }; + VecPt3d segmentPoints = {{958452.39285714, 498553.35714286, 7531.0}, + {958450.41714286, 498552.80285714, 7531.0}}; BSHP sorter = BSHP(new GmMultiPolyIntersectionSorterTerse()); - BSHP mpi = GmMultiPolyIntersector::New(xmGrid->GetLocations(), cellPolys, sorter, 0); + BSHP mpi = + GmMultiPolyIntersector::New(xmGrid->GetLocations(), cellPolys, sorter, 0); VecInt polyIds; VecDbl tValues; VecPt3d points; - mpi->TraverseLineSegment(segmentPoints[0].x, segmentPoints[0].y, segmentPoints[1].x, segmentPoints[1].y, polyIds, tValues, points); + mpi->TraverseLineSegment(segmentPoints[0].x, segmentPoints[0].y, segmentPoints[1].x, + segmentPoints[1].y, polyIds, tValues, points); VecInt expectedPolyIds = {1836, -1}; - VecPt3d expectedPoints = { - {958452.39285714, 498553.35714286, 0.0}, - {958450.41714286, 498552.80285714, 0.0} - }; + VecPt3d expectedPoints = {{958452.39285714, 498553.35714286, 0.0}, + {958450.41714286, 498552.80285714, 0.0}}; VecDbl expectedTValues = {0.0, 1.0}; TS_ASSERT_EQUALS(expectedPolyIds, polyIds); TS_ASSERT_DELTA_VECPT3D(expectedPoints, points, 1.0e-3); @@ -2550,10 +2515,11 @@ void GmMultiPolyIntersector2IntermediateTests::testPointsNearEdgePoints() {958436.79571429, 498597.85, 7531.0}, {958434.84428571, 498597.37, 7531.0}, }; - mpi->TraverseLineSegment(segmentPoints[0].x, segmentPoints[0].y, segmentPoints[1].x, segmentPoints[1].y, polyIds, tValues, points); + mpi->TraverseLineSegment(segmentPoints[0].x, segmentPoints[0].y, segmentPoints[1].x, + segmentPoints[1].y, polyIds, tValues, points); expectedPolyIds = {2387, -1}; expectedPoints = { - {958436.79571429, 498597.85,0.0}, + {958436.79571429, 498597.85, 0.0}, {958434.84428571, 498597.37, 0.0}, }; TS_ASSERT_EQUALS(expectedPolyIds, polyIds); @@ -2561,6 +2527,109 @@ void GmMultiPolyIntersector2IntermediateTests::testPointsNearEdgePoints() TS_ASSERT_DELTA_VEC(expectedTValues, tValues, 1.0e-5); } // GmMultiPolyIntersector2IntermediateTests::testPointsNearEdgePoints +//------------------------------------------------------------------------------ +/// \brief GmMultiPolyIntersectionSorterTerse normally produces pairs of +/// records everywhere except the start and end nodes. This case makes +/// it produce a record in the middle of the sequence with no matching +/// one, which the intersector used to respond to by associating a node +/// with a cell it wasn't even attached to. +/// +/// The interesting detail about this case is that the left side of the +/// rightmost segment intersects the element at a shallow angle. The +/// point of entry is far enough from the node that it doesn't +/// immediately snap to the node, but the point of exit is close enough +/// that it *does* snap to the node. At the scale where you can see the +/// whole area of the arc, it looks like it just intersects the cell at +/// its node, but if you zoom way in, you can see it really just cuts +/// across the very tip of the cell. +/// +/// Initial attempts at trimming the geometry to a smaller case failed, +/// and it turned out the reason why is that this bug depends on whether +/// a point of intersection is close enough to the node to snap to it. +/// It seems there are two different tolerances at play: One is affected +/// by the extents of the grid, and the other by the extents of the +/// smallest cell in the grid. If too many small cells are deleted, or +/// cells around the boundary are removed, the issue goes away because +/// the tolerance is larger. +/// +/// When opening the grid file in XMS, it may appear to be empty. This +/// is just because there are some cells far away from the area of +/// interest to keep the tolerances right. +//------------------------------------------------------------------------------ +void GmMultiPolyIntersector2IntermediateTests::testBug15785() +{ + std::string testFilesPath(XMS_TEST_PATH); + std::string inFile = testFilesPath + "bug-0015785.xmugrid"; + std::shared_ptr xmGrid = XmReadUGridFromAsciiFile(inFile); + TS_REQUIRE_NOT_NULL(xmGrid); + + xms::VecInt2d cellPolys; + cellPolys.assign(xmGrid->GetCellCount(), xms::VecInt()); + for (int cellIdx = 0; cellIdx < xmGrid->GetCellCount(); ++cellIdx) + { + xms::VecInt& polygon = cellPolys[cellIdx]; + polygon.reserve(xmGrid->GetCellEdgeCount(cellIdx)); + for (int edgeIdx = 0; edgeIdx < xmGrid->GetCellEdgeCount(cellIdx); ++edgeIdx) + { + // int idx1, idx2; + xms::XmEdge pairIdx = xmGrid->GetCellEdge(cellIdx, edgeIdx); + polygon.push_back(pairIdx.GetFirst()); + } + } + + VecPt3d segmentPoints = {{-38517.0778, 254331.439, 720.32778}, + {-38520.7629, 254331.003, 722.184779}}; + + BSHP sorter = + BSHP(new GmMultiPolyIntersectionSorterTerse()); + BSHP mpi = + GmMultiPolyIntersector::New(xmGrid->GetLocations(), cellPolys, sorter, 0); + VecInt polyIds; + VecDbl tValues; + VecPt3d points; + mpi->TraverseLineSegment(segmentPoints[0].x, segmentPoints[0].y, segmentPoints[1].x, + segmentPoints[1].y, polyIds, tValues, points); + + VecInt expectedPolyIds = {27, 25, 26, 24, -1}; + VecPt3d expectedPoints = { + {-38517.077799999999, 254331.43900000001}, {-38518.518340051305, 254331.26856352275}, + {-38519.841556860658, 254331.11200806187}, {-38519.841600000000, 254331.11199999999}, + {-38520.762900000002, 254331.00300000000}, + }; + VecDbl expectedTValues = {-0.0000000000000000, 0.39090935152470824, 0.74998150949921127, + 0.74999330957784383, 1.0000000000000000}; + TS_ASSERT_EQUALS(expectedPolyIds, polyIds); + TS_ASSERT_DELTA_VECPT3D(expectedPoints, points, 1.0e-3); + TS_ASSERT_DELTA_VEC(expectedTValues, tValues, 1.0e-5); + + segmentPoints = { + {-38520.7629, 254331.003, 722.184779}, + {-38524.7269, 254330.534, 724.372841}, + }; + mpi->TraverseLineSegment(segmentPoints[0].x, segmentPoints[0].y, segmentPoints[1].x, + segmentPoints[1].y, polyIds, tValues, points); + expectedPolyIds = {21, 22, 18, 16, 12, 13, 10, 9, -1}; + expectedPoints = { + {-38520.762900000002, 254331.00300000000, 0.0000000000000000}, + {-38521.555688419845, 254330.90920137012, 0.0000000000000000}, + {-38521.555765061355, 254330.90919230226, 0.0000000000000000}, + {-38522.348534325967, 254330.81539593873, 0.0000000000000000}, + {-38523.141298664712, 254330.72160015799, 0.0000000000000000}, + {-38523.934039003907, 254330.62780721675, 0.0000000000000000}, + {-38523.934171293266, 254330.62779156497, 0.0000000000000000}, + {-38523.934302363145, 254330.62777605746, 0.0000000000000000}, + {-38524.726900000001, 254330.53400000001, 0.0000000000000000}, + }; + expectedTValues = { + -0.0000000000000000, 0.19999707866855634, 0.20001641305637374, + 0.40000865942642250, 0.59999966314584430, 0.79998461248892860, + 0.80001798518236922, 0.80005105023759870, 1.0000000000000000, + }; + TS_ASSERT_EQUALS(expectedPolyIds, polyIds); + TS_ASSERT_DELTA_VECPT3D(expectedPoints, points, 1.0e-3); + TS_ASSERT_DELTA_VEC(expectedTValues, tValues, 1.0e-5); +} // GmMultiPolyIntersector2IntermediateTests::testBug15785 + //} // namespace xms #endif \ No newline at end of file diff --git a/xmsgrid/geometry/GmMultiPolyIntersector.t.h b/xmsgrid/geometry/GmMultiPolyIntersector.t.h index c851d2ac..038cd4ef 100644 --- a/xmsgrid/geometry/GmMultiPolyIntersector.t.h +++ b/xmsgrid/geometry/GmMultiPolyIntersector.t.h @@ -93,6 +93,7 @@ class GmMultiPolyIntersector2IntermediateTests : public CxxTest::TestSuite void testBug13273(); void testPointOnPolygonVertex(); void testPointsNearEdgePoints(); + void testBug15785(); }; //} // namespace xms From fb0ca5deaefd123e9a1ba4fc3e2ff9d5e138f00a Mon Sep 17 00:00:00 2001 From: Thomas Christiansen Date: Tue, 22 Jul 2025 14:39:36 -0600 Subject: [PATCH 4/4] Update a test baseline. --- xmsgrid/geometry/GmPtSearch.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/xmsgrid/geometry/GmPtSearch.cpp b/xmsgrid/geometry/GmPtSearch.cpp index 134e5148..b76d3f1b 100644 --- a/xmsgrid/geometry/GmPtSearch.cpp +++ b/xmsgrid/geometry/GmPtSearch.cpp @@ -1109,8 +1109,10 @@ void PtSearchUnitTests::testPtsInBox() base = {0, 1, 4}; TS_ASSERT_EQUALS_VEC(base, n); + // The first point must be less than or equal to the second one in both + // coordinates or the search won't find anything. p.PtsInBoxInRtree(Pt3d(1.0, 0.0), Pt3d(0.0, -1.0), n); - base = {0, 1, 4}; + base = {}; TS_ASSERT_EQUALS_VEC(base, n); } // PtSearchUnitTests::testPtsInBox