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