Skip to content

Fix various meterdist bugs #14

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
patrickbr merged 4 commits into from
Aug 16, 2026
Merged

Fix various meterdist bugs #14

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
75b1710
add a failing tests between short linestring and short polygon which …
patrickbr Aug 16, 2026
a88f341
use distFunc in outer polygon ring distance computation for small lin…
patrickbr Aug 16, 2026
cad665d
internal conversion to double in computation of projection progressio…
patrickbr Aug 16, 2026
e514160
add 1FFF0FFF2 constant which is the CORRECT de9im matrix between two …
patrickbr Aug 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()
1 change: 1 addition & 0 deletions geo/DE9IMatrix.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand Down
84 changes: 43 additions & 41 deletions geo/Geo.tpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1235,10 +1235,10 @@ std::pair<double, bool> withinDist(const Point<T>& p, const XSortedRing<T>& 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;
}
}

Expand Down Expand Up @@ -1356,8 +1356,9 @@ double withinDist(const Point<T>& p, const XSortedLine<T>& line, double maxDist,
distFunc(p, line.rawLine().front().p,
std::numeric_limits<double>::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;

Expand Down Expand Up @@ -1395,8 +1396,8 @@ double withinDist(const Point<T>& p, const XSortedLine<T>& 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;
}
}

Expand Down Expand Up @@ -3619,15 +3620,17 @@ double withinDist(const LineSegment<T>& ls1, const LineSegment<T>& 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;
}

Expand Down Expand Up @@ -3663,15 +3666,17 @@ double dist(const LineSegment<T>& ls1, const LineSegment<T>& 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;
}

Expand Down Expand Up @@ -4223,7 +4228,7 @@ double withinDist(const Polygon<T>& poly, const Line<T>& 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));
Expand Down Expand Up @@ -5143,7 +5148,8 @@ double distToSegment(T lax, T lay, T lbx, T lby, T px, T py, DF&& distFunc) {
return distFunc(Point<T>{px, py}, Point<T>{lax, lay},
std::numeric_limits<double>::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<T>{px, py}, Point<T>{lax, lay},
Expand All @@ -5154,8 +5160,8 @@ double distToSegment(T lax, T lay, T lbx, T lby, T px, T py, DF&& distFunc) {
}

return distFunc(Point<T>{px, py},
Point<T>{static_cast<T>(lax + t * (lbx - lax)),
static_cast<T>(lay + t * (lby - lay))},
Point<T>{static_cast<T>(lax * 1.0 + t * (lbx - lax)),
static_cast<T>(lay * 1.0 + t * (lby - lay))},
std::numeric_limits<double>::max());
}

Expand Down Expand Up @@ -6454,9 +6460,8 @@ double withinDist(const std::vector<XSortedTuple<T>>& 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<T>::max() * 1.0, pad.xPadding);
Expand Down Expand Up @@ -6557,8 +6562,8 @@ double withinDist(const std::vector<XSortedTuple<T>>& 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;
}

Expand Down Expand Up @@ -6649,8 +6654,8 @@ double withinDist(const std::vector<XSortedTuple<T>>& 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
Expand Down Expand Up @@ -6994,8 +6999,7 @@ std::tuple<double, double, bool> probeDistanceUpperBound(
}
}

return {upperBound, euclideanUpperBound,
stepA == 1 && stepB == 1 && !pruned};
return {upperBound, euclideanUpperBound, stepA == 1 && stepB == 1 && !pruned};
}

// _____________________________________________________________________________
Expand Down Expand Up @@ -7085,10 +7089,10 @@ Padding splitPadding(double padding, const Box<T>& boxA, const Box<T>& boxB) {
LineSegment<T>{Point<T>{0, boxB.getLowerLeft().getY()},
Point<T>{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))};
}

// _____________________________________________________________________________
Expand Down Expand Up @@ -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<double>::infinity();
double max2 = 0;
Expand All @@ -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;
}
Expand Down
118 changes: 111 additions & 7 deletions tests/GeoTestDist.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,50 +16,56 @@ 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<char>(f)), {});
}

struct LargeTestGeoms {
// unsorted variants
MultiPolygon<double> germany, spain;
Polygon<double> saimaa;
Polygon<double> vauban;
Collection<double> flixbus;

// xsorted variants
XSortedMultiPolygon<double> germanyX, spainX, saimaaX;
XSortedMultiPolygon<double> germanyX, spainX, saimaaX, vaubanX;
XSortedCollection<double> flixbusX;

// web mercator variants for the meter distance tests
MultiPolygon<double> germanyM, spainM;
Polygon<double> saimaaM;
Polygon<double> vaubanM;
Collection<double> flixbusM;

XSortedMultiPolygon<double> germanyMX, spainMX, saimaaMX;
XSortedMultiPolygon<double> germanyMX, spainMX, saimaaMX, vaubanMX;
XSortedCollection<double> 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<char>(f)), {});
}

LargeTestGeoms()
: germany(multiPolygonFromWKT<double>(readTestDataset("germany.tsv"))),
spain(multiPolygonFromWKT<double>(readTestDataset("spain.tsv"))),
saimaa(polygonFromWKT<double>(readTestDataset("saimaa.tsv"))),
vauban(polygonFromWKT<double>(readTestDataset("vauban.tsv"))),
flixbus(collectionFromWKT<double>(readTestDataset("flixbus.tsv"))),
germanyX(germany),
spainX(spain),
saimaaX(saimaa),
vaubanX(vauban),
flixbusX(flixbus),
germanyM(multiPolygonFromWKTProj<double>(readTestDataset("germany.tsv"),
util::geo::projectToWebMerc<double>)),
spainM(multiPolygonFromWKTProj<double>(readTestDataset("spain.tsv"),
util::geo::projectToWebMerc<double>)),
saimaaM(polygonFromWKTProj<double>(readTestDataset("saimaa.tsv"),
util::geo::projectToWebMerc<double>)),
vaubanM(polygonFromWKTProj<double>(readTestDataset("vauban.tsv"),
util::geo::projectToWebMerc<double>)),
flixbusM(collectionFromWKTProj<double>(readTestDataset("flixbus.tsv"),
util::geo::projectToWebMerc<double>)),
germanyMX(germanyM),
spainMX(spainM),
saimaaMX(saimaaM),
vaubanMX(vaubanM),
flixbusMX(flixbusM) {}
};

Expand Down Expand Up @@ -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<double>("LINESTRING(7.8824970 48.0228303,7.8823288 48.0227874,7.8820604 48.0227417,7.8819946 48.0227305)", util::geo::projectToWebMerc<double>);
auto lineX = XSortedLine<double>(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));
}

// _____________________________________________________________________________
Expand Down Expand Up @@ -683,6 +717,9 @@ static void testDistOther() {
auto point = pointFromWKT<double>("POINT(4.5 4.5)");
auto point2 = pointFromWKT<double>("POINT(11 11)");

auto lineFreiburgHbf = lineFromWKT<double>("");
auto polygonFreiburg = polygonFromWKT<double>("");

// web mercator copies, for the meter distance assertions
auto polyWithInnerM = polygonFromWKTProj<double>(
"POLYGON((0 0, 10 0, 10 10, 0 10, 0 0), (4 4, 5 4, 5 5, 4 5, 4 4))",
Expand Down Expand Up @@ -945,6 +982,73 @@ static void testDistLimitedPrecision() {
auto point_b = pointFromWKT<bool>("POINT(1 1)");
auto point2_b = pointFromWKT<bool>("POINT(0 0)");
TEST(util::geo::dist(point_b, point2_b), ==, approx(sqrt(2)));

auto germanyCoarse = polygonFromWKTProj<int32_t>(
"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<double>(p, CRS84);
return Point<int32_t>{static_cast<int32_t>(proj.getX() * 10),
static_cast<int32_t>(proj.getY() * 10)};
});
auto londonCoarse = polygonFromWKTProj<int32_t>(
"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<double>(p, CRS84);
return Point<int32_t>{static_cast<int32_t>(proj.getX() * 10),
static_cast<int32_t>(proj.getY() * 10)};
});

auto germanyCoarseRaw = polygonFromWKTProj<double>(
"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<double>);
auto londonCoarseRaw = polygonFromWKTProj<double>(
"POLYGON((-0.1198608 51.5027451,-0.1197395 51.5027354,-0.1194922 "
"51.5039381,-0.1196135 51.5039478,-0.1198608 51.5027451))",
util::geo::projectToWebMerc<double>);

TEST(util::geo::webMercMeterDist(germanyCoarseRaw, londonCoarseRaw), ==,
approx(426521.22769));

TEST(
util::geo::withinDist(
germanyCoarse, londonCoarse, 426521.0 + 10,
[](double euDistUp, double distUp, Box<int32_t> boxa,
Box<int32_t> 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<int32_t> a, const Point<int32_t> 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));
}

// _____________________________________________________________________________
Expand Down
Loading
Loading