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/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"); diff --git a/geo/Geo.tpp b/geo/Geo.tpp index 023b106..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; } @@ -4223,7 +4228,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)); @@ -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 4333da1..96ffd07 100644 --- a/tests/GeoTestDist.cpp +++ b/tests/GeoTestDist.cpp @@ -16,38 +16,41 @@ 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 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; - 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"))), + 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)), @@ -55,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) {} }; @@ -433,6 +439,34 @@ 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 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::webMercMeterDist(line, g.vaubanM), ==, + approx(6449.59555)); } // _____________________________________________________________________________ @@ -683,6 +717,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))", @@ -945,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)); } // _____________________________________________________________________________ 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))