Skip to content

Commit

Permalink
GEOSClipByRect: preserve Z coordinate (#1095)
Browse files Browse the repository at this point in the history 
For [multi]linestrings/linerings, Z coordinate is interpolated
from that of the nearest existing vertices.

For [multi]polygons, Z coordinate of points created on the boundary
is interpolated as above.  For interior points, however, it is populated
using the ElevationModel.
  • Loading branch information
@vadimcn
vadimcn authored May 23, 2024
1 parent 1945569 commit a13c2aa
Show file tree
Hide file tree
Showing 2 changed files with 121 additions and 24 deletions.
54 changes: 35 additions & 19 deletions src/operation/intersection/RectangleIntersection.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
#include <geos/operation/intersection/RectangleIntersection.h>
#include <geos/operation/intersection/Rectangle.h>
#include <geos/operation/intersection/RectangleIntersectionBuilder.h>
#include <geos/operation/overlayng/ElevationModel.h>
#include <geos/operation/predicate/RectangleIntersects.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/geom/CoordinateSequence.h>
Expand All @@ -34,6 +35,7 @@

using geos::operation::intersection::Rectangle;
using geos::operation::intersection::RectangleIntersectionBuilder;
using geos::operation::overlayng::ElevationModel;
using namespace geos::geom;
using namespace geos::algorithm;
namespace geos {
Expand Down Expand Up @@ -62,15 +64,18 @@ different(double x1, double y1, double x2, double y2)

inline
void
clip_one_edge(double& x1, double& y1, double x2, double y2, double limit)
clip_one_edge(double& x1, double& y1, double& z1, double x2, double y2, double z2, double limit)
{
if(x2 == limit) {
y1 = y2;
x1 = x2;
z1 = z2;
}

if(x1 != x2) {
y1 += (y2 - y1) * (limit - x1) / (x2 - x1);
double fraction = (limit - x1) / (x2 - x1);
y1 += (y2 - y1) * fraction;
z1 += (z2 - z1) * fraction;
x1 = limit;
}
}
Expand All @@ -86,22 +91,22 @@ clip_one_edge(double& x1, double& y1, double x2, double y2, double limit)
*/

void
clip_to_edges(double& x1, double& y1,
double x2, double y2,
clip_to_edges(double& x1, double& y1, double& z1,
double x2, double y2, double z2,
const Rectangle& rect)
{
if(x1 < rect.xmin()) {
clip_one_edge(x1, y1, x2, y2, rect.xmin());
clip_one_edge(x1, y1, z1, x2, y2, z2, rect.xmin());
}
else if(x1 > rect.xmax()) {
clip_one_edge(x1, y1, x2, y2, rect.xmax());
clip_one_edge(x1, y1, z1, x2, y2, z2, rect.xmax());
}

if(y1 < rect.ymin()) {
clip_one_edge(y1, x1, y2, x2, rect.ymin());
clip_one_edge(y1, x1, z1, y2, x2, z2, rect.ymin());
}
else if(y1 > rect.ymax()) {
clip_one_edge(y1, x1, y2, x2, rect.ymax());
clip_one_edge(y1, x1, z1, y2, x2, z2, rect.ymax());
}
}

Expand Down Expand Up @@ -153,6 +158,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,

double x0 = 0;
double y0 = 0;
double z0 = 0;
bool add_start = false;

// Start iterating
Expand All @@ -164,6 +170,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,

double x = cs[i].x;
double y = cs[i].y;
double z = cs[i].z;
Rectangle::Position pos = rect.position(x, y);

if(pos == Rectangle::Outside) {
Expand Down Expand Up @@ -199,12 +206,14 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
// Establish new position
x = cs[i].x;
y = cs[i].y;
z = cs[i].z;
pos = rect.position(x, y);

// Handle all possible cases
x0 = cs[i - 1].x;
y0 = cs[i - 1].y;
clip_to_edges(x0, y0, x, y, rect);
z0 = cs[i - 1].z;
clip_to_edges(x0, y0, z0, x, y, z, rect);

if(pos == Rectangle::Inside) {
add_start = true; // x0,y0 must have clipped the rectangle
Expand All @@ -218,7 +227,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
// which will then enter the Outside section.

// Clip the other end too
clip_to_edges(x, y, x0, y0, rect);
clip_to_edges(x, y, z, x0, y0, z0, rect);

Rectangle::Position prev_pos = rect.position(x0, y0);
pos = rect.position(x, y);
Expand All @@ -229,8 +238,8 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
!Rectangle::onSameEdge(prev_pos, pos) // discard if travels along edge
) {
auto coords = detail::make_unique<geom::CoordinateSequence>(2u);
coords->setAt(Coordinate(x0, y0), 0);
coords->setAt(Coordinate(x, y), 1);
coords->setAt(Coordinate(x0, y0, z0), 0);
coords->setAt(Coordinate(x, y, z), 1);
auto line = _gf->createLineString(std::move(coords));
parts.add(line.release());
}
Expand Down Expand Up @@ -267,6 +276,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
while(!go_outside && ++i < n) {
x = cs[i].x;
y = cs[i].y;
z = cs[i].z;

Rectangle::Position prev_pos = pos;
pos = rect.position(x, y);
Expand All @@ -278,7 +288,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
go_outside = true;

// Clip the outside point to edges
clip_to_edges(x, y, cs[i - 1].x, cs[i - 1].y, rect);
clip_to_edges(x, y, z, cs[i - 1].x, cs[i - 1].y, cs[i - 1].z, rect);
pos = rect.position(x, y);

// Does the line exit through the inside of the box?
Expand All @@ -291,15 +301,15 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
if(start_index < i - 1 || add_start || through_box) {
auto coords = detail::make_unique<CoordinateSequence>();
if(add_start) {
coords->add(Coordinate(x0, y0));
coords->add(Coordinate(x0, y0, z0));
add_start = false;
}
//line->addSubLineString(&g, start_index, i-1);
coords->add(cs.begin() + static_cast<long>(start_index),
cs.begin() + static_cast<long>(i));

if(through_box) {
coords->add(Coordinate(x, y));
coords->add(Coordinate(x, y, z));
}

auto line = _gf->createLineString(std::move(coords));
Expand All @@ -316,7 +326,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
//geom::LineString * line = new geom::LineString();
if(add_start) {
//line->addPoint(x0,y0);
coords->add(Coordinate(x0, y0));
coords->add(Coordinate(x0, y0, z0));
add_start = false;
}
//line->addSubLineString(&g, start_index, i-1);
Expand Down Expand Up @@ -349,7 +359,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
//geom::LineString * line = new geom::LineString();
if(add_start) {
//line->addPoint(x0,y0);
coords->add(Coordinate(x0, y0));
coords->add(Coordinate(x0, y0, z0));
add_start = false;
}
//line->addSubLineString(&g, start_index, i-1);
Expand Down Expand Up @@ -521,7 +531,6 @@ RectangleIntersection::clip_polygon_to_polygons(const geom::Polygon* g,

parts.reconnectPolygons(rect);
parts.release(toParts);

}

/**
Expand Down Expand Up @@ -682,7 +691,14 @@ std::unique_ptr<geom::Geometry>
RectangleIntersection::clip(const geom::Geometry& g, const Rectangle& rect)
{
RectangleIntersection ri(g, rect);
return ri.clip();
std::unique_ptr<geom::Geometry> result = ri.clip();

if (g.hasZ()) {
std::unique_ptr<ElevationModel> elevModel = ElevationModel::create(g);
elevModel->populateZ(*result);
}

return result;
}

std::unique_ptr<geom::Geometry>
Expand Down
91 changes: 86 additions & 5 deletions tests/unit/capi/GEOSClipByRectTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,27 @@ struct test_capigeosclipbyrect_data : public capitest::utility {
ensure(equal);
}

void
checkClipByRectIdentical(
const char* wktInput,
double xmin, double ymin,
double xmax, double ymax,
const char* wktExpected)
{
input_ = fromWKT(wktInput);
result_ = GEOSClipByRect(input_, xmin, ymin, xmax, ymax);
expected_ = fromWKT(wktExpected);

bool equal = GEOSEqualsIdentical(result_, expected_) == 1;
if (!equal) {
wkt_ = GEOSWKTWriter_write(wktw_, expected_);
std::printf("EXP: %s\n", wkt_);
GEOSFree(wkt_);
wkt_ = GEOSWKTWriter_write(wktw_, result_);
std::printf("OBT: %s\n", wkt_);
}
ensure(equal);
}
};

typedef test_group<test_capigeosclipbyrect_data> group;
Expand Down Expand Up @@ -70,7 +91,7 @@ template<>
template<>
void object::test<2>()
{
checkClipByRect(
checkClipByRectIdentical(
"POINT(15 15)",
10, 10, 20, 20,
"POINT(15 15)"
Expand Down Expand Up @@ -106,7 +127,7 @@ template<>
template<>
void object::test<5>()
{
checkClipByRect(
checkClipByRectIdentical(
"LINESTRING(15 15, 16 15)",
10, 10, 20, 20,
"LINESTRING(15 15, 16 15)"
Expand All @@ -130,7 +151,7 @@ template<>
template<>
void object::test<7>()
{
checkClipByRect(
checkClipByRectIdentical(
"LINESTRING(10 5, 25 20)",
10, 10, 20, 20,
"LINESTRING (15 10, 20 15)"
Expand Down Expand Up @@ -191,7 +212,7 @@ template<>
void object::test<12>()
{
const char* wkt = "POLYGON((1 1, 1 30, 30 30, 30 1, 1 1),(10 10, 20 10, 20 20, 10 20, 10 10))";
checkClipByRect(
checkClipByRectIdentical(
wkt,
0, 0, 40, 40,
wkt);
Expand All @@ -202,7 +223,7 @@ template<>
template<>
void object::test<13>()
{
checkClipByRect(
checkClipByRectIdentical(
"POLYGON((0 0, 0 30, 30 30, 30 0, 0 0),(10 10, 20 10, 20 20, 10 20, 10 10))",
5, 5, 15, 15,
"POLYGON ((5 5, 5 15, 10 15, 10 10, 15 10, 15 5, 5 5))"
Expand Down Expand Up @@ -282,6 +303,66 @@ void object::test<17>()
ensure(result_ == nullptr);
}

/// Point Z inside
template<>
template<>
void object::test<18>()
{
checkClipByRectIdentical(
"POINT Z (15 15 100)",
10, 10, 20, 20,
"POINT Z (15 15 100)"
);
}

/// Line Z inside
template<>
template<>
void object::test<19>()
{
checkClipByRectIdentical(
"LINESTRING Z (15 15 0, 16 15 100)",
10, 10, 20, 20,
"LINESTRING Z (15 15 0, 16 15 100)"
);
}

/// Line Z splitting rectangle
template<>
template<>
void object::test<20>()
{
checkClipByRectIdentical(
"LINESTRING Z (0 15 0, 100 15 100)",
10, 10, 20, 20,
"LINESTRING Z (10 15 10, 20 15 20)"
);
}

/// Polygon Z overlapping rectangle
template<>
template<>
void object::test<21>()
{
checkClipByRectIdentical(
"POLYGON Z ((0 0 100, 0 30 100, 30 30 100, 30 0 100, 0 0 100),(10 10 100, 20 10 100, 20 20 100, 10 20 100, 10 10 100))",
5, 5, 15, 15,
"POLYGON Z ((5 5 100, 5 15 100, 10 15 100, 10 10 100, 15 10 100, 15 5 100, 5 5 100))"
);
}

/// Polygon Z enclosing rectangle
template<>
template<>
void object::test<22>()
{
checkClipByRectIdentical(
"POLYGON Z ((0 0 100, 0 30 100, 30 30 100, 30 0 100, 0 0 100))",
5, 5, 15, 15,
"POLYGON Z ((5 5 100, 5 15 100, 15 15 100, 15 5 100, 5 5 100))"
);
}


} // namespace tut

0 comments on commit a13c2aa

Please sign in to comment.