GEOS
3.3.5
|
00001 /********************************************************************** 00002 * $Id: BinaryOp.h 3606 2012-04-10 15:02:23Z strk $ 00003 * 00004 * GEOS - Geometry Engine Open Source 00005 * http://geos.refractions.net 00006 * 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/geom/Geometry.h> 00053 #include <geos/geom/PrecisionModel.h> 00054 #include <geos/precision/CommonBitsRemover.h> 00055 #include <geos/precision/SimpleGeometryPrecisionReducer.h> 00056 00057 #include <geos/operation/overlay/snap/GeometrySnapper.h> 00058 00059 #include <geos/simplify/TopologyPreservingSimplifier.h> 00060 #include <geos/operation/valid/IsValidOp.h> 00061 #include <geos/operation/valid/TopologyValidationError.h> 00062 #include <geos/util/TopologyException.h> 00063 #include <geos/util.h> 00064 00065 #include <memory> // for auto_ptr 00066 00067 //#define GEOS_DEBUG_BINARYOP 1 00068 00069 #ifdef GEOS_DEBUG_BINARYOP 00070 # include <iostream> 00071 # include <iomanip> 00072 # include <sstream> 00073 #endif 00074 00075 00076 /* 00077 * Always try original input first 00078 */ 00079 #ifndef USE_ORIGINAL_INPUT 00080 # define USE_ORIGINAL_INPUT 1 00081 #endif 00082 00083 /* 00084 * Define this to use PrecisionReduction policy 00085 * in an attempt at by-passing binary operation 00086 * robustness problems (handles TopologyExceptions) 00087 */ 00088 #ifndef USE_PRECISION_REDUCTION_POLICY 00089 //# define USE_PRECISION_REDUCTION_POLICY 1 00090 #endif 00091 00092 /* 00093 * Define this to use TopologyPreserving simplification policy 00094 * in an attempt at by-passing binary operation 00095 * robustness problems (handles TopologyExceptions) 00096 */ 00097 #ifndef USE_TP_SIMPLIFY_POLICY 00098 //# define USE_TP_SIMPLIFY_POLICY 1 00099 #endif 00100 00101 /* 00102 * Use common bits removal policy. 00103 * If enabled, this would be tried /before/ 00104 * Geometry snapping. 00105 */ 00106 #ifndef USE_COMMONBITS_POLICY 00107 # define USE_COMMONBITS_POLICY 1 00108 #endif 00109 00110 /* 00111 * Check validity of operation performed 00112 * by common bits removal policy. 00113 * 00114 * This matches what EnhancedPrecisionOp does in JTS 00115 * and fixes 5 tests of invalid outputs in our testsuite 00116 * (stmlf-cases-20061020-invalid-output.xml) 00117 * and breaks 1 test (robustness-invalid-output.xml) so much 00118 * to prevent a result. 00119 * 00120 */ 00121 #define GEOS_CHECK_COMMONBITS_VALIDITY 1 00122 00123 /* 00124 * Use snapping policy 00125 */ 00126 #ifndef USE_SNAPPING_POLICY 00127 # define USE_SNAPPING_POLICY 1 00128 #endif 00129 00130 namespace geos { 00131 namespace geom { // geos::geom 00132 00133 inline bool 00134 check_valid(const Geometry& g, const std::string& label) 00135 { 00136 operation::valid::IsValidOp ivo(&g); 00137 if ( ! ivo.isValid() ) 00138 { 00139 #ifdef GEOS_DEBUG_BINARYOP 00140 using operation::valid::TopologyValidationError; 00141 TopologyValidationError* err = ivo.getValidationError(); 00142 std::cerr << label << " is INVALID: " 00143 << err->toString() 00144 << " (" << std::setprecision(20) 00145 << err->getCoordinate() << ")" 00146 << std::endl; 00147 #else 00148 (void)label; 00149 #endif 00150 return false; 00151 } 00152 return true; 00153 } 00154 00155 /* A single component may become a multi component */ 00156 inline std::auto_ptr<Geometry> 00157 fix_self_intersections(std::auto_ptr<Geometry> g, const std::string& label) 00158 { 00159 // TODO: check for presence of self-intersections first ? 00160 return g->Union(); 00161 } 00162 00163 00169 template <class BinOp> 00170 std::auto_ptr<Geometry> 00171 SnapOp(const Geometry* g0, const Geometry *g1, BinOp _Op) 00172 { 00173 typedef std::auto_ptr<Geometry> GeomPtr; 00174 00175 #define CBR_BEFORE_SNAPPING 1 00176 00177 //using geos::precision::GeometrySnapper; 00178 using geos::operation::overlay::snap::GeometrySnapper; 00179 00180 // Snap tolerance must be computed on the original 00181 // (not commonbits-removed) geoms 00182 double snapTolerance = GeometrySnapper::computeOverlaySnapTolerance(*g0, *g1); 00183 #if GEOS_DEBUG_BINARYOP 00184 std::cerr<< std::setprecision(20) << "Computed snap tolerance: "<<snapTolerance<<std::endl; 00185 #endif 00186 00187 00188 #if CBR_BEFORE_SNAPPING 00189 // Compute common bits 00190 geos::precision::CommonBitsRemover cbr; 00191 cbr.add(g0); cbr.add(g1); 00192 #if GEOS_DEBUG_BINARYOP 00193 std::cerr<<"Computed common bits: "<<cbr.getCommonCoordinate()<<std::endl; 00194 #endif 00195 00196 // Now remove common bits 00197 GeomPtr rG0( cbr.removeCommonBits(g0->clone()) ); 00198 GeomPtr rG1( cbr.removeCommonBits(g1->clone()) ); 00199 00200 #if GEOS_DEBUG_BINARYOP 00201 check_valid(*rG0, "CBR: removed-bits geom 0"); 00202 check_valid(*rG1, "CBR: removed-bits geom 1"); 00203 #endif 00204 00205 const Geometry& operand0 = *rG0; 00206 const Geometry& operand1 = *rG1; 00207 #else // don't CBR before snapping 00208 const Geometry& operand0 = *g0; 00209 const Geometry& operand1 = *g1; 00210 #endif 00211 00212 00213 GeometrySnapper snapper0( operand0 ); 00214 GeomPtr snapG0( snapper0.snapTo(operand1, snapTolerance) ); 00215 snapG0 = fix_self_intersections(snapG0, "SNAP: snapped geom 0"); 00216 00217 // NOTE: second geom is snapped on the snapped first one 00218 GeometrySnapper snapper1( operand1 ); 00219 GeomPtr snapG1( snapper1.snapTo(*snapG0, snapTolerance) ); 00220 snapG1 = fix_self_intersections(snapG1, "SNAP: snapped geom 1"); 00221 00222 // Run the binary op 00223 GeomPtr result( _Op(snapG0.get(), snapG1.get()) ); 00224 00225 #if GEOS_DEBUG_BINARYOP 00226 check_valid(*result, "SNAP: result (before common-bits addition"); 00227 #endif 00228 00229 #if CBR_BEFORE_SNAPPING 00230 // Add common bits back in 00231 cbr.addCommonBits( result.get() ); 00232 result = fix_self_intersections(result, "SNAP: result (after common-bits addition)"); 00233 #endif 00234 00235 #if GEOS_DEBUG_BINARYOP 00236 check_valid(*result, "SNAP: result (after common-bits addition"); 00237 #endif 00238 00239 return result; 00240 } 00241 00242 template <class BinOp> 00243 std::auto_ptr<Geometry> 00244 BinaryOp(const Geometry* g0, const Geometry *g1, BinOp _Op) 00245 { 00246 typedef std::auto_ptr<Geometry> GeomPtr; 00247 00248 GeomPtr ret; 00249 geos::util::TopologyException origException; 00250 00251 #ifdef USE_ORIGINAL_INPUT 00252 // Try with original input 00253 try 00254 { 00255 #if GEOS_DEBUG_BINARYOP 00256 std::cerr << "Trying with original input." << std::endl; 00257 #endif 00258 ret.reset(_Op(g0, g1)); 00259 return ret; 00260 } 00261 catch (const geos::util::TopologyException& ex) 00262 { 00263 origException=ex; 00264 #if GEOS_DEBUG_BINARYOP 00265 std::cerr << "Original exception: " << ex.what() << std::endl; 00266 #endif 00267 } 00268 00269 //#if GEOS_DEBUG_BINARYOP 00270 // Should we just give up here ? 00271 check_valid(*g0, "Input geom 0"); 00272 check_valid(*g1, "Input geom 1"); 00273 //#endif 00274 00275 #endif // USE_ORIGINAL_INPUT 00276 00277 00278 #ifdef USE_COMMONBITS_POLICY 00279 // Try removing common bits (possibly obsoleted by snapping below) 00280 // 00281 // NOTE: this policy was _later_ implemented 00282 // in JTS as EnhancedPrecisionOp 00283 // TODO: consider using the now-ported EnhancedPrecisionOp 00284 // here too 00285 // 00286 try 00287 { 00288 GeomPtr rG0; 00289 GeomPtr rG1; 00290 precision::CommonBitsRemover cbr; 00291 00292 #if GEOS_DEBUG_BINARYOP 00293 std::cerr << "Trying with Common Bits Remover (CBR)" << std::endl; 00294 #endif 00295 00296 cbr.add(g0); 00297 cbr.add(g1); 00298 00299 rG0.reset( cbr.removeCommonBits(g0->clone()) ); 00300 rG1.reset( cbr.removeCommonBits(g1->clone()) ); 00301 00302 #if GEOS_DEBUG_BINARYOP 00303 check_valid(*rG0, "CBR: geom 0 (after common-bits removal)"); 00304 check_valid(*rG1, "CBR: geom 1 (after common-bits removal)"); 00305 #endif 00306 00307 ret.reset( _Op(rG0.get(), rG1.get()) ); 00308 00309 #if GEOS_DEBUG_BINARYOP 00310 check_valid(*ret, "CBR: result (before common-bits addition)"); 00311 #endif 00312 00313 cbr.addCommonBits( ret.get() ); 00314 00315 #if GEOS_DEBUG_BINARYOP 00316 check_valid(*ret, "CBR: result (after common-bits addition)"); 00317 #endif 00318 00319 // Common-bits removal policy could introduce self-intersections 00320 ret = fix_self_intersections(ret, "CBR: result (after common-bits addition)"); 00321 00322 #if GEOS_DEBUG_BINARYOP 00323 check_valid(*ret, "CBR: result (after common-bits addition and fix_self_intersections)"); 00324 #endif 00325 00326 #if GEOS_CHECK_COMMONBITS_VALIDITY 00327 // check that result is a valid geometry after the 00328 // reshift to orginal precision (see EnhancedPrecisionOp) 00329 using operation::valid::IsValidOp; 00330 using operation::valid::TopologyValidationError; 00331 IsValidOp ivo(ret.get()); 00332 if ( ! ivo.isValid() ) 00333 { 00334 TopologyValidationError* e = ivo.getValidationError(); 00335 throw geos::util::TopologyException( 00336 "Result of overlay became invalid " 00337 "after re-addin common bits of operand " 00338 "coordinates: " + e->toString(), 00339 e->getCoordinate()); 00340 } 00341 #endif // GEOS_CHECK_COMMONBITS_VALIDITY 00342 00343 return ret; 00344 } 00345 catch (const geos::util::TopologyException& ex) 00346 { 00347 ::geos::ignore_unused_variable_warning(ex); 00348 #if GEOS_DEBUG_BINARYOP 00349 std::cerr << "CBR: " << ex.what() << std::endl; 00350 #endif 00351 } 00352 #endif 00353 00354 // Try with snapping 00355 // 00356 // TODO: possible optimization would be reusing the 00357 // already common-bit-removed inputs and just 00358 // apply geometry snapping, whereas the current 00359 // SnapOp function does both. 00360 // { 00361 #if USE_SNAPPING_POLICY 00362 00363 #if GEOS_DEBUG_BINARYOP 00364 std::cerr << "Trying with snapping " << std::endl; 00365 #endif 00366 00367 try { 00368 ret = SnapOp(g0, g1, _Op); 00369 #if GEOS_DEBUG_BINARYOP 00370 std::cerr << "SnapOp succeeded" << std::endl; 00371 #endif 00372 return ret; 00373 00374 } 00375 catch (const geos::util::TopologyException& ex) 00376 { 00377 ::geos::ignore_unused_variable_warning(ex); 00378 #if GEOS_DEBUG_BINARYOP 00379 std::cerr << "SNAP: " << ex.what() << std::endl; 00380 #endif 00381 } 00382 00383 #endif // USE_SNAPPING_POLICY } 00384 00385 00386 00387 // { 00388 #if USE_PRECISION_REDUCTION_POLICY 00389 00390 00391 // Try reducing precision 00392 try 00393 { 00394 int maxPrecision=25; 00395 00396 for (int precision=maxPrecision; precision; --precision) 00397 { 00398 std::auto_ptr<PrecisionModel> pm(new PrecisionModel(precision)); 00399 #if GEOS_DEBUG_BINARYOP 00400 std::cerr << "Trying with precision " << precision << std::endl; 00401 #endif 00402 00403 precision::SimpleGeometryPrecisionReducer reducer( pm.get() ); 00404 GeomPtr rG0( reducer.reduce(g0) ); 00405 GeomPtr rG1( reducer.reduce(g1) ); 00406 00407 try 00408 { 00409 ret.reset( _Op(rG0.get(), rG1.get()) ); 00410 return ret; 00411 } 00412 catch (const geos::util::TopologyException& ex) 00413 { 00414 if ( precision == 1 ) throw ex; 00415 #if GEOS_DEBUG_BINARYOP 00416 std::cerr << "Reduced with precision (" << precision << "): " 00417 << ex.what() << std::endl; 00418 #endif 00419 } 00420 00421 } 00422 00423 } 00424 catch (const geos::util::TopologyException& ex) 00425 { 00426 #if GEOS_DEBUG_BINARYOP 00427 std::cerr << "Reduced: " << ex.what() << std::endl; 00428 #endif 00429 } 00430 00431 #endif 00432 // USE_PRECISION_REDUCTION_POLICY } 00433 00434 // { 00435 #if USE_TP_SIMPLIFY_POLICY 00436 00437 // Try simplifying 00438 try 00439 { 00440 00441 double maxTolerance = 0.04; 00442 double minTolerance = 0.01; 00443 double tolStep = 0.01; 00444 00445 for (double tol = minTolerance; tol <= maxTolerance; tol += tolStep) 00446 { 00447 #if GEOS_DEBUG_BINARYOP 00448 std::cerr << "Trying simplifying with tolerance " << tol << std::endl; 00449 #endif 00450 00451 GeomPtr rG0( simplify::TopologyPreservingSimplifier::simplify(g0, tol) ); 00452 GeomPtr rG1( simplify::TopologyPreservingSimplifier::simplify(g1, tol) ); 00453 00454 try 00455 { 00456 ret.reset( _Op(rG0.get(), rG1.get()) ); 00457 return ret; 00458 } 00459 catch (const geos::util::TopologyException& ex) 00460 { 00461 if ( tol >= maxTolerance ) throw ex; 00462 #if GEOS_DEBUG_BINARYOP 00463 std::cerr << "Simplified with tolerance (" << tol << "): " 00464 << ex.what() << std::endl; 00465 #endif 00466 } 00467 00468 } 00469 00470 return ret; 00471 00472 } 00473 catch (const geos::util::TopologyException& ex) 00474 { 00475 #if GEOS_DEBUG_BINARYOP 00476 std::cerr << "Simplified: " << ex.what() << std::endl; 00477 #endif 00478 } 00479 00480 #endif 00481 // USE_TP_SIMPLIFY_POLICY } 00482 00483 throw origException; 00484 } 00485 00486 00487 } // namespace geos::geom 00488 } // namespace geos 00489 00490 #endif // GEOS_GEOM_BINARYOP_H