From 853e633f469a3cd21be177cc7baf3a39623be514 Mon Sep 17 00:00:00 2001 From: Patrick Brosi Date: Fri, 21 Feb 2025 12:36:43 +0100 Subject: [PATCH] first rough implementation of --within-distance, already works for points --- src/spatialjoin/SpatialJoinMain.cpp | 10 ++++ src/spatialjoin/Sweeper.cpp | 86 ++++++++++++++++++++++++----- src/spatialjoin/Sweeper.h | 6 ++ src/spatialjoin/tests/TestMain.cpp | 18 +++--- 4 files changed, 97 insertions(+), 23 deletions(-) diff --git a/src/spatialjoin/SpatialJoinMain.cpp b/src/spatialjoin/SpatialJoinMain.cpp index 8d3cc03..78d7aff 100755 --- a/src/spatialjoin/SpatialJoinMain.cpp +++ b/src/spatialjoin/SpatialJoinMain.cpp @@ -70,6 +70,8 @@ void printHelp(int argc, char** argv) { << "separator between crossing geometry IDs\n" << std::setw(41) << " --suffix (default: '\\n')" << "suffix added at the beginning of every relation\n\n" + << std::setw(41) << " --within-distance (default: '')" + << "if set to non-negative value, only compute for each object the objects within the given distance\n\n" << std::setfill(' ') << std::left << "Geometric computation:\n" << std::setw(41) << " --no-box-ids" << "disable box id criteria for contains/covers/intersect computation\n" @@ -121,6 +123,7 @@ int main(int argc, char** argv) { std::string overlaps = " overlaps "; std::string crosses = " crosses "; std::string suffix = "\n"; + double withinDist = -1; bool useBoxIds = true; bool useArea = true; @@ -171,6 +174,8 @@ int main(int argc, char** argv) { state = 13; } else if (cur == "--cache-max-size") { state = 14; + } else if (cur == "--within-distance") { + state = 15; } else if (cur == "--no-box-ids") { useBoxIds = false; } else if (cur == "--no-surface-area") { @@ -248,6 +253,10 @@ int main(int argc, char** argv) { geomCacheMaxSize = atoi(cur.c_str()); state = 0; break; + case 15: + withinDist = atof(cur.c_str()); + state = 0; + break; } } @@ -277,6 +286,7 @@ int main(int argc, char** argv) { useFastSweepSkip, useInnerOuter, noGeometryChecks, + withinDist, {}, [](const std::string& s) { LOGTO(INFO, std::cerr) << s; }, [](const std::string& s) { std::cerr << s; }, diff --git a/src/spatialjoin/Sweeper.cpp b/src/spatialjoin/Sweeper.cpp index 9c9a9e9..460dde7 100644 --- a/src/spatialjoin/Sweeper.cpp +++ b/src/spatialjoin/Sweeper.cpp @@ -5,8 +5,8 @@ #include #include -#include #include +#include #include #include #include @@ -31,6 +31,7 @@ using sj::innerouter::Mode; using util::writeAll; using util::geo::area; using util::geo::getBoundingBox; +using util::geo::FPoint; using util::geo::I32Box; using util::geo::I32Line; using util::geo::I32MultiLine; @@ -167,7 +168,7 @@ void Sweeper::add(const I32Polygon& poly, const std::string& gidR, size_t subid, std::string gid = (side ? ("B" + gidR) : ("A" + gidR)); WriteCand cur; - const auto& box = getBoundingBox(poly); + const auto& box = getPaddedBoundingBox(poly); I32XSortedPolygon spoly(poly); if (spoly.empty()) return; @@ -188,7 +189,7 @@ void Sweeper::add(const I32Polygon& poly, const std::string& gidR, size_t subid, I32Box box45; if (_cfg.useDiagBox) { auto polyR = util::geo::rotateSinCos(poly, sin45, cos45, I32Point(0, 0)); - box45 = getBoundingBox(polyR); + box45 = getPaddedBoundingBox(polyR); } cur.subid = subid; @@ -320,7 +321,7 @@ void Sweeper::add(const I32Line& line, const std::string& gidR, size_t subid, WriteCand cur; - const auto& box = getBoundingBox(line); + const auto& box = getPaddedBoundingBox(line); BoxIdList boxIds; std::map cutouts; @@ -337,7 +338,7 @@ void Sweeper::add(const I32Line& line, const std::string& gidR, size_t subid, I32Box box45; if (_cfg.useDiagBox) { auto polyR = util::geo::rotateSinCos(line, sin45, cos45, I32Point(0, 0)); - box45 = getBoundingBox(polyR); + box45 = getPaddedBoundingBox(polyR); } cur.subid = subid; @@ -443,25 +444,27 @@ void Sweeper::add(const I32Point& point, const std::string& gidR, size_t subid, cur.subid = subid; cur.gid = gid; + const auto& box = getPaddedBoundingBox(point); + auto pointR = util::geo::rotateSinCos(point, sin45, cos45, I32Point(0, 0)); cur.boxvalIn = {0, // placeholder, will be overwritten later on - point.getY(), - point.getY(), - point.getX(), + box.getLowerLeft().getY(), + box.getUpperRight().getY(), + box.getLowerLeft().getX(), false, POINT, 0, - util::geo::getBoundingBox(pointR), + getPaddedBoundingBox(pointR), side, false}; cur.boxvalOut = {0, // placeholder, will be overwritten later on - point.getY(), - point.getY(), - point.getX(), + box.getLowerLeft().getY(), + box.getUpperRight().getY(), + box.getUpperRight().getX(), true, POINT, 0, - util::geo::getBoundingBox(pointR), + getPaddedBoundingBox(pointR), side, false}; @@ -867,7 +870,16 @@ void Sweeper::flush() { for (size_t side = 0; side < 2; side++) { for (size_t i = 0; i < _multiIds[side].size(); i++) { - diskAdd({i, 1, 0, _multiLeftX[side][i] - 1, false, POINT, 0.0, {}, side}); + diskAdd({i, + 1, + 0, + _multiLeftX[side][i] - 1, + false, + POINT, + 0.0, + {}, + side, + false}); } } @@ -1773,6 +1785,13 @@ void Sweeper::writeRel(size_t t, const std::string& a, const std::string& b, _stats[t].timeWrite += TOOK(ts); } +// ____________________________________________________________________________ +void Sweeper::writeDist(size_t t, const std::string& a, size_t aSub, + const std::string& b, size_t bSub, double dist) { + writeRel(t, a, b, "\t" + std::to_string(dist) + "\t"); + writeRel(t, b, a, "\t" + std::to_string(dist) + "\t"); +} + // ____________________________________________________________________________ void Sweeper::writeIntersect(size_t t, const std::string& a, const std::string& b) { @@ -1850,6 +1869,26 @@ void Sweeper::doCheck(const BoxVal cur, const SweepVal sv, size_t t) { // every 10000 checks, update our position if (_checks[t] % 10000 == 0) _atomicCurX[t] = _curX[t]; + if (_cfg.withinDist >= 0) { + if (cur.type == POINT && sv.type == POINT) { + auto p1 = util::geo::centroid(cur.b45); + auto p2 = util::geo::centroid(sv.b45); + p1 = util::geo::rotateSinCos(p1, -sin45, -cos45, I32Point(0, 0)); + p2 = util::geo::rotateSinCos(p2, -sin45, -cos45, I32Point(0, 0)); + + auto dist = util::geo::webMercMeterDist(FPoint{(p1.getX() * 1.0) / (PREC * 1.0), (p1.getY() * 1.0) / (PREC * 1.0)}, FPoint{(p2.getX() * 1.0) / (PREC * 1.0), (p2.getY() * 1.0) / (PREC * 1.0)}); + + if (dist <= _cfg.withinDist) { + auto a = _pointCache.get(cur.id, cur.large ? -1 : t); + auto b = _pointCache.get(sv.id, sv.large ? -1 : t); + + writeDist(t, a->id, a->subId, b->id, b->subId, dist); + } + } + + return; + } + if ((cur.type == SIMPLE_POLYGON || cur.type == POLYGON) && (sv.type == SIMPLE_POLYGON || sv.type == POLYGON)) { auto ts = TIME(); @@ -3135,6 +3174,25 @@ std::string Sweeper::intToBase126(uint64_t id) { return ret; } +// _____________________________________________________________________________ +template +I32Box Sweeper::getPaddedBoundingBox(const G& geom) const { + auto bbox = getBoundingBox(geom); + + if (_cfg.withinDist >= 0) { + double xScaleFactor = + std::min(util::geo::webMercDistFactor(I32Point{bbox.getLowerLeft().getX() / PREC, bbox.getLowerLeft().getY() / PREC}), + util::geo::webMercDistFactor(I32Point{bbox.getUpperRight().getX() / PREC, bbox.getUpperRight().getY() / PREC})); + + uint32_t padX = (_cfg.withinDist / xScaleFactor) * PREC; + uint32_t padY = _cfg.withinDist * PREC; + + return {{bbox.getLowerLeft().getX() - padX, bbox.getLowerLeft().getY() - padY}, {bbox.getUpperRight().getX() + padX, bbox.getUpperRight().getY() + padY}}; + } + + return bbox; +} + // _____________________________________________________________________________ uint64_t Sweeper::base126ToInt(const std::string& id) { uint64_t ret = 0; diff --git a/src/spatialjoin/Sweeper.h b/src/spatialjoin/Sweeper.h index ed5401f..fd00a2d 100644 --- a/src/spatialjoin/Sweeper.h +++ b/src/spatialjoin/Sweeper.h @@ -137,6 +137,7 @@ struct SweeperCfg { bool useFastSweepSkip; bool useInnerOuter; bool noGeometryChecks; + double withinDist; std::function writeRelCb; std::function logCb; @@ -363,6 +364,9 @@ class Sweeper { bool refRelated(const std::string& a, const std::string& b) const; + template + util::geo::I32Box getPaddedBoundingBox(const G& geom) const; + void diskAdd(const BoxVal& bv); void multiOut(size_t t, const std::string& gid); @@ -381,6 +385,8 @@ class Sweeper { size_t bSub); void writeEquals(size_t t, const std::string& a, size_t aSub, const std::string& b, size_t bSub); + void writeDist(size_t t, const std::string& a, size_t aSub, + const std::string& b, size_t bSub, double dist); void writeTouches(size_t t, const std::string& a, size_t aSub, const std::string& b, size_t bSub); void writeNotTouches(size_t t, const std::string& a, size_t aSub, diff --git a/src/spatialjoin/tests/TestMain.cpp b/src/spatialjoin/tests/TestMain.cpp index cb5f52e..8407639 100644 --- a/src/spatialjoin/tests/TestMain.cpp +++ b/src/spatialjoin/tests/TestMain.cpp @@ -78,63 +78,63 @@ int main(int, char**) { NUM_THREADS, NUM_THREADS, 1000, "$", " intersects ", " contains ", " covers ", " touches ", " equals ", " overlaps ", " crosses ", "$\n", false, false, false, - false, false, false, false, false, + false, false, false, false, false, -1, {}, {}, {}, {}}; sj::SweeperCfg all{ NUM_THREADS, NUM_THREADS, 1000, "$", " intersects ", " contains ", " covers ", " touches ", " equals ", " overlaps ", " crosses ", "$\n", true, true, true, - true, true, true, true, false, + true, true, true, true, false, -1, {}, {}, {}, {}}; sj::SweeperCfg noSurfaceArea{ NUM_THREADS, NUM_THREADS, 1000, "$", " intersects ", " contains ", " covers ", " touches ", " equals ", " overlaps ", " crosses ", "$\n", true, false, true, - true, true, true, true, false, + true, true, true, true, false, -1, {}, {}, {}, {}}; sj::SweeperCfg noBoxIds{ NUM_THREADS, NUM_THREADS, 1000, "$", " intersects ", " contains ", " covers ", " touches ", " equals ", " overlaps ", " crosses ", "$\n", false, true, true, - true, true, true, true, false, + true, true, true, true, false, -1, {}, {}, {}, {}}; sj::SweeperCfg noObb{ NUM_THREADS, NUM_THREADS, 1000, "$", " intersects ", " contains ", " covers ", " touches ", " equals ", " overlaps ", " crosses ", "$\n", true, true, false, - true, true, true, true, false, + true, true, true, true, false, -1, {}, {}, {}, {}}; sj::SweeperCfg noCutouts{ NUM_THREADS, NUM_THREADS, 1000, "$", " intersects ", " contains ", " covers ", " touches ", " equals ", " overlaps ", " crosses ", "$\n", true, true, true, - false, true, true, true, false, + false, true, true, true, false, -1, {}, {}, {}, {}}; sj::SweeperCfg noDiagBox{ NUM_THREADS, NUM_THREADS, 1000, "$", " intersects ", " contains ", " covers ", " touches ", " equals ", " overlaps ", " crosses ", "$\n", true, true, true, - true, false, true, true, false, + true, false, true, true, false, -1, {}, {}, {}, {}}; sj::SweeperCfg noFastSweep{ NUM_THREADS, NUM_THREADS, 1000, "$", " intersects ", " contains ", " covers ", " touches ", " equals ", " overlaps ", " crosses ", "$\n", true, true, true, - true, true, false, true, false, + true, true, false, true, false, -1, {}, {}, {}, {}}; sj::SweeperCfg noInnerOuter{ NUM_THREADS, NUM_THREADS, 1000, "$", " intersects ", " contains ", " covers ", " touches ", " equals ", " overlaps ", " crosses ", "$\n", true, true, true, - true, true, true, false, false, + true, true, true, false, false, -1, {}, {}, {}, {}}; std::vector cfgs{baseline, all, noSurfaceArea,