GEOS  3.3.3
BinaryOp.h
1 /**********************************************************************
2  * $Id: BinaryOp.h 3552 2011-12-15 14:18:38Z strk $
3  *
4  * GEOS - Geometry Engine Open Source
5  * http://geos.refractions.net
6  *
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/geom/Geometry.h>
53 #include <geos/geom/PrecisionModel.h>
54 #include <geos/precision/CommonBitsRemover.h>
55 #include <geos/precision/SimpleGeometryPrecisionReducer.h>
56 
57 #include <geos/operation/overlay/snap/GeometrySnapper.h>
58 
59 #include <geos/simplify/TopologyPreservingSimplifier.h>
60 #include <geos/operation/valid/IsValidOp.h>
61 #include <geos/operation/valid/TopologyValidationError.h>
62 #include <geos/util/TopologyException.h>
63 #include <geos/util.h>
64 
65 #include <memory> // for auto_ptr
66 
67 //#define GEOS_DEBUG_BINARYOP 1
68 
69 #ifdef GEOS_DEBUG_BINARYOP
70 # include <iostream>
71 # include <iomanip>
72 # include <sstream>
73 #endif
74 
75 
76 /*
77  * Always try original input first
78  */
79 #ifndef USE_ORIGINAL_INPUT
80 # define USE_ORIGINAL_INPUT 1
81 #endif
82 
83 /*
84  * Define this to use PrecisionReduction policy
85  * in an attempt at by-passing binary operation
86  * robustness problems (handles TopologyExceptions)
87  */
88 #ifndef USE_PRECISION_REDUCTION_POLICY
89 //# define USE_PRECISION_REDUCTION_POLICY 1
90 #endif
91 
92 /*
93  * Define this to use TopologyPreserving simplification policy
94  * in an attempt at by-passing binary operation
95  * robustness problems (handles TopologyExceptions)
96  */
97 #ifndef USE_TP_SIMPLIFY_POLICY
98 //# define USE_TP_SIMPLIFY_POLICY 1
99 #endif
100 
101 /*
102  * Use common bits removal policy.
103  * If enabled, this would be tried /before/
104  * Geometry snapping.
105  */
106 #ifndef USE_COMMONBITS_POLICY
107 # define USE_COMMONBITS_POLICY 1
108 #endif
109 
110 /*
111  * Check validity of operation performed
112  * by common bits removal policy.
113  *
114  * This matches what EnhancedPrecisionOp does in JTS
115  * and fixes 5 tests of invalid outputs in our testsuite
116  * (stmlf-cases-20061020-invalid-output.xml)
117  * and breaks 1 test (robustness-invalid-output.xml) so much
118  * to prevent a result.
119  *
120  */
121 #define GEOS_CHECK_COMMONBITS_VALIDITY 1
122 
123 /*
124  * Use snapping policy
125  */
126 #ifndef USE_SNAPPING_POLICY
127 # define USE_SNAPPING_POLICY 1
128 #endif
129 
130 namespace geos {
131 namespace geom { // geos::geom
132 
133 inline bool
134 check_valid(const Geometry& g, const std::string& label)
135 {
136  operation::valid::IsValidOp ivo(&g);
137  if ( ! ivo.isValid() )
138  {
139 #ifdef GEOS_DEBUG_BINARYOP
140  using operation::valid::TopologyValidationError;
141  TopologyValidationError* err = ivo.getValidationError();
142  std::cerr << label << " is INVALID: "
143  << err->toString()
144  << " (" << std::setprecision(20)
145  << err->getCoordinate() << ")"
146  << std::endl;
147 #else
148  (void)label;
149 #endif
150  return false;
151  }
152  return true;
153 }
154 
155 /* A single component may become a multi component */
156 inline std::auto_ptr<Geometry>
157 fix_snap_collapses(std::auto_ptr<Geometry> g, const std::string& label)
158 {
159 
160  // Areal geometries may become self-intersecting on snapping
161  if ( g->getGeometryTypeId() == GEOS_POLYGON ||
162  g->getGeometryTypeId() == GEOS_MULTIPOLYGON )
163  {
164 
165  // TODO: use only ConsistentAreaTester
166  if ( ! check_valid(*g, label) ) {
167 #if GEOS_DEBUG_BINARYOP
168  std::cerr << label << ": self-unioning" << std::endl;
169 #endif
170  g = g->Union();
171 #if GEOS_DEBUG_BINARYOP
172  std::stringstream ss;
173  ss << label << " self-unioned";
174  check_valid(*g, ss.str());
175 #endif
176  }
177 
178  }
179 
180  // TODO: check linear collapses ? (!isSimple)
181 
182  return g;
183 }
184 
185 
191 template <class BinOp>
192 std::auto_ptr<Geometry>
193 SnapOp(const Geometry* g0, const Geometry *g1, BinOp _Op)
194 {
195  typedef std::auto_ptr<Geometry> GeomPtr;
196 
197 #define CBR_BEFORE_SNAPPING 1
198 
199  //using geos::precision::GeometrySnapper;
201 
202  // Snap tolerance must be computed on the original
203  // (not commonbits-removed) geoms
204  double snapTolerance = GeometrySnapper::computeOverlaySnapTolerance(*g0, *g1);
205 #if GEOS_DEBUG_BINARYOP
206  std::cerr<< std::setprecision(20) << "Computed snap tolerance: "<<snapTolerance<<std::endl;
207 #endif
208 
209 
210 #if CBR_BEFORE_SNAPPING
211  // Compute common bits
213  cbr.add(g0); cbr.add(g1);
214 #if GEOS_DEBUG_BINARYOP
215  std::cerr<<"Computed common bits: "<<cbr.getCommonCoordinate()<<std::endl;
216 #endif
217 
218  // Now remove common bits
219  GeomPtr rG0( cbr.removeCommonBits(g0->clone()) );
220  GeomPtr rG1( cbr.removeCommonBits(g1->clone()) );
221 
222 #if GEOS_DEBUG_BINARYOP
223  check_valid(*rG0, "CBR: removed-bits geom 0");
224  check_valid(*rG1, "CBR: removed-bits geom 1");
225 #endif
226 
227  const Geometry& operand0 = *rG0;
228  const Geometry& operand1 = *rG1;
229 #else // don't CBR before snapping
230  const Geometry& operand0 = *g0;
231  const Geometry& operand1 = *g1;
232 #endif
233 
234 
235  GeometrySnapper snapper0( operand0 );
236  GeomPtr snapG0( snapper0.snapTo(operand1, snapTolerance) );
237  snapG0 = fix_snap_collapses(snapG0, "SNAP: snapped geom 0");
238 
239  // NOTE: second geom is snapped on the snapped first one
240  GeometrySnapper snapper1( operand1 );
241  GeomPtr snapG1( snapper1.snapTo(*snapG0, snapTolerance) );
242  snapG1 = fix_snap_collapses(snapG1, "SNAP: snapped geom 1");
243 
244 
245  // Run the binary op
246  GeomPtr result( _Op(snapG0.get(), snapG1.get()) );
247 
248 #if GEOS_DEBUG_BINARYOP
249  check_valid(*result, "SNAP: result (before common-bits addition");
250 #endif
251 
252 #if CBR_BEFORE_SNAPPING
253  // Add common bits back in
254  cbr.addCommonBits( result.get() );
255 #endif
256 
257 #if GEOS_DEBUG_BINARYOP
258  check_valid(*result, "SNAP: result (after common-bits addition");
259 #endif
260 
261  return result;
262 }
263 
264 template <class BinOp>
265 std::auto_ptr<Geometry>
266 BinaryOp(const Geometry* g0, const Geometry *g1, BinOp _Op)
267 {
268  typedef std::auto_ptr<Geometry> GeomPtr;
269 
270  GeomPtr ret;
271  geos::util::TopologyException origException;
272 
273 #ifdef USE_ORIGINAL_INPUT
274  // Try with original input
275  try
276  {
277 #if GEOS_DEBUG_BINARYOP
278  std::cerr << "Trying with original input." << std::endl;
279 #endif
280  ret.reset(_Op(g0, g1));
281  return ret;
282  }
283  catch (const geos::util::TopologyException& ex)
284  {
285  origException=ex;
286 #if GEOS_DEBUG_BINARYOP
287  std::cerr << "Original exception: " << ex.what() << std::endl;
288 #endif
289  }
290 #endif // USE_ORIGINAL_INPUT
291 
292 
293 #ifdef USE_COMMONBITS_POLICY
294  // Try removing common bits (possibly obsoleted by snapping below)
295  //
296  // NOTE: this policy was _later_ implemented
297  // in JTS as EnhancedPrecisionOp
298  // TODO: consider using the now-ported EnhancedPrecisionOp
299  // here too
300  //
301  try
302  {
303  GeomPtr rG0;
304  GeomPtr rG1;
305  precision::CommonBitsRemover cbr;
306 
307 #if GEOS_DEBUG_BINARYOP
308  std::cerr << "Trying with Common Bits Remover (CBR)" << std::endl;
309 #endif
310 
311  cbr.add(g0);
312  cbr.add(g1);
313 
314  rG0.reset( cbr.removeCommonBits(g0->clone()) );
315  rG1.reset( cbr.removeCommonBits(g1->clone()) );
316 
317 #if GEOS_DEBUG_BINARYOP
318  check_valid(*rG0, "CBR: geom 0 (after common-bits removal)");
319  check_valid(*rG1, "CBR: geom 1 (after common-bits removal)");
320 #endif
321 
322  ret.reset( _Op(rG0.get(), rG1.get()) );
323 
324 #if GEOS_DEBUG_BINARYOP
325  check_valid(*ret, "CBR: result (before common-bits addition)");
326 #endif
327 
328  cbr.addCommonBits( ret.get() );
329 
330 #if GEOS_DEBUG_BINARYOP
331  check_valid(*ret, "CBR: result (before common-bits addition)");
332 #endif
333 
334 #if GEOS_CHECK_COMMONBITS_VALIDITY
335  // check that result is a valid geometry after the
336  // reshift to orginal precision (see EnhancedPrecisionOp)
337  using operation::valid::IsValidOp;
338  using operation::valid::TopologyValidationError;
339  IsValidOp ivo(ret.get());
340  if ( ! ivo.isValid() )
341  {
342  TopologyValidationError* e = ivo.getValidationError();
344  "Result of overlay became invalid "
345  "after re-addin common bits of operand "
346  "coordinates: " + e->toString(),
347  e->getCoordinate());
348  }
349 #endif // GEOS_CHECK_COMMONBITS_VALIDITY
350 
351  return ret;
352  }
353  catch (const geos::util::TopologyException& ex)
354  {
355  ::geos::ignore_unused_variable_warning(ex);
356 #if GEOS_DEBUG_BINARYOP
357  std::cerr << "CBR: " << ex.what() << std::endl;
358 #endif
359  }
360 #endif
361 
362  // Try with snapping
363  //
364  // TODO: possible optimization would be reusing the
365  // already common-bit-removed inputs and just
366  // apply geometry snapping, whereas the current
367  // SnapOp function does both.
368 // {
369 #if USE_SNAPPING_POLICY
370 
371 #if GEOS_DEBUG_BINARYOP
372  std::cerr << "Trying with snapping " << std::endl;
373 #endif
374 
375  try {
376  ret = SnapOp(g0, g1, _Op);
377 #if GEOS_DEBUG_BINARYOP
378  std::cerr << "SnapOp succeeded" << std::endl;
379 #endif
380  return ret;
381 
382  }
383  catch (const geos::util::TopologyException& ex)
384  {
385  ::geos::ignore_unused_variable_warning(ex);
386 #if GEOS_DEBUG_BINARYOP
387  std::cerr << "SNAP: " << ex.what() << std::endl;
388 #endif
389  }
390 
391 #endif // USE_SNAPPING_POLICY }
392 
393 
394 
395 // {
396 #if USE_PRECISION_REDUCTION_POLICY
397 
398 
399  // Try reducing precision
400  try
401  {
402  int maxPrecision=25;
403 
404  for (int precision=maxPrecision; precision; --precision)
405  {
406  std::auto_ptr<PrecisionModel> pm(new PrecisionModel(precision));
407 #if GEOS_DEBUG_BINARYOP
408  std::cerr << "Trying with precision " << precision << std::endl;
409 #endif
410 
411  precision::SimpleGeometryPrecisionReducer reducer( pm.get() );
412  GeomPtr rG0( reducer.reduce(g0) );
413  GeomPtr rG1( reducer.reduce(g1) );
414 
415  try
416  {
417  ret.reset( _Op(rG0.get(), rG1.get()) );
418  return ret;
419  }
420  catch (const geos::util::TopologyException& ex)
421  {
422  if ( precision == 1 ) throw ex;
423 #if GEOS_DEBUG_BINARYOP
424  std::cerr << "Reduced with precision (" << precision << "): "
425  << ex.what() << std::endl;
426 #endif
427  }
428 
429  }
430 
431  }
432  catch (const geos::util::TopologyException& ex)
433  {
434 #if GEOS_DEBUG_BINARYOP
435  std::cerr << "Reduced: " << ex.what() << std::endl;
436 #endif
437  }
438 
439 #endif
440 // USE_PRECISION_REDUCTION_POLICY }
441 
442 // {
443 #if USE_TP_SIMPLIFY_POLICY
444 
445  // Try simplifying
446  try
447  {
448 
449  double maxTolerance = 0.04;
450  double minTolerance = 0.01;
451  double tolStep = 0.01;
452 
453  for (double tol = minTolerance; tol <= maxTolerance; tol += tolStep)
454  {
455 #if GEOS_DEBUG_BINARYOP
456  std::cerr << "Trying simplifying with tolerance " << tol << std::endl;
457 #endif
458 
459  GeomPtr rG0( simplify::TopologyPreservingSimplifier::simplify(g0, tol) );
460  GeomPtr rG1( simplify::TopologyPreservingSimplifier::simplify(g1, tol) );
461 
462  try
463  {
464  ret.reset( _Op(rG0.get(), rG1.get()) );
465  return ret;
466  }
467  catch (const geos::util::TopologyException& ex)
468  {
469  if ( tol >= maxTolerance ) throw ex;
470 #if GEOS_DEBUG_BINARYOP
471  std::cerr << "Simplified with tolerance (" << tol << "): "
472  << ex.what() << std::endl;
473 #endif
474  }
475 
476  }
477 
478  return ret;
479 
480  }
481  catch (const geos::util::TopologyException& ex)
482  {
483 #if GEOS_DEBUG_BINARYOP
484  std::cerr << "Simplified: " << ex.what() << std::endl;
485 #endif
486  }
487 
488 #endif
489 // USE_TP_SIMPLIFY_POLICY }
490 
491  throw origException;
492 }
493 
494 
495 } // namespace geos::geom
496 } // namespace geos
497 
498 #endif // GEOS_GEOM_BINARYOP_H