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   ::geos::ignore_unused_variable_warning(label);
00189 #ifdef GEOS_DEBUG_BINARYOP
00190         std::cerr << label << " fix_self_intersection (UnaryUnion)" << std::endl;
00191 #endif
00192 
00193   // Only multi-components can be fixed by UnaryUnion
00194   if ( ! dynamic_cast<const GeometryCollection*>(g.get()) ) return g;
00195 
00196   using operation::valid::IsValidOp;
00197 
00198   IsValidOp ivo(g.get());
00199 
00200   // Polygon is valid, nothing to do
00201   if ( ivo.isValid() ) return g;
00202 
00203   // Not all invalidities can be fixed by this code
00204 
00205   using operation::valid::TopologyValidationError;
00206   TopologyValidationError* err = ivo.getValidationError();
00207   switch ( err->getErrorType() ) {
00208     case TopologyValidationError::eRingSelfIntersection:
00209     case TopologyValidationError::eTooFewPoints: // collapsed lines
00210 #ifdef GEOS_DEBUG_BINARYOP
00211             std::cerr << label << " ATTEMPT_TO_FIX: " << err->getErrorType() << ": " << *g << std::endl;
00212 #endif
00213       g = g->Union();
00214 #ifdef GEOS_DEBUG_BINARYOP
00215             std::cerr << label << " ATTEMPT_TO_FIX succeeded.. " << std::endl;
00216 #endif
00217       return g;
00218     case TopologyValidationError::eSelfIntersection:
00219       // this one is within a single component, won't be fixed
00220     default:
00221 #ifdef GEOS_DEBUG_BINARYOP
00222             std::cerr << label << " invalidity is: " << err->getErrorType() << std::endl;
00223 #endif
00224       return g;
00225   }
00226 }
00227 
00228 
00234 template <class BinOp>
00235 std::auto_ptr<Geometry>
00236 SnapOp(const Geometry* g0, const Geometry *g1, BinOp _Op)
00237 {
00238         typedef std::auto_ptr<Geometry> GeomPtr;
00239 
00240 #define CBR_BEFORE_SNAPPING 1
00241 
00242         //using geos::precision::GeometrySnapper;
00243         using geos::operation::overlay::snap::GeometrySnapper;
00244 
00245         // Snap tolerance must be computed on the original
00246         // (not commonbits-removed) geoms
00247         double snapTolerance = GeometrySnapper::computeOverlaySnapTolerance(*g0, *g1);
00248 #if GEOS_DEBUG_BINARYOP
00249         std::cerr<< std::setprecision(20) << "Computed snap tolerance: "<<snapTolerance<<std::endl;
00250 #endif
00251 
00252 
00253 #if CBR_BEFORE_SNAPPING
00254         // Compute common bits
00255         geos::precision::CommonBitsRemover cbr;
00256         cbr.add(g0); cbr.add(g1);
00257 #if GEOS_DEBUG_BINARYOP
00258         std::cerr<<"Computed common bits: "<<cbr.getCommonCoordinate()<<std::endl;
00259 #endif
00260 
00261         // Now remove common bits
00262         GeomPtr rG0( cbr.removeCommonBits(g0->clone()) );
00263         GeomPtr rG1( cbr.removeCommonBits(g1->clone()) );
00264 
00265 #if GEOS_DEBUG_BINARYOP
00266         check_valid(*rG0, "CBR: removed-bits geom 0");
00267         check_valid(*rG1, "CBR: removed-bits geom 1");
00268 #endif
00269 
00270         const Geometry& operand0 = *rG0;
00271         const Geometry& operand1 = *rG1;
00272 #else // don't CBR before snapping
00273         const Geometry& operand0 = *g0;
00274         const Geometry& operand1 = *g1;
00275 #endif
00276 
00277 
00278         GeometrySnapper snapper0( operand0 );
00279         GeomPtr snapG0( snapper0.snapTo(operand1, snapTolerance) );
00280         //snapG0 = fix_self_intersections(snapG0, "SNAP: snapped geom 0");
00281 
00282         // NOTE: second geom is snapped on the snapped first one
00283         GeometrySnapper snapper1( operand1 );
00284         GeomPtr snapG1( snapper1.snapTo(*snapG0, snapTolerance) );
00285         //snapG1 = fix_self_intersections(snapG1, "SNAP: snapped geom 1");
00286 
00287         // Run the binary op
00288         GeomPtr result( _Op(snapG0.get(), snapG1.get()) );
00289 
00290 #if GEOS_DEBUG_BINARYOP
00291         check_valid(*result, "SNAP: result (before common-bits addition");
00292 #endif
00293 
00294 #if CBR_BEFORE_SNAPPING
00295         // Add common bits back in
00296         cbr.addCommonBits( result.get() );
00297         //result = fix_self_intersections(result, "SNAP: result (after common-bits addition)");
00298 
00299   check_valid(*result, "CBR: result (after common-bits addition)", true);
00300 
00301 #endif
00302 
00303         return result;
00304 }
00305 
00306 template <class BinOp>
00307 std::auto_ptr<Geometry>
00308 BinaryOp(const Geometry* g0, const Geometry *g1, BinOp _Op)
00309 {
00310         typedef std::auto_ptr<Geometry> GeomPtr;
00311 
00312         GeomPtr ret;
00313         geos::util::TopologyException origException;
00314 
00315 #ifdef USE_ORIGINAL_INPUT
00316         // Try with original input
00317         try
00318         {
00319 #if GEOS_DEBUG_BINARYOP
00320                 std::cerr << "Trying with original input." << std::endl;
00321 #endif
00322                 ret.reset(_Op(g0, g1));
00323                 return ret;
00324         }
00325         catch (const geos::util::TopologyException& ex)
00326         {
00327                 origException=ex;
00328 #if GEOS_DEBUG_BINARYOP
00329                 std::cerr << "Original exception: " << ex.what() << std::endl;
00330 #endif
00331         }
00332 
00333   check_valid(*g0, "Input geom 0", true, true);
00334   check_valid(*g1, "Input geom 1", true, true);
00335 
00336 #if GEOS_DEBUG_BINARYOP
00337   // Should we just give up here ?
00338   check_valid(*g0, "Input geom 0");
00339   check_valid(*g1, "Input geom 1");
00340 #endif
00341 
00342 #endif // USE_ORIGINAL_INPUT
00343 
00344 #ifdef USE_COMMONBITS_POLICY
00345         // Try removing common bits (possibly obsoleted by snapping below)
00346         //
00347         // NOTE: this policy was _later_ implemented 
00348         //       in JTS as EnhancedPrecisionOp
00349         // TODO: consider using the now-ported EnhancedPrecisionOp
00350         //       here too
00351         // 
00352         try
00353         {
00354                 GeomPtr rG0;
00355                 GeomPtr rG1;
00356                 precision::CommonBitsRemover cbr;
00357 
00358 #if GEOS_DEBUG_BINARYOP
00359                 std::cerr << "Trying with Common Bits Remover (CBR)" << std::endl;
00360 #endif
00361 
00362                 cbr.add(g0);
00363                 cbr.add(g1);
00364 
00365                 rG0.reset( cbr.removeCommonBits(g0->clone()) );
00366                 rG1.reset( cbr.removeCommonBits(g1->clone()) );
00367 
00368 #if GEOS_DEBUG_BINARYOP
00369                 check_valid(*rG0, "CBR: geom 0 (after common-bits removal)");
00370                 check_valid(*rG1, "CBR: geom 1 (after common-bits removal)");
00371 #endif
00372 
00373                 ret.reset( _Op(rG0.get(), rG1.get()) );
00374 
00375 #if GEOS_DEBUG_BINARYOP
00376                 check_valid(*ret, "CBR: result (before common-bits addition)");
00377 #endif
00378 
00379                 cbr.addCommonBits( ret.get() ); 
00380 
00381                 check_valid(*ret, "CBR: result (after common-bits addition)", true);
00382 
00383 #if GEOS_CHECK_COMMONBITS_VALIDITY
00384                 // check that result is a valid geometry after the
00385                 // reshift to orginal precision (see EnhancedPrecisionOp)
00386                 using operation::valid::IsValidOp;
00387                 using operation::valid::TopologyValidationError;
00388                 IsValidOp ivo(ret.get());
00389                 if ( ! ivo.isValid() )
00390                 {
00391                         TopologyValidationError* e = ivo.getValidationError();
00392                         throw geos::util::TopologyException(
00393                                 "Result of overlay became invalid "
00394                                 "after re-addin common bits of operand "
00395                                 "coordinates: " + e->toString(),
00396                                 e->getCoordinate());
00397                 }
00398 #endif // GEOS_CHECK_COMMONBITS_VALIDITY
00399 
00400                 return ret;
00401         }
00402         catch (const geos::util::TopologyException& ex)
00403         {
00404         ::geos::ignore_unused_variable_warning(ex);
00405 #if GEOS_DEBUG_BINARYOP
00406                 std::cerr << "CBR: " << ex.what() << std::endl;
00407 #endif
00408         }
00409 #endif
00410 
00411         // Try with snapping
00412         //
00413         // TODO: possible optimization would be reusing the
00414         //       already common-bit-removed inputs and just
00415         //       apply geometry snapping, whereas the current
00416         //       SnapOp function does both.
00417 // {
00418 #if USE_SNAPPING_POLICY
00419 
00420 #if GEOS_DEBUG_BINARYOP
00421         std::cerr << "Trying with snapping " << std::endl;
00422 #endif
00423 
00424         try {
00425                 ret = SnapOp(g0, g1, _Op);
00426 #if GEOS_DEBUG_BINARYOP
00427         std::cerr << "SnapOp succeeded" << std::endl;
00428 #endif
00429                 return ret;
00430                 
00431         }
00432         catch (const geos::util::TopologyException& ex)
00433         {
00434         ::geos::ignore_unused_variable_warning(ex);
00435 #if GEOS_DEBUG_BINARYOP
00436                 std::cerr << "SNAP: " << ex.what() << std::endl;
00437 #endif
00438         }
00439 
00440 #endif // USE_SNAPPING_POLICY }
00441 
00442 // {
00443 #if USE_PRECISION_REDUCTION_POLICY
00444 
00445 
00446         // Try reducing precision
00447         try
00448         {
00449                 long unsigned int g0scale = 
00450             static_cast<long unsigned int>(g0->getFactory()->getPrecisionModel()->getScale());
00451                 long unsigned int g1scale = 
00452             static_cast<long unsigned int>(g1->getFactory()->getPrecisionModel()->getScale());
00453 
00454 #if GEOS_DEBUG_BINARYOP
00455                 std::cerr << "Original input scales are: "
00456               << g0scale
00457               << " and " 
00458               << g1scale
00459               << std::endl;
00460 #endif
00461 
00462                 double maxScale = 1e16;
00463 
00464     // Don't use a scale bigger than the input one
00465     if ( g0scale && g0scale < maxScale ) maxScale = g0scale;
00466     if ( g1scale && g1scale < maxScale ) maxScale = g1scale;
00467 
00468 
00469                 for (double scale=maxScale; scale >= 1; scale /= 10)
00470                 {
00471                         PrecisionModel pm(scale);
00472                         GeometryFactory gf(&pm);
00473 #if GEOS_DEBUG_BINARYOP
00474                         std::cerr << "Trying with scale " << scale << std::endl;
00475 #endif
00476 
00477                         precision::GeometryPrecisionReducer reducer( gf );
00478                         GeomPtr rG0( reducer.reduce(*g0) );
00479                         GeomPtr rG1( reducer.reduce(*g1) );
00480 
00481                         try
00482                         {
00483                                 ret.reset( _Op(rG0.get(), rG1.get()) );
00484         // restore original precision (least precision between inputs)
00485         if ( g0->getFactory()->getPrecisionModel()->compareTo( g1->getFactory()->getPrecisionModel() ) < 0 ) {
00486           ret.reset( g0->getFactory()->createGeometry(ret.get()) );
00487         }
00488         else {
00489           ret.reset( g1->getFactory()->createGeometry(ret.get()) );
00490         }
00491                                 return ret;
00492                         }
00493                         catch (const geos::util::TopologyException& ex)
00494                         {
00495                                 if ( scale == 1 ) throw ex;
00496 #if GEOS_DEBUG_BINARYOP
00497                                 std::cerr << "Reduced with scale (" << scale << "): "
00498                                           << ex.what() << std::endl;
00499 #endif
00500                         }
00501 
00502                 }
00503 
00504         }
00505         catch (const geos::util::TopologyException& ex)
00506         {
00507 #if GEOS_DEBUG_BINARYOP
00508                 std::cerr << "Reduced: " << ex.what() << std::endl;
00509 #endif
00510         ::geos::ignore_unused_variable_warning(ex);
00511         }
00512 
00513 #endif
00514 // USE_PRECISION_REDUCTION_POLICY }
00515 
00516 
00517 
00518 
00519 
00520 // {
00521 #if USE_TP_SIMPLIFY_POLICY 
00522 
00523         // Try simplifying
00524         try
00525         {
00526 
00527                 double maxTolerance = 0.04;
00528                 double minTolerance = 0.01;
00529                 double tolStep = 0.01;
00530 
00531                 for (double tol = minTolerance; tol <= maxTolerance; tol += tolStep)
00532                 {
00533 #if GEOS_DEBUG_BINARYOP
00534                         std::cerr << "Trying simplifying with tolerance " << tol << std::endl;
00535 #endif
00536 
00537                         GeomPtr rG0( simplify::TopologyPreservingSimplifier::simplify(g0, tol) );
00538                         GeomPtr rG1( simplify::TopologyPreservingSimplifier::simplify(g1, tol) );
00539 
00540                         try
00541                         {
00542                                 ret.reset( _Op(rG0.get(), rG1.get()) );
00543                                 return ret;
00544                         }
00545                         catch (const geos::util::TopologyException& ex)
00546                         {
00547                                 if ( tol >= maxTolerance ) throw ex;
00548 #if GEOS_DEBUG_BINARYOP
00549                                 std::cerr << "Simplified with tolerance (" << tol << "): "
00550                                           << ex.what() << std::endl;
00551 #endif
00552                         }
00553 
00554                 }
00555 
00556                 return ret;
00557 
00558         }
00559         catch (const geos::util::TopologyException& ex)
00560         {
00561 #if GEOS_DEBUG_BINARYOP
00562                 std::cerr << "Simplified: " << ex.what() << std::endl;
00563 #endif
00564         }
00565 
00566 #endif
00567 // USE_TP_SIMPLIFY_POLICY }
00568 
00569         throw origException;
00570 }
00571 
00572 
00573 } // namespace geos::geom
00574 } // namespace geos
00575 
00576 #endif // GEOS_GEOM_BINARYOP_H

Generated on 30 Dec 2015 for GEOS by  doxygen 1.4.7