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