GEOS 3.8.1
BinaryOp.h
1/**********************************************************************
2 *
3 * GEOS - Geometry Engine Open Source
4 * http://geos.osgeo.org
5 *
6 * Copyright (C) 2013 Sandro Santilli <strk@kbt.io>
7 * Copyright (C) 2006 Refractions Research Inc.
8 *
9 * This is free software; you can redistribute and/or modify it under
10 * the terms of the GNU Lesser General Public Licence as published
11 * by the Free Software Foundation.
12 * See the COPYING file for more information.
13 *
14 **********************************************************************
15 *
16 * Last port: ORIGINAL WORK
17 *
18 **********************************************************************
19 *
20 * This file provides a single templated function, taking two
21 * const Geometry pointers, applying a binary operator to them
22 * and returning a result Geometry in an unique_ptr<>.
23 *
24 * The binary operator is expected to take two const Geometry pointers
25 * and return a newly allocated Geometry pointer, possibly throwing
26 * a TopologyException to signal it couldn't succeed due to robustness
27 * issues.
28 *
29 * This function will catch TopologyExceptions and try again with
30 * slightly modified versions of the input. The following heuristic
31 * is used:
32 *
33 * - Try with original input.
34 * - Try removing common bits from input coordinate values
35 * - Try snaping input geometries to each other
36 * - Try snaping input coordinates to a increasing grid (size from 1/25 to 1)
37 * - Try simplifiying input with increasing tolerance (from 0.01 to 0.04)
38 *
39 * If none of the step succeeds the original exception is thrown.
40 *
41 * Note that you can skip Grid snapping, Geometry snapping and Simplify policies
42 * by a compile-time define when building geos.
43 * See USE_TP_SIMPLIFY_POLICY, USE_PRECISION_REDUCTION_POLICY and
44 * USE_SNAPPING_POLICY macros below.
45 *
46 *
47 **********************************************************************/
48
49#ifndef GEOS_GEOM_BINARYOP_H
50#define GEOS_GEOM_BINARYOP_H
51
52#include <geos/algorithm/BoundaryNodeRule.h>
53#include <geos/geom/Geometry.h>
54#include <geos/geom/GeometryCollection.h>
55#include <geos/geom/Polygon.h>
56#include <geos/geom/PrecisionModel.h>
57#include <geos/geom/GeometryFactory.h>
58#include <geos/precision/CommonBitsRemover.h>
59#include <geos/precision/SimpleGeometryPrecisionReducer.h>
60#include <geos/precision/GeometryPrecisionReducer.h>
61
62#include <geos/operation/overlay/snap/GeometrySnapper.h>
63
64#include <geos/simplify/TopologyPreservingSimplifier.h>
65#include <geos/operation/IsSimpleOp.h>
66#include <geos/operation/valid/IsValidOp.h>
67#include <geos/operation/valid/TopologyValidationError.h>
68#include <geos/util/TopologyException.h>
69#include <geos/util.h>
70
71#include <memory> // for unique_ptr
72
73//#define GEOS_DEBUG_BINARYOP 1
74#define GEOS_DEBUG_BINARYOP_PRINT_INVALID 1
75
76#ifdef GEOS_DEBUG_BINARYOP
77# include <iostream>
78# include <iomanip>
79# include <sstream>
80#endif
81
82
83/*
84 * Always try original input first
85 */
86#ifndef USE_ORIGINAL_INPUT
87# define USE_ORIGINAL_INPUT 1
88#endif
89
90/*
91 * Check validity of operation between original geometries
92 */
93#define GEOS_CHECK_ORIGINAL_RESULT_VALIDITY 0
94
95
96/*
97 * Define this to use PrecisionReduction policy
98 * in an attempt at by-passing binary operation
99 * robustness problems (handles TopologyExceptions)
100 */
101#ifndef USE_PRECISION_REDUCTION_POLICY
102# define USE_PRECISION_REDUCTION_POLICY 1
103#endif
104
105/*
106 * Check validity of operation performed
107 * by precision reduction policy.
108 *
109 * Precision reduction policy reduces precision of inputs
110 * and restores it in the result. The restore phase may
111 * introduce invalidities.
112 *
113 */
114#define GEOS_CHECK_PRECISION_REDUCTION_VALIDITY 0
115
116/*
117 * Define this to use TopologyPreserving simplification policy
118 * in an attempt at by-passing binary operation
119 * robustness problems (handles TopologyExceptions)
120 */
121#ifndef USE_TP_SIMPLIFY_POLICY
122//# define USE_TP_SIMPLIFY_POLICY 1
123#endif
124
125/*
126 * Use common bits removal policy.
127 * If enabled, this would be tried /before/
128 * Geometry snapping.
129 */
130#ifndef USE_COMMONBITS_POLICY
131# define USE_COMMONBITS_POLICY 1
132#endif
133
134/*
135 * Check validity of operation performed
136 * by common bits removal policy.
137 *
138 * This matches what EnhancedPrecisionOp does in JTS
139 * and fixes 5 tests of invalid outputs in our testsuite
140 * (stmlf-cases-20061020-invalid-output.xml)
141 * and breaks 1 test (robustness-invalid-output.xml) so much
142 * to prevent a result.
143 *
144 */
145#define GEOS_CHECK_COMMONBITS_VALIDITY 1
146
147/*
148 * Use snapping policy
149 */
150#ifndef USE_SNAPPING_POLICY
151# define USE_SNAPPING_POLICY 1
152#endif
153
154/* Remove common bits before snapping */
155#ifndef CBR_BEFORE_SNAPPING
156# define CBR_BEFORE_SNAPPING 1
157#endif
158
159
160/*
161 * Check validity of result from SnapOp
162 */
163#define GEOS_CHECK_SNAPPINGOP_VALIDITY 0
164
165
166namespace geos {
167namespace geom { // geos::geom
168
169inline bool
170check_valid(const Geometry& g, const std::string& label, bool doThrow = false, bool validOnly = false)
171{
172 if(g.isLineal()) {
173 if(! validOnly) {
174 operation::IsSimpleOp sop(g, algorithm::BoundaryNodeRule::getBoundaryEndPoint());
175 if(! sop.isSimple()) {
176 if(doThrow) {
177 throw geos::util::TopologyException(
178 label + " is not simple");
179 }
180 return false;
181 }
182 }
183 }
184 else {
185 operation::valid::IsValidOp ivo(&g);
186 if(! ivo.isValid()) {
187 using operation::valid::TopologyValidationError;
188 TopologyValidationError* err = ivo.getValidationError();
189#ifdef GEOS_DEBUG_BINARYOP
190 std::cerr << label << " is INVALID: "
191 << err->toString()
192 << " (" << std::setprecision(20)
193 << err->getCoordinate() << ")"
194 << std::endl
195#ifdef GEOS_DEBUG_BINARYOP_PRINT_INVALID
196 << "<A>" << std::endl
197 << g.toString()
198 << std::endl
199 << "</A>" << std::endl
200#endif
201 ;
202#endif
203 if(doThrow) {
204 throw geos::util::TopologyException(
205 label + " is invalid: " + err->toString(),
206 err->getCoordinate());
207 }
208 return false;
209 }
210 }
211 return true;
212}
213
214/*
215 * Attempt to fix noding of multilines and
216 * self-intersection of multipolygons
217 *
218 * May return the input untouched.
219 */
220inline std::unique_ptr<Geometry>
221fix_self_intersections(std::unique_ptr<Geometry> g, const std::string& label)
222{
223 ::geos::ignore_unused_variable_warning(label);
224#ifdef GEOS_DEBUG_BINARYOP
225 std::cerr << label << " fix_self_intersection (UnaryUnion)" << std::endl;
226#endif
227
228 // Only multi-components can be fixed by UnaryUnion
229 if(! dynamic_cast<const GeometryCollection*>(g.get())) {
230 return g;
231 }
232
233 using operation::valid::IsValidOp;
234
235 IsValidOp ivo(g.get());
236
237 // Polygon is valid, nothing to do
238 if(ivo.isValid()) {
239 return g;
240 }
241
242 // Not all invalidities can be fixed by this code
243
244 using operation::valid::TopologyValidationError;
245 TopologyValidationError* err = ivo.getValidationError();
246 switch(err->getErrorType()) {
247 case TopologyValidationError::eRingSelfIntersection:
248 case TopologyValidationError::eTooFewPoints: // collapsed lines
249#ifdef GEOS_DEBUG_BINARYOP
250 std::cerr << label << " ATTEMPT_TO_FIX: " << err->getErrorType() << ": " << *g << std::endl;
251#endif
252 g = g->Union();
253#ifdef GEOS_DEBUG_BINARYOP
254 std::cerr << label << " ATTEMPT_TO_FIX succeeded.. " << std::endl;
255#endif
256 return g;
257 case TopologyValidationError::eSelfIntersection:
258 // this one is within a single component, won't be fixed
259 default:
260#ifdef GEOS_DEBUG_BINARYOP
261 std::cerr << label << " invalidity is: " << err->getErrorType() << std::endl;
262#endif
263 return g;
264 }
265}
266
267
273template <class BinOp>
274std::unique_ptr<Geometry>
275SnapOp(const Geometry* g0, const Geometry* g1, BinOp _Op)
276{
277 typedef std::unique_ptr<Geometry> GeomPtr;
278
279 //using geos::precision::GeometrySnapper;
280 using geos::operation::overlay::snap::GeometrySnapper;
281
282 // Snap tolerance must be computed on the original
283 // (not commonbits-removed) geoms
284 double snapTolerance = GeometrySnapper::computeOverlaySnapTolerance(*g0, *g1);
285#if GEOS_DEBUG_BINARYOP
286 std::cerr << std::setprecision(20) << "Computed snap tolerance: " << snapTolerance << std::endl;
287#endif
288
289
290#if CBR_BEFORE_SNAPPING
291 // Compute common bits
292 geos::precision::CommonBitsRemover cbr;
293 cbr.add(g0);
294 cbr.add(g1);
295#if GEOS_DEBUG_BINARYOP
296 std::cerr << "Computed common bits: " << cbr.getCommonCoordinate() << std::endl;
297#endif
298
299 // Now remove common bits
300 GeomPtr rG0 = g0->clone();
301 cbr.removeCommonBits(rG0.get());
302 GeomPtr rG1 = g1->clone();
303 cbr.removeCommonBits(rG1.get());
304
305#if GEOS_DEBUG_BINARYOP
306 check_valid(*rG0, "CBR: removed-bits geom 0");
307 check_valid(*rG1, "CBR: removed-bits geom 1");
308#endif
309
310 const Geometry& operand0 = *rG0;
311 const Geometry& operand1 = *rG1;
312#else // don't CBR before snapping
313 const Geometry& operand0 = *g0;
314 const Geometry& operand1 = *g1;
315#endif
316
317
318 GeometrySnapper snapper0(operand0);
319 GeomPtr snapG0(snapper0.snapTo(operand1, snapTolerance));
320 //snapG0 = fix_self_intersections(snapG0, "SNAP: snapped geom 0");
321
322 // NOTE: second geom is snapped on the snapped first one
323 GeometrySnapper snapper1(operand1);
324 GeomPtr snapG1(snapper1.snapTo(*snapG0, snapTolerance));
325 //snapG1 = fix_self_intersections(snapG1, "SNAP: snapped geom 1");
326
327 // Run the binary op
328 GeomPtr result(_Op(snapG0.get(), snapG1.get()));
329
330#if GEOS_DEBUG_BINARYOP
331 check_valid(*result, "SNAP: result (before common-bits addition");
332#endif
333
334#if CBR_BEFORE_SNAPPING
335 // Add common bits back in
336 cbr.addCommonBits(result.get());
337 //result = fix_self_intersections(result, "SNAP: result (after common-bits addition)");
338
339 check_valid(*result, "CBR: result (after common-bits addition)", true);
340
341#endif
342
343 return result;
344}
345
346template <class BinOp>
347std::unique_ptr<Geometry>
348BinaryOp(const Geometry* g0, const Geometry* g1, BinOp _Op)
349{
350 typedef std::unique_ptr<Geometry> GeomPtr;
351
352 GeomPtr ret;
353 geos::util::TopologyException origException;
354
355#ifdef USE_ORIGINAL_INPUT
356 // Try with original input
357 try {
358#if GEOS_DEBUG_BINARYOP
359 std::cerr << "Trying with original input." << std::endl;
360#endif
361 ret.reset(_Op(g0, g1));
362
363#if GEOS_CHECK_ORIGINAL_RESULT_VALIDITY
364 check_valid(*ret, "Overlay result between original inputs", true, true);
365#endif
366
367#if GEOS_DEBUG_BINARYOP
368 std::cerr << "Attempt with original input succeeded" << std::endl;
369#endif
370 return ret;
371 }
372 catch(const geos::util::TopologyException& ex) {
373 origException = ex;
374#if GEOS_DEBUG_BINARYOP
375 std::cerr << "Original exception: " << ex.what() << std::endl;
376#endif
377 }
378#endif // USE_ORIGINAL_INPUT
379
380 check_valid(*g0, "Input geom 0", true, true);
381 check_valid(*g1, "Input geom 1", true, true);
382
383#if USE_COMMONBITS_POLICY
384 // Try removing common bits (possibly obsoleted by snapping below)
385 //
386 // NOTE: this policy was _later_ implemented
387 // in JTS as EnhancedPrecisionOp
388 // TODO: consider using the now-ported EnhancedPrecisionOp
389 // here too
390 //
391 try {
392 GeomPtr rG0;
393 GeomPtr rG1;
394 precision::CommonBitsRemover cbr;
395
396#if GEOS_DEBUG_BINARYOP
397 std::cerr << "Trying with Common Bits Remover (CBR)" << std::endl;
398#endif
399
400 cbr.add(g0);
401 cbr.add(g1);
402
403 rG0 = g0->clone();
404 cbr.removeCommonBits(rG0.get());
405
406 rG1 = g1->clone();
407 cbr.removeCommonBits(rG1.get());
408
409#if GEOS_DEBUG_BINARYOP
410 check_valid(*rG0, "CBR: geom 0 (after common-bits removal)");
411 check_valid(*rG1, "CBR: geom 1 (after common-bits removal)");
412#endif
413
414 ret.reset(_Op(rG0.get(), rG1.get()));
415
416#if GEOS_DEBUG_BINARYOP
417 check_valid(*ret, "CBR: result (before common-bits addition)");
418#endif
419
420 cbr.addCommonBits(ret.get());
421
422#if GEOS_CHECK_COMMONBITS_VALIDITY
423 check_valid(*ret, "CBR: result (after common-bits addition)", true);
424#endif
425
426#if GEOS_DEBUG_BINARYOP
427 std::cerr << "Attempt with CBR succeeded" << std::endl;
428#endif
429
430 return ret;
431 }
432 catch(const geos::util::TopologyException& ex) {
433 ::geos::ignore_unused_variable_warning(ex);
434#if GEOS_DEBUG_BINARYOP
435 std::cerr << "CBR: " << ex.what() << std::endl;
436#endif
437 }
438#endif
439
440 // Try with snapping
441 //
442 // TODO: possible optimization would be reusing the
443 // already common-bit-removed inputs and just
444 // apply geometry snapping, whereas the current
445 // SnapOp function does both.
446// {
447#if USE_SNAPPING_POLICY
448
449#if GEOS_DEBUG_BINARYOP
450 std::cerr << "Trying with snapping " << std::endl;
451#endif
452
453 try {
454 ret = SnapOp(g0, g1, _Op);
455#if GEOS_CHECK_SNAPPINGOP_VALIDITY
456 check_valid(*ret, "SNAP: result", true, true);
457#endif
458#if GEOS_DEBUG_BINARYOP
459 std::cerr << "SnapOp succeeded" << std::endl;
460#endif
461 return ret;
462
463 }
464 catch(const geos::util::TopologyException& ex) {
465 ::geos::ignore_unused_variable_warning(ex);
466#if GEOS_DEBUG_BINARYOP
467 std::cerr << "SNAP: " << ex.what() << std::endl;
468#endif
469 }
470
471#endif // USE_SNAPPING_POLICY }
472
473// {
474#if USE_PRECISION_REDUCTION_POLICY
475
476
477 // Try reducing precision
478 try {
479 long unsigned int g0scale =
480 static_cast<long unsigned int>(g0->getFactory()->getPrecisionModel()->getScale());
481 long unsigned int g1scale =
482 static_cast<long unsigned int>(g1->getFactory()->getPrecisionModel()->getScale());
483
484#if GEOS_DEBUG_BINARYOP
485 std::cerr << "Original input scales are: "
486 << g0scale
487 << " and "
488 << g1scale
489 << std::endl;
490#endif
491
492 double maxScale = 1e16;
493
494 // Don't use a scale bigger than the input one
495 if(g0scale && static_cast<double>(g0scale) < maxScale) {
496 maxScale = static_cast<double>(g0scale);
497 }
498 if(g1scale && static_cast<double>(g1scale) < maxScale) {
499 maxScale = static_cast<double>(g1scale);
500 }
501
502
503 for(double scale = maxScale; scale >= 1; scale /= 10) {
504 PrecisionModel pm(scale);
505 GeometryFactory::Ptr gf = GeometryFactory::create(&pm);
506#if GEOS_DEBUG_BINARYOP
507 std::cerr << "Trying with scale " << scale << std::endl;
508#endif
509
510 precision::GeometryPrecisionReducer reducer(*gf);
511 GeomPtr rG0(reducer.reduce(*g0));
512 GeomPtr rG1(reducer.reduce(*g1));
513
514#if GEOS_DEBUG_BINARYOP
515 check_valid(*rG0, "PR: geom 0 (after precision reduction)");
516 check_valid(*rG1, "PR: geom 1 (after precision reduction)");
517#endif
518
519 try {
520 ret.reset(_Op(rG0.get(), rG1.get()));
521 // restore original precision (least precision between inputs)
522 if(g0->getFactory()->getPrecisionModel()->compareTo(g1->getFactory()->getPrecisionModel()) < 0) {
523 ret.reset(g0->getFactory()->createGeometry(ret.get()));
524 }
525 else {
526 ret.reset(g1->getFactory()->createGeometry(ret.get()));
527 }
528
529#if GEOS_CHECK_PRECISION_REDUCTION_VALIDITY
530 check_valid(*ret, "PR: result (after restore of original precision)", true);
531#endif
532
533#if GEOS_DEBUG_BINARYOP
534 std::cerr << "Attempt with scale " << scale << " succeeded" << std::endl;
535#endif
536 return ret;
537 }
538 catch(const geos::util::TopologyException& ex) {
539#if GEOS_DEBUG_BINARYOP
540 std::cerr << "Reduced with scale (" << scale << "): "
541 << ex.what() << std::endl;
542#endif
543 if(scale == 1) {
544 throw ex;
545 }
546 }
547
548 }
549
550 }
551 catch(const geos::util::TopologyException& ex) {
552#if GEOS_DEBUG_BINARYOP
553 std::cerr << "Reduced: " << ex.what() << std::endl;
554#endif
555 ::geos::ignore_unused_variable_warning(ex);
556 }
557
558#endif
559// USE_PRECISION_REDUCTION_POLICY }
560
561
562
563
564
565// {
566#if USE_TP_SIMPLIFY_POLICY
567
568 // Try simplifying
569 try {
570
571 double maxTolerance = 0.04;
572 double minTolerance = 0.01;
573 double tolStep = 0.01;
574
575 for(double tol = minTolerance; tol <= maxTolerance; tol += tolStep) {
576#if GEOS_DEBUG_BINARYOP
577 std::cerr << "Trying simplifying with tolerance " << tol << std::endl;
578#endif
579
580 GeomPtr rG0(simplify::TopologyPreservingSimplifier::simplify(g0, tol));
581 GeomPtr rG1(simplify::TopologyPreservingSimplifier::simplify(g1, tol));
582
583 try {
584 ret.reset(_Op(rG0.get(), rG1.get()));
585 return ret;
586 }
587 catch(const geos::util::TopologyException& ex) {
588 if(tol >= maxTolerance) {
589 throw ex;
590 }
591#if GEOS_DEBUG_BINARYOP
592 std::cerr << "Simplified with tolerance (" << tol << "): "
593 << ex.what() << std::endl;
594#endif
595 }
596
597 }
598
599 return ret;
600
601 }
602 catch(const geos::util::TopologyException& ex) {
603#if GEOS_DEBUG_BINARYOP
604 std::cerr << "Simplified: " << ex.what() << std::endl;
605#endif
606 }
607
608#endif
609// USE_TP_SIMPLIFY_POLICY }
610
611#if GEOS_DEBUG_BINARYOP
612 std::cerr << "No attempts worked to union " << std::endl;
613 std::cerr << "Input geometries:" << std::endl
614 << "<A>" << std::endl
615 << g0->toString() << std::endl
616 << "</A>" << std::endl
617 << "<B>" << std::endl
618 << g1->toString() << std::endl
619 << "</B>" << std::endl;
620#endif
621
622 throw origException;
623}
624
625
626} // namespace geos::geom
627} // namespace geos
628
629#endif // GEOS_GEOM_BINARYOP_H
geos::algorithm::BoundaryNodeRule::getBoundaryEndPoint
static const BoundaryNodeRule & getBoundaryEndPoint()
The Endpoint Boundary Node Rule.
geos::geom::GeometryFactory::create
static GeometryFactory::Ptr create()
Constructs a GeometryFactory that generates Geometries having a floating PrecisionModel and a spatial...
geos::geom::Geometry
Basic implementation of Geometry, constructed and destructed by GeometryFactory.
Definition: Geometry.h:188
geos::geom::Geometry::clone
virtual std::unique_ptr< Geometry > clone() const =0
Make a deep-copy of this Geometry.
geos::operation::overlay::snap::GeometrySnapper
Snaps the vertices and segments of a geom::Geometry to another Geometry's vertices.
Definition: GeometrySnapper.h:58
geos::precision::CommonBitsRemover
Allow computing and removing common mantissa bits from one or more Geometries.
Definition: CommonBitsRemover.h:40
geos::precision::CommonBitsRemover::getCommonCoordinate
geom::Coordinate & getCommonCoordinate()
geos::precision::CommonBitsRemover::addCommonBits
geom::Geometry * addCommonBits(geom::Geometry *geom)
Adds the common coordinate bits back into a Geometry. The coordinates of the Geometry are changed.
geos::precision::CommonBitsRemover::removeCommonBits
void removeCommonBits(geom::Geometry *geom)
Removes the common coordinate bits from a Geometry. The coordinates of the Geometry are changed.
geos::precision::CommonBitsRemover::add
void add(const geom::Geometry *geom)
geos::util::TopologyException
Indicates an invalid or inconsistent topological situation encountered during processing.
Definition: TopologyException.h:35
geos::geom::SnapOp
std::unique_ptr< Geometry > SnapOp(const Geometry *g0, const Geometry *g1, BinOp _Op)
Apply a binary operation to the given geometries after snapping them to each other after common-bits ...
Definition: BinaryOp.h:275
geos
Basic namespace for all GEOS functionalities.
Definition: IndexedNestedRingTester.h:25
Generated by doxygen 1.9.5