Main Page | Namespace List | Class Hierarchy | Class List | File List | Namespace Members | Class Members | Related Pages

BinaryOp.h

00001 /**********************************************************************
00002  * $Id: BinaryOp.h 3207 2011-02-13 19:22:24Z mloskot $
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 #endif
00073 
00074 
00075 /*
00076  * Always try original input first
00077  */
00078 #ifndef USE_ORIGINAL_INPUT
00079 # define USE_ORIGINAL_INPUT 1
00080 #endif
00081 
00082 /*
00083  * Define this to use PrecisionReduction policy
00084  * in an attempt at by-passing binary operation
00085  * robustness problems (handles TopologyExceptions)
00086  */
00087 #ifndef USE_PRECISION_REDUCTION_POLICY
00088 //# define USE_PRECISION_REDUCTION_POLICY 1
00089 #endif
00090 
00091 /*
00092  * Define this to use TopologyPreserving simplification policy
00093  * in an attempt at by-passing binary operation
00094  * robustness problems (handles TopologyExceptions)
00095  */
00096 #ifndef USE_TP_SIMPLIFY_POLICY 
00097 //# define USE_TP_SIMPLIFY_POLICY 1
00098 #endif
00099 
00100 /*
00101  * Use common bits removal policy.
00102  * If enabled, this would be tried /before/
00103  * Geometry snapping.
00104  */
00105 #ifndef USE_COMMONBITS_POLICY
00106 # define USE_COMMONBITS_POLICY 1
00107 #endif
00108 
00109 /*
00110  * Check validity of operation performed
00111  * by common bits removal policy.
00112  *
00113  * This matches what EnhancedPrecisionOp does in JTS
00114  * and fixes 5 tests of invalid outputs in our testsuite
00115  * (stmlf-cases-20061020-invalid-output.xml)
00116  * and breaks 1 test (robustness-invalid-output.xml) so much
00117  * to prevent a result.
00118  *
00119  */
00120 #define GEOS_CHECK_COMMONBITS_VALIDITY 1
00121 
00122 /*
00123  * Use snapping policy
00124  */
00125 #ifndef USE_SNAPPING_POLICY
00126 # define USE_SNAPPING_POLICY 1
00127 #endif
00128 
00129 namespace geos {
00130 namespace geom { // geos::geom
00131 
00132 inline bool
00133 check_valid(const Geometry& g, const std::string& label)
00134 {
00135         operation::valid::IsValidOp ivo(&g);
00136         if ( ! ivo.isValid() )
00137         {
00138 #ifdef GEOS_DEBUG_BINARYOP
00139                 using operation::valid::TopologyValidationError;
00140                 TopologyValidationError* err = ivo.getValidationError();
00141                 std::cerr << label << " is INVALID: "
00142                         << err->toString()
00143                         << " (" << std::setprecision(20)
00144                         << err->getCoordinate() << ")"
00145                         << std::endl;
00146 #else
00147         (void)label;
00148 #endif
00149                 return false;
00150         } 
00151         return true;
00152 }
00153 
00159 template <class BinOp>
00160 std::auto_ptr<Geometry>
00161 SnapOp(const Geometry* g0, const Geometry *g1, BinOp _Op)
00162 {
00163         typedef std::auto_ptr<Geometry> GeomPtr;
00164 
00165 #define CBR_BEFORE_SNAPPING 1
00166 
00167         //using geos::precision::GeometrySnapper;
00168         using geos::operation::overlay::snap::GeometrySnapper;
00169 
00170         // Snap tolerance must be computed on the original
00171         // (not commonbits-removed) geoms
00172         double snapTolerance = GeometrySnapper::computeOverlaySnapTolerance(*g0, *g1);
00173 #if GEOS_DEBUG_BINARYOP
00174         std::cerr<< std::setprecision(20) << "Computed snap tolerance: "<<snapTolerance<<std::endl;
00175 #endif
00176 
00177 
00178 #if CBR_BEFORE_SNAPPING
00179         // Compute common bits
00180         geos::precision::CommonBitsRemover cbr;
00181         cbr.add(g0); cbr.add(g1);
00182 #if GEOS_DEBUG_BINARYOP
00183         std::cerr<<"Computed common bits: "<<cbr.getCommonCoordinate()<<std::endl;
00184 #endif
00185 
00186         // Now remove common bits
00187         GeomPtr rG0( cbr.removeCommonBits(g0->clone()) );
00188         GeomPtr rG1( cbr.removeCommonBits(g1->clone()) );
00189 
00190 #if GEOS_DEBUG_BINARYOP
00191         check_valid(*rG0, "CBR: removed-bits geom 0");
00192         check_valid(*rG1, "CBR: removed-bits geom 1");
00193 #endif
00194 
00195         const Geometry& operand0 = *rG0;
00196         const Geometry& operand1 = *rG1;
00197 #else // don't CBR before snapping
00198         const Geometry& operand0 = *g0
00199         const Geometry& operand1 = *g1
00200 #endif
00201 
00202 
00203         GeometrySnapper snapper0( operand0 );
00204         GeomPtr snapG0( snapper0.snapTo(operand1, snapTolerance) );
00205 
00206         // NOTE: second geom is snapped on the snapped first one
00207         GeometrySnapper snapper1( operand1 );
00208         GeomPtr snapG1( snapper1.snapTo(*snapG0, snapTolerance) );
00209 
00210 #if GEOS_DEBUG_BINARYOP
00211         check_valid(*snapG0, "SNAP: snapped geom 0");
00212         check_valid(*snapG1, "SNAP: snapped geom 1");
00213 #endif
00214 
00215         // Run the binary op
00216         GeomPtr result( _Op(snapG0.get(), snapG1.get()) );
00217 
00218 #if GEOS_DEBUG_BINARYOP
00219         check_valid(*result, "SNAP: result (before common-bits addition");
00220 #endif
00221 
00222 #if CBR_BEFORE_SNAPPING
00223         // Add common bits back in
00224         cbr.addCommonBits( result.get() );
00225 #endif
00226 
00227 #if GEOS_DEBUG_BINARYOP
00228         check_valid(*result, "SNAP: result (after common-bits addition");
00229 #endif
00230 
00231         return result;
00232 }
00233 
00234 template <class BinOp>
00235 std::auto_ptr<Geometry>
00236 BinaryOp(const Geometry* g0, const Geometry *g1, BinOp _Op)
00237 {
00238         typedef std::auto_ptr<Geometry> GeomPtr;
00239 
00240         GeomPtr ret;
00241         geos::util::TopologyException origException;
00242 
00243 #ifdef USE_ORIGINAL_INPUT
00244         // Try with original input
00245         try
00246         {
00247 #if GEOS_DEBUG_BINARYOP
00248                 std::cerr << "Trying with original input." << std::endl;
00249 #endif
00250                 ret.reset(_Op(g0, g1));
00251                 return ret;
00252         }
00253         catch (const geos::util::TopologyException& ex)
00254         {
00255                 origException=ex;
00256 #if GEOS_DEBUG_BINARYOP
00257                 std::cerr << "Original exception: " << ex.what() << std::endl;
00258 #endif
00259         }
00260 #endif // USE_ORIGINAL_INPUT
00261 
00262 
00263 #ifdef USE_COMMONBITS_POLICY
00264         // Try removing common bits (possibly obsoleted by snapping below)
00265         //
00266         // NOTE: this policy was _later_ implemented 
00267         //       in JTS as EnhancedPrecisionOp
00268         // TODO: consider using the now-ported EnhancedPrecisionOp
00269         //       here too
00270         // 
00271         try
00272         {
00273                 GeomPtr rG0;
00274                 GeomPtr rG1;
00275                 precision::CommonBitsRemover cbr;
00276 
00277 #if GEOS_DEBUG_BINARYOP
00278                 std::cerr << "Trying with Common Bits Remover (CBR)" << std::endl;
00279 #endif
00280 
00281                 cbr.add(g0);
00282                 cbr.add(g1);
00283 
00284                 rG0.reset( cbr.removeCommonBits(g0->clone()) );
00285                 rG1.reset( cbr.removeCommonBits(g1->clone()) );
00286 
00287 #if GEOS_DEBUG_BINARYOP
00288                 check_valid(*rG0, "CBR: geom 0 (after common-bits removal)");
00289                 check_valid(*rG1, "CBR: geom 1 (after common-bits removal)");
00290 #endif
00291 
00292                 ret.reset( _Op(rG0.get(), rG1.get()) );
00293 
00294 #if GEOS_DEBUG_BINARYOP
00295                 check_valid(*ret, "CBR: result (before common-bits addition)");
00296 #endif
00297 
00298                 cbr.addCommonBits( ret.get() ); 
00299 
00300 #if GEOS_DEBUG_BINARYOP
00301                 check_valid(*ret, "CBR: result (before common-bits addition)");
00302 #endif
00303 
00304 #if GEOS_CHECK_COMMONBITS_VALIDITY
00305                 // check that result is a valid geometry after the
00306                 // reshift to orginal precision (see EnhancedPrecisionOp)
00307                 using operation::valid::IsValidOp;
00308                 using operation::valid::TopologyValidationError;
00309                 IsValidOp ivo(ret.get());
00310                 if ( ! ivo.isValid() )
00311                 {
00312                         TopologyValidationError* e = ivo.getValidationError();
00313                         throw geos::util::TopologyException(
00314                                 "Result of overlay became invalid "
00315                                 "after re-addin common bits of operand "
00316                                 "coordinates: " + e->toString(),
00317                                 e->getCoordinate());
00318                 }
00319 #endif // GEOS_CHECK_COMMONBITS_VALIDITY
00320 
00321                 return ret;
00322         }
00323         catch (const geos::util::TopologyException& ex)
00324         {
00325         ::geos::ignore_unused_variable_warning(ex);
00326 #if GEOS_DEBUG_BINARYOP
00327                 std::cerr << "CBR: " << ex.what() << std::endl;
00328 #endif
00329         }
00330 #endif
00331 
00332         // Try with snapping
00333         //
00334         // TODO: possible optimization would be reusing the
00335         //       already common-bit-removed inputs and just
00336         //       apply geometry snapping, whereas the current
00337         //       SnapOp function does both.
00338 // {
00339 #if USE_SNAPPING_POLICY
00340 
00341 #if GEOS_DEBUG_BINARYOP
00342         std::cerr << "Trying with snapping " << std::endl;
00343 #endif
00344 
00345         try {
00346                 ret = SnapOp(g0, g1, _Op);
00347 #if GEOS_DEBUG_BINARYOP
00348         std::cerr << "SnapOp succeeded" << std::endl;
00349 #endif
00350                 return ret;
00351                 
00352         }
00353         catch (const geos::util::TopologyException& ex)
00354         {
00355         ::geos::ignore_unused_variable_warning(ex);
00356 #if GEOS_DEBUG_BINARYOP
00357                 std::cerr << "SNAP: " << ex.what() << std::endl;
00358 #endif
00359         }
00360 
00361 #endif // USE_SNAPPING_POLICY }
00362 
00363 
00364 
00365 // {
00366 #if USE_PRECISION_REDUCTION_POLICY
00367 
00368 
00369         // Try reducing precision
00370         try
00371         {
00372                 int maxPrecision=25;
00373 
00374                 for (int precision=maxPrecision; precision; --precision)
00375                 {
00376                         std::auto_ptr<PrecisionModel> pm(new PrecisionModel(precision));
00377 #if GEOS_DEBUG_BINARYOP
00378                         std::cerr << "Trying with precision " << precision << std::endl;
00379 #endif
00380 
00381                         precision::SimpleGeometryPrecisionReducer reducer( pm.get() );
00382                         GeomPtr rG0( reducer.reduce(g0) );
00383                         GeomPtr rG1( reducer.reduce(g1) );
00384 
00385                         try
00386                         {
00387                                 ret.reset( _Op(rG0.get(), rG1.get()) );
00388                                 return ret;
00389                         }
00390                         catch (const geos::util::TopologyException& ex)
00391                         {
00392                                 if ( precision == 1 ) throw ex;
00393 #if GEOS_DEBUG_BINARYOP
00394                                 std::cerr << "Reduced with precision (" << precision << "): "
00395                                           << ex.what() << std::endl;
00396 #endif
00397                         }
00398 
00399                 }
00400 
00401         }
00402         catch (const geos::util::TopologyException& ex)
00403         {
00404 #if GEOS_DEBUG_BINARYOP
00405                 std::cerr << "Reduced: " << ex.what() << std::endl;
00406 #endif
00407         }
00408 
00409 #endif
00410 // USE_PRECISION_REDUCTION_POLICY }
00411 
00412 // {
00413 #if USE_TP_SIMPLIFY_POLICY 
00414 
00415         // Try simplifying
00416         try
00417         {
00418 
00419                 double maxTolerance = 0.04;
00420                 double minTolerance = 0.01;
00421                 double tolStep = 0.01;
00422 
00423                 for (double tol = minTolerance; tol <= maxTolerance; tol += tolStep)
00424                 {
00425 #if GEOS_DEBUG_BINARYOP
00426                         std::cerr << "Trying simplifying with tolerance " << tol << std::endl;
00427 #endif
00428 
00429                         GeomPtr rG0( simplify::TopologyPreservingSimplifier::simplify(g0, tol) );
00430                         GeomPtr rG1( simplify::TopologyPreservingSimplifier::simplify(g1, tol) );
00431 
00432                         try
00433                         {
00434                                 ret.reset( _Op(rG0.get(), rG1.get()) );
00435                                 return ret;
00436                         }
00437                         catch (const geos::util::TopologyException& ex)
00438                         {
00439                                 if ( tol >= maxTolerance ) throw ex;
00440 #if GEOS_DEBUG_BINARYOP
00441                                 std::cerr << "Simplified with tolerance (" << tol << "): "
00442                                           << ex.what() << std::endl;
00443 #endif
00444                         }
00445 
00446                 }
00447 
00448                 return ret;
00449 
00450         }
00451         catch (const geos::util::TopologyException& ex)
00452         {
00453 #if GEOS_DEBUG_BINARYOP
00454                 std::cerr << "Simplified: " << ex.what() << std::endl;
00455 #endif
00456         }
00457 
00458 #endif
00459 // USE_TP_SIMPLIFY_POLICY }
00460 
00461         throw origException;
00462 }
00463 
00464 
00465 } // namespace geos::geom
00466 } // namespace geos
00467 
00468 #endif // GEOS_GEOM_BINARYOP_H

Generated on Sun Aug 21 23:21:01 2011 for GEOS by  doxygen 1.3.9.1