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