GEOS  3.8.0dev
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 
166 namespace geos {
167 namespace geom { // geos::geom
168 
169 inline bool
170 check_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) {
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) {
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  */
220 inline std::unique_ptr<Geometry>
221 fix_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 
273 template <class BinOp>
274 std::unique_ptr<Geometry>
275 SnapOp(const Geometry* g0, const Geometry* g1, BinOp _Op)
276 {
277  typedef std::unique_ptr<Geometry> GeomPtr;
278 
279  //using geos::precision::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
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 
346 template <class BinOp>
347 std::unique_ptr<Geometry>
348 BinaryOp(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
Allow computing and removing common mantissa bits from one or more Geometries.
Definition: CommonBitsRemover.h:40
void removeCommonBits(geom::Geometry *geom)
Removes the common coordinate bits from a Geometry. The coordinates of the Geometry are changed...
Basic implementation of Geometry, constructed and destructed by GeometryFactory.
Definition: Geometry.h:188
geom::Coordinate & getCommonCoordinate()
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
Basic namespace for all GEOS functionalities.
Definition: IndexedNestedRingTester.h:25
Snaps the vertices and segments of a geom::Geometry to another Geometry's vertices.
Definition: GeometrySnapper.h:58
virtual std::unique_ptr< Geometry > clone() const =0
Make a deep-copy of this Geometry.
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...