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/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  if(doThrow) {
179  label + " is not simple");
180  }
181  return false;
182  }
183  }
184  }
185  else {
186  operation::valid::IsValidOp ivo(&g);
187  if(! ivo.isValid()) {
188  using operation::valid::TopologyValidationError;
189  TopologyValidationError* err = ivo.getValidationError();
190 #ifdef GEOS_DEBUG_BINARYOP
191  std::cerr << label << " is INVALID: "
192  << err->toString()
193  << " (" << std::setprecision(20)
194  << err->getCoordinate() << ")"
195  << std::endl
196 #ifdef GEOS_DEBUG_BINARYOP_PRINT_INVALID
197  << "<A>" << std::endl
198  << g.toString()
199  << std::endl
200  << "</A>" << std::endl
201 #endif
202  ;
203 #endif
204  if(doThrow) {
206  label + " is invalid: " + err->toString(),
207  err->getCoordinate());
208  }
209  return false;
210  }
211  }
212  return true;
213 }
214 
215 /*
216  * Attempt to fix noding of multilines and
217  * self-intersection of multipolygons
218  *
219  * May return the input untouched.
220  */
221 inline std::unique_ptr<Geometry>
222 fix_self_intersections(std::unique_ptr<Geometry> g, const std::string& label)
223 {
224  ::geos::ignore_unused_variable_warning(label);
225 #ifdef GEOS_DEBUG_BINARYOP
226  std::cerr << label << " fix_self_intersection (UnaryUnion)" << std::endl;
227 #endif
228 
229  // Only multi-components can be fixed by UnaryUnion
230  if(! dynamic_cast<const GeometryCollection*>(g.get())) {
231  return g;
232  }
233 
234  using operation::valid::IsValidOp;
235 
236  IsValidOp ivo(g.get());
237 
238  // Polygon is valid, nothing to do
239  if(ivo.isValid()) {
240  return g;
241  }
242 
243  // Not all invalidities can be fixed by this code
244 
245  using operation::valid::TopologyValidationError;
246  TopologyValidationError* err = ivo.getValidationError();
247  switch(err->getErrorType()) {
248  case TopologyValidationError::eRingSelfIntersection:
249  case TopologyValidationError::eTooFewPoints: // collapsed lines
250 #ifdef GEOS_DEBUG_BINARYOP
251  std::cerr << label << " ATTEMPT_TO_FIX: " << err->getErrorType() << ": " << *g << std::endl;
252 #endif
253  g = g->Union();
254 #ifdef GEOS_DEBUG_BINARYOP
255  std::cerr << label << " ATTEMPT_TO_FIX succeeded.. " << std::endl;
256 #endif
257  return g;
258  case TopologyValidationError::eSelfIntersection:
259  // this one is within a single component, won't be fixed
260  default:
261 #ifdef GEOS_DEBUG_BINARYOP
262  std::cerr << label << " invalidity is: " << err->getErrorType() << std::endl;
263 #endif
264  return g;
265  }
266 }
267 
268 
274 template <class BinOp>
275 std::unique_ptr<Geometry>
276 SnapOp(const Geometry* g0, const Geometry* g1, BinOp _Op)
277 {
278  typedef std::unique_ptr<Geometry> GeomPtr;
279 
280  //using geos::precision::GeometrySnapper;
282 
283  // Snap tolerance must be computed on the original
284  // (not commonbits-removed) geoms
285  double snapTolerance = GeometrySnapper::computeOverlaySnapTolerance(*g0, *g1);
286 #if GEOS_DEBUG_BINARYOP
287  std::cerr << std::setprecision(20) << "Computed snap tolerance: " << snapTolerance << std::endl;
288 #endif
289 
290 
291 #if CBR_BEFORE_SNAPPING
292  // Compute common bits
294  cbr.add(g0);
295  cbr.add(g1);
296 #if GEOS_DEBUG_BINARYOP
297  std::cerr << "Computed common bits: " << cbr.getCommonCoordinate() << std::endl;
298 #endif
299 
300  // Now remove common bits
301  GeomPtr rG0(cbr.removeCommonBits(g0->clone()));
302  GeomPtr rG1(cbr.removeCommonBits(g1->clone()));
303 
304 #if GEOS_DEBUG_BINARYOP
305  check_valid(*rG0, "CBR: removed-bits geom 0");
306  check_valid(*rG1, "CBR: removed-bits geom 1");
307 #endif
308 
309  const Geometry& operand0 = *rG0;
310  const Geometry& operand1 = *rG1;
311 #else // don't CBR before snapping
312  const Geometry& operand0 = *g0;
313  const Geometry& operand1 = *g1;
314 #endif
315 
316 
317  GeometrySnapper snapper0(operand0);
318  GeomPtr snapG0(snapper0.snapTo(operand1, snapTolerance));
319  //snapG0 = fix_self_intersections(snapG0, "SNAP: snapped geom 0");
320 
321  // NOTE: second geom is snapped on the snapped first one
322  GeometrySnapper snapper1(operand1);
323  GeomPtr snapG1(snapper1.snapTo(*snapG0, snapTolerance));
324  //snapG1 = fix_self_intersections(snapG1, "SNAP: snapped geom 1");
325 
326  // Run the binary op
327  GeomPtr result(_Op(snapG0.get(), snapG1.get()));
328 
329 #if GEOS_DEBUG_BINARYOP
330  check_valid(*result, "SNAP: result (before common-bits addition");
331 #endif
332 
333 #if CBR_BEFORE_SNAPPING
334  // Add common bits back in
335  cbr.addCommonBits(result.get());
336  //result = fix_self_intersections(result, "SNAP: result (after common-bits addition)");
337 
338  check_valid(*result, "CBR: result (after common-bits addition)", true);
339 
340 #endif
341 
342  return result;
343 }
344 
345 template <class BinOp>
346 std::unique_ptr<Geometry>
347 BinaryOp(const Geometry* g0, const Geometry* g1, BinOp _Op)
348 {
349  typedef std::unique_ptr<Geometry> GeomPtr;
350 
351  GeomPtr ret;
352  geos::util::TopologyException origException;
353 
354 #ifdef USE_ORIGINAL_INPUT
355  // Try with original input
356  try {
357 #if GEOS_DEBUG_BINARYOP
358  std::cerr << "Trying with original input." << std::endl;
359 #endif
360  ret.reset(_Op(g0, g1));
361 
362 #if GEOS_CHECK_ORIGINAL_RESULT_VALIDITY
363  check_valid(*ret, "Overlay result between original inputs", true, true);
364 #endif
365 
366 #if GEOS_DEBUG_BINARYOP
367  std::cerr << "Attempt with original input succeeded" << std::endl;
368 #endif
369  return ret;
370  }
371  catch(const geos::util::TopologyException& ex) {
372  origException = ex;
373 #if GEOS_DEBUG_BINARYOP
374  std::cerr << "Original exception: " << ex.what() << std::endl;
375 #endif
376  }
377 #endif // USE_ORIGINAL_INPUT
378 
379  check_valid(*g0, "Input geom 0", true, true);
380  check_valid(*g1, "Input geom 1", true, true);
381 
382 #if USE_COMMONBITS_POLICY
383  // Try removing common bits (possibly obsoleted by snapping below)
384  //
385  // NOTE: this policy was _later_ implemented
386  // in JTS as EnhancedPrecisionOp
387  // TODO: consider using the now-ported EnhancedPrecisionOp
388  // here too
389  //
390  try {
391  GeomPtr rG0;
392  GeomPtr rG1;
393  precision::CommonBitsRemover cbr;
394 
395 #if GEOS_DEBUG_BINARYOP
396  std::cerr << "Trying with Common Bits Remover (CBR)" << std::endl;
397 #endif
398 
399  cbr.add(g0);
400  cbr.add(g1);
401 
402  rG0.reset(cbr.removeCommonBits(g0->clone()));
403  rG1.reset(cbr.removeCommonBits(g1->clone()));
404 
405 #if GEOS_DEBUG_BINARYOP
406  check_valid(*rG0, "CBR: geom 0 (after common-bits removal)");
407  check_valid(*rG1, "CBR: geom 1 (after common-bits removal)");
408 #endif
409 
410  ret.reset(_Op(rG0.get(), rG1.get()));
411 
412 #if GEOS_DEBUG_BINARYOP
413  check_valid(*ret, "CBR: result (before common-bits addition)");
414 #endif
415 
416  cbr.addCommonBits(ret.get());
417 
418 #if GEOS_CHECK_COMMONBITS_VALIDITY
419  check_valid(*ret, "CBR: result (after common-bits addition)", true);
420 #endif
421 
422 #if GEOS_DEBUG_BINARYOP
423  std::cerr << "Attempt with CBR succeeded" << std::endl;
424 #endif
425 
426  return ret;
427  }
428  catch(const geos::util::TopologyException& ex) {
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  ::geos::ignore_unused_variable_warning(ex);
462 #if GEOS_DEBUG_BINARYOP
463  std::cerr << "SNAP: " << ex.what() << std::endl;
464 #endif
465  }
466 
467 #endif // USE_SNAPPING_POLICY }
468 
469 // {
470 #if USE_PRECISION_REDUCTION_POLICY
471 
472 
473  // Try reducing precision
474  try {
475  long unsigned int g0scale =
476  static_cast<long unsigned int>(g0->getFactory()->getPrecisionModel()->getScale());
477  long unsigned int g1scale =
478  static_cast<long unsigned int>(g1->getFactory()->getPrecisionModel()->getScale());
479 
480 #if GEOS_DEBUG_BINARYOP
481  std::cerr << "Original input scales are: "
482  << g0scale
483  << " and "
484  << g1scale
485  << std::endl;
486 #endif
487 
488  double maxScale = 1e16;
489 
490  // Don't use a scale bigger than the input one
491  if(g0scale && static_cast<double>(g0scale) < maxScale) {
492  maxScale = static_cast<double>(g0scale);
493  }
494  if(g1scale && static_cast<double>(g1scale) < maxScale) {
495  maxScale = static_cast<double>(g1scale);
496  }
497 
498 
499  for(double scale = maxScale; scale >= 1; scale /= 10) {
500  PrecisionModel pm(scale);
501  GeometryFactory::Ptr gf = GeometryFactory::create(&pm);
502 #if GEOS_DEBUG_BINARYOP
503  std::cerr << "Trying with scale " << scale << std::endl;
504 #endif
505 
506  precision::GeometryPrecisionReducer reducer(*gf);
507  GeomPtr rG0(reducer.reduce(*g0));
508  GeomPtr rG1(reducer.reduce(*g1));
509 
510 #if GEOS_DEBUG_BINARYOP
511  check_valid(*rG0, "PR: geom 0 (after precision reduction)");
512  check_valid(*rG1, "PR: geom 1 (after precision reduction)");
513 #endif
514 
515  try {
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 #if GEOS_DEBUG_BINARYOP
536  std::cerr << "Reduced with scale (" << scale << "): "
537  << ex.what() << std::endl;
538 #endif
539  if(scale == 1) {
540  throw ex;
541  }
542  }
543 
544  }
545 
546  }
547  catch(const geos::util::TopologyException& ex) {
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  double maxTolerance = 0.04;
568  double minTolerance = 0.01;
569  double tolStep = 0.01;
570 
571  for(double tol = minTolerance; tol <= maxTolerance; tol += tolStep) {
572 #if GEOS_DEBUG_BINARYOP
573  std::cerr << "Trying simplifying with tolerance " << tol << std::endl;
574 #endif
575 
576  GeomPtr rG0(simplify::TopologyPreservingSimplifier::simplify(g0, tol));
577  GeomPtr rG1(simplify::TopologyPreservingSimplifier::simplify(g1, tol));
578 
579  try {
580  ret.reset(_Op(rG0.get(), rG1.get()));
581  return ret;
582  }
583  catch(const geos::util::TopologyException& ex) {
584  if(tol >= maxTolerance) {
585  throw ex;
586  }
587 #if GEOS_DEBUG_BINARYOP
588  std::cerr << "Simplified with tolerance (" << tol << "): "
589  << ex.what() << std::endl;
590 #endif
591  }
592 
593  }
594 
595  return ret;
596 
597  }
598  catch(const geos::util::TopologyException& ex) {
599 #if GEOS_DEBUG_BINARYOP
600  std::cerr << "Simplified: " << ex.what() << std::endl;
601 #endif
602  }
603 
604 #endif
605 // USE_TP_SIMPLIFY_POLICY }
606 
607 #if GEOS_DEBUG_BINARYOP
608  std::cerr << "No attempts worked to union " << std::endl;
609  std::cerr << "Input geometries:" << std::endl
610  << "<A>" << std::endl
611  << g0->toString() << std::endl
612  << "</A>" << std::endl
613  << "<B>" << std::endl
614  << g1->toString() << std::endl
615  << "</B>" << std::endl;
616 #endif
617 
618  throw origException;
619 }
620 
621 
622 } // namespace geos::geom
623 } // namespace geos
624 
625 #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:187
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:276
Basic namespace for all GEOS functionalities.
Definition: IndexedNestedRingTester.h:25
Snaps the vertices and segments of a Geometry to another Geometry'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...