GEOS
3.3.6
|
00001 /********************************************************************** 00002 * $Id: BinaryOp.h 3713 2012-09-10 08:05:28Z 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 #ifdef GEOS_DEBUG_BINARYOP 00160 std::cerr << label << " fix_self_intersection (UnaryUnion)" << std::endl; 00161 #endif 00162 // TODO: check for the sole self-intersections case ? 00163 if ( ! check_valid(*g, label) ) return g->Union(); 00164 else return std::auto_ptr<Geometry>(g->clone()); 00165 } 00166 00167 00173 template <class BinOp> 00174 std::auto_ptr<Geometry> 00175 SnapOp(const Geometry* g0, const Geometry *g1, BinOp _Op) 00176 { 00177 typedef std::auto_ptr<Geometry> GeomPtr; 00178 00179 #define CBR_BEFORE_SNAPPING 1 00180 00181 //using geos::precision::GeometrySnapper; 00182 using geos::operation::overlay::snap::GeometrySnapper; 00183 00184 // Snap tolerance must be computed on the original 00185 // (not commonbits-removed) geoms 00186 double snapTolerance = GeometrySnapper::computeOverlaySnapTolerance(*g0, *g1); 00187 #if GEOS_DEBUG_BINARYOP 00188 std::cerr<< std::setprecision(20) << "Computed snap tolerance: "<<snapTolerance<<std::endl; 00189 #endif 00190 00191 00192 #if CBR_BEFORE_SNAPPING 00193 // Compute common bits 00194 geos::precision::CommonBitsRemover cbr; 00195 cbr.add(g0); cbr.add(g1); 00196 #if GEOS_DEBUG_BINARYOP 00197 std::cerr<<"Computed common bits: "<<cbr.getCommonCoordinate()<<std::endl; 00198 #endif 00199 00200 // Now remove common bits 00201 GeomPtr rG0( cbr.removeCommonBits(g0->clone()) ); 00202 GeomPtr rG1( cbr.removeCommonBits(g1->clone()) ); 00203 00204 #if GEOS_DEBUG_BINARYOP 00205 check_valid(*rG0, "CBR: removed-bits geom 0"); 00206 check_valid(*rG1, "CBR: removed-bits geom 1"); 00207 #endif 00208 00209 const Geometry& operand0 = *rG0; 00210 const Geometry& operand1 = *rG1; 00211 #else // don't CBR before snapping 00212 const Geometry& operand0 = *g0; 00213 const Geometry& operand1 = *g1; 00214 #endif 00215 00216 00217 GeometrySnapper snapper0( operand0 ); 00218 GeomPtr snapG0( snapper0.snapTo(operand1, snapTolerance) ); 00219 snapG0 = fix_self_intersections(snapG0, "SNAP: snapped geom 0"); 00220 00221 // NOTE: second geom is snapped on the snapped first one 00222 GeometrySnapper snapper1( operand1 ); 00223 GeomPtr snapG1( snapper1.snapTo(*snapG0, snapTolerance) ); 00224 snapG1 = fix_self_intersections(snapG1, "SNAP: snapped geom 1"); 00225 00226 // Run the binary op 00227 GeomPtr result( _Op(snapG0.get(), snapG1.get()) ); 00228 00229 #if GEOS_DEBUG_BINARYOP 00230 check_valid(*result, "SNAP: result (before common-bits addition"); 00231 #endif 00232 00233 #if CBR_BEFORE_SNAPPING 00234 // Add common bits back in 00235 cbr.addCommonBits( result.get() ); 00236 result = fix_self_intersections(result, "SNAP: result (after common-bits addition)"); 00237 #endif 00238 00239 #if GEOS_DEBUG_BINARYOP 00240 check_valid(*result, "SNAP: result (after common-bits addition"); 00241 #endif 00242 00243 return result; 00244 } 00245 00246 template <class BinOp> 00247 std::auto_ptr<Geometry> 00248 BinaryOp(const Geometry* g0, const Geometry *g1, BinOp _Op) 00249 { 00250 typedef std::auto_ptr<Geometry> GeomPtr; 00251 00252 GeomPtr ret; 00253 geos::util::TopologyException origException; 00254 00255 #ifdef USE_ORIGINAL_INPUT 00256 // Try with original input 00257 try 00258 { 00259 #if GEOS_DEBUG_BINARYOP 00260 std::cerr << "Trying with original input." << std::endl; 00261 #endif 00262 ret.reset(_Op(g0, g1)); 00263 return ret; 00264 } 00265 catch (const geos::util::TopologyException& ex) 00266 { 00267 origException=ex; 00268 #if GEOS_DEBUG_BINARYOP 00269 std::cerr << "Original exception: " << ex.what() << std::endl; 00270 #endif 00271 } 00272 00273 //#if GEOS_DEBUG_BINARYOP 00274 // Should we just give up here ? 00275 check_valid(*g0, "Input geom 0"); 00276 check_valid(*g1, "Input geom 1"); 00277 //#endif 00278 00279 #endif // USE_ORIGINAL_INPUT 00280 00281 00282 #ifdef USE_COMMONBITS_POLICY 00283 // Try removing common bits (possibly obsoleted by snapping below) 00284 // 00285 // NOTE: this policy was _later_ implemented 00286 // in JTS as EnhancedPrecisionOp 00287 // TODO: consider using the now-ported EnhancedPrecisionOp 00288 // here too 00289 // 00290 try 00291 { 00292 GeomPtr rG0; 00293 GeomPtr rG1; 00294 precision::CommonBitsRemover cbr; 00295 00296 #if GEOS_DEBUG_BINARYOP 00297 std::cerr << "Trying with Common Bits Remover (CBR)" << std::endl; 00298 #endif 00299 00300 cbr.add(g0); 00301 cbr.add(g1); 00302 00303 rG0.reset( cbr.removeCommonBits(g0->clone()) ); 00304 rG1.reset( cbr.removeCommonBits(g1->clone()) ); 00305 00306 #if GEOS_DEBUG_BINARYOP 00307 check_valid(*rG0, "CBR: geom 0 (after common-bits removal)"); 00308 check_valid(*rG1, "CBR: geom 1 (after common-bits removal)"); 00309 #endif 00310 00311 ret.reset( _Op(rG0.get(), rG1.get()) ); 00312 00313 #if GEOS_DEBUG_BINARYOP 00314 check_valid(*ret, "CBR: result (before common-bits addition)"); 00315 #endif 00316 00317 cbr.addCommonBits( ret.get() ); 00318 00319 #if GEOS_DEBUG_BINARYOP 00320 check_valid(*ret, "CBR: result (after common-bits addition)"); 00321 #endif 00322 00323 // Common-bits removal policy could introduce self-intersections 00324 ret = fix_self_intersections(ret, "CBR: result (after common-bits addition)"); 00325 00326 #if GEOS_DEBUG_BINARYOP 00327 check_valid(*ret, "CBR: result (after common-bits addition and fix_self_intersections)"); 00328 #endif 00329 00330 #if GEOS_CHECK_COMMONBITS_VALIDITY 00331 // check that result is a valid geometry after the 00332 // reshift to orginal precision (see EnhancedPrecisionOp) 00333 using operation::valid::IsValidOp; 00334 using operation::valid::TopologyValidationError; 00335 IsValidOp ivo(ret.get()); 00336 if ( ! ivo.isValid() ) 00337 { 00338 TopologyValidationError* e = ivo.getValidationError(); 00339 throw geos::util::TopologyException( 00340 "Result of overlay became invalid " 00341 "after re-addin common bits of operand " 00342 "coordinates: " + e->toString(), 00343 e->getCoordinate()); 00344 } 00345 #endif // GEOS_CHECK_COMMONBITS_VALIDITY 00346 00347 return ret; 00348 } 00349 catch (const geos::util::TopologyException& ex) 00350 { 00351 ::geos::ignore_unused_variable_warning(ex); 00352 #if GEOS_DEBUG_BINARYOP 00353 std::cerr << "CBR: " << ex.what() << std::endl; 00354 #endif 00355 } 00356 #endif 00357 00358 // Try with snapping 00359 // 00360 // TODO: possible optimization would be reusing the 00361 // already common-bit-removed inputs and just 00362 // apply geometry snapping, whereas the current 00363 // SnapOp function does both. 00364 // { 00365 #if USE_SNAPPING_POLICY 00366 00367 #if GEOS_DEBUG_BINARYOP 00368 std::cerr << "Trying with snapping " << std::endl; 00369 #endif 00370 00371 try { 00372 ret = SnapOp(g0, g1, _Op); 00373 #if GEOS_DEBUG_BINARYOP 00374 std::cerr << "SnapOp succeeded" << std::endl; 00375 #endif 00376 return ret; 00377 00378 } 00379 catch (const geos::util::TopologyException& ex) 00380 { 00381 ::geos::ignore_unused_variable_warning(ex); 00382 #if GEOS_DEBUG_BINARYOP 00383 std::cerr << "SNAP: " << ex.what() << std::endl; 00384 #endif 00385 } 00386 00387 #endif // USE_SNAPPING_POLICY } 00388 00389 00390 00391 // { 00392 #if USE_PRECISION_REDUCTION_POLICY 00393 00394 00395 // Try reducing precision 00396 try 00397 { 00398 int maxPrecision=25; 00399 00400 for (int precision=maxPrecision; precision; --precision) 00401 { 00402 std::auto_ptr<PrecisionModel> pm(new PrecisionModel(precision)); 00403 #if GEOS_DEBUG_BINARYOP 00404 std::cerr << "Trying with precision " << precision << std::endl; 00405 #endif 00406 00407 precision::SimpleGeometryPrecisionReducer reducer( pm.get() ); 00408 GeomPtr rG0( reducer.reduce(g0) ); 00409 GeomPtr rG1( reducer.reduce(g1) ); 00410 00411 try 00412 { 00413 ret.reset( _Op(rG0.get(), rG1.get()) ); 00414 return ret; 00415 } 00416 catch (const geos::util::TopologyException& ex) 00417 { 00418 if ( precision == 1 ) throw ex; 00419 #if GEOS_DEBUG_BINARYOP 00420 std::cerr << "Reduced with precision (" << precision << "): " 00421 << ex.what() << std::endl; 00422 #endif 00423 } 00424 00425 } 00426 00427 } 00428 catch (const geos::util::TopologyException& ex) 00429 { 00430 #if GEOS_DEBUG_BINARYOP 00431 std::cerr << "Reduced: " << ex.what() << std::endl; 00432 #endif 00433 } 00434 00435 #endif 00436 // USE_PRECISION_REDUCTION_POLICY } 00437 00438 // { 00439 #if USE_TP_SIMPLIFY_POLICY 00440 00441 // Try simplifying 00442 try 00443 { 00444 00445 double maxTolerance = 0.04; 00446 double minTolerance = 0.01; 00447 double tolStep = 0.01; 00448 00449 for (double tol = minTolerance; tol <= maxTolerance; tol += tolStep) 00450 { 00451 #if GEOS_DEBUG_BINARYOP 00452 std::cerr << "Trying simplifying with tolerance " << tol << std::endl; 00453 #endif 00454 00455 GeomPtr rG0( simplify::TopologyPreservingSimplifier::simplify(g0, tol) ); 00456 GeomPtr rG1( simplify::TopologyPreservingSimplifier::simplify(g1, tol) ); 00457 00458 try 00459 { 00460 ret.reset( _Op(rG0.get(), rG1.get()) ); 00461 return ret; 00462 } 00463 catch (const geos::util::TopologyException& ex) 00464 { 00465 if ( tol >= maxTolerance ) throw ex; 00466 #if GEOS_DEBUG_BINARYOP 00467 std::cerr << "Simplified with tolerance (" << tol << "): " 00468 << ex.what() << std::endl; 00469 #endif 00470 } 00471 00472 } 00473 00474 return ret; 00475 00476 } 00477 catch (const geos::util::TopologyException& ex) 00478 { 00479 #if GEOS_DEBUG_BINARYOP 00480 std::cerr << "Simplified: " << ex.what() << std::endl; 00481 #endif 00482 } 00483 00484 #endif 00485 // USE_TP_SIMPLIFY_POLICY } 00486 00487 throw origException; 00488 } 00489 00490 00491 } // namespace geos::geom 00492 } // namespace geos 00493 00494 #endif // GEOS_GEOM_BINARYOP_H