Skip to content

Commit

Permalink
Improve overlay robustness
Browse files Browse the repository at this point in the history 
 - Validate CBR results before accepting them
 - Enable overlay node validator even for FIXED precision
 - Enable geometry-reduction policy
 - Bail out on exception from overlay if any input is invalid

Fixes bug libgeos#459

git-svn-id: http://svn.osgeo.org/geos/trunk@3853 5242fede-7e19-0410-aef8-94bd7d2200fb
  • Loading branch information
Sandro Santilli committed Jul 31, 2013
1 parent e8723b7 commit f38b4ef
Show file tree
Hide file tree
Showing 10 changed files with 184 additions and 91 deletions.
3 changes: 3 additions & 0 deletions NEWS
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ Changes in 3.4.0
- io::Writer::reserve() method
- CAPI: GEOSNearestPoints
- Add --cclibs, --static-clibs and --static-cclibs to geos-config (#497)
- Early bail out of overlay exception if input is invalid

- C++ API changes:
- New noding::GeometryNoder class
Expand All @@ -23,6 +24,8 @@ Changes in 3.4.0
- NodedSegmentString takes ownership of CoordinateSenuence now
- io::Writer's toString() returns by const ref, write() takes a const ref
- Bug fixes / improvements
- Improve Overlay robustness by reducing input precision on topology
exception and by refusing to accept unnoded output (#459)
- Improve Buffer robustness by reducing input precision on topology
exception (#605)
- Fixed Linear Referencing API to handle MultiLineStrings consistently
Expand Down
130 changes: 85 additions & 45 deletions include/geos/geom/BinaryOp.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,16 +49,21 @@
#ifndef GEOS_GEOM_BINARYOP_H
#define GEOS_GEOM_BINARYOP_H

#include <geos/algorithm/BoundaryNodeRule.h>
#include <geos/geom/Geometry.h>
#include <geos/geom/GeometryCollection.h>
#include <geos/geom/Polygon.h>
#include <geos/geom/Lineal.h>
#include <geos/geom/PrecisionModel.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/precision/CommonBitsRemover.h>
#include <geos/precision/SimpleGeometryPrecisionReducer.h>
#include <geos/precision/GeometryPrecisionReducer.h>

#include <geos/operation/overlay/snap/GeometrySnapper.h>

#include <geos/simplify/TopologyPreservingSimplifier.h>
#include <geos/operation/IsSimpleOp.h>
#include <geos/operation/valid/IsValidOp.h>
#include <geos/operation/valid/TopologyValidationError.h>
#include <geos/util/TopologyException.h>
Expand Down Expand Up @@ -88,7 +93,7 @@
* robustness problems (handles TopologyExceptions)
*/
#ifndef USE_PRECISION_REDUCTION_POLICY
//# define USE_PRECISION_REDUCTION_POLICY 1
# define USE_PRECISION_REDUCTION_POLICY 1
#endif

/*
Expand Down Expand Up @@ -133,24 +138,41 @@ namespace geos {
namespace geom { // geos::geom

inline bool
check_valid(const Geometry& g, const std::string& label)
check_valid(const Geometry& g, const std::string& label, bool doThrow=false, bool validOnly=false)
{
operation::valid::IsValidOp ivo(&g);
if ( ! ivo.isValid() )
{
if ( dynamic_cast<const Lineal*>(&g) ) {
if ( ! validOnly ) {
operation::IsSimpleOp sop(g, algorithm::BoundaryNodeRule::getBoundaryEndPoint());
if ( ! sop.isSimple() )
{
if ( doThrow ) {
throw geos::util::TopologyException(
label + " is not simple");
}
return false;
}
}
} else {
operation::valid::IsValidOp ivo(&g);
if ( ! ivo.isValid() )
{
using operation::valid::TopologyValidationError;
TopologyValidationError* err = ivo.getValidationError();
#ifdef GEOS_DEBUG_BINARYOP
using operation::valid::TopologyValidationError;
TopologyValidationError* err = ivo.getValidationError();
std::cerr << label << " is INVALID: "
<< err->toString()
<< " (" << std::setprecision(20)
<< err->getCoordinate() << ")"
<< std::endl;
#else
(void)label;
#endif
return false;
}
std::cerr << label << " is INVALID: "
<< err->toString()
<< " (" << std::setprecision(20)
<< err->getCoordinate() << ")"
<< std::endl;
#endif
if ( doThrow ) {
throw geos::util::TopologyException(
label + " is invalid: " + err->toString(),
err->getCoordinate());
}
return false;
}
}
return true;
}

Expand Down Expand Up @@ -254,12 +276,12 @@ SnapOp(const Geometry* g0, const Geometry *g1, BinOp _Op)

GeometrySnapper snapper0( operand0 );
GeomPtr snapG0( snapper0.snapTo(operand1, snapTolerance) );
snapG0 = fix_self_intersections(snapG0, "SNAP: snapped geom 0");
//snapG0 = fix_self_intersections(snapG0, "SNAP: snapped geom 0");

// NOTE: second geom is snapped on the snapped first one
GeometrySnapper snapper1( operand1 );
GeomPtr snapG1( snapper1.snapTo(*snapG0, snapTolerance) );
snapG1 = fix_self_intersections(snapG1, "SNAP: snapped geom 1");
//snapG1 = fix_self_intersections(snapG1, "SNAP: snapped geom 1");

// Run the binary op
GeomPtr result( _Op(snapG0.get(), snapG1.get()) );
Expand All @@ -271,11 +293,10 @@ SnapOp(const Geometry* g0, const Geometry *g1, BinOp _Op)
#if CBR_BEFORE_SNAPPING
// Add common bits back in
cbr.addCommonBits( result.get() );
result = fix_self_intersections(result, "SNAP: result (after common-bits addition)");
#endif
//result = fix_self_intersections(result, "SNAP: result (after common-bits addition)");

check_valid(*result, "CBR: result (after common-bits addition)", true);

#if GEOS_DEBUG_BINARYOP
check_valid(*result, "SNAP: result (after common-bits addition");
#endif

return result;
Expand Down Expand Up @@ -308,6 +329,9 @@ BinaryOp(const Geometry* g0, const Geometry *g1, BinOp _Op)
#endif
}

check_valid(*g0, "Input geom 0", true, true);
check_valid(*g1, "Input geom 1", true, true);

#if GEOS_DEBUG_BINARYOP
// Should we just give up here ?
check_valid(*g0, "Input geom 0");
Expand All @@ -316,7 +340,6 @@ BinaryOp(const Geometry* g0, const Geometry *g1, BinOp _Op)

#endif // USE_ORIGINAL_INPUT


#ifdef USE_COMMONBITS_POLICY
// Try removing common bits (possibly obsoleted by snapping below)
//
Expand Down Expand Up @@ -354,16 +377,7 @@ BinaryOp(const Geometry* g0, const Geometry *g1, BinOp _Op)

cbr.addCommonBits( ret.get() );

#if GEOS_DEBUG_BINARYOP
check_valid(*ret, "CBR: result (after common-bits addition)");
#endif

// Common-bits removal policy could introduce self-intersections
ret = fix_self_intersections(ret, "CBR: result (after common-bits addition)");

#if GEOS_DEBUG_BINARYOP
check_valid(*ret, "CBR: result (after common-bits addition and fix_self_intersections)");
#endif
check_valid(*ret, "CBR: result (after common-bits addition)", true);

#if GEOS_CHECK_COMMONBITS_VALIDITY
// check that result is a valid geometry after the
Expand Down Expand Up @@ -424,38 +438,60 @@ BinaryOp(const Geometry* g0, const Geometry *g1, BinOp _Op)

#endif // USE_SNAPPING_POLICY }



// {
#if USE_PRECISION_REDUCTION_POLICY


// Try reducing precision
try
{
int maxPrecision=25;
long unsigned int g0scale = g0->getFactory()->getPrecisionModel()->getScale();
long unsigned int g1scale = g1->getFactory()->getPrecisionModel()->getScale();

#if GEOS_DEBUG_BINARYOP
std::cerr << "Original input scales are: "
<< g0scale
<< " and "
<< g1scale
<< std::endl;
#endif

long unsigned int maxScale = 1e16;

for (int precision=maxPrecision; precision; --precision)
// Don't use a scale bigger than the input one
if ( g0scale && g0scale < maxScale ) maxScale = g0scale;
if ( g1scale && g1scale < maxScale ) maxScale = g1scale;


for (long unsigned int scale=maxScale; scale >= 1; scale /= 10)
{
std::auto_ptr<PrecisionModel> pm(new PrecisionModel(precision));
PrecisionModel pm(scale);
GeometryFactory gf(&pm);
#if GEOS_DEBUG_BINARYOP
std::cerr << "Trying with precision " << precision << std::endl;
std::cerr << "Trying with scale " << scale << std::endl;
#endif

precision::SimpleGeometryPrecisionReducer reducer( pm.get() );
GeomPtr rG0( reducer.reduce(g0) );
GeomPtr rG1( reducer.reduce(g1) );
precision::GeometryPrecisionReducer reducer( gf );
GeomPtr rG0( reducer.reduce(*g0) );
GeomPtr rG1( reducer.reduce(*g1) );

try
{
ret.reset( _Op(rG0.get(), rG1.get()) );
// restore original precision (least precision between inputs)
if ( g0->getFactory()->getPrecisionModel()->compareTo( g1->getFactory()->getPrecisionModel() ) < 0 ) {
ret.reset( g0->getFactory()->createGeometry(ret.get()) );
}
else {
ret.reset( g1->getFactory()->createGeometry(ret.get()) );
}
return ret;
}
catch (const geos::util::TopologyException& ex)
{
if ( precision == 1 ) throw ex;
if ( scale == 1 ) throw ex;
#if GEOS_DEBUG_BINARYOP
std::cerr << "Reduced with precision (" << precision << "): "
std::cerr << "Reduced with scale (" << scale << "): "
<< ex.what() << std::endl;
#endif
}
Expand All @@ -473,6 +509,10 @@ BinaryOp(const Geometry* g0, const Geometry *g1, BinOp _Op)
#endif
// USE_PRECISION_REDUCTION_POLICY }





// {
#if USE_TP_SIMPLIFY_POLICY

Expand Down
53 changes: 46 additions & 7 deletions include/geos/operation/overlay/OverlayOp.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
*
***********************************************************************
*
* Last port: operation/overlay/OverlayOp.java rev. 1.31 (JTS-1.10)
* Last port: operation/overlay/OverlayOp.java r567 (JTS-1.12+)
*
**********************************************************************/

Expand Down Expand Up @@ -59,7 +59,7 @@ namespace geos {
namespace operation { // geos::operation
namespace overlay { // geos::operation::overlay

/// Computes the overlay of two Geometry.
/// Computes the geometric overlay of two Geometry.
//
/// The overlay can be used to determine any
/// boolean combination of the geometries.
Expand All @@ -74,17 +74,42 @@ class GEOS_DLL OverlayOp: public GeometryGraphOperation {
/// the resultants of the overlay.
///
enum OpCode {
opINTERSECTION=1,
opUNION,
opDIFFERENCE,
opSYMDIFFERENCE
/// The code for the Intersection overlay operation.
opINTERSECTION = 1,
/// The code for the Union overlay operation.
opUNION = 2,
/// The code for the Difference overlay operation.
opDIFFERENCE = 3,
/// The code for the Symmetric Difference overlay operation.
opSYMDIFFERENCE = 4
};

/**
* Computes an overlay operation for the given geometry arguments.
*
* @param geom0 the first geometry argument
* @param geom1 the second geometry argument
* @param opCode the code for the desired overlay operation
* @return the result of the overlay operation
* @throws TopologyException if a robustness problem is encountered
*/
static geom::Geometry* overlayOp(const geom::Geometry *geom0,
const geom::Geometry *geom1,
OpCode opCode);
//throw(TopologyException *);

/**
* Tests whether a point with a given topological {@link Label}
* relative to two geometries is contained in
* the result of overlaying the geometries using
* a given overlay operation.
*
* The method handles arguments of {@link Location#NONE} correctly
*
* @param label the topological label of the point
* @param opCode the code for the overlay operation to test
* @return true if the label locations correspond to the overlayOpCode
*/
static bool isResultOfOp(const geomgraph::Label& label, OpCode opCode);

/// This method will handle arguments of Location.NULL correctly
Expand All @@ -102,9 +127,23 @@ class GEOS_DLL OverlayOp: public GeometryGraphOperation {

virtual ~OverlayOp(); // FIXME: virtual ?

geom::Geometry* getResultGeometry(OpCode funcCode);
/**
* Gets the result of the overlay for a given overlay operation.
*
* Note: this method can be called once only.
*
* @param overlayOpCode the overlay operation to perform
* @return the compute result geometry
* @throws TopologyException if a robustness problem is encountered
*/
geom::Geometry* getResultGeometry(OpCode overlayOpCode);
// throw(TopologyException *);

/**
* Gets the graph constructed to compute the overlay.
*
* @return the overlay graph
*/
geomgraph::PlanarGraph& getGraph() { return graph; }

/** \brief
Expand Down
Loading

0 comments on commit f38b4ef

Please sign in to comment.