Skip to content

Commit

Permalink
efficient --within-distance for point/polygon and point/linestring
Browse files Browse the repository at this point in the history
  • Loading branch information
@patrickbr
patrickbr committed Feb 25, 2025
1 parent 853e633 commit d22cf72
Show file tree
Hide file tree
Showing 3 changed files with 227 additions and 28 deletions.
243 changes: 216 additions & 27 deletions src/spatialjoin/Sweeper.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,8 @@ using sj::boxids::packBoxIds;
using sj::innerouter::Mode;
using util::writeAll;
using util::geo::area;
using util::geo::getBoundingBox;
using util::geo::FPoint;
using util::geo::getBoundingBox;
using util::geo::I32Box;
using util::geo::I32Line;
using util::geo::I32MultiLine;
Expand Down Expand Up @@ -337,8 +337,8 @@ 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 = getPaddedBoundingBox(polyR);
auto lineR = util::geo::rotateSinCos(line, sin45, cos45, I32Point(0, 0));
box45 = getPaddedBoundingBox(lineR);
}

cur.subid = subid;
Expand Down Expand Up @@ -1787,7 +1787,7 @@ void Sweeper::writeRel(size_t t, const std::string& a, const std::string& b,

// ____________________________________________________________________________
void Sweeper::writeDist(size_t t, const std::string& a, size_t aSub,
const std::string& b, size_t bSub, double dist) {
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");
}
Expand Down Expand Up @@ -1860,7 +1860,7 @@ void Sweeper::selfCheck(const BoxVal cur, size_t t) {
}

// ____________________________________________________________________________
void Sweeper::doCheck(const BoxVal cur, const SweepVal sv, size_t t) {
void Sweeper::doDistCheck(const BoxVal cur, const SweepVal sv, size_t t) {
_checks[t]++;
_curX[t] = cur.val;

Expand All @@ -1869,26 +1869,128 @@ 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));
if (false && 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)});
auto dist = meterDist(p1, p2);

if (dist <= _cfg.withinDist) {
auto a = _pointCache.get(cur.id, cur.large ? -1 : t);
auto b = _pointCache.get(sv.id, sv.large ? -1 : t);
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);
}
writeDist(t, a->id, a->subId, b->id, b->subId, dist);
}
} else if (cur.type == POINT &&
(sv.type == POLYGON || sv.type == SIMPLE_POLYGON)) {
auto p = util::geo::centroid(cur.b45);
p = util::geo::rotateSinCos(p, -sin45, cos45, I32Point(0, 0));

const Area* a;
std::shared_ptr<Area> asp;
Area al;

if (sv.type == SIMPLE_POLYGON) {
auto p = _simpleAreaCache.get(sv.id, sv.large ? -1 : t);
al = areaFromSimpleArea(p.get());
a = &al;
} else {
asp = _areaCache.get(sv.id, sv.large ? -1 : t);
a = asp.get();
}

double dist = distCheck(p, a, t);

if (dist <= _cfg.withinDist) {
auto b = _pointCache.get(cur.id, cur.large ? -1 : t);
writeDist(t, a->id, a->subId, b->id, b->subId, dist);
}
} else if ((cur.type == POLYGON || cur.type == SIMPLE_POLYGON) &&
sv.type == POINT) {
auto p = util::geo::centroid(sv.b45);
p = util::geo::rotateSinCos(p, -sin45, cos45, I32Point(0, 0));

const Area* a;
std::shared_ptr<Area> asp;
Area al;

if (cur.type == SIMPLE_POLYGON) {
auto p = _simpleAreaCache.get(cur.id, cur.large ? -1 : t);
al = areaFromSimpleArea(p.get());
a = &al;
} else {
asp = _areaCache.get(cur.id, cur.large ? -1 : t);
a = asp.get();
}

double dist = distCheck(p, a, t);

if (dist <= _cfg.withinDist) {
auto b = _pointCache.get(sv.id, sv.large ? -1 : t);
writeDist(t, a->id, a->subId, b->id, b->subId, dist);
}
} else if ((cur.type == SIMPLE_LINE || cur.type == LINE) && sv.type == POINT) {
auto p = util::geo::centroid(sv.b45);
p = util::geo::rotateSinCos(p, -sin45, cos45, I32Point(0, 0));

double dist = std::numeric_limits<double>::max();

if (cur.type == SIMPLE_LINE) {
auto b = _simpleLineCache.get(cur.id, cur.large ? -1 : t);
dist = distCheck(p, b.get(), t);

if (dist <= _cfg.withinDist) {
auto a = _pointCache.get(sv.id, sv.large ? -1 : t);
writeDist(t, a->id, a->subId, b->id, 0, dist);
}
} else {
auto b = _lineCache.get(cur.id, cur.large ? -1 : t);
dist = distCheck(p, b.get(), t);

if (dist <= _cfg.withinDist) {
auto a = _pointCache.get(sv.id, sv.large ? -1 : t);
writeDist(t, a->id, a->subId, b->id, b->subId, dist);
}
}
} else if ((sv.type == SIMPLE_LINE || sv.type == LINE) && cur.type == POINT) {
auto p = util::geo::centroid(cur.b45);
p = util::geo::rotateSinCos(p, -sin45, cos45, I32Point(0, 0));

double dist = std::numeric_limits<double>::max();

if (sv.type == SIMPLE_LINE) {
auto b = _simpleLineCache.get(sv.id, sv.large ? -1 : t);
dist = distCheck(p, b.get(), t);

if (dist <= _cfg.withinDist) {
auto a = _pointCache.get(cur.id, cur.large ? -1 : t);
writeDist(t, a->id, a->subId, b->id, 0, dist);
}
} else {
auto b = _lineCache.get(sv.id, sv.large ? -1 : t);
dist = distCheck(p, b.get(), t);

return;
if (dist <= _cfg.withinDist) {
auto a = _pointCache.get(cur.id, cur.large ? -1 : t);
writeDist(t, a->id, a->subId, b->id, b->subId, dist);
}
}
}

}

// ____________________________________________________________________________
void Sweeper::doCheck(const BoxVal cur, const SweepVal sv, size_t t) {
_checks[t]++;
_curX[t] = cur.val;

if (cur.type == sv.type && cur.id == sv.id) return selfCheck(cur, t);

// every 10000 checks, update our position
if (_checks[t] % 10000 == 0) _atomicCurX[t] = _curX[t];

if ((cur.type == SIMPLE_POLYGON || cur.type == POLYGON) &&
(sv.type == SIMPLE_POLYGON || sv.type == POLYGON)) {
auto ts = TIME();
Expand Down Expand Up @@ -2756,10 +2858,15 @@ void Sweeper::processQueue(size_t t) {
JobBatch batch;
while ((batch = _jobs.get()).size()) {
for (const auto& job : batch) {
if (job.multiOut.empty())
doCheck(job.boxVal, job.sweepVal, t);
else
if (job.multiOut.empty()) {
if (_cfg.withinDist >= 0) {
doDistCheck(job.boxVal, job.sweepVal, t);
} else {
doCheck(job.boxVal, job.sweepVal, t);
}
} else {
multiOut(t, job.multiOut);
}
}
}
} catch (const std::runtime_error& e) {
Expand Down Expand Up @@ -3174,20 +3281,40 @@ std::string Sweeper::intToBase126(uint64_t id) {
return ret;
}

// _____________________________________________________________________________
double Sweeper::getMaxScaleFactor(const I32Box& bbox) const {
double invXScaleFactor = 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}));

return 1.0 / invXScaleFactor;
}

// _____________________________________________________________________________
double Sweeper::getMaxScaleFactor(const I32Point& p) const {
return 1.0 / util::geo::webMercDistFactor(
I32Point{p.getX() / PREC, p.getY() / PREC});
}

// _____________________________________________________________________________
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}));
auto bbox = getBoundingBox(geom);
double xScaleFactor = getMaxScaleFactor(bbox);

uint32_t padX = (_cfg.withinDist / xScaleFactor) * PREC;
uint32_t padY = _cfg.withinDist * PREC;
uint32_t padX = ((_cfg.withinDist / 2) * xScaleFactor) * PREC;
uint32_t padY = (_cfg.withinDist / 2) * PREC;

return {{bbox.getLowerLeft().getX() - padX, bbox.getLowerLeft().getY() - padY}, {bbox.getUpperRight().getX() + padX, bbox.getUpperRight().getY() + padY}};
return {
{bbox.getLowerLeft().getX() - padX, bbox.getLowerLeft().getY() - padY},
{bbox.getUpperRight().getX() + padX,
bbox.getUpperRight().getY() + padY}};
}

return bbox;
Expand All @@ -3203,3 +3330,65 @@ uint64_t Sweeper::base126ToInt(const std::string& id) {
}
return ret;
}

// _____________________________________________________________________________
double Sweeper::meterDist(const I32Point& p1, const I32Point& p2) {
return 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)});
}

// _____________________________________________________________________________
double Sweeper::distCheck(const I32Point& a, const Area* b, size_t t) const {
if (_cfg.useBoxIds) {
auto ts = TIME();
auto r = boxIdIsect({{1, 0}, {getBoxId(a), 0}}, b->boxIds);
_stats[t].timeBoxIdIsectAreaPoint += TOOK(ts);

// all boxes of a are fully contained in b, we are contained
if (r.first) return 0;
}

auto ts = TIME();
double scaleFactor =
std::max(getMaxScaleFactor(a), getMaxScaleFactor(b->box)) * PREC;

auto dist = util::geo::withinDist<int32_t>(
a, b->geom, _cfg.withinDist * scaleFactor, _cfg.withinDist,
&Sweeper::meterDist);
_stats[t].timeFullGeoCheckAreaPoint += TOOK(ts);
_stats[t].fullGeoChecksAreaPoint++;

return dist;
}

// _____________________________________________________________________________
double Sweeper::distCheck(const I32Point& a, const Line* b, size_t t) const {
auto ts = TIME();
double scaleFactor =
std::max(getMaxScaleFactor(a), getMaxScaleFactor(b->box)) * PREC;

auto dist = util::geo::withinDist<int32_t>(
a, b->geom, _cfg.withinDist * scaleFactor, _cfg.withinDist,
&Sweeper::meterDist);
_stats[t].timeFullGeoCheckLinePoint += TOOK(ts);
_stats[t].fullGeoChecksLinePoint++;

return dist;
}

// _____________________________________________________________________________
double Sweeper::distCheck(const I32Point& a, const SimpleLine* b, size_t t) const {
auto ts = TIME();

auto p2 = projectOn(b->a, a,
b->b);

auto dist = Sweeper::meterDist(a, p2);

_stats[t].timeFullGeoCheckLinePoint += TOOK(ts);
_stats[t].fullGeoChecksLinePoint++;

return dist;
}
10 changes: 10 additions & 0 deletions src/spatialjoin/Sweeper.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -362,11 +362,18 @@ class Sweeper {
std::tuple<bool, bool> check(const util::geo::I32Point& a, const Line* b,
size_t t) const;

double distCheck(const util::geo::I32Point& a, const Area* b, size_t t) const;
double distCheck(const util::geo::I32Point& a, const Line* b, size_t t) const;
double distCheck(const util::geo::I32Point& a, const SimpleLine* b, size_t t) const;

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

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

double getMaxScaleFactor(const util::geo::I32Box& geom) const;
double getMaxScaleFactor(const util::geo::I32Point& geom) const;

void diskAdd(const BoxVal& bv);

void multiOut(size_t t, const std::string& gid);
Expand Down Expand Up @@ -403,6 +410,7 @@ class Sweeper {
const std::string& b, size_t bSub);

void doCheck(BoxVal cur, SweepVal sv, size_t t);
void doDistCheck(BoxVal cur, SweepVal sv, size_t t);
void selfCheck(BoxVal cur, size_t t);

void processQueue(size_t t);
Expand All @@ -416,6 +424,8 @@ class Sweeper {
bool notTouches(const std::string& a, const std::string& b);
bool notCrosses(const std::string& a, const std::string& b);

static double meterDist(const util::geo::I32Point& p1, const util::geo::I32Point& p2);

void fillBatch(JobBatch* batch,
const sj::IntervalIdx<int32_t, SweepVal>* actives,
const BoxVal* cur) const;
Expand Down
2 changes: 1 addition & 1 deletion src/util
Submodule util updated 2 files
+168 −20 geo/Geo.h
+1 −0 geo/Polygon.h

0 comments on commit d22cf72

Please sign in to comment.