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