BinaryOp.h

00001 /**********************************************************************
00002  *
00003  * GEOS - Geometry Engine Open Source
00004  * http://geos.osgeo.org
00005  *
00006  * Copyright (C) 2013 Sandro Santilli <strk@keybit.net>
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/algorithm/BoundaryNodeRule.h>
00053 #include <geos/geom/Geometry.h>
00054 #include <geos/geom/GeometryCollection.h>
00055 #include <geos/geom/Polygon.h>
00056 #include <geos/geom/Lineal.h>
00057 #include <geos/geom/PrecisionModel.h>
00058 #include <geos/geom/GeometryFactory.h>
00059 #include <geos/precision/CommonBitsRemover.h>
00060 #include <geos/precision/SimpleGeometryPrecisionReducer.h>
00061 #include <geos/precision/GeometryPrecisionReducer.h>
00062 
00063 #include <geos/operation/overlay/snap/GeometrySnapper.h>
00064 
00065 #include <geos/simplify/TopologyPreservingSimplifier.h>
00066 #include <geos/operation/IsSimpleOp.h>
00067 #include <geos/operation/valid/IsValidOp.h>
00068 #include <geos/operation/valid/TopologyValidationError.h>
00069 #include <geos/util/TopologyException.h>
00070 #include <geos/util.h>
00071 
00072 #include <memory> // for auto_ptr
00073 
00074 //#define GEOS_DEBUG_BINARYOP 1
00075 
00076 #ifdef GEOS_DEBUG_BINARYOP
00077 # include <iostream>
00078 # include <iomanip>
00079 # include <sstream>
00080 #endif
00081 
00082 
00083 /*
00084  * Always try original input first
00085  */
00086 #ifndef USE_ORIGINAL_INPUT
00087 # define USE_ORIGINAL_INPUT 1
00088 #endif
00089 
00090 /*
00091  * Define this to use PrecisionReduction policy
00092  * in an attempt at by-passing binary operation
00093  * robustness problems (handles TopologyExceptions)
00094  */
00095 #ifndef USE_PRECISION_REDUCTION_POLICY
00096 # define USE_PRECISION_REDUCTION_POLICY 1
00097 #endif
00098 
00099 /*
00100  * Define this to use TopologyPreserving simplification policy
00101  * in an attempt at by-passing binary operation
00102  * robustness problems (handles TopologyExceptions)
00103  */
00104 #ifndef USE_TP_SIMPLIFY_POLICY 
00105 //# define USE_TP_SIMPLIFY_POLICY 1
00106 #endif
00107 
00108 /*
00109  * Use common bits removal policy.
00110  * If enabled, this would be tried /before/
00111  * Geometry snapping.
00112  */
00113 #ifndef USE_COMMONBITS_POLICY
00114 # define USE_COMMONBITS_POLICY 1
00115 #endif
00116 
00117 /*
00118  * Check validity of operation performed
00119  * by common bits removal policy.
00120  *
00121  * This matches what EnhancedPrecisionOp does in JTS
00122  * and fixes 5 tests of invalid outputs in our testsuite
00123  * (stmlf-cases-20061020-invalid-output.xml)
00124  * and breaks 1 test (robustness-invalid-output.xml) so much
00125  * to prevent a result.
00126  *
00127  */
00128 #define GEOS_CHECK_COMMONBITS_VALIDITY 1
00129 
00130 /*
00131  * Use snapping policy
00132  */
00133 #ifndef USE_SNAPPING_POLICY
00134 # define USE_SNAPPING_POLICY 1
00135 #endif
00136 
00137 namespace geos {
00138 namespace geom { // geos::geom
00139 
00140 inline bool
00141 check_valid(const Geometry& g, const std::string& label, bool doThrow=false, bool validOnly=false)
00142 {
00143   if ( dynamic_cast<const Lineal*>(&g) ) {
00144     if ( ! validOnly ) {
00145       operation::IsSimpleOp sop(g, algorithm::BoundaryNodeRule::getBoundaryEndPoint());
00146       if ( ! sop.isSimple() )
00147       {
00148         if ( doThrow ) {
00149           throw geos::util::TopologyException(
00150             label + " is not simple");
00151         }
00152         return false;
00153       }
00154     }
00155   } else {
00156     operation::valid::IsValidOp ivo(&g);
00157     if ( ! ivo.isValid() )
00158     {
00159       using operation::valid::TopologyValidationError;
00160       TopologyValidationError* err = ivo.getValidationError();
00161 #ifdef GEOS_DEBUG_BINARYOP
00162       std::cerr << label << " is INVALID: "
00163         << err->toString()
00164         << " (" << std::setprecision(20)
00165         << err->getCoordinate() << ")"
00166         << std::endl;
00167 #endif
00168       if ( doThrow ) {
00169         throw geos::util::TopologyException(
00170           label + " is invalid: " + err->toString(),
00171                 err->getCoordinate());
00172       }
00173       return false;
00174     } 
00175   }
00176         return true;
00177 }
00178 
00179 /*
00180  * Attempt to fix noding of multilines and
00181  * self-intersection of multipolygons 
00182  *
00183  * May return the input untouched.
00184  */
00185 inline std::auto_ptr<Geometry>
00186 fix_self_intersections(std::auto_ptr<Geometry> g, const std::string& label)
00187 {
00188 #ifdef GEOS_DEBUG_BINARYOP
00189         std::cerr << label << " fix_self_intersection (UnaryUnion)" << std::endl;
00190 #endif
00191 
00192   // Only multi-components can be fixed by UnaryUnion
00193   if ( ! dynamic_cast<const GeometryCollection*>(g.get()) ) return g;
00194 
00195   using operation::valid::IsValidOp;
00196 
00197   IsValidOp ivo(g.get());
00198 
00199   // Polygon is valid, nothing to do
00200   if ( ivo.isValid() ) return g;
00201 
00202   // Not all invalidities can be fixed by this code
00203 
00204   using operation::valid::TopologyValidationError;
00205   TopologyValidationError* err = ivo.getValidationError();
00206   switch ( err->getErrorType() ) {
00207     case TopologyValidationError::eRingSelfIntersection:
00208     case TopologyValidationError::eTooFewPoints: // collapsed lines
00209 #ifdef GEOS_DEBUG_BINARYOP
00210             std::cerr << label << " ATTEMPT_TO_FIX: " << err->getErrorType() << ": " << *g << std::endl;
00211 #endif
00212       g = g->Union();
00213 #ifdef GEOS_DEBUG_BINARYOP
00214             std::cerr << label << " ATTEMPT_TO_FIX succeeded.. " << std::endl;
00215 #endif
00216       return g;
00217     case TopologyValidationError::eSelfIntersection:
00218       // this one is within a single component, won't be fixed
00219     default:
00220 #ifdef GEOS_DEBUG_BINARYOP
00221             std::cerr << label << " invalidity is: " << err->getErrorType() << std::endl;
00222 #endif
00223       return g;
00224   }
00225 }
00226 
00227 
00233 template <class BinOp>
00234 std::auto_ptr<Geometry>
00235 SnapOp(const Geometry* g0, const Geometry *g1, BinOp _Op)
00236 {
00237         typedef std::auto_ptr<Geometry> GeomPtr;
00238 
00239 #define CBR_BEFORE_SNAPPING 1
00240 
00241         //using geos::precision::GeometrySnapper;
00242         using geos::operation::overlay::snap::GeometrySnapper;
00243 
00244         // Snap tolerance must be computed on the original
00245         // (not commonbits-removed) geoms
00246         double snapTolerance = GeometrySnapper::computeOverlaySnapTolerance(*g0, *g1);
00247 #if GEOS_DEBUG_BINARYOP
00248         std::cerr<< std::setprecision(20) << "Computed snap tolerance: "<<snapTolerance<<std::endl;
00249 #endif
00250 
00251 
00252 #if CBR_BEFORE_SNAPPING
00253         // Compute common bits
00254         geos::precision::CommonBitsRemover cbr;
00255         cbr.add(g0); cbr.add(g1);
00256 #if GEOS_DEBUG_BINARYOP
00257         std::cerr<<"Computed common bits: "<<cbr.getCommonCoordinate()<<std::endl;
00258 #endif
00259 
00260         // Now remove common bits
00261         GeomPtr rG0( cbr.removeCommonBits(g0->clone()) );
00262         GeomPtr rG1( cbr.removeCommonBits(g1->clone()) );
00263 
00264 #if GEOS_DEBUG_BINARYOP
00265         check_valid(*rG0, "CBR: removed-bits geom 0");
00266         check_valid(*rG1, "CBR: removed-bits geom 1");
00267 #endif
00268 
00269         const Geometry& operand0 = *rG0;
00270         const Geometry& operand1 = *rG1;
00271 #else // don't CBR before snapping
00272         const Geometry& operand0 = *g0;
00273         const Geometry& operand1 = *g1;
00274 #endif
00275 
00276 
00277         GeometrySnapper snapper0( operand0 );
00278         GeomPtr snapG0( snapper0.snapTo(operand1, snapTolerance) );
00279         //snapG0 = fix_self_intersections(snapG0, "SNAP: snapped geom 0");
00280 
00281         // NOTE: second geom is snapped on the snapped first one
00282         GeometrySnapper snapper1( operand1 );
00283         GeomPtr snapG1( snapper1.snapTo(*snapG0, snapTolerance) );
00284         //snapG1 = fix_self_intersections(snapG1, "SNAP: snapped geom 1");
00285 
00286         // Run the binary op
00287         GeomPtr result( _Op(snapG0.get(), snapG1.get()) );
00288 
00289 #if GEOS_DEBUG_BINARYOP
00290         check_valid(*result, "SNAP: result (before common-bits addition");
00291 #endif
00292 
00293 #if CBR_BEFORE_SNAPPING
00294         // Add common bits back in
00295         cbr.addCommonBits( result.get() );
00296         //result = fix_self_intersections(result, "SNAP: result (after common-bits addition)");
00297 
00298   check_valid(*result, "CBR: result (after common-bits addition)", true);
00299 
00300 #endif
00301 
00302         return result;
00303 }
00304 
00305 template <class BinOp>
00306 std::auto_ptr<Geometry>
00307 BinaryOp(const Geometry* g0, const Geometry *g1, BinOp _Op)
00308 {
00309         typedef std::auto_ptr<Geometry> GeomPtr;
00310 
00311         GeomPtr ret;
00312         geos::util::TopologyException origException;
00313 
00314 #ifdef USE_ORIGINAL_INPUT
00315         // Try with original input
00316         try
00317         {
00318 #if GEOS_DEBUG_BINARYOP
00319                 std::cerr << "Trying with original input." << std::endl;
00320 #endif
00321                 ret.reset(_Op(g0, g1));
00322                 return ret;
00323         }
00324         catch (const geos::util::TopologyException& ex)
00325         {
00326                 origException=ex;
00327 #if GEOS_DEBUG_BINARYOP
00328                 std::cerr << "Original exception: " << ex.what() << std::endl;
00329 #endif
00330         }
00331 
00332   check_valid(*g0, "Input geom 0", true, true);
00333   check_valid(*g1, "Input geom 1", true, true);
00334 
00335 #if GEOS_DEBUG_BINARYOP
00336   // Should we just give up here ?
00337   check_valid(*g0, "Input geom 0");
00338   check_valid(*g1, "Input geom 1");
00339 #endif
00340 
00341 #endif // USE_ORIGINAL_INPUT
00342 
00343 #ifdef USE_COMMONBITS_POLICY
00344         // Try removing common bits (possibly obsoleted by snapping below)
00345         //
00346         // NOTE: this policy was _later_ implemented 
00347         //       in JTS as EnhancedPrecisionOp
00348         // TODO: consider using the now-ported EnhancedPrecisionOp
00349         //       here too
00350         // 
00351         try
00352         {
00353                 GeomPtr rG0;
00354                 GeomPtr rG1;
00355                 precision::CommonBitsRemover cbr;
00356 
00357 #if GEOS_DEBUG_BINARYOP
00358                 std::cerr << "Trying with Common Bits Remover (CBR)" << std::endl;
00359 #endif
00360 
00361                 cbr.add(g0);
00362                 cbr.add(g1);
00363 
00364                 rG0.reset( cbr.removeCommonBits(g0->clone()) );
00365                 rG1.reset( cbr.removeCommonBits(g1->clone()) );
00366 
00367 #if GEOS_DEBUG_BINARYOP
00368                 check_valid(*rG0, "CBR: geom 0 (after common-bits removal)");
00369                 check_valid(*rG1, "CBR: geom 1 (after common-bits removal)");
00370 #endif
00371 
00372                 ret.reset( _Op(rG0.get(), rG1.get()) );
00373 
00374 #if GEOS_DEBUG_BINARYOP
00375                 check_valid(*ret, "CBR: result (before common-bits addition)");
00376 #endif
00377 
00378                 cbr.addCommonBits( ret.get() ); 
00379 
00380                 check_valid(*ret, "CBR: result (after common-bits addition)", true);
00381 
00382 #if GEOS_CHECK_COMMONBITS_VALIDITY
00383                 // check that result is a valid geometry after the
00384                 // reshift to orginal precision (see EnhancedPrecisionOp)
00385                 using operation::valid::IsValidOp;
00386                 using operation::valid::TopologyValidationError;
00387                 IsValidOp ivo(ret.get());
00388                 if ( ! ivo.isValid() )
00389                 {
00390                         TopologyValidationError* e = ivo.getValidationError();
00391                         throw geos::util::TopologyException(
00392                                 "Result of overlay became invalid "
00393                                 "after re-addin common bits of operand "
00394                                 "coordinates: " + e->toString(),
00395                                 e->getCoordinate());
00396                 }
00397 #endif // GEOS_CHECK_COMMONBITS_VALIDITY
00398 
00399                 return ret;
00400         }
00401         catch (const geos::util::TopologyException& ex)
00402         {
00403         ::geos::ignore_unused_variable_warning(ex);
00404 #if GEOS_DEBUG_BINARYOP
00405                 std::cerr << "CBR: " << ex.what() << std::endl;
00406 #endif
00407         }
00408 #endif
00409 
00410         // Try with snapping
00411         //
00412         // TODO: possible optimization would be reusing the
00413         //       already common-bit-removed inputs and just
00414         //       apply geometry snapping, whereas the current
00415         //       SnapOp function does both.
00416 // {
00417 #if USE_SNAPPING_POLICY
00418 
00419 #if GEOS_DEBUG_BINARYOP
00420         std::cerr << "Trying with snapping " << std::endl;
00421 #endif
00422 
00423         try {
00424                 ret = SnapOp(g0, g1, _Op);
00425 #if GEOS_DEBUG_BINARYOP
00426         std::cerr << "SnapOp succeeded" << std::endl;
00427 #endif
00428                 return ret;
00429                 
00430         }
00431         catch (const geos::util::TopologyException& ex)
00432         {
00433         ::geos::ignore_unused_variable_warning(ex);
00434 #if GEOS_DEBUG_BINARYOP
00435                 std::cerr << "SNAP: " << ex.what() << std::endl;
00436 #endif
00437         }
00438 
00439 #endif // USE_SNAPPING_POLICY }
00440 
00441 // {
00442 #if USE_PRECISION_REDUCTION_POLICY
00443 
00444 
00445         // Try reducing precision
00446         try
00447         {
00448                 long unsigned int g0scale = g0->getFactory()->getPrecisionModel()->getScale();
00449                 long unsigned int g1scale = g1->getFactory()->getPrecisionModel()->getScale();
00450 
00451 #if GEOS_DEBUG_BINARYOP
00452                 std::cerr << "Original input scales are: "
00453               << g0scale
00454               << " and " 
00455               << g1scale
00456               << std::endl;
00457 #endif
00458 
00459                 double maxScale = 1e16;
00460 
00461     // Don't use a scale bigger than the input one
00462     if ( g0scale && g0scale < maxScale ) maxScale = g0scale;
00463     if ( g1scale && g1scale < maxScale ) maxScale = g1scale;
00464 
00465 
00466                 for (double scale=maxScale; scale >= 1; scale /= 10)
00467                 {
00468                         PrecisionModel pm(scale);
00469                         GeometryFactory gf(&pm);
00470 #if GEOS_DEBUG_BINARYOP
00471                         std::cerr << "Trying with scale " << scale << std::endl;
00472 #endif
00473 
00474                         precision::GeometryPrecisionReducer reducer( gf );
00475                         GeomPtr rG0( reducer.reduce(*g0) );
00476                         GeomPtr rG1( reducer.reduce(*g1) );
00477 
00478                         try
00479                         {
00480                                 ret.reset( _Op(rG0.get(), rG1.get()) );
00481         // restore original precision (least precision between inputs)
00482         if ( g0->getFactory()->getPrecisionModel()->compareTo( g1->getFactory()->getPrecisionModel() ) < 0 ) {
00483           ret.reset( g0->getFactory()->createGeometry(ret.get()) );
00484         }
00485         else {
00486           ret.reset( g1->getFactory()->createGeometry(ret.get()) );
00487         }
00488                                 return ret;
00489                         }
00490                         catch (const geos::util::TopologyException& ex)
00491                         {
00492                                 if ( scale == 1 ) throw ex;
00493 #if GEOS_DEBUG_BINARYOP
00494                                 std::cerr << "Reduced with scale (" << scale << "): "
00495                                           << ex.what() << std::endl;
00496 #endif
00497                         }
00498 
00499                 }
00500 
00501         }
00502         catch (const geos::util::TopologyException& ex)
00503         {
00504 #if GEOS_DEBUG_BINARYOP
00505                 std::cerr << "Reduced: " << ex.what() << std::endl;
00506 #endif
00507         }
00508 
00509 #endif
00510 // USE_PRECISION_REDUCTION_POLICY }
00511 
00512 
00513 
00514 
00515 
00516 // {
00517 #if USE_TP_SIMPLIFY_POLICY 
00518 
00519         // Try simplifying
00520         try
00521         {
00522 
00523                 double maxTolerance = 0.04;
00524                 double minTolerance = 0.01;
00525                 double tolStep = 0.01;
00526 
00527                 for (double tol = minTolerance; tol <= maxTolerance; tol += tolStep)
00528                 {
00529 #if GEOS_DEBUG_BINARYOP
00530                         std::cerr << "Trying simplifying with tolerance " << tol << std::endl;
00531 #endif
00532 
00533                         GeomPtr rG0( simplify::TopologyPreservingSimplifier::simplify(g0, tol) );
00534                         GeomPtr rG1( simplify::TopologyPreservingSimplifier::simplify(g1, tol) );
00535 
00536                         try
00537                         {
00538                                 ret.reset( _Op(rG0.get(), rG1.get()) );
00539                                 return ret;
00540                         }
00541                         catch (const geos::util::TopologyException& ex)
00542                         {
00543                                 if ( tol >= maxTolerance ) throw ex;
00544 #if GEOS_DEBUG_BINARYOP
00545                                 std::cerr << "Simplified with tolerance (" << tol << "): "
00546                                           << ex.what() << std::endl;
00547 #endif
00548                         }
00549 
00550                 }
00551 
00552                 return ret;
00553 
00554         }
00555         catch (const geos::util::TopologyException& ex)
00556         {
00557 #if GEOS_DEBUG_BINARYOP
00558                 std::cerr << "Simplified: " << ex.what() << std::endl;
00559 #endif
00560         }
00561 
00562 #endif
00563 // USE_TP_SIMPLIFY_POLICY }
00564 
00565         throw origException;
00566 }
00567 
00568 
00569 } // namespace geos::geom
00570 } // namespace geos
00571 
00572 #endif // GEOS_GEOM_BINARYOP_H

Generated on 23 Sep 2013 for GEOS by  doxygen 1.4.7