GEOS
3.6.2
|
00001 /********************************************************************** 00002 * 00003 * GEOS - Geometry Engine Open Source 00004 * http://geos.osgeo.org 00005 * 00006 * Copyright (C) 2013 Sandro Santilli <strk@keybit.net> 00007 * Copyright (C) 2006 Refractions Research Inc. 00008 * 00009 * This is free software; you can redistribute and/or modify it under 00010 * the terms of the GNU Lesser General Public Licence as published 00011 * by the Free Software Foundation. 00012 * See the COPYING file for more information. 00013 * 00014 ********************************************************************** 00015 * 00016 * Last port: ORIGINAL WORK 00017 * 00018 ********************************************************************** 00019 * 00020 * This file provides a single templated function, taking two 00021 * const Geometry pointers, applying a binary operator to them 00022 * and returning a result Geometry in an auto_ptr<>. 00023 * 00024 * The binary operator is expected to take two const Geometry pointers 00025 * and return a newly allocated Geometry pointer, possibly throwing 00026 * a TopologyException to signal it couldn't succeed due to robustness 00027 * issues. 00028 * 00029 * This function will catch TopologyExceptions and try again with 00030 * slightly modified versions of the input. The following heuristic 00031 * is used: 00032 * 00033 * - Try with original input. 00034 * - Try removing common bits from input coordinate values 00035 * - Try snaping input geometries to each other 00036 * - Try snaping input coordinates to a increasing grid (size from 1/25 to 1) 00037 * - Try simplifiying input with increasing tolerance (from 0.01 to 0.04) 00038 * 00039 * If none of the step succeeds the original exception is thrown. 00040 * 00041 * Note that you can skip Grid snapping, Geometry snapping and Simplify policies 00042 * by a compile-time define when building geos. 00043 * See USE_TP_SIMPLIFY_POLICY, USE_PRECISION_REDUCTION_POLICY and 00044 * USE_SNAPPING_POLICY macros below. 00045 * 00046 * 00047 **********************************************************************/ 00048 00049 #ifndef GEOS_GEOM_BINARYOP_H 00050 #define GEOS_GEOM_BINARYOP_H 00051 00052 #include <geos/algorithm/BoundaryNodeRule.h> 00053 #include <geos/geom/Geometry.h> 00054 #include <geos/geom/GeometryCollection.h> 00055 #include <geos/geom/Polygon.h> 00056 #include <geos/geom/Lineal.h> 00057 #include <geos/geom/PrecisionModel.h> 00058 #include <geos/geom/GeometryFactory.h> 00059 #include <geos/precision/CommonBitsRemover.h> 00060 #include <geos/precision/SimpleGeometryPrecisionReducer.h> 00061 #include <geos/precision/GeometryPrecisionReducer.h> 00062 00063 #include <geos/operation/overlay/snap/GeometrySnapper.h> 00064 00065 #include <geos/simplify/TopologyPreservingSimplifier.h> 00066 #include <geos/operation/IsSimpleOp.h> 00067 #include <geos/operation/valid/IsValidOp.h> 00068 #include <geos/operation/valid/TopologyValidationError.h> 00069 #include <geos/util/TopologyException.h> 00070 #include <geos/util.h> 00071 00072 #include <memory> // for auto_ptr 00073 00074 //#define GEOS_DEBUG_BINARYOP 1 00075 00076 #ifdef GEOS_DEBUG_BINARYOP 00077 # include <iostream> 00078 # include <iomanip> 00079 # include <sstream> 00080 #endif 00081 00082 00083 /* 00084 * Always try original input first 00085 */ 00086 #ifndef USE_ORIGINAL_INPUT 00087 # define USE_ORIGINAL_INPUT 1 00088 #endif 00089 00090 /* 00091 * Define this to use PrecisionReduction policy 00092 * in an attempt at by-passing binary operation 00093 * robustness problems (handles TopologyExceptions) 00094 */ 00095 #ifndef USE_PRECISION_REDUCTION_POLICY 00096 # define USE_PRECISION_REDUCTION_POLICY 1 00097 #endif 00098 00099 /* 00100 * Define this to use TopologyPreserving simplification policy 00101 * in an attempt at by-passing binary operation 00102 * robustness problems (handles TopologyExceptions) 00103 */ 00104 #ifndef USE_TP_SIMPLIFY_POLICY 00105 //# define USE_TP_SIMPLIFY_POLICY 1 00106 #endif 00107 00108 /* 00109 * Use common bits removal policy. 00110 * If enabled, this would be tried /before/ 00111 * Geometry snapping. 00112 */ 00113 #ifndef USE_COMMONBITS_POLICY 00114 # define USE_COMMONBITS_POLICY 1 00115 #endif 00116 00117 /* 00118 * Check validity of operation performed 00119 * by common bits removal policy. 00120 * 00121 * This matches what EnhancedPrecisionOp does in JTS 00122 * and fixes 5 tests of invalid outputs in our testsuite 00123 * (stmlf-cases-20061020-invalid-output.xml) 00124 * and breaks 1 test (robustness-invalid-output.xml) so much 00125 * to prevent a result. 00126 * 00127 */ 00128 #define GEOS_CHECK_COMMONBITS_VALIDITY 1 00129 00130 /* 00131 * Use snapping policy 00132 */ 00133 #ifndef USE_SNAPPING_POLICY 00134 # define USE_SNAPPING_POLICY 1 00135 #endif 00136 00137 namespace geos { 00138 namespace geom { // geos::geom 00139 00140 inline bool 00141 check_valid(const Geometry& g, const std::string& label, bool doThrow=false, bool validOnly=false) 00142 { 00143 if ( dynamic_cast<const Lineal*>(&g) ) { 00144 if ( ! validOnly ) { 00145 operation::IsSimpleOp sop(g, algorithm::BoundaryNodeRule::getBoundaryEndPoint()); 00146 if ( ! sop.isSimple() ) 00147 { 00148 if ( doThrow ) { 00149 throw geos::util::TopologyException( 00150 label + " is not simple"); 00151 } 00152 return false; 00153 } 00154 } 00155 } else { 00156 operation::valid::IsValidOp ivo(&g); 00157 if ( ! ivo.isValid() ) 00158 { 00159 using operation::valid::TopologyValidationError; 00160 TopologyValidationError* err = ivo.getValidationError(); 00161 #ifdef GEOS_DEBUG_BINARYOP 00162 std::cerr << label << " is INVALID: " 00163 << err->toString() 00164 << " (" << std::setprecision(20) 00165 << err->getCoordinate() << ")" 00166 << std::endl; 00167 #endif 00168 if ( doThrow ) { 00169 throw geos::util::TopologyException( 00170 label + " is invalid: " + err->toString(), 00171 err->getCoordinate()); 00172 } 00173 return false; 00174 } 00175 } 00176 return true; 00177 } 00178 00179 /* 00180 * Attempt to fix noding of multilines and 00181 * self-intersection of multipolygons 00182 * 00183 * May return the input untouched. 00184 */ 00185 inline std::auto_ptr<Geometry> 00186 fix_self_intersections(std::auto_ptr<Geometry> g, const std::string& label) 00187 { 00188 ::geos::ignore_unused_variable_warning(label); 00189 #ifdef GEOS_DEBUG_BINARYOP 00190 std::cerr << label << " fix_self_intersection (UnaryUnion)" << std::endl; 00191 #endif 00192 00193 // Only multi-components can be fixed by UnaryUnion 00194 if ( ! dynamic_cast<const GeometryCollection*>(g.get()) ) return g; 00195 00196 using operation::valid::IsValidOp; 00197 00198 IsValidOp ivo(g.get()); 00199 00200 // Polygon is valid, nothing to do 00201 if ( ivo.isValid() ) return g; 00202 00203 // Not all invalidities can be fixed by this code 00204 00205 using operation::valid::TopologyValidationError; 00206 TopologyValidationError* err = ivo.getValidationError(); 00207 switch ( err->getErrorType() ) { 00208 case TopologyValidationError::eRingSelfIntersection: 00209 case TopologyValidationError::eTooFewPoints: // collapsed lines 00210 #ifdef GEOS_DEBUG_BINARYOP 00211 std::cerr << label << " ATTEMPT_TO_FIX: " << err->getErrorType() << ": " << *g << std::endl; 00212 #endif 00213 g = g->Union(); 00214 #ifdef GEOS_DEBUG_BINARYOP 00215 std::cerr << label << " ATTEMPT_TO_FIX succeeded.. " << std::endl; 00216 #endif 00217 return g; 00218 case TopologyValidationError::eSelfIntersection: 00219 // this one is within a single component, won't be fixed 00220 default: 00221 #ifdef GEOS_DEBUG_BINARYOP 00222 std::cerr << label << " invalidity is: " << err->getErrorType() << std::endl; 00223 #endif 00224 return g; 00225 } 00226 } 00227 00228 00234 template <class BinOp> 00235 std::auto_ptr<Geometry> 00236 SnapOp(const Geometry* g0, const Geometry *g1, BinOp _Op) 00237 { 00238 typedef std::auto_ptr<Geometry> GeomPtr; 00239 00240 #define CBR_BEFORE_SNAPPING 1 00241 00242 //using geos::precision::GeometrySnapper; 00243 using geos::operation::overlay::snap::GeometrySnapper; 00244 00245 // Snap tolerance must be computed on the original 00246 // (not commonbits-removed) geoms 00247 double snapTolerance = GeometrySnapper::computeOverlaySnapTolerance(*g0, *g1); 00248 #if GEOS_DEBUG_BINARYOP 00249 std::cerr<< std::setprecision(20) << "Computed snap tolerance: "<<snapTolerance<<std::endl; 00250 #endif 00251 00252 00253 #if CBR_BEFORE_SNAPPING 00254 // Compute common bits 00255 geos::precision::CommonBitsRemover cbr; 00256 cbr.add(g0); cbr.add(g1); 00257 #if GEOS_DEBUG_BINARYOP 00258 std::cerr<<"Computed common bits: "<<cbr.getCommonCoordinate()<<std::endl; 00259 #endif 00260 00261 // Now remove common bits 00262 GeomPtr rG0( cbr.removeCommonBits(g0->clone()) ); 00263 GeomPtr rG1( cbr.removeCommonBits(g1->clone()) ); 00264 00265 #if GEOS_DEBUG_BINARYOP 00266 check_valid(*rG0, "CBR: removed-bits geom 0"); 00267 check_valid(*rG1, "CBR: removed-bits geom 1"); 00268 #endif 00269 00270 const Geometry& operand0 = *rG0; 00271 const Geometry& operand1 = *rG1; 00272 #else // don't CBR before snapping 00273 const Geometry& operand0 = *g0; 00274 const Geometry& operand1 = *g1; 00275 #endif 00276 00277 00278 GeometrySnapper snapper0( operand0 ); 00279 GeomPtr snapG0( snapper0.snapTo(operand1, snapTolerance) ); 00280 //snapG0 = fix_self_intersections(snapG0, "SNAP: snapped geom 0"); 00281 00282 // NOTE: second geom is snapped on the snapped first one 00283 GeometrySnapper snapper1( operand1 ); 00284 GeomPtr snapG1( snapper1.snapTo(*snapG0, snapTolerance) ); 00285 //snapG1 = fix_self_intersections(snapG1, "SNAP: snapped geom 1"); 00286 00287 // Run the binary op 00288 GeomPtr result( _Op(snapG0.get(), snapG1.get()) ); 00289 00290 #if GEOS_DEBUG_BINARYOP 00291 check_valid(*result, "SNAP: result (before common-bits addition"); 00292 #endif 00293 00294 #if CBR_BEFORE_SNAPPING 00295 // Add common bits back in 00296 cbr.addCommonBits( result.get() ); 00297 //result = fix_self_intersections(result, "SNAP: result (after common-bits addition)"); 00298 00299 check_valid(*result, "CBR: result (after common-bits addition)", true); 00300 00301 #endif 00302 00303 return result; 00304 } 00305 00306 template <class BinOp> 00307 std::auto_ptr<Geometry> 00308 BinaryOp(const Geometry* g0, const Geometry *g1, BinOp _Op) 00309 { 00310 typedef std::auto_ptr<Geometry> GeomPtr; 00311 00312 GeomPtr ret; 00313 geos::util::TopologyException origException; 00314 00315 #ifdef USE_ORIGINAL_INPUT 00316 // Try with original input 00317 try 00318 { 00319 #if GEOS_DEBUG_BINARYOP 00320 std::cerr << "Trying with original input." << std::endl; 00321 #endif 00322 ret.reset(_Op(g0, g1)); 00323 return ret; 00324 } 00325 catch (const geos::util::TopologyException& ex) 00326 { 00327 origException=ex; 00328 #if GEOS_DEBUG_BINARYOP 00329 std::cerr << "Original exception: " << ex.what() << std::endl; 00330 #endif 00331 } 00332 00333 check_valid(*g0, "Input geom 0", true, true); 00334 check_valid(*g1, "Input geom 1", true, true); 00335 00336 #if GEOS_DEBUG_BINARYOP 00337 // Should we just give up here ? 00338 check_valid(*g0, "Input geom 0"); 00339 check_valid(*g1, "Input geom 1"); 00340 #endif 00341 00342 #endif // USE_ORIGINAL_INPUT 00343 00344 #ifdef USE_COMMONBITS_POLICY 00345 // Try removing common bits (possibly obsoleted by snapping below) 00346 // 00347 // NOTE: this policy was _later_ implemented 00348 // in JTS as EnhancedPrecisionOp 00349 // TODO: consider using the now-ported EnhancedPrecisionOp 00350 // here too 00351 // 00352 try 00353 { 00354 GeomPtr rG0; 00355 GeomPtr rG1; 00356 precision::CommonBitsRemover cbr; 00357 00358 #if GEOS_DEBUG_BINARYOP 00359 std::cerr << "Trying with Common Bits Remover (CBR)" << std::endl; 00360 #endif 00361 00362 cbr.add(g0); 00363 cbr.add(g1); 00364 00365 rG0.reset( cbr.removeCommonBits(g0->clone()) ); 00366 rG1.reset( cbr.removeCommonBits(g1->clone()) ); 00367 00368 #if GEOS_DEBUG_BINARYOP 00369 check_valid(*rG0, "CBR: geom 0 (after common-bits removal)"); 00370 check_valid(*rG1, "CBR: geom 1 (after common-bits removal)"); 00371 #endif 00372 00373 ret.reset( _Op(rG0.get(), rG1.get()) ); 00374 00375 #if GEOS_DEBUG_BINARYOP 00376 check_valid(*ret, "CBR: result (before common-bits addition)"); 00377 #endif 00378 00379 cbr.addCommonBits( ret.get() ); 00380 00381 check_valid(*ret, "CBR: result (after common-bits addition)", true); 00382 00383 #if GEOS_CHECK_COMMONBITS_VALIDITY 00384 // check that result is a valid geometry after the 00385 // reshift to orginal precision (see EnhancedPrecisionOp) 00386 using operation::valid::IsValidOp; 00387 using operation::valid::TopologyValidationError; 00388 IsValidOp ivo(ret.get()); 00389 if ( ! ivo.isValid() ) 00390 { 00391 TopologyValidationError* e = ivo.getValidationError(); 00392 throw geos::util::TopologyException( 00393 "Result of overlay became invalid " 00394 "after re-addin common bits of operand " 00395 "coordinates: " + e->toString(), 00396 e->getCoordinate()); 00397 } 00398 #endif // GEOS_CHECK_COMMONBITS_VALIDITY 00399 00400 return ret; 00401 } 00402 catch (const geos::util::TopologyException& ex) 00403 { 00404 ::geos::ignore_unused_variable_warning(ex); 00405 #if GEOS_DEBUG_BINARYOP 00406 std::cerr << "CBR: " << ex.what() << std::endl; 00407 #endif 00408 } 00409 #endif 00410 00411 // Try with snapping 00412 // 00413 // TODO: possible optimization would be reusing the 00414 // already common-bit-removed inputs and just 00415 // apply geometry snapping, whereas the current 00416 // SnapOp function does both. 00417 // { 00418 #if USE_SNAPPING_POLICY 00419 00420 #if GEOS_DEBUG_BINARYOP 00421 std::cerr << "Trying with snapping " << std::endl; 00422 #endif 00423 00424 try { 00425 ret = SnapOp(g0, g1, _Op); 00426 #if GEOS_DEBUG_BINARYOP 00427 std::cerr << "SnapOp succeeded" << std::endl; 00428 #endif 00429 return ret; 00430 00431 } 00432 catch (const geos::util::TopologyException& ex) 00433 { 00434 ::geos::ignore_unused_variable_warning(ex); 00435 #if GEOS_DEBUG_BINARYOP 00436 std::cerr << "SNAP: " << ex.what() << std::endl; 00437 #endif 00438 } 00439 00440 #endif // USE_SNAPPING_POLICY } 00441 00442 // { 00443 #if USE_PRECISION_REDUCTION_POLICY 00444 00445 00446 // Try reducing precision 00447 try 00448 { 00449 long unsigned int g0scale = 00450 static_cast<long unsigned int>(g0->getFactory()->getPrecisionModel()->getScale()); 00451 long unsigned int g1scale = 00452 static_cast<long unsigned int>(g1->getFactory()->getPrecisionModel()->getScale()); 00453 00454 #if GEOS_DEBUG_BINARYOP 00455 std::cerr << "Original input scales are: " 00456 << g0scale 00457 << " and " 00458 << g1scale 00459 << std::endl; 00460 #endif 00461 00462 double maxScale = 1e16; 00463 00464 // Don't use a scale bigger than the input one 00465 if ( g0scale && g0scale < maxScale ) maxScale = g0scale; 00466 if ( g1scale && g1scale < maxScale ) maxScale = g1scale; 00467 00468 00469 for (double scale=maxScale; scale >= 1; scale /= 10) 00470 { 00471 PrecisionModel pm(scale); 00472 GeometryFactory::unique_ptr gf = GeometryFactory::create(&pm); 00473 #if GEOS_DEBUG_BINARYOP 00474 std::cerr << "Trying with scale " << scale << std::endl; 00475 #endif 00476 00477 precision::GeometryPrecisionReducer reducer( *gf ); 00478 GeomPtr rG0( reducer.reduce(*g0) ); 00479 GeomPtr rG1( reducer.reduce(*g1) ); 00480 00481 try 00482 { 00483 ret.reset( _Op(rG0.get(), rG1.get()) ); 00484 // restore original precision (least precision between inputs) 00485 if ( g0->getFactory()->getPrecisionModel()->compareTo( g1->getFactory()->getPrecisionModel() ) < 0 ) { 00486 ret.reset( g0->getFactory()->createGeometry(ret.get()) ); 00487 } 00488 else { 00489 ret.reset( g1->getFactory()->createGeometry(ret.get()) ); 00490 } 00491 return ret; 00492 } 00493 catch (const geos::util::TopologyException& ex) 00494 { 00495 if ( scale == 1 ) throw ex; 00496 #if GEOS_DEBUG_BINARYOP 00497 std::cerr << "Reduced with scale (" << scale << "): " 00498 << ex.what() << std::endl; 00499 #endif 00500 } 00501 00502 } 00503 00504 } 00505 catch (const geos::util::TopologyException& ex) 00506 { 00507 #if GEOS_DEBUG_BINARYOP 00508 std::cerr << "Reduced: " << ex.what() << std::endl; 00509 #endif 00510 ::geos::ignore_unused_variable_warning(ex); 00511 } 00512 00513 #endif 00514 // USE_PRECISION_REDUCTION_POLICY } 00515 00516 00517 00518 00519 00520 // { 00521 #if USE_TP_SIMPLIFY_POLICY 00522 00523 // Try simplifying 00524 try 00525 { 00526 00527 double maxTolerance = 0.04; 00528 double minTolerance = 0.01; 00529 double tolStep = 0.01; 00530 00531 for (double tol = minTolerance; tol <= maxTolerance; tol += tolStep) 00532 { 00533 #if GEOS_DEBUG_BINARYOP 00534 std::cerr << "Trying simplifying with tolerance " << tol << std::endl; 00535 #endif 00536 00537 GeomPtr rG0( simplify::TopologyPreservingSimplifier::simplify(g0, tol) ); 00538 GeomPtr rG1( simplify::TopologyPreservingSimplifier::simplify(g1, tol) ); 00539 00540 try 00541 { 00542 ret.reset( _Op(rG0.get(), rG1.get()) ); 00543 return ret; 00544 } 00545 catch (const geos::util::TopologyException& ex) 00546 { 00547 if ( tol >= maxTolerance ) throw ex; 00548 #if GEOS_DEBUG_BINARYOP 00549 std::cerr << "Simplified with tolerance (" << tol << "): " 00550 << ex.what() << std::endl; 00551 #endif 00552 } 00553 00554 } 00555 00556 return ret; 00557 00558 } 00559 catch (const geos::util::TopologyException& ex) 00560 { 00561 #if GEOS_DEBUG_BINARYOP 00562 std::cerr << "Simplified: " << ex.what() << std::endl; 00563 #endif 00564 } 00565 00566 #endif 00567 // USE_TP_SIMPLIFY_POLICY } 00568 00569 throw origException; 00570 } 00571 00572 00573 } // namespace geos::geom 00574 } // namespace geos 00575 00576 #endif // GEOS_GEOM_BINARYOP_H