BinaryOp.h

00001 /**********************************************************************
00002  * $Id: BinaryOp.h 3750 2013-01-17 12:04:00Z strk $
00003  *
00004  * GEOS - Geometry Engine Open Source
00005  * http://geos.refractions.net
00006  *
00007  * Copyright (C) 2013 Sandro Santilli <strk@keybit.net>
00008  * Copyright (C) 2006 Refractions Research Inc.
00009  *
00010  * This is free software; you can redistribute and/or modify it under
00011  * the terms of the GNU Lesser General Public Licence as published
00012  * by the Free Software Foundation. 
00013  * See the COPYING file for more information.
00014  *
00015  **********************************************************************
00016  *
00017  * Last port: ORIGINAL WORK
00018  *
00019  **********************************************************************
00020  *
00021  * This file provides a single templated function, taking two
00022  * const Geometry pointers, applying a binary operator to them
00023  * and returning a result Geometry in an auto_ptr<>.
00024  *
00025  * The binary operator is expected to take two const Geometry pointers
00026  * and return a newly allocated Geometry pointer, possibly throwing
00027  * a TopologyException to signal it couldn't succeed due to robustness
00028  * issues.
00029  *
00030  * This function will catch TopologyExceptions and try again with
00031  * slightly modified versions of the input. The following heuristic
00032  * is used:
00033  *
00034  *      - Try with original input.
00035  *      - Try removing common bits from input coordinate values
00036  *      - Try snaping input geometries to each other
00037  *      - Try snaping input coordinates to a increasing grid (size from 1/25 to 1)
00038  *      - Try simplifiying input with increasing tolerance (from 0.01 to 0.04)
00039  *
00040  * If none of the step succeeds the original exception is thrown.
00041  *
00042  * Note that you can skip Grid snapping, Geometry snapping and Simplify policies
00043  * by a compile-time define when building geos.
00044  * See USE_TP_SIMPLIFY_POLICY, USE_PRECISION_REDUCTION_POLICY and
00045  * USE_SNAPPING_POLICY macros below.
00046  *
00047  *
00048  **********************************************************************/
00049 
00050 #ifndef GEOS_GEOM_BINARYOP_H
00051 #define GEOS_GEOM_BINARYOP_H
00052 
00053 #include <geos/geom/Geometry.h>
00054 #include <geos/geom/GeometryCollection.h>
00055 #include <geos/geom/Polygon.h>
00056 #include <geos/geom/PrecisionModel.h>
00057 #include <geos/precision/CommonBitsRemover.h>
00058 #include <geos/precision/SimpleGeometryPrecisionReducer.h>
00059 
00060 #include <geos/operation/overlay/snap/GeometrySnapper.h>
00061 
00062 #include <geos/simplify/TopologyPreservingSimplifier.h>
00063 #include <geos/operation/valid/IsValidOp.h>
00064 #include <geos/operation/valid/TopologyValidationError.h>
00065 #include <geos/util/TopologyException.h>
00066 #include <geos/util.h>
00067 
00068 #include <memory> // for auto_ptr
00069 
00070 //#define GEOS_DEBUG_BINARYOP 1
00071 
00072 #ifdef GEOS_DEBUG_BINARYOP
00073 # include <iostream>
00074 # include <iomanip>
00075 # include <sstream>
00076 #endif
00077 
00078 
00079 /*
00080  * Always try original input first
00081  */
00082 #ifndef USE_ORIGINAL_INPUT
00083 # define USE_ORIGINAL_INPUT 1
00084 #endif
00085 
00086 /*
00087  * Define this to use PrecisionReduction policy
00088  * in an attempt at by-passing binary operation
00089  * robustness problems (handles TopologyExceptions)
00090  */
00091 #ifndef USE_PRECISION_REDUCTION_POLICY
00092 //# define USE_PRECISION_REDUCTION_POLICY 1
00093 #endif
00094 
00095 /*
00096  * Define this to use TopologyPreserving simplification policy
00097  * in an attempt at by-passing binary operation
00098  * robustness problems (handles TopologyExceptions)
00099  */
00100 #ifndef USE_TP_SIMPLIFY_POLICY 
00101 //# define USE_TP_SIMPLIFY_POLICY 1
00102 #endif
00103 
00104 /*
00105  * Use common bits removal policy.
00106  * If enabled, this would be tried /before/
00107  * Geometry snapping.
00108  */
00109 #ifndef USE_COMMONBITS_POLICY
00110 # define USE_COMMONBITS_POLICY 1
00111 #endif
00112 
00113 /*
00114  * Check validity of operation performed
00115  * by common bits removal policy.
00116  *
00117  * This matches what EnhancedPrecisionOp does in JTS
00118  * and fixes 5 tests of invalid outputs in our testsuite
00119  * (stmlf-cases-20061020-invalid-output.xml)
00120  * and breaks 1 test (robustness-invalid-output.xml) so much
00121  * to prevent a result.
00122  *
00123  */
00124 #define GEOS_CHECK_COMMONBITS_VALIDITY 1
00125 
00126 /*
00127  * Use snapping policy
00128  */
00129 #ifndef USE_SNAPPING_POLICY
00130 # define USE_SNAPPING_POLICY 1
00131 #endif
00132 
00133 namespace geos {
00134 namespace geom { // geos::geom
00135 
00136 inline bool
00137 check_valid(const Geometry& g, const std::string& label)
00138 {
00139         operation::valid::IsValidOp ivo(&g);
00140         if ( ! ivo.isValid() )
00141         {
00142 #ifdef GEOS_DEBUG_BINARYOP
00143                 using operation::valid::TopologyValidationError;
00144                 TopologyValidationError* err = ivo.getValidationError();
00145                 std::cerr << label << " is INVALID: "
00146                         << err->toString()
00147                         << " (" << std::setprecision(20)
00148                         << err->getCoordinate() << ")"
00149                         << std::endl;
00150 #else
00151         (void)label;
00152 #endif
00153                 return false;
00154         } 
00155         return true;
00156 }
00157 
00158 /*
00159  * Attempt to fix noding of multilines and
00160  * self-intersection of multipolygons 
00161  *
00162  * May return the input untouched.
00163  */
00164 inline std::auto_ptr<Geometry>
00165 fix_self_intersections(std::auto_ptr<Geometry> g, const std::string& label)
00166 {
00167 #ifdef GEOS_DEBUG_BINARYOP
00168         std::cerr << label << " fix_self_intersection (UnaryUnion)" << std::endl;
00169 #endif
00170 
00171   // Only multi-components can be fixed by UnaryUnion
00172   if ( ! dynamic_cast<const GeometryCollection*>(g.get()) ) return g;
00173 
00174   using operation::valid::IsValidOp;
00175 
00176   IsValidOp ivo(g.get());
00177 
00178   // Polygon is valid, nothing to do
00179   if ( ivo.isValid() ) return g;
00180 
00181   // Not all invalidities can be fixed by this code
00182 
00183   using operation::valid::TopologyValidationError;
00184   TopologyValidationError* err = ivo.getValidationError();
00185   switch ( err->getErrorType() ) {
00186     case TopologyValidationError::eRingSelfIntersection:
00187     case TopologyValidationError::eTooFewPoints: // collapsed lines
00188 #ifdef GEOS_DEBUG_BINARYOP
00189             std::cerr << label << " ATTEMPT_TO_FIX: " << err->getErrorType() << ": " << *g << std::endl;
00190 #endif
00191       g = g->Union();
00192 #ifdef GEOS_DEBUG_BINARYOP
00193             std::cerr << label << " ATTEMPT_TO_FIX succeeded.. " << std::endl;
00194 #endif
00195       return g;
00196     case TopologyValidationError::eSelfIntersection:
00197       // this one is within a single component, won't be fixed
00198     default:
00199 #ifdef GEOS_DEBUG_BINARYOP
00200             std::cerr << label << " invalidity is: " << err->getErrorType() << std::endl;
00201 #endif
00202       return g;
00203   }
00204 }
00205 
00206 
00212 template <class BinOp>
00213 std::auto_ptr<Geometry>
00214 SnapOp(const Geometry* g0, const Geometry *g1, BinOp _Op)
00215 {
00216         typedef std::auto_ptr<Geometry> GeomPtr;
00217 
00218 #define CBR_BEFORE_SNAPPING 1
00219 
00220         //using geos::precision::GeometrySnapper;
00221         using geos::operation::overlay::snap::GeometrySnapper;
00222 
00223         // Snap tolerance must be computed on the original
00224         // (not commonbits-removed) geoms
00225         double snapTolerance = GeometrySnapper::computeOverlaySnapTolerance(*g0, *g1);
00226 #if GEOS_DEBUG_BINARYOP
00227         std::cerr<< std::setprecision(20) << "Computed snap tolerance: "<<snapTolerance<<std::endl;
00228 #endif
00229 
00230 
00231 #if CBR_BEFORE_SNAPPING
00232         // Compute common bits
00233         geos::precision::CommonBitsRemover cbr;
00234         cbr.add(g0); cbr.add(g1);
00235 #if GEOS_DEBUG_BINARYOP
00236         std::cerr<<"Computed common bits: "<<cbr.getCommonCoordinate()<<std::endl;
00237 #endif
00238 
00239         // Now remove common bits
00240         GeomPtr rG0( cbr.removeCommonBits(g0->clone()) );
00241         GeomPtr rG1( cbr.removeCommonBits(g1->clone()) );
00242 
00243 #if GEOS_DEBUG_BINARYOP
00244         check_valid(*rG0, "CBR: removed-bits geom 0");
00245         check_valid(*rG1, "CBR: removed-bits geom 1");
00246 #endif
00247 
00248         const Geometry& operand0 = *rG0;
00249         const Geometry& operand1 = *rG1;
00250 #else // don't CBR before snapping
00251         const Geometry& operand0 = *g0;
00252         const Geometry& operand1 = *g1;
00253 #endif
00254 
00255 
00256         GeometrySnapper snapper0( operand0 );
00257         GeomPtr snapG0( snapper0.snapTo(operand1, snapTolerance) );
00258         snapG0 = fix_self_intersections(snapG0, "SNAP: snapped geom 0");
00259 
00260         // NOTE: second geom is snapped on the snapped first one
00261         GeometrySnapper snapper1( operand1 );
00262         GeomPtr snapG1( snapper1.snapTo(*snapG0, snapTolerance) );
00263         snapG1 = fix_self_intersections(snapG1, "SNAP: snapped geom 1");
00264 
00265         // Run the binary op
00266         GeomPtr result( _Op(snapG0.get(), snapG1.get()) );
00267 
00268 #if GEOS_DEBUG_BINARYOP
00269         check_valid(*result, "SNAP: result (before common-bits addition");
00270 #endif
00271 
00272 #if CBR_BEFORE_SNAPPING
00273         // Add common bits back in
00274         cbr.addCommonBits( result.get() );
00275         result = fix_self_intersections(result, "SNAP: result (after common-bits addition)");
00276 #endif
00277 
00278 #if GEOS_DEBUG_BINARYOP
00279         check_valid(*result, "SNAP: result (after common-bits addition");
00280 #endif
00281 
00282         return result;
00283 }
00284 
00285 template <class BinOp>
00286 std::auto_ptr<Geometry>
00287 BinaryOp(const Geometry* g0, const Geometry *g1, BinOp _Op)
00288 {
00289         typedef std::auto_ptr<Geometry> GeomPtr;
00290 
00291         GeomPtr ret;
00292         geos::util::TopologyException origException;
00293 
00294 #ifdef USE_ORIGINAL_INPUT
00295         // Try with original input
00296         try
00297         {
00298 #if GEOS_DEBUG_BINARYOP
00299                 std::cerr << "Trying with original input." << std::endl;
00300 #endif
00301                 ret.reset(_Op(g0, g1));
00302                 return ret;
00303         }
00304         catch (const geos::util::TopologyException& ex)
00305         {
00306                 origException=ex;
00307 #if GEOS_DEBUG_BINARYOP
00308                 std::cerr << "Original exception: " << ex.what() << std::endl;
00309 #endif
00310         }
00311 
00312 #if GEOS_DEBUG_BINARYOP
00313   // Should we just give up here ?
00314   check_valid(*g0, "Input geom 0");
00315   check_valid(*g1, "Input geom 1");
00316 #endif
00317 
00318 #endif // USE_ORIGINAL_INPUT
00319 
00320 
00321 #ifdef USE_COMMONBITS_POLICY
00322         // Try removing common bits (possibly obsoleted by snapping below)
00323         //
00324         // NOTE: this policy was _later_ implemented 
00325         //       in JTS as EnhancedPrecisionOp
00326         // TODO: consider using the now-ported EnhancedPrecisionOp
00327         //       here too
00328         // 
00329         try
00330         {
00331                 GeomPtr rG0;
00332                 GeomPtr rG1;
00333                 precision::CommonBitsRemover cbr;
00334 
00335 #if GEOS_DEBUG_BINARYOP
00336                 std::cerr << "Trying with Common Bits Remover (CBR)" << std::endl;
00337 #endif
00338 
00339                 cbr.add(g0);
00340                 cbr.add(g1);
00341 
00342                 rG0.reset( cbr.removeCommonBits(g0->clone()) );
00343                 rG1.reset( cbr.removeCommonBits(g1->clone()) );
00344 
00345 #if GEOS_DEBUG_BINARYOP
00346                 check_valid(*rG0, "CBR: geom 0 (after common-bits removal)");
00347                 check_valid(*rG1, "CBR: geom 1 (after common-bits removal)");
00348 #endif
00349 
00350                 ret.reset( _Op(rG0.get(), rG1.get()) );
00351 
00352 #if GEOS_DEBUG_BINARYOP
00353                 check_valid(*ret, "CBR: result (before common-bits addition)");
00354 #endif
00355 
00356                 cbr.addCommonBits( ret.get() ); 
00357 
00358 #if GEOS_DEBUG_BINARYOP
00359                 check_valid(*ret, "CBR: result (after common-bits addition)");
00360 #endif
00361 
00362                 // Common-bits removal policy could introduce self-intersections
00363                 ret = fix_self_intersections(ret, "CBR: result (after common-bits addition)");
00364 
00365 #if GEOS_DEBUG_BINARYOP
00366                 check_valid(*ret, "CBR: result (after common-bits addition and fix_self_intersections)");
00367 #endif
00368 
00369 #if GEOS_CHECK_COMMONBITS_VALIDITY
00370                 // check that result is a valid geometry after the
00371                 // reshift to orginal precision (see EnhancedPrecisionOp)
00372                 using operation::valid::IsValidOp;
00373                 using operation::valid::TopologyValidationError;
00374                 IsValidOp ivo(ret.get());
00375                 if ( ! ivo.isValid() )
00376                 {
00377                         TopologyValidationError* e = ivo.getValidationError();
00378                         throw geos::util::TopologyException(
00379                                 "Result of overlay became invalid "
00380                                 "after re-addin common bits of operand "
00381                                 "coordinates: " + e->toString(),
00382                                 e->getCoordinate());
00383                 }
00384 #endif // GEOS_CHECK_COMMONBITS_VALIDITY
00385 
00386                 return ret;
00387         }
00388         catch (const geos::util::TopologyException& ex)
00389         {
00390         ::geos::ignore_unused_variable_warning(ex);
00391 #if GEOS_DEBUG_BINARYOP
00392                 std::cerr << "CBR: " << ex.what() << std::endl;
00393 #endif
00394         }
00395 #endif
00396 
00397         // Try with snapping
00398         //
00399         // TODO: possible optimization would be reusing the
00400         //       already common-bit-removed inputs and just
00401         //       apply geometry snapping, whereas the current
00402         //       SnapOp function does both.
00403 // {
00404 #if USE_SNAPPING_POLICY
00405 
00406 #if GEOS_DEBUG_BINARYOP
00407         std::cerr << "Trying with snapping " << std::endl;
00408 #endif
00409 
00410         try {
00411                 ret = SnapOp(g0, g1, _Op);
00412 #if GEOS_DEBUG_BINARYOP
00413         std::cerr << "SnapOp succeeded" << std::endl;
00414 #endif
00415                 return ret;
00416                 
00417         }
00418         catch (const geos::util::TopologyException& ex)
00419         {
00420         ::geos::ignore_unused_variable_warning(ex);
00421 #if GEOS_DEBUG_BINARYOP
00422                 std::cerr << "SNAP: " << ex.what() << std::endl;
00423 #endif
00424         }
00425 
00426 #endif // USE_SNAPPING_POLICY }
00427 
00428 
00429 
00430 // {
00431 #if USE_PRECISION_REDUCTION_POLICY
00432 
00433 
00434         // Try reducing precision
00435         try
00436         {
00437                 int maxPrecision=25;
00438 
00439                 for (int precision=maxPrecision; precision; --precision)
00440                 {
00441                         std::auto_ptr<PrecisionModel> pm(new PrecisionModel(precision));
00442 #if GEOS_DEBUG_BINARYOP
00443                         std::cerr << "Trying with precision " << precision << std::endl;
00444 #endif
00445 
00446                         precision::SimpleGeometryPrecisionReducer reducer( pm.get() );
00447                         GeomPtr rG0( reducer.reduce(g0) );
00448                         GeomPtr rG1( reducer.reduce(g1) );
00449 
00450                         try
00451                         {
00452                                 ret.reset( _Op(rG0.get(), rG1.get()) );
00453                                 return ret;
00454                         }
00455                         catch (const geos::util::TopologyException& ex)
00456                         {
00457                                 if ( precision == 1 ) throw ex;
00458 #if GEOS_DEBUG_BINARYOP
00459                                 std::cerr << "Reduced with precision (" << precision << "): "
00460                                           << ex.what() << std::endl;
00461 #endif
00462                         }
00463 
00464                 }
00465 
00466         }
00467         catch (const geos::util::TopologyException& ex)
00468         {
00469 #if GEOS_DEBUG_BINARYOP
00470                 std::cerr << "Reduced: " << ex.what() << std::endl;
00471 #endif
00472         }
00473 
00474 #endif
00475 // USE_PRECISION_REDUCTION_POLICY }
00476 
00477 // {
00478 #if USE_TP_SIMPLIFY_POLICY 
00479 
00480         // Try simplifying
00481         try
00482         {
00483 
00484                 double maxTolerance = 0.04;
00485                 double minTolerance = 0.01;
00486                 double tolStep = 0.01;
00487 
00488                 for (double tol = minTolerance; tol <= maxTolerance; tol += tolStep)
00489                 {
00490 #if GEOS_DEBUG_BINARYOP
00491                         std::cerr << "Trying simplifying with tolerance " << tol << std::endl;
00492 #endif
00493 
00494                         GeomPtr rG0( simplify::TopologyPreservingSimplifier::simplify(g0, tol) );
00495                         GeomPtr rG1( simplify::TopologyPreservingSimplifier::simplify(g1, tol) );
00496 
00497                         try
00498                         {
00499                                 ret.reset( _Op(rG0.get(), rG1.get()) );
00500                                 return ret;
00501                         }
00502                         catch (const geos::util::TopologyException& ex)
00503                         {
00504                                 if ( tol >= maxTolerance ) throw ex;
00505 #if GEOS_DEBUG_BINARYOP
00506                                 std::cerr << "Simplified with tolerance (" << tol << "): "
00507                                           << ex.what() << std::endl;
00508 #endif
00509                         }
00510 
00511                 }
00512 
00513                 return ret;
00514 
00515         }
00516         catch (const geos::util::TopologyException& ex)
00517         {
00518 #if GEOS_DEBUG_BINARYOP
00519                 std::cerr << "Simplified: " << ex.what() << std::endl;
00520 #endif
00521         }
00522 
00523 #endif
00524 // USE_TP_SIMPLIFY_POLICY }
00525 
00526         throw origException;
00527 }
00528 
00529 
00530 } // namespace geos::geom
00531 } // namespace geos
00532 
00533 #endif // GEOS_GEOM_BINARYOP_H

Generated on Thu Mar 14 16:37:09 2013 for GEOS by  doxygen 1.4.7