GEOS  3.3.2
BinaryOp.h
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