diff --git a/Makefile b/Makefile index b651ec64..57227d4c 100644 --- a/Makefile +++ b/Makefile @@ -157,6 +157,7 @@ test: \ test_attribute_store \ test_config_validator \ test_deque_map \ + test_geom \ test_helpers \ test_options_parser \ test_pbf_reader \ @@ -196,6 +197,11 @@ test_deque_map: \ test/deque_map.test.o $(CXX) $(CXXFLAGS) -o test.deque_map $^ $(INC) $(LIB) $(LDFLAGS) && ./test.deque_map +test_geom: \ + src/geom.o \ + test/geom.test.o + $(CXX) $(CXXFLAGS) -o test.geom $^ $(INC) $(LIB) $(LDFLAGS) && ./test.geom + test_helpers: \ src/helpers.o \ src/external/libdeflate/lib/adler32.o \ diff --git a/include/geom.h b/include/geom.h index a599da42..0a35df52 100644 --- a/include/geom.h +++ b/include/geom.h @@ -81,8 +81,10 @@ void make_valid(MultiPolygon &mp); // Attempt to repair an invalid areal geometry in place: dissolve-based // make_valid first (preserves area), then a zero-width buffer as a last // resort. Returns true if mp is valid afterwards; on failure mp is left as the -// best-effort input so callers never regress. -bool repair_multi_polygon(MultiPolygon &mp, bool knownInvalid = false); +// best-effort input so callers never regress. If wholeValid is given, it is set +// to true only if mp as a whole is known to pass geom::is_valid afterwards (a true +// result alone can leave individually valid polygons that overlap each other). +bool repair_multi_polygon(MultiPolygon &mp, bool knownInvalid = false, bool *wholeValid = nullptr); void union_many(std::vector &mps); @@ -91,6 +93,10 @@ char bit_code(Point const &p, Box const &bbox); // fast_clip returns true if the geometry was modified (false means identity). bool fast_clip(Ring &points, Box const &bbox); bool fast_clip(MultiPolygon &mp, Box const &bbox); +// Exact clip of a valid multipolygon to bbox that keeps it valid (fast_clip can +// leave zero-width bridges along the box edge). Returns false if it could not +// clip robustly; `out` is then unspecified and the caller must use geom::intersection. +bool rect_clip(const MultiPolygon &mp, Box const &bbox, MultiPolygon &out); // Cheap conservative validity tests, run before geom::is_valid. They return // true only when the geometry is certainly valid (is_valid would return true diff --git a/include/osm_mem_tiles.h b/include/osm_mem_tiles.h index e8d7ed9b..1966bc4b 100644 --- a/include/osm_mem_tiles.h +++ b/include/osm_mem_tiles.h @@ -61,7 +61,8 @@ class OsmMemTiles : public TileDataSource { // CorrectGeometry found them Valid (or Reversed: valid after reversing the ring) // on the same lon/latp coordinates. bool polygonKnownValid(NodeID objectID) const override { - return objectID >= OSM_THRESHOLD && (IS_WAY(objectID) || IS_WAY_REV(objectID)); + if (objectID < OSM_THRESHOLD) return TileDataSource::polygonKnownValid(objectID); + return IS_WAY(objectID) || IS_WAY_REV(objectID); } const NodeStore& nodeStore; diff --git a/include/tile_data.h b/include/tile_data.h index 5bc9fa67..3e8d20e9 100644 --- a/include/tile_data.h +++ b/include/tile_data.h @@ -298,6 +298,8 @@ class TileDataSource { std::vector linestringStores; std::vector multilinestringStores; std::vector multipolygonStores; + // Per multipolygonStores entry: stored geometry is known to pass geom::is_valid + std::vector> multipolygonValid; ClipCache multiPolygonClipCache; ClipCache multiLinestringClipCache; @@ -359,6 +361,7 @@ class TileDataSource { pointStores[0].push_back(Point(0,0)); linestringStores[0].push_back(linestring_t()); multipolygonStores[0].push_back(multi_polygon_t()); + multipolygonValid[0].push_back(false); multilinestringStores[0].push_back(multi_linestring_t()); } void reportSize() const; @@ -380,7 +383,11 @@ class TileDataSource { virtual void populateMultiPolygon(MultiPolygon& dst, NodeID objectID); // True if the polygon populateMultiPolygon returns for objectID is known to // have passed geom::is_valid on exactly the same coordinates. - virtual bool polygonKnownValid(NodeID objectID) const { return false; } + virtual bool polygonKnownValid(NodeID objectID) const { + const auto& valid = multipolygonValid[getShard(objectID)]; + const auto offset = getId(objectID); + return offset < valid.size() && valid[offset]; + } unsigned int getIndexZoom() const { return indexZoom; } inline size_t getId(NodeID id) const { @@ -415,7 +422,8 @@ class TileDataSource { return shard.at(offset); } - NodeID storeMultiPolygon(const MultiPolygon& src); + // knownValid: src has passed geom::is_valid (see polygonKnownValid) + NodeID storeMultiPolygon(const MultiPolygon& src, bool knownValid = false); multi_polygon_t const &retrieveMultiPolygon(NodeID id) const { const auto& shardId = getShard(id); diff --git a/src/geom.cpp b/src/geom.cpp index 50f54705..867a84f9 100644 --- a/src/geom.cpp +++ b/src/geom.cpp @@ -484,11 +484,15 @@ static bool repair_one_polygon(const Polygon &p, double minArea, MultiPolygon &o return false; } -bool repair_multi_polygon(MultiPolygon &mp, bool knownInvalid) +bool repair_multi_polygon(MultiPolygon &mp, bool knownInvalid, bool *wholeValid) { + if (wholeValid) *wholeValid = false; if (!knownInvalid) { const bool pre = multiPolygonDefinitelyValid(mp); - if (pre || geom::is_valid(mp)) return true; + if (pre || geom::is_valid(mp)) { + if (wholeValid) *wholeValid = true; + return true; + } } // Repair PER POLYGON, area-preserving. Running make_valid/buffer on the whole @@ -527,6 +531,8 @@ bool repair_multi_polygon(MultiPolygon &mp, bool knownInvalid) } mp = std::move(out); + // A single valid polygon is a valid multipolygon + if (wholeValid) *wholeValid = allValid && mp.size() == 1; return allValid; } @@ -699,3 +705,323 @@ bool fast_clip(MultiPolygon &mp, Box const &bbox) { mp.end()); return modified || mp.size() != sizeBefore; } + +// --------------- +// Exact clipping of a valid polygon to an axis-aligned box. +// +// Sutherland-Hodgman (fast_clip) is exact for convex input but leaves zero-width +// "bridges" along the box edge wherever a concave polygon leaves and re-enters +// the box; Boost reports those as self-intersections. This clipper instead cuts +// every ring into the chains that lie inside the box and reconnects them along +// the box boundary (Weiler-Atherton specialised to a rectangle), so a valid input +// gives a valid output. Whenever the input does not fit its assumptions (numerical +// trouble, a hole it cannot place) it gives up and returns false, and the caller +// falls back to Boost. + +namespace { + +struct RectClipper { + double x0, y0, x1, y1, w, h, perimeter; + + explicit RectClipper(const Box &box) : + x0(box.min_corner().x()), y0(box.min_corner().y()), + x1(box.max_corner().x()), y1(box.max_corner().y()), + w(x1 - x0), h(y1 - y0), perimeter(2 * (x1 - x0) + 2 * (y1 - y0)) {} + + bool inside(const Point &p) const { + return p.x() >= x0 && p.x() <= x1 && p.y() >= y0 && p.y() <= y1; + } + bool strictlyInside(const Point &p) const { + return p.x() > x0 && p.x() < x1 && p.y() > y0 && p.y() < y1; + } + // Segment a->b, both inside the box, runs through its interior (not along an edge) + bool throughInterior(const Point &a, const Point &b) const { + return strictlyInside(Point((a.x() + b.x()) / 2, (a.y() + b.y()) / 2)); + } + // Segment b->c, b on the boundary, starts off into the interior + bool entersInterior(const Point &b, const Point &c) const { + double t0, t1; + Point pa, pb; + return clipSegment(b, c, t0, t1, pa, pb) && t0 == 0 && t1 > 0 && throughInterior(b, pb); + } + // A ring p->t->q touches the boundary at t (not a corner) from inside. True if the + // polygon interior (right of the ring) covers the boundary on both sides of t, + // i.e. the edge direction lies in the clockwise sweep from t->q to t->p. + bool pinchAt(const Point &p, const Point &t, const Point &q) const { + if ((t.x() == x0 || t.x() == x1) && (t.y() == y0 || t.y() == y1)) return false; + const double ex = (t.y() == y0 || t.y() == y1) ? 1 : 0, ey = 1 - ex; + const double qx = q.x() - t.x(), qy = q.y() - t.y(); + auto clockwiseFromQ = [&](double vx, double vy) { + const double a = -std::atan2(qx * vy - qy * vx, qx * vx + qy * vy); + return a < 0 ? a + 2 * M_PI : a; + }; + return clockwiseFromQ(ex, ey) < clockwiseFromQ(p.x() - t.x(), p.y() - t.y()); + } + + // Position along the box boundary, clockwise from the top-left corner. + // Returns false if p is not exactly on the boundary. + bool boundaryPos(const Point &p, double &s) const { + if (p.y() == y1) s = p.x() - x0; + else if (p.x() == x1) s = w + (y1 - p.y()); + else if (p.y() == y0) s = w + h + (x1 - p.x()); + else if (p.x() == x0) s = 2 * w + h + (p.y() - y0); + else return false; + return true; + } + + // Liang-Barsky: the part [t0,t1] of segment a->b inside the closed box, with + // its end points snapped exactly onto the box edge they were clipped against. + bool clipSegment(const Point &a, const Point &b, double &t0, double &t1, Point &pa, Point &pb) const { + const double dx = b.x() - a.x(), dy = b.y() - a.y(); + const double p[4] = { -dx, dx, -dy, dy }; + const double q[4] = { a.x() - x0, x1 - a.x(), a.y() - y0, y1 - a.y() }; + int e0 = -1, e1 = -1; + t0 = 0; t1 = 1; + for (int k = 0; k < 4; k++) { + if (p[k] == 0) { + if (q[k] < 0) return false; + continue; + } + const double r = q[k] / p[k]; + if (p[k] < 0) { if (r > t0) { t0 = r; e0 = k; } } + else { if (r < t1) { t1 = r; e1 = k; } } + } + if (t0 > t1) return false; + pa = e0 < 0 ? a : onEdge(a, dx, dy, t0, e0); + pb = e1 < 0 ? b : onEdge(a, dx, dy, t1, e1); + return true; + } + + Point onEdge(const Point &a, double dx, double dy, double t, int edge) const { + switch (edge) { + case 0: return Point(x0, std::min(std::max(a.y() + t * dy, y0), y1)); + case 1: return Point(x1, std::min(std::max(a.y() + t * dy, y0), y1)); + case 2: return Point(std::min(std::max(a.x() + t * dx, x0), x1), y0); + default: return Point(std::min(std::max(a.x() + t * dx, x0), x1), y1); + } + } + + // Corners in clockwise boundary order, starting top-left + Point corner(int i) const { + switch (i) { + case 0: return Point(x0, y1); + case 1: return Point(x1, y1); + case 2: return Point(x1, y0); + default: return Point(x0, y0); + } + } + double cornerPos(int i) const { + switch (i) { + case 0: return 0; + case 1: return w; + case 2: return w + h; + default: return 2 * w + h; + } + } +}; + +enum class RingClip { Inside, Outside, Cut, Fail }; + +struct RectChain { + std::vector points; + double entry, exit; + bool used = false; + bool touchExit = false; // ends at a vertex touching the boundary from inside; the next chain starts there +}; + +// Cut a closed ring into the chains inside the box. Inside: every vertex is +// inside (so is the ring). Outside: no part of the ring enters the box. +RingClip clipRing(const RectClipper &rc, const Ring &ring, std::vector &chains) { + if (ring.size() < 4) return RingClip::Fail; + const std::size_t n = ring.size() - 1; // closed: ring[n] == ring[0] + std::size_t start = n; + for (std::size_t i = 0; i < n; i++) + if (!rc.inside(ring[i])) { start = i; break; } + if (start == n) return RingClip::Inside; + + // Start at a vertex outside the box, so every chain begins with an entry. + std::vector current; + bool building = false, any = false; + auto finish = [&]() -> bool { + RectChain c; + if (current.size() < 2 || !rc.boundaryPos(current.front(), c.entry) || !rc.boundaryPos(current.back(), c.exit)) return false; + c.points.swap(current); + chains.push_back(std::move(c)); + current.clear(); + building = false; + any = true; + return true; + }; + for (std::size_t k = 0; k < n; k++) { + const Point &a = ring[(start + k) % n]; + const Point &b = ring[(start + k + 1) % n]; + double t0, t1; + Point pa, pb; + const bool hit = rc.clipSegment(a, b, t0, t1, pa, pb) && t1 > t0; + if (!hit || (building && t0 > 0)) { + if (building && !finish()) return RingClip::Fail; + if (!hit) continue; + } + if (!building) { + current.push_back(pa); + building = true; + } + if (!geom::equals(pb, current.back())) current.push_back(pb); + if (t1 < 1) { + if (!finish()) return RingClip::Fail; + } else if (current.size() >= 2 && !rc.strictlyInside(b) && rc.throughInterior(current[current.size() - 2], b) && + rc.entersInterior(b, ring[(start + k + 2) % n]) && + rc.pinchAt(current[current.size() - 2], b, ring[(start + k + 2) % n])) { + // b touches the boundary from inside (reached through the interior, left + // into the interior) and the polygon lies along the boundary on both sides + // of it. Kept in one chain, the boundary walk would run through b and make + // the ring touch itself there: split the chain at b, so b is an exit and + // an entry, with the polygon pinched in two there. + if (!finish()) return RingClip::Fail; + chains.back().touchExit = true; + current.push_back(b); + building = true; + } + } + // The walk ends at the outside start vertex, so the last chain must have exited + if (building) return RingClip::Fail; + return any ? RingClip::Cut : RingClip::Outside; +} + +// Crossing-number test: 1 inside, -1 outside, 0 on the boundary. +int pointInRing(const Point &p, const Ring &ring) { + bool in = false; + for (std::size_t i = 0, j = ring.size() - 1; i < ring.size(); j = i++) { + const Point &a = ring[i], &b = ring[j]; + if (geom::equals(a, p)) return 0; + if ((a.y() > p.y()) != (b.y() > p.y())) { + const double x = a.x() + (p.y() - a.y()) * (b.x() - a.x()) / (b.y() - a.y()); + if (x == p.x()) return 0; + if (x > p.x()) in = !in; + } + } + return in ? 1 : -1; +} + +// Where a ring that does not cross `outer` lies relative to it (0: unsure). +int ringInRing(const Ring &ring, const Ring &outer) { + for (const Point &p : ring) { + const int r = pointInRing(p, outer); + if (r != 0) return r; + } + return 0; +} + +bool clipPolygon(const RectClipper &rc, const Polygon &poly, MultiPolygon &out) { + std::vector chains; + const RingClip outerClip = clipRing(rc, poly.outer(), chains); + if (outerClip == RingClip::Fail) return false; + if (outerClip == RingClip::Inside) { + out.push_back(poly); + return true; + } + + std::vector insideHoles, outsideHoles; + for (const Ring &inner : poly.inners()) { + switch (clipRing(rc, inner, chains)) { + case RingClip::Fail: return false; + case RingClip::Inside: insideHoles.push_back(&inner); break; + case RingClip::Outside: outsideHoles.push_back(&inner); break; + case RingClip::Cut: break; + } + } + + const std::size_t firstOuter = out.size(); + if (chains.empty()) { + // No ring crosses the box: it is either inside the polygon or disjoint from it + const Point centre((rc.x0 + rc.x1) / 2, (rc.y0 + rc.y1) / 2); + const int inOuter = pointInRing(centre, poly.outer()); + if (inOuter == 0) return false; + bool inside = inOuter > 0; + for (const Ring *hole : outsideHoles) { + const int inHole = pointInRing(centre, *hole); + if (inHole == 0) return false; + if (inHole > 0) inside = false; + } + if (!inside) return insideHoles.empty(); + Polygon boxPoly; + for (int i = 0; i < 4; i++) boxPoly.outer().push_back(rc.corner(i)); + boxPoly.outer().push_back(rc.corner(0)); + out.push_back(std::move(boxPoly)); + } else { + // Polygon interiors are on the right of every ring (clockwise outers, + // anticlockwise holes), so from each exit the boundary walk goes clockwise + // to the next entry. + std::vector> entries; + entries.reserve(chains.size()); + for (std::size_t i = 0; i < chains.size(); i++) entries.emplace_back(chains[i].entry, i); + std::sort(entries.begin(), entries.end()); + + for (std::size_t first = 0; first < chains.size(); first++) { + if (chains[first].used) continue; + Ring ring; + std::size_t c = first; + while (true) { + RectChain &chain = chains[c]; + chain.used = true; + for (const Point &p : chain.points) + if (ring.empty() || !geom::equals(p, ring.back())) ring.push_back(p); + + // From a touch exit, the entry at the same point is the far side of the + // touch: walk on to the next one + auto it = chain.touchExit + ? std::upper_bound(entries.begin(), entries.end(), std::make_pair(chain.exit, std::numeric_limits::max())) + : std::lower_bound(entries.begin(), entries.end(), std::make_pair(chain.exit, std::size_t(0))); + if (it == entries.end()) it = entries.begin(); + const std::size_t next = it->second; + if (chains[next].used && next != first) return false; + + double span = chains[next].entry - chain.exit; + if (span < 0 || (span == 0 && chain.touchExit)) span += rc.perimeter; + std::pair passed[4]; + int numPassed = 0; + for (int k = 0; k < 4; k++) { + double d = rc.cornerPos(k) - chain.exit; + if (d < 0) d += rc.perimeter; + if (d > 0 && d < span) passed[numPassed++] = { d, k }; + } + std::sort(passed, passed + numPassed); + for (int k = 0; k < numPassed; k++) ring.push_back(rc.corner(passed[k].second)); + + if (next == first) break; + c = next; + } + if (!geom::equals(ring.front(), ring.back())) ring.push_back(ring.front()); + if (ring.size() < 4 || geom::area(ring) <= 0) return false; + Polygon p; + p.outer() = std::move(ring); + out.push_back(std::move(p)); + } + } + + // Holes wholly inside the box go to whichever new outer contains them + for (const Ring *hole : insideHoles) { + bool placed = false; + for (std::size_t i = firstOuter; i < out.size() && !placed; i++) { + const int r = ringInRing(*hole, out[i].outer()); + if (r == 0) return false; + if (r > 0) { + out[i].inners().push_back(*hole); + placed = true; + } + } + if (!placed) return false; + } + return true; +} + +} // namespace + +bool rect_clip(const MultiPolygon &mp, Box const &bbox, MultiPolygon &out) { + out.clear(); + const RectClipper rc(bbox); + if (!(rc.w > 0 && rc.h > 0)) return false; + for (const Polygon &poly : mp) + if (!clipPolygon(rc, poly, out)) return false; + return true; +} diff --git a/src/osm_lua_processing.cpp b/src/osm_lua_processing.cpp index c719eaac..ba34304e 100644 --- a/src/osm_lua_processing.cpp +++ b/src/osm_lua_processing.cpp @@ -864,8 +864,13 @@ void OsmLuaProcessing::Layer(const string &layerName, bool area) { if (isRelation) { try { MultiPolygon mp = multiPolygonCached(); - if(CorrectGeometry(mp) == CorrectGeometryResult::Invalid) return; - NodeID id = osmMemTiles.storeMultiPolygon(mp); + auto correctionResult = CorrectGeometry(mp); + if(correctionResult == CorrectGeometryResult::Invalid) return; + // make_valid usually succeeds; checking that once here saves + // re-validating the relation in every tile it is clipped to. + const bool knownValid = correctionResult == CorrectGeometryResult::Valid || + PrecheckValid(mp) || geom::is_valid(mp); + NodeID id = osmMemTiles.storeMultiPolygon(mp, knownValid); OutputObject oo(geomType, layers.layerMap[layerName], id, 0, layerMinZoom); outputs.push_back(std::make_pair(std::move(oo), attributes)); } catch(std::out_of_range &err) { @@ -893,7 +898,7 @@ void OsmLuaProcessing::Layer(const string &layerName, bool area) { id = USE_WAY_STORE_REVERSED | originalOsmID; wayEmitted = true; } else - id = osmMemTiles.storeMultiPolygon(mp); + id = osmMemTiles.storeMultiPolygon(mp, correctionResult != CorrectGeometryResult::Corrected); OutputObject oo(geomType, layers.layerMap[layerName], id, 0, layerMinZoom); outputs.push_back(std::make_pair(std::move(oo), attributes)); } diff --git a/src/tile_data.cpp b/src/tile_data.cpp index 3fff009a..89b367e0 100644 --- a/src/tile_data.cpp +++ b/src/tile_data.cpp @@ -34,8 +34,8 @@ static const char* validityFailureName(geom::validity_failure_type failure) { // Thin wrapper around the shared repair_multi_polygon() that adds per-object // verbose diagnostics. See geom.cpp for the dissolve + zero-width buffer logic. -static bool repairMultiPolygon(MultiPolygon &mp, NodeID objectID) { - bool ok = repair_multi_polygon(mp); +static bool repairMultiPolygon(MultiPolygon &mp, NodeID objectID, bool *wholeValid = nullptr) { + bool ok = repair_multi_polygon(mp, false, wholeValid); if (!ok && verbose) std::cerr << ("multipolygon repair failed for object " + std::to_string(objectID) + "\n"); return ok; @@ -60,6 +60,7 @@ TileDataSource::TileDataSource(size_t threadNum, unsigned int indexZoom, bool in linestringStores(threadNum), multilinestringStores(threadNum), multipolygonStores(threadNum), + multipolygonValid(threadNum), multiPolygonClipCache(ClipCache(threadNum, indexZoom)), multiLinestringClipCache(ClipCache(threadNum, indexZoom)) { @@ -398,6 +399,20 @@ Geometry TileDataSource::buildWayGeometry(OutputGeometryType const geomType, geom::correct(mp); geom::validity_failure_type failure = geom::validity_failure_type::no_failure; bool valid = multiPolygonDefinitelyValid(mp); + if (!valid && inputKnownValid) { + // fast_clip leaves zero-width bridges along the box edge where a concave + // polygon leaves and re-enters the box. Clip the valid input exactly + // instead, which is much cheaper than the is_valid/intersection/repair + // path below; rect_clip bails out on anything it can't handle. + MultiPolygon original; + if (cachedClip == nullptr) populateMultiPolygon(original, objectID); + MultiPolygon clipped; + if (rect_clip(cachedClip == nullptr ? original : input, box, clipped) && + (multiPolygonDefinitelyValid(clipped) || geom::is_valid(clipped))) { + multiPolygonClipCache.add(bbox, objectID, clipped, true); + return clipped; + } + } if (!valid) valid = geom::is_valid(mp,failure); if (!valid) { @@ -432,14 +447,16 @@ Geometry TileDataSource::buildWayGeometry(OutputGeometryType const geomType, // The intersection result can itself still be invalid for very complex // multipolygons (e.g. large reservoirs), which previously produced // dropped or holey tiles. Repair it before returning. - repairMultiPolygon(output, objectID); - multiPolygonClipCache.add(bbox, objectID, output); + bool outputValid; + repairMultiPolygon(output, objectID, &outputValid); + // Children of this tile can then take the fast paths for valid input + multiPolygonClipCache.add(bbox, objectID, output, outputValid); return output; } else if (!valid) { // occasionally also wrong_topological_dimension, disconnected_interior: // defects geom::correct cannot mend. Repair mp in place; on failure it // is left unchanged so behaviour never regresses. - repairMultiPolygon(mp, objectID); + repairMultiPolygon(mp, objectID, &valid); } } @@ -658,7 +675,7 @@ NodeID TileDataSource::storeLinestring(const Linestring& src) { return rv; } -NodeID TileDataSource::storeMultiPolygon(const MultiPolygon& src) { +NodeID TileDataSource::storeMultiPolygon(const MultiPolygon& src, bool knownValid) { const auto& store = multipolygonStore.get(this); multi_polygon_t dst; @@ -676,6 +693,8 @@ NodeID TileDataSource::storeMultiPolygon(const MultiPolygon& src) { NodeID offset = store.second->size(); store.second->emplace_back(std::move(dst)); + // The store is leased to this thread, so its validity flags are too + multipolygonValid[store.first].push_back(knownValid); NodeID rv = (store.first << (TILE_DATA_ID_SIZE - shardBits)) + offset; return rv; } diff --git a/test/geom.test.cpp b/test/geom.test.cpp new file mode 100644 index 00000000..8464f04a --- /dev/null +++ b/test/geom.test.cpp @@ -0,0 +1,84 @@ +#include +#include +#include "external/minunit.h" +#include "geom.h" + +static MultiPolygon fromWkt(const char *wkt) { + MultiPolygon mp; + geom::read_wkt(wkt, mp); + geom::correct(mp); + return mp; +} + +// rect_clip must give a valid result of the expected area +static void checkClip(const char *wkt, const Box &box, std::size_t expectedPolygons, double expectedArea) { + MultiPolygon mp = fromWkt(wkt); + mu_check(geom::is_valid(mp)); + + MultiPolygon clipped; + mu_check(rect_clip(mp, box, clipped)); + mu_check(geom::is_valid(clipped)); + mu_check(clipped.size() == expectedPolygons); + mu_check(std::abs(geom::area(clipped) - expectedArea) < 1e-9); +} + +MU_TEST(test_rect_clip_concave) { + // A U shape whose arms leave the box through its top: Sutherland-Hodgman + // joins the arms with a zero-width bridge along the box edge + const char *u = "MULTIPOLYGON(((0 0,0 10,2 10,2 2,8 2,8 10,10 10,10 0,0 0)))"; + checkClip(u, Box(Point(-1, -1), Point(11, 5)), 1, 32); + // Cutting the base off leaves the two arms as separate polygons + checkClip(u, Box(Point(-1, 4), Point(11, 12)), 2, 24); +} + +MU_TEST(test_rect_clip_holes) { + const char *holes = "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0),(1 1,3 1,3 3,1 3,1 1),(6 4,9 4,9 6,6 6,6 4)))"; + // One hole inside the box, one cut by its right edge + checkClip(holes, Box(Point(-1, -1), Point(7, 7)), 1, 43); + // Box wholly inside the polygon: the box itself + checkClip(holes, Box(Point(4, 7), Point(5, 8)), 1, 1); + // Box wholly inside a hole: nothing + checkClip(holes, Box(Point(1.5, 1.5), Point(2.5, 2.5)), 0, 0); +} + +MU_TEST(test_rect_clip_trivial) { + const char *square = "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))"; + // Disjoint + checkClip(square, Box(Point(20, 20), Point(30, 30)), 0, 0); + // Polygon wholly inside the box is returned as it is + MultiPolygon mp = fromWkt(square), clipped; + mu_check(rect_clip(mp, Box(Point(-1, -1), Point(11, 11)), clipped)); + mu_check(geom::equals(mp, clipped)); +} + +MU_TEST(test_rect_clip_vertex_on_edge) { + // The ring crosses the box edge at a (duplicated) vertex exactly on it. + // (Boost's intersection returns nothing for the first box.) + const char *poly = "MULTIPOLYGON(((0 0,0 4,0 4,0 10,10 10,10 6,12 4,10 2,10 0,0 0)))"; + checkClip(poly, Box(Point(-1, 4), Point(11, 11)), 1, 61.5); + checkClip(poly, Box(Point(-1, -1), Point(11, 11)), 1, 103); +} + +MU_TEST(test_rect_clip_touch_from_inside) { + // A notch from the top comes down to touch the box's bottom edge at (5 0), then + // goes back up. The box walk must not run through the touch: two polygons. + const char *notch = "MULTIPOLYGON(((-5 -5,-5 8,4 8,5 0,6 8,15 8,15 -5,-5 -5)))"; + checkClip(notch, Box(Point(0, 0), Point(10, 10)), 2, 72); + // The same touching the left edge + const char *sideNotch = "MULTIPOLYGON(((-5 -5,-5 15,8 15,8 6,0 5,8 4,8 -5,-5 -5)))"; + checkClip(sideNotch, Box(Point(0, 0), Point(10, 10)), 2, 72); +} + +MU_TEST_SUITE(test_suite_geom) { + MU_RUN_TEST(test_rect_clip_concave); + MU_RUN_TEST(test_rect_clip_holes); + MU_RUN_TEST(test_rect_clip_trivial); + MU_RUN_TEST(test_rect_clip_vertex_on_edge); + MU_RUN_TEST(test_rect_clip_touch_from_inside); +} + +int main() { + MU_RUN_SUITE(test_suite_geom); + MU_REPORT(); + return MU_EXIT_CODE; +}