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