Skip to content

Commit

Permalink
Add GEOSPreparedContainsXY
Browse files Browse the repository at this point in the history
  • Loading branch information
@dbaston
dbaston committed Sep 9, 2022
1 parent e3169bd commit c812592
Show file tree
Hide file tree
Showing 5 changed files with 107 additions and 43 deletions.
104 changes: 62 additions & 42 deletions benchmarks/capi/GEOSPreparedContainsPerfTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,80 +32,100 @@
class GEOSPreparedContainsPerfTest {

public:
void test(const GEOSGeometry* g, std::size_t num_points) {
void test(const std::vector<GEOSGeometry*>& geoms, std::size_t num_points) {
using namespace geos::geom;

double xmin, xmax, ymin, ymax;
geos::util::Profile sw("GEOSPreparedContains");
geos::util::Profile sw_xy("GEOSPreparedContainsXY");

GEOSGeom_getXMin(g, &xmin);
GEOSGeom_getXMax(g, &xmax);
GEOSGeom_getYMin(g, &ymin);
GEOSGeom_getYMax(g, &ymax);
std::size_t hits = 0;
std::size_t hits_xy = 0;

std::default_random_engine e(12345);
std::uniform_real_distribution<> xdist(xmin, xmax);
std::uniform_real_distribution<> ydist(ymin, ymax);
for (const auto& g: geoms) {
double xmin, xmax, ymin, ymax;

std::vector<Coordinate> coords(num_points);
std::generate(coords.begin(), coords.end(), [&xdist, &ydist, &e]() {
return Coordinate(xdist(e), ydist(e));
});
GEOSGeom_getXMin(g, &xmin);
GEOSGeom_getXMax(g, &xmax);
GEOSGeom_getYMin(g, &ymin);
GEOSGeom_getYMax(g, &ymax);

geos::util::Profile sw("GEOSPreparedContains");
sw.start();

std::size_t hits = 0;
auto prep = GEOSPrepare(g);
for (const auto& c : coords) {
auto pt = GEOSGeom_createPointFromXY(c.x, c.y);
std::default_random_engine e(12345);
std::uniform_real_distribution<> xdist(xmin, xmax);
std::uniform_real_distribution<> ydist(ymin, ymax);

if (GEOSPreparedContains(prep, pt)) {
hits++;
}
std::vector<Coordinate> coords(num_points);
std::generate(coords.begin(), coords.end(), [&xdist, &ydist, &e]() {
return Coordinate(xdist(e), ydist(e));
});

GEOSGeom_destroy(pt);
}

GEOSPreparedGeom_destroy(prep);
sw.start();
auto prep = GEOSPrepare(g);
for (const auto& c : coords) {
auto pt = GEOSGeom_createPointFromXY(c.x, c.y);
if (GEOSPreparedContains(prep, pt)) {
hits++;
}

sw.stop();
GEOSGeom_destroy(pt);
}
GEOSPreparedGeom_destroy(prep);
sw.stop();

sw_xy.start();
prep = GEOSPrepare(g);
for (const auto& c : coords) {
if (GEOSPreparedContainsXY(prep, c.x, c.y)) {
hits_xy++;
}
}
GEOSPreparedGeom_destroy(prep);
sw_xy.stop();
}

std::cout << sw.name << ": " << hits << " hits from " << num_points << " points in " << sw.getTotFormatted() << std::endl;

std::cout << sw_xy.name << ": " << hits_xy << " hits from " << num_points << " points in " << sw_xy.getTotFormatted() << std::endl;
}
};

int main(int argc, char** argv) {
if (argc != 3) {
std::cout << "perf_geospreparedcontins performs a specified number of point-in-polygon tests" << std::endl;
std::cout << "on randomly generated points from the bounding box of a single geometry provided" << std::endl;
std::cout << "perf_geospreparedcontians performs a specified number of point-in-polygon tests" << std::endl;
std::cout << "on randomly generated points from the bounding box of a each geometry provided" << std::endl;
std::cout << "in a file as WKT." << std::endl;
std::cout << std::endl;
std::cout << "Usage: perf_geospreparedcontins [wktfile] [n]" << std::endl;
std::cout << "Usage: perf_geospreparedcontains [wktfile] [n]" << std::endl;
return 0;
}

initGEOS(nullptr, nullptr);

GEOSPreparedContainsPerfTest tester;

std::size_t n = static_cast<std::size_t>(std::atoi(argv[2]));
std::cout << "Performing " << n << " point-in-polygon tests." << std::endl;

std::string fname{argv[1]};
std::cout << "Reading shape from " << fname << std::endl;
std::cout << "Reading shapes from " << fname << std::endl;

std::vector<GEOSGeometry*> geoms;

std::ifstream f(fname);
std::stringstream buff;
buff << f.rdbuf();
std::string line;
long i = 0;
while(std::getline(f, line)) {
if (line != "") {
geoms.push_back(GEOSGeomFromWKT(line.c_str()));
i++;
}
}
f.close();

std::string wkt = buff.str();
buff.clear();

initGEOS(nullptr, nullptr);
GEOSGeometry* g = GEOSGeomFromWKT(wkt.c_str());
wkt.clear();
std::cout << "Read " << geoms.size() << " geometries." << std::endl;

tester.test(g, n);
tester.test(geoms, n);

GEOSGeom_destroy(g);
for (auto& g : geoms) {
GEOSGeom_destroy(g);
}
}
6 changes: 6 additions & 0 deletions capi/geos_c.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1449,6 +1449,12 @@ extern "C" {
return GEOSPreparedContains_r(handle, pg1, g2);
}

char
GEOSPreparedContainsXY(const geos::geom::prep::PreparedGeometry* pg1, double x, double y)
{
return GEOSPreparedContainsXY_r(handle, pg1, x, y);
}

char
GEOSPreparedContainsProperly(const geos::geom::prep::PreparedGeometry* pg1, const Geometry* g2)
{
Expand Down
11 changes: 11 additions & 0 deletions capi/geos_c.h.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1111,6 +1111,12 @@ extern char GEOS_DLL GEOSPreparedContains_r(
const GEOSPreparedGeometry* pg1,
const GEOSGeometry* g2);

extern char GEOS_DLL GEOSPreparedContainsXY_r(
GEOSContextHandle_t handle,
const GEOSPreparedGeometry* pg1,
double x,
double y);

/** \see GEOSPreparedContainsProperly */
extern char GEOS_DLL GEOSPreparedContainsProperly_r(
GEOSContextHandle_t handle,
Expand Down Expand Up @@ -4284,6 +4290,11 @@ extern char GEOS_DLL GEOSPreparedContains(
const GEOSPreparedGeometry* pg1,
const GEOSGeometry* g2);

extern char GEOS_DLL GEOSPreparedContainsXY(
const GEOSPreparedGeometry* pg1,
double x,
double y);

/**
* Using a \ref GEOSPreparedGeometry do a high performance
* calculation of whether the provided geometry is contained properly.
Expand Down
16 changes: 15 additions & 1 deletion capi/geos_ts_c.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -207,6 +207,8 @@ typedef struct GEOSContextHandle_HS {
uint8_t WKBOutputDims;
int WKBByteOrder;
int initialized;
std::unique_ptr<geos::geom::Point> point2d;
geos::geom::Coordinate* point2dc;

GEOSContextHandle_HS()
:
Expand All @@ -216,10 +218,13 @@ typedef struct GEOSContextHandle_HS {
noticeData(nullptr),
errorMessageOld(nullptr),
errorMessageNew(nullptr),
errorData(nullptr)
errorData(nullptr),
point2d(nullptr)
{
memset(msgBuffer, 0, sizeof(msgBuffer));
geomFactory = GeometryFactory::getDefaultInstance();
point2d.reset(geomFactory->createPoint({1, 1}));
point2dc = (geos::geom::Coordinate*) point2d->getCoordinate();
WKBOutputDims = 2;
WKBByteOrder = getMachineByteOrder();
setNoticeHandler(nullptr);
Expand Down Expand Up @@ -3395,6 +3400,15 @@ extern "C" {
});
}

char
GEOSPreparedContainsXY_r(GEOSContextHandle_t extHandle,
const geos::geom::prep::PreparedGeometry* pg, double x, double y)
{
extHandle->point2d->setXY(x, y);

return GEOSPreparedContains_r(extHandle, pg, extHandle->point2d.get());
}

char
GEOSPreparedContainsProperly_r(GEOSContextHandle_t extHandle,
const geos::geom::prep::PreparedGeometry* pg, const Geometry* g)
Expand Down
13 changes: 13 additions & 0 deletions tests/unit/capi/GEOSPreparedGeometryTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -422,6 +422,19 @@ void object::test<14>
ensure_equals(ret, 2);
}

template<>
template<>
void object::test<15>
()
{
geom1_ = GEOSGeomFromWKT("POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))") ;
prepGeom1_ = GEOSPrepare(geom1_);

ensure_equals(GEOSPreparedContainsXY(prepGeom1_, 0.5, 0.5), 1);
ensure_equals(GEOSPreparedContainsXY(prepGeom1_, 1.5, 0.5), 0);
ensure_equals(GEOSPreparedContainsXY(prepGeom1_, 0.75, 0.5), 1);
}


} // namespace tut

0 comments on commit c812592

Please sign in to comment.