Skip to content

Commit

Permalink
first rough implementation of --within-distance, already works for po…
Browse files Browse the repository at this point in the history 
…ints
  • Loading branch information
@patrickbr
patrickbr committed Feb 21, 2025
1 parent 51e2321 commit 853e633
Show file tree
Hide file tree
Showing 4 changed files with 97 additions and 23 deletions.
10 changes: 10 additions & 0 deletions src/spatialjoin/SpatialJoinMain.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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") {
Expand Down Expand Up @@ -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;
}
}

Expand Down Expand Up @@ -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; },
Expand Down
86 changes: 72 additions & 14 deletions src/spatialjoin/Sweeper.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,8 @@

#include <algorithm>
#include <cassert>
#include <cmath>
#include <climits>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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<int32_t, size_t> cutouts;

Expand All @@ -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;
Expand Down Expand Up @@ -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};

Expand Down Expand Up @@ -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});
}
}

Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -3135,6 +3174,25 @@ std::string Sweeper::intToBase126(uint64_t id) {
return ret;
}

// _____________________________________________________________________________
template <typename G>
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;
Expand Down
6 changes: 6 additions & 0 deletions src/spatialjoin/Sweeper.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -137,6 +137,7 @@ struct SweeperCfg {
bool useFastSweepSkip;
bool useInnerOuter;
bool noGeometryChecks;
double withinDist;
std::function<void(size_t t, const char* a, const char* b, const char* pred)>
writeRelCb;
std::function<void(const std::string&)> logCb;
Expand Down Expand Up @@ -363,6 +364,9 @@ class Sweeper {

bool refRelated(const std::string& a, const std::string& b) const;

template <typename G>
util::geo::I32Box getPaddedBoundingBox(const G& geom) const;

void diskAdd(const BoxVal& bv);

void multiOut(size_t t, const std::string& gid);
Expand All @@ -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,
Expand Down
18 changes: 9 additions & 9 deletions src/spatialjoin/tests/TestMain.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<sj::SweeperCfg> cfgs{baseline, all, noSurfaceArea,
Expand Down

0 comments on commit 853e633

Please sign in to comment.