From 75b1710fd8980d12501b2bff937a687c0024954c Mon Sep 17 00:00:00 2001 From: Patrick Brosi Date: Sun, 16 Aug 2026 14:00:34 +0200 Subject: [PATCH 1/4] add a failing tests between short linestring and short polygon which incorrectly returns the euclidean distance in web mercator instead of the meter distsance, see #12 --- CMakeLists.txt | 2 ++ tests/GeoTestDist.cpp | 24 +++++++++++++++++++----- 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0f0beb0..a6907ed 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -90,6 +90,8 @@ if (ZLIB_FOUND) target_link_libraries(pb_util_xml ${ZLIB_LIBRARIES}) endif(ZLIB_FOUND) +target_link_libraries(pb_util -pthread) + if(CMAKE_TESTING_ENABLED) add_subdirectory(tests) endif() diff --git a/tests/GeoTestDist.cpp b/tests/GeoTestDist.cpp index 4333da1..f995812 100644 --- a/tests/GeoTestDist.cpp +++ b/tests/GeoTestDist.cpp @@ -16,11 +16,16 @@ using namespace util; using namespace util::geo; // _____________________________________________________________________________ +std::string readTestDataset(const std::string& name) { + std::ifstream f(std::string(TEST_DATASETS) + "/" + name, std::ios::binary); + return std::string((std::istreambuf_iterator(f)), {}); +} struct LargeTestGeoms { // unsorted variants MultiPolygon germany, spain; Polygon saimaa; + Polygon vaubaun; Collection flixbus; // xsorted variants @@ -35,15 +40,11 @@ struct LargeTestGeoms { XSortedMultiPolygon germanyMX, spainMX, saimaaMX; XSortedCollection flixbusMX; - std::string readTestDataset(const std::string& name) { - std::ifstream f(std::string(TEST_DATASETS) + "/" + name, std::ios::binary); - return std::string((std::istreambuf_iterator(f)), {}); - } - LargeTestGeoms() : germany(multiPolygonFromWKT(readTestDataset("germany.tsv"))), spain(multiPolygonFromWKT(readTestDataset("spain.tsv"))), saimaa(polygonFromWKT(readTestDataset("saimaa.tsv"))), + vaubaun(polygonFromWKT(readTestDataset("vauban.tsv"))), flixbus(collectionFromWKT(readTestDataset("flixbus.tsv"))), germanyX(germany), spainX(spain), @@ -433,6 +434,16 @@ static void testDistComplexGeoms(const LargeTestGeoms& g) { TEST(util::geo::dist(g.spain, g.flixbus), ==, approx(7.00409)); TEST(util::geo::webMercMeterDist(g.spainM, g.flixbusM), ==, approx(703461.25144)); + + auto vauban = polygonFromWKT(readTestDataset("vauban.tsv")); + auto line = lineFromWKT("LINESTRING(7.8824970 48.0228303,7.8823288 48.0227874,7.8820604 48.0227417,7.8819946 48.0227305)"); + + TEST(util::geo::withinDist(vauban, line, 10), ==, approx(0.06998)); + TEST(util::geo::withinDist(vauban, line, 0.06998), ==, + approx(0.06998)); + TEST(util::geo::dist(vauban, line), ==, approx(0.06998)); + TEST(util::geo::webMercMeterDist(vauban, line), !=, + approx(0.06998)); } // _____________________________________________________________________________ @@ -683,6 +694,9 @@ static void testDistOther() { auto point = pointFromWKT("POINT(4.5 4.5)"); auto point2 = pointFromWKT("POINT(11 11)"); + auto lineFreiburgHbf = lineFromWKT(""); + auto polygonFreiburg = polygonFromWKT(""); + // web mercator copies, for the meter distance assertions auto polyWithInnerM = polygonFromWKTProj( "POLYGON((0 0, 10 0, 10 10, 0 10, 0 0), (4 4, 5 4, 5 5, 4 5, 4 4))", From a88f34103234aae0d400117caad3efaf67e97ae7 Mon Sep 17 00:00:00 2001 From: Patrick Brosi Date: Sun, 16 Aug 2026 14:37:46 +0200 Subject: [PATCH 2/4] use distFunc in outer polygon ring distance computation for small line/polygon pairs, fixes #12 --- geo/Geo.tpp | 2 +- tests/GeoTestDist.cpp | 47 +++++++++++++++++++++++++++++---------- tests/datasets/vauban.tsv | 1 + 3 files changed, 37 insertions(+), 13 deletions(-) create mode 100644 tests/datasets/vauban.tsv diff --git a/geo/Geo.tpp b/geo/Geo.tpp index 023b106..3afc92b 100644 --- a/geo/Geo.tpp +++ b/geo/Geo.tpp @@ -4223,7 +4223,7 @@ double withinDist(const Polygon& poly, const Line& l, double maxDist, } if (intersects(l, poly)) return 0; - double d = dist(l, poly.getOuter()); + double d = dist(l, poly.getOuter(), distFunc); for (const auto& inner : poly.getInners()) { d = std::min(d, dist(l, inner, distFunc)); diff --git a/tests/GeoTestDist.cpp b/tests/GeoTestDist.cpp index f995812..f128905 100644 --- a/tests/GeoTestDist.cpp +++ b/tests/GeoTestDist.cpp @@ -25,30 +25,32 @@ struct LargeTestGeoms { // unsorted variants MultiPolygon germany, spain; Polygon saimaa; - Polygon vaubaun; + Polygon vauban; Collection flixbus; // xsorted variants - XSortedMultiPolygon germanyX, spainX, saimaaX; + XSortedMultiPolygon germanyX, spainX, saimaaX, vaubanX; XSortedCollection flixbusX; // web mercator variants for the meter distance tests MultiPolygon germanyM, spainM; Polygon saimaaM; + Polygon vaubanM; Collection flixbusM; - XSortedMultiPolygon germanyMX, spainMX, saimaaMX; + XSortedMultiPolygon germanyMX, spainMX, saimaaMX, vaubanMX; XSortedCollection flixbusMX; LargeTestGeoms() : germany(multiPolygonFromWKT(readTestDataset("germany.tsv"))), spain(multiPolygonFromWKT(readTestDataset("spain.tsv"))), saimaa(polygonFromWKT(readTestDataset("saimaa.tsv"))), - vaubaun(polygonFromWKT(readTestDataset("vauban.tsv"))), + vauban(polygonFromWKT(readTestDataset("vauban.tsv"))), flixbus(collectionFromWKT(readTestDataset("flixbus.tsv"))), germanyX(germany), spainX(spain), saimaaX(saimaa), + vaubanX(vauban), flixbusX(flixbus), germanyM(multiPolygonFromWKTProj(readTestDataset("germany.tsv"), util::geo::projectToWebMerc)), @@ -56,11 +58,14 @@ struct LargeTestGeoms { util::geo::projectToWebMerc)), saimaaM(polygonFromWKTProj(readTestDataset("saimaa.tsv"), util::geo::projectToWebMerc)), + vaubanM(polygonFromWKTProj(readTestDataset("vauban.tsv"), + util::geo::projectToWebMerc)), flixbusM(collectionFromWKTProj(readTestDataset("flixbus.tsv"), util::geo::projectToWebMerc)), germanyMX(germanyM), spainMX(spainM), saimaaMX(saimaaM), + vaubanMX(vaubanM), flixbusMX(flixbusM) {} }; @@ -435,15 +440,33 @@ static void testDistComplexGeoms(const LargeTestGeoms& g) { TEST(util::geo::webMercMeterDist(g.spainM, g.flixbusM), ==, approx(703461.25144)); - auto vauban = polygonFromWKT(readTestDataset("vauban.tsv")); - auto line = lineFromWKT("LINESTRING(7.8824970 48.0228303,7.8823288 48.0227874,7.8820604 48.0227417,7.8819946 48.0227305)"); + auto line = lineFromWKTProj("LINESTRING(7.8824970 48.0228303,7.8823288 48.0227874,7.8820604 48.0227417,7.8819946 48.0227305)", util::geo::projectToWebMerc); + auto lineX = XSortedLine(line); + + TEST(util::geo::withinDist(g.vaubanM, line, 10), ==, approx(9638.74057)); + TEST(util::geo::withinDist(g.vaubanM, line, 9638.74057), ==, + approx(9638.74057)); + TEST(util::geo::dist(g.vaubanM, line), ==, approx(9638.74057)); + TEST(util::geo::webMercMeterDist(g.vaubanM, line), !=, + approx(util::geo::dist(g.vaubanM, line))); + + TEST(util::geo::webMercMeterDist(g.vaubanM, line), ==, + util::geo::webMercMeterDist(line, g.vaubanM)); + + TEST(util::geo::webMercMeterDist(g.vaubanM, line), ==, + util::geo::webMercMeterDist(lineX, g.vaubanMX)); + + TEST(util::geo::webMercMeterDist(line, g.vaubanM), ==, + util::geo::webMercMeterDist(lineX, g.vaubanMX)); + + TEST(util::geo::webMercMeterDist(line, g.vaubanM), ==, + util::geo::webMercMeterDist(g.vaubanMX, lineX)); + + TEST(util::geo::webMercMeterDist(g.vaubanM, line), ==, + approx(6449.59555)); - TEST(util::geo::withinDist(vauban, line, 10), ==, approx(0.06998)); - TEST(util::geo::withinDist(vauban, line, 0.06998), ==, - approx(0.06998)); - TEST(util::geo::dist(vauban, line), ==, approx(0.06998)); - TEST(util::geo::webMercMeterDist(vauban, line), !=, - approx(0.06998)); + TEST(util::geo::webMercMeterDist(line, g.vaubanM), ==, + approx(6449.59555)); } // _____________________________________________________________________________ diff --git a/tests/datasets/vauban.tsv b/tests/datasets/vauban.tsv new file mode 100644 index 0000000..4f47b0a --- /dev/null +++ b/tests/datasets/vauban.tsv @@ -0,0 +1 @@ +POLYGON((7.8162398 47.9763992,7.8164091 47.9771168,7.8184260 47.9773586,7.8188708 47.9774137,7.8203575 47.9776000,7.8215326 47.9778054,7.8218164 47.9778508,7.8222874 47.9779510,7.8232282 47.9781804,7.8236506 47.9783205,7.8238930 47.9781692,7.8243199 47.9779988,7.8249092 47.9778342,7.8260388 47.9776750,7.8263735 47.9776291,7.8264594 47.9776118,7.8265967 47.9775774,7.8272282 47.9773588,7.8280287 47.9770847,7.8282490 47.9770455,7.8289259 47.9769504,7.8288224 47.9767615,7.8287569 47.9765447,7.8287244 47.9763899,7.8287254 47.9761625,7.8287071 47.9756249,7.8287547 47.9756224,7.8288581 47.9756171,7.8288765 47.9756179,7.8294394 47.9755304,7.8297272 47.9754895,7.8299600 47.9755356,7.8301187 47.9755427,7.8302953 47.9755230,7.8304103 47.9754422,7.8304890 47.9752820,7.8305973 47.9749780,7.8307303 47.9748233,7.8309394 47.9746820,7.8309789 47.9746628,7.8310331 47.9746403,7.8309907 47.9745980,7.8308425 47.9744499,7.8307335 47.9743494,7.8305891 47.9741868,7.8304897 47.9740748,7.8304368 47.9739872,7.8304150 47.9739281,7.8301687 47.9732326,7.8300145 47.9727803,7.8299752 47.9727856,7.8286519 47.9729133,7.8284741 47.9729338,7.8284577 47.9729356,7.8283686 47.9729608,7.8282577 47.9729886,7.8279156 47.9731000,7.8275250 47.9732172,7.8273224 47.9729975,7.8270519 47.9726753,7.8269310 47.9725313,7.8267142 47.9722826,7.8265342 47.9721603,7.8263582 47.9724352,7.8254674 47.9727822,7.8253557 47.9727983,7.8247968 47.9730695,7.8246720 47.9731070,7.8243281 47.9732846,7.8240726 47.9733767,7.8237091 47.9734688,7.8234438 47.9735214,7.8230999 47.9736004,7.8227659 47.9736661,7.8223827 47.9737714,7.8220192 47.9739293,7.8216753 47.9740213,7.8212528 47.9741200,7.8208755 47.9742592,7.8206099 47.9743112,7.8202504 47.9743382,7.8201648 47.9743592,7.8199448 47.9744699,7.8197331 47.9745807,7.8196349 47.9746596,7.8194384 47.9747649,7.8192222 47.9748833,7.8189569 47.9749951,7.8186818 47.9751332,7.8186425 47.9751924,7.8184460 47.9753832,7.8183478 47.9755279,7.8182293 47.9756539,7.8180716 47.9757415,7.8177380 47.9758758,7.8174084 47.9759909,7.8167310 47.9762228,7.8163872 47.9763599,7.8162398 47.9763992)) From cad665d033ffacd353af913cd5a664e5eef4d0cc Mon Sep 17 00:00:00 2001 From: Patrick Brosi Date: Sun, 16 Aug 2026 16:23:37 +0200 Subject: [PATCH 3/4] internal conversion to double in computation of projection progression, bake fail from #13 into test --- geo/Geo.tpp | 82 ++++++++++++++++++++++--------------------- tests/GeoTestDist.cpp | 67 +++++++++++++++++++++++++++++++++++ 2 files changed, 109 insertions(+), 40 deletions(-) diff --git a/geo/Geo.tpp b/geo/Geo.tpp index 3afc92b..1fff158 100644 --- a/geo/Geo.tpp +++ b/geo/Geo.tpp @@ -1235,10 +1235,10 @@ std::pair withinDist(const Point& p, const XSortedRing& ph, } if (euclideanDist < euclideanDistUpperBound) { euclideanDistUpperBound = euclideanDist; - padding = paddingFunc(euclideanDistUpperBound, minDist, - getBoundingBox(p), ph.boundingBox()); - xPadding = splitPadding(padding, getBoundingBox(p), ph.boundingBox()) - .xPadding; + padding = paddingFunc(euclideanDistUpperBound, minDist, getBoundingBox(p), + ph.boundingBox()); + xPadding = + splitPadding(padding, getBoundingBox(p), ph.boundingBox()).xPadding; } } @@ -1356,8 +1356,9 @@ double withinDist(const Point& p, const XSortedLine& line, double maxDist, distFunc(p, line.rawLine().front().p, std::numeric_limits::max())); - auto padding = paddingFunc(euclideanDistUpperBound, std::min(maxDist, minDist), - getBoundingBox(p), line.boundingBox()); + auto padding = + paddingFunc(euclideanDistUpperBound, std::min(maxDist, minDist), + getBoundingBox(p), line.boundingBox()); auto xPadding = splitPadding(padding, getBoundingBox(p), line.boundingBox()).xPadding; @@ -1395,8 +1396,8 @@ double withinDist(const Point& p, const XSortedLine& line, double maxDist, euclideanDistUpperBound = euclideanDist; padding = paddingFunc(euclideanDistUpperBound, std::min(maxDist, minDist), getBoundingBox(p), line.boundingBox()); - xPadding = splitPadding(padding, getBoundingBox(p), line.boundingBox()) - .xPadding; + xPadding = + splitPadding(padding, getBoundingBox(p), line.boundingBox()).xPadding; } } @@ -3619,15 +3620,17 @@ double withinDist(const LineSegment& ls1, const LineSegment& ls2, double d = distToSegment(ls2.first.getX(), ls2.first.getY(), ls2.second.getX(), ls2.second.getY(), ls1.first.getX(), ls1.first.getY(), distFunc); - d = std::min(d, distToSegment(ls2.first.getX(), ls2.first.getY(), - ls2.second.getX(), ls2.second.getY(), - ls1.second.getX(), ls1.second.getY(), distFunc)); + d = std::min( + d, distToSegment(ls2.first.getX(), ls2.first.getY(), ls2.second.getX(), + ls2.second.getY(), ls1.second.getX(), ls1.second.getY(), + distFunc)); d = std::min(d, distToSegment(ls1.first.getX(), ls1.first.getY(), ls1.second.getX(), ls1.second.getY(), ls2.first.getX(), ls2.first.getY(), distFunc)); - d = std::min(d, distToSegment(ls1.first.getX(), ls1.first.getY(), - ls1.second.getX(), ls1.second.getY(), - ls2.second.getX(), ls2.second.getY(), distFunc)); + d = std::min( + d, distToSegment(ls1.first.getX(), ls1.first.getY(), ls1.second.getX(), + ls1.second.getY(), ls2.second.getX(), ls2.second.getY(), + distFunc)); return d; } @@ -3663,15 +3666,17 @@ double dist(const LineSegment& ls1, const LineSegment& ls2, double d = distToSegment(ls2.first.getX(), ls2.first.getY(), ls2.second.getX(), ls2.second.getY(), ls1.first.getX(), ls1.first.getY(), distFunc); - d = std::min(d, distToSegment(ls2.first.getX(), ls2.first.getY(), - ls2.second.getX(), ls2.second.getY(), - ls1.second.getX(), ls1.second.getY(), distFunc)); + d = std::min( + d, distToSegment(ls2.first.getX(), ls2.first.getY(), ls2.second.getX(), + ls2.second.getY(), ls1.second.getX(), ls1.second.getY(), + distFunc)); d = std::min(d, distToSegment(ls1.first.getX(), ls1.first.getY(), ls1.second.getX(), ls1.second.getY(), ls2.first.getX(), ls2.first.getY(), distFunc)); - d = std::min(d, distToSegment(ls1.first.getX(), ls1.first.getY(), - ls1.second.getX(), ls1.second.getY(), - ls2.second.getX(), ls2.second.getY(), distFunc)); + d = std::min( + d, distToSegment(ls1.first.getX(), ls1.first.getY(), ls1.second.getX(), + ls1.second.getY(), ls2.second.getX(), ls2.second.getY(), + distFunc)); return d; } @@ -5143,7 +5148,8 @@ double distToSegment(T lax, T lay, T lbx, T lby, T px, T py, DF&& distFunc) { return distFunc(Point{px, py}, Point{lax, lay}, std::numeric_limits::max()); - double t = ((px - lax) * (lbx - lax) + (py - lay) * (lby - lay)) / d; + double t = + ((px - lax) * 1.0 * (lbx - lax) + (py - lay) * 1.0 * (lby - lay)) / d; if (t < 0) { return distFunc(Point{px, py}, Point{lax, lay}, @@ -5154,8 +5160,8 @@ double distToSegment(T lax, T lay, T lbx, T lby, T px, T py, DF&& distFunc) { } return distFunc(Point{px, py}, - Point{static_cast(lax + t * (lbx - lax)), - static_cast(lay + t * (lby - lay))}, + Point{static_cast(lax * 1.0 + t * (lbx - lax)), + static_cast(lay * 1.0 + t * (lby - lay))}, std::numeric_limits::max()); } @@ -6454,9 +6460,8 @@ double withinDist(const std::vector>& ls1, size_t k = 0; // position in OUT ls2 size_t ls2OutSize = 0; - double padding = - paddingFunc(euclideanDistUpperBound, std::min(minDist, maxDist), boxA, - boxB); + double padding = paddingFunc(euclideanDistUpperBound, + std::min(minDist, maxDist), boxA, boxB); const auto pad = splitPadding(padding, boxA, boxB); T xPadding = std::min(std::numeric_limits::max() * 1.0, pad.xPadding); @@ -6557,8 +6562,8 @@ double withinDist(const std::vector>& ls1, ls1seg); if (processActives(activesB, ls1seg, euclideanDistUpperBound, minDist, - maxDist, padding, xPadding, yPadding, - box, boxB, boxA, segs, paddingFunc, distFunc)) + maxDist, padding, xPadding, yPadding, box, boxB, + boxA, segs, paddingFunc, distFunc)) return minDist; } @@ -6649,8 +6654,8 @@ double withinDist(const std::vector>& ls1, ls2OutSeg); if (processActives(activesA, ls2OutSeg, euclideanDistUpperBound, minDist, - maxDist, padding, xPadding, yPadding, - box, boxA, boxB, segs, paddingFunc, distFunc)) + maxDist, padding, xPadding, yPadding, box, boxA, boxB, + segs, paddingFunc, distFunc)) return minDist; // advance to next OUT @@ -6994,8 +6999,7 @@ std::tuple probeDistanceUpperBound( } } - return {upperBound, euclideanUpperBound, - stepA == 1 && stepB == 1 && !pruned}; + return {upperBound, euclideanUpperBound, stepA == 1 && stepB == 1 && !pruned}; } // _____________________________________________________________________________ @@ -7085,10 +7089,10 @@ Padding splitPadding(double padding, const Box& boxA, const Box& boxB) { LineSegment{Point{0, boxB.getLowerLeft().getY()}, Point{0, boxB.getUpperRight().getY()}}); - return {sqrt(std::max(0.0, padding * padding - - minEuclideanYDist * minEuclideanYDist)), - sqrt(std::max(0.0, padding * padding - - minEuclideanXDist * minEuclideanXDist))}; + return {sqrt(std::max( + 0.0, padding * padding - minEuclideanYDist * minEuclideanYDist)), + sqrt(std::max( + 0.0, padding * padding - minEuclideanXDist * minEuclideanXDist))}; } // _____________________________________________________________________________ @@ -7282,8 +7286,7 @@ double webMercMeterDistLocalSearchPadding(double euclideanDistanceUpperBound, auto boxBStar = util::geo::intersection(paddedA, boxB); // may be empty! - if (boxBStar.isNull()) - return factorNew2 * euclideanDistanceUpperBound; + if (boxBStar.isNull()) return factorNew2 * euclideanDistanceUpperBound; double min2 = std::numeric_limits::infinity(); double max2 = 0; @@ -7305,8 +7308,7 @@ double webMercMeterDistLocalSearchPadding(double euclideanDistanceUpperBound, double factorNew3 = max2 / min2; - if (factorNew2 < factorNew3) - return factorNew2 * euclideanDistanceUpperBound; + if (factorNew2 < factorNew3) return factorNew2 * euclideanDistanceUpperBound; return factorNew3 * euclideanDistanceUpperBound; } diff --git a/tests/GeoTestDist.cpp b/tests/GeoTestDist.cpp index f128905..96ffd07 100644 --- a/tests/GeoTestDist.cpp +++ b/tests/GeoTestDist.cpp @@ -982,6 +982,73 @@ static void testDistLimitedPrecision() { auto point_b = pointFromWKT("POINT(1 1)"); auto point2_b = pointFromWKT("POINT(0 0)"); TEST(util::geo::dist(point_b, point2_b), ==, approx(sqrt(2))); + + auto germanyCoarse = polygonFromWKTProj( + "POLYGON((7.20369317867016 53.62121249029073, 9.335040870259194 " + "54.77156944262062, 13.97127141588071 53.7058383745324, " + "14.77327338230339 51.01654754091759, 11.916828022441791 " + "50.36932046223437, 13.674640551587391 48.68663848319227, " + "12.773761630400273 47.74969625921073, 7.58917 47.59002, 8.03916 " + "49.01783, 6.50056816701192 49.535220384133375, 6.0391423781112 " + "51.804566644690524, 7.20369317867016 53.62121249029073))", + [](const DPoint& p, CRSType) { + auto proj = util::geo::projectToWebMerc(p, CRS84); + return Point{static_cast(proj.getX() * 10), + static_cast(proj.getY() * 10)}; + }); + auto londonCoarse = polygonFromWKTProj( + "POLYGON((-0.1198608 51.5027451,-0.1197395 51.5027354,-0.1194922 " + "51.5039381,-0.1196135 51.5039478,-0.1198608 51.5027451))", + [](const DPoint& p, CRSType) { + auto proj = util::geo::projectToWebMerc(p, CRS84); + return Point{static_cast(proj.getX() * 10), + static_cast(proj.getY() * 10)}; + }); + + auto germanyCoarseRaw = polygonFromWKTProj( + "POLYGON((7.20369317867016 53.62121249029073, 9.335040870259194 " + "54.77156944262062, 13.97127141588071 53.7058383745324, " + "14.77327338230339 51.01654754091759, 11.916828022441791 " + "50.36932046223437, 13.674640551587391 48.68663848319227, " + "12.773761630400273 47.74969625921073, 7.58917 47.59002, 8.03916 " + "49.01783, 6.50056816701192 49.535220384133375, 6.0391423781112 " + "51.804566644690524, 7.20369317867016 53.62121249029073))", + util::geo::projectToWebMerc); + auto londonCoarseRaw = polygonFromWKTProj( + "POLYGON((-0.1198608 51.5027451,-0.1197395 51.5027354,-0.1194922 " + "51.5039381,-0.1196135 51.5039478,-0.1198608 51.5027451))", + util::geo::projectToWebMerc); + + TEST(util::geo::webMercMeterDist(germanyCoarseRaw, londonCoarseRaw), ==, + approx(426521.22769)); + + TEST( + util::geo::withinDist( + germanyCoarse, londonCoarse, 426521.0 + 10, + [](double euDistUp, double distUp, Box boxa, + Box boxb) -> double { + euDistUp = euDistUp / 10.0; + DBox boxAD{{(boxa.getLowerLeft().getX() * 1.0) / 10.0, + (boxa.getLowerLeft().getY() * 1.0) / 10.0}, + {(boxa.getUpperRight().getX() * 1.0) / 10.0, + (boxa.getUpperRight().getY() * 1.0) / 10.0}}; + DBox boxBD{{(boxb.getLowerLeft().getX() * 1.0) / 10.0, + (boxb.getLowerLeft().getY() * 1.0) / 10.0}, + {(boxb.getUpperRight().getX() * 1.0) / 10.0, + (boxb.getUpperRight().getY() * 1.0) / 10.0}}; + return webMercMeterDistLocalSearchPadding(euDistUp, distUp, boxAD, + boxBD) * + 10.0; + }, + 426521 * 1.05, + [](const Point a, const Point b, double) -> double { + DPoint aReal{(a.getX() * 1.0) / 10.0, (a.getY() * 1.0) / 10.0}; + DPoint bReal{(b.getX() * 1.0) / 10.0, (b.getY() * 1.0) / 10.0}; + return haversineWebMerc(aReal, bReal); + }), + // NOTE: difference because of precision to only 10 cm because of coarse + // projection + ==, approx(426521.18896)); } // _____________________________________________________________________________ From e5141609290b8b1d026d5961d3ff81510d875797 Mon Sep 17 00:00:00 2001 From: Patrick Brosi Date: Sun, 16 Aug 2026 18:01:53 +0200 Subject: [PATCH 4/4] add 1FFF0FFF2 constant which is the CORRECT de9im matrix between two equivalent linestrings --- geo/DE9IMatrix.h | 1 + 1 file changed, 1 insertion(+) diff --git a/geo/DE9IMatrix.h b/geo/DE9IMatrix.h index 309c594..03b20d7 100644 --- a/geo/DE9IMatrix.h +++ b/geo/DE9IMatrix.h @@ -352,6 +352,7 @@ static CONSTEXPR DE9IMatrix M2FFF1FFF2("2FFF1FFF2"); static CONSTEXPR DE9IMatrix M2FF1FF212("2FF1FF212"); static CONSTEXPR DE9IMatrix M1FF0FF212("1FF0FF212"); static CONSTEXPR DE9IMatrix M10FF0FFF2("10FF0FFF2"); +static CONSTEXPR DE9IMatrix M1FFF0FFF2("1FFF0FFF2"); static CONSTEXPR DE9IMatrix MFF1FF0212("FF1FF0212"); static CONSTEXPR DE9IMatrix M2F2FFF2F2("2F2FFF2F2");