00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
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>
00066
00067
00068
00069 #ifdef GEOS_DEBUG_BINARYOP
00070 # include <iostream>
00071 # include <iomanip>
00072 # include <sstream>
00073 #endif
00074
00075
00076
00077
00078
00079 #ifndef USE_ORIGINAL_INPUT
00080 # define USE_ORIGINAL_INPUT 1
00081 #endif
00082
00083
00084
00085
00086
00087
00088 #ifndef USE_PRECISION_REDUCTION_POLICY
00089
00090 #endif
00091
00092
00093
00094
00095
00096
00097 #ifndef USE_TP_SIMPLIFY_POLICY
00098
00099 #endif
00100
00101
00102
00103
00104
00105
00106 #ifndef USE_COMMONBITS_POLICY
00107 # define USE_COMMONBITS_POLICY 1
00108 #endif
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121 #define GEOS_CHECK_COMMONBITS_VALIDITY 1
00122
00123
00124
00125
00126 #ifndef USE_SNAPPING_POLICY
00127 # define USE_SNAPPING_POLICY 1
00128 #endif
00129
00130 namespace geos {
00131 namespace 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
00156 inline std::auto_ptr<Geometry>
00157 fix_snap_collapses(std::auto_ptr<Geometry> g, const std::string& label)
00158 {
00159
00160
00161 if ( g->getGeometryTypeId() == GEOS_POLYGON ||
00162 g->getGeometryTypeId() == GEOS_MULTIPOLYGON )
00163 {
00164
00165
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
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
00200 using geos::operation::overlay::snap::GeometrySnapper;
00201
00202
00203
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
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
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
00240 GeometrySnapper snapper1( operand1 );
00241 GeomPtr snapG1( snapper1.snapTo(*snapG0, snapTolerance) );
00242 snapG1 = fix_snap_collapses(snapG1, "SNAP: snapped geom 1");
00243
00244
00245
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
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
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
00295
00296
00297
00298
00299
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
00336
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
00363
00364
00365
00366
00367
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
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
00441
00442
00443 #if USE_TP_SIMPLIFY_POLICY
00444
00445
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
00490
00491 throw origException;
00492 }
00493
00494
00495 }
00496 }
00497
00498 #endif // GEOS_GEOM_BINARYOP_H