libgeos++-dev 3.6.2-1build2 / usr / include / geos / geom / BinaryOp.h

/**********************************************************************
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.osgeo.org
 *
 * Copyright (C) 2013 Sandro Santilli <strk@keybit.net>
 * Copyright (C) 2006 Refractions Research Inc.
 *
 * This is free software; you can redistribute and/or modify it under
 * the terms of the GNU Lesser General Public Licence as published
 * by the Free Software Foundation. 
 * See the COPYING file for more information.
 *
 **********************************************************************
 *
 * Last port: ORIGINAL WORK
 *
 **********************************************************************
 *
 * This file provides a single templated function, taking two
 * const Geometry pointers, applying a binary operator to them
 * and returning a result Geometry in an auto_ptr<>.
 *
 * The binary operator is expected to take two const Geometry pointers
 * and return a newly allocated Geometry pointer, possibly throwing
 * a TopologyException to signal it couldn't succeed due to robustness
 * issues.
 *
 * This function will catch TopologyExceptions and try again with
 * slightly modified versions of the input. The following heuristic
 * is used:
 *
 *	- Try with original input.
 *	- Try removing common bits from input coordinate values
 *	- Try snaping input geometries to each other
 *	- Try snaping input coordinates to a increasing grid (size from 1/25 to 1)
 *	- Try simplifiying input with increasing tolerance (from 0.01 to 0.04)
 *
 * If none of the step succeeds the original exception is thrown.
 *
 * Note that you can skip Grid snapping, Geometry snapping and Simplify policies
 * by a compile-time define when building geos.
 * See USE_TP_SIMPLIFY_POLICY, USE_PRECISION_REDUCTION_POLICY and
 * USE_SNAPPING_POLICY macros below.
 *
 *
 **********************************************************************/

#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>
#include <geos/util.h>

#include <memory> // for auto_ptr

//#define GEOS_DEBUG_BINARYOP 1

#ifdef GEOS_DEBUG_BINARYOP
# include <iostream>
# include <iomanip>
# include <sstream>
#endif


/*
 * Always try original input first
 */
#ifndef USE_ORIGINAL_INPUT
# define USE_ORIGINAL_INPUT 1
#endif

/*
 * Define this to use PrecisionReduction policy
 * in an attempt at by-passing binary operation
 * robustness problems (handles TopologyExceptions)
 */
#ifndef USE_PRECISION_REDUCTION_POLICY
# define USE_PRECISION_REDUCTION_POLICY 1
#endif

/*
 * Define this to use TopologyPreserving simplification policy
 * in an attempt at by-passing binary operation
 * robustness problems (handles TopologyExceptions)
 */
#ifndef USE_TP_SIMPLIFY_POLICY 
//# define USE_TP_SIMPLIFY_POLICY 1
#endif

/*
 * Use common bits removal policy.
 * If enabled, this would be tried /before/
 * Geometry snapping.
 */
#ifndef USE_COMMONBITS_POLICY
# define USE_COMMONBITS_POLICY 1
#endif

/*
 * Check validity of operation performed
 * by common bits removal policy.
 *
 * This matches what EnhancedPrecisionOp does in JTS
 * and fixes 5 tests of invalid outputs in our testsuite
 * (stmlf-cases-20061020-invalid-output.xml)
 * and breaks 1 test (robustness-invalid-output.xml) so much
 * to prevent a result.
 *
 */
#define GEOS_CHECK_COMMONBITS_VALIDITY 1

/*
 * Use snapping policy
 */
#ifndef USE_SNAPPING_POLICY
# define USE_SNAPPING_POLICY 1
#endif

namespace geos {
namespace geom { // geos::geom

inline bool
check_valid(const Geometry& g, const std::string& label, bool doThrow=false, bool validOnly=false)
{
  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
      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;
}

/*
 * Attempt to fix noding of multilines and
 * self-intersection of multipolygons 
 *
 * May return the input untouched.
 */
inline std::auto_ptr<Geometry>
fix_self_intersections(std::auto_ptr<Geometry> g, const std::string& label)
{
  ::geos::ignore_unused_variable_warning(label);
#ifdef GEOS_DEBUG_BINARYOP
	std::cerr << label << " fix_self_intersection (UnaryUnion)" << std::endl;
#endif

  // Only multi-components can be fixed by UnaryUnion
  if ( ! dynamic_cast<const GeometryCollection*>(g.get()) ) return g;

  using operation::valid::IsValidOp;

  IsValidOp ivo(g.get());

  // Polygon is valid, nothing to do
  if ( ivo.isValid() ) return g;

  // Not all invalidities can be fixed by this code

  using operation::valid::TopologyValidationError;
  TopologyValidationError* err = ivo.getValidationError();
  switch ( err->getErrorType() ) {
    case TopologyValidationError::eRingSelfIntersection:
    case TopologyValidationError::eTooFewPoints: // collapsed lines
#ifdef GEOS_DEBUG_BINARYOP
	    std::cerr << label << " ATTEMPT_TO_FIX: " << err->getErrorType() << ": " << *g << std::endl;
#endif
      g = g->Union();
#ifdef GEOS_DEBUG_BINARYOP
	    std::cerr << label << " ATTEMPT_TO_FIX succeeded.. " << std::endl;
#endif
      return g;
    case TopologyValidationError::eSelfIntersection:
      // this one is within a single component, won't be fixed
    default:
#ifdef GEOS_DEBUG_BINARYOP
	    std::cerr << label << " invalidity is: " << err->getErrorType() << std::endl;
#endif
      return g;
  }
}


/// \brief
/// Apply a binary operation to the given geometries
/// after snapping them to each other after common-bits
/// removal.
///
template <class BinOp>
std::auto_ptr<Geometry>
SnapOp(const Geometry* g0, const Geometry *g1, BinOp _Op)
{
	typedef std::auto_ptr<Geometry> GeomPtr;

#define CBR_BEFORE_SNAPPING 1

	//using geos::precision::GeometrySnapper;
	using geos::operation::overlay::snap::GeometrySnapper;

	// Snap tolerance must be computed on the original
	// (not commonbits-removed) geoms
	double snapTolerance = GeometrySnapper::computeOverlaySnapTolerance(*g0, *g1);
#if GEOS_DEBUG_BINARYOP
	std::cerr<< std::setprecision(20) << "Computed snap tolerance: "<<snapTolerance<<std::endl;
#endif


#if CBR_BEFORE_SNAPPING
	// Compute common bits
	geos::precision::CommonBitsRemover cbr;
	cbr.add(g0); cbr.add(g1);
#if GEOS_DEBUG_BINARYOP
	std::cerr<<"Computed common bits: "<<cbr.getCommonCoordinate()<<std::endl;
#endif

	// Now remove common bits
	GeomPtr rG0( cbr.removeCommonBits(g0->clone()) );
	GeomPtr rG1( cbr.removeCommonBits(g1->clone()) );

#if GEOS_DEBUG_BINARYOP
	check_valid(*rG0, "CBR: removed-bits geom 0");
	check_valid(*rG1, "CBR: removed-bits geom 1");
#endif

	const Geometry& operand0 = *rG0;
	const Geometry& operand1 = *rG1;
#else // don't CBR before snapping
	const Geometry& operand0 = *g0;
	const Geometry& operand1 = *g1;
#endif


	GeometrySnapper snapper0( operand0 );
	GeomPtr snapG0( snapper0.snapTo(operand1, snapTolerance) );
	//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");

	// Run the binary op
	GeomPtr result( _Op(snapG0.get(), snapG1.get()) );

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

#if CBR_BEFORE_SNAPPING
	// Add common bits back in
	cbr.addCommonBits( result.get() );
	//result = fix_self_intersections(result, "SNAP: result (after common-bits addition)");

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

#endif

	return result;
}

template <class BinOp>
std::auto_ptr<Geometry>
BinaryOp(const Geometry* g0, const Geometry *g1, BinOp _Op)
{
	typedef std::auto_ptr<Geometry> GeomPtr;

	GeomPtr ret;
	geos::util::TopologyException origException;

#ifdef USE_ORIGINAL_INPUT
	// Try with original input
	try
	{
#if GEOS_DEBUG_BINARYOP
		std::cerr << "Trying with original input." << std::endl;
#endif
		ret.reset(_Op(g0, g1));
		return ret;
	}
	catch (const geos::util::TopologyException& ex)
	{
		origException=ex;
#if GEOS_DEBUG_BINARYOP
		std::cerr << "Original exception: " << ex.what() << std::endl;
#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");
  check_valid(*g1, "Input geom 1");
#endif

#endif // USE_ORIGINAL_INPUT

#ifdef USE_COMMONBITS_POLICY
	// Try removing common bits (possibly obsoleted by snapping below)
	//
	// NOTE: this policy was _later_ implemented 
	//       in JTS as EnhancedPrecisionOp
	// TODO: consider using the now-ported EnhancedPrecisionOp
	//       here too
	// 
	try
	{
		GeomPtr rG0;
		GeomPtr rG1;
		precision::CommonBitsRemover cbr;

#if GEOS_DEBUG_BINARYOP
		std::cerr << "Trying with Common Bits Remover (CBR)" << std::endl;
#endif

		cbr.add(g0);
		cbr.add(g1);

		rG0.reset( cbr.removeCommonBits(g0->clone()) );
		rG1.reset( cbr.removeCommonBits(g1->clone()) );

#if GEOS_DEBUG_BINARYOP
		check_valid(*rG0, "CBR: geom 0 (after common-bits removal)");
		check_valid(*rG1, "CBR: geom 1 (after common-bits removal)");
#endif

		ret.reset( _Op(rG0.get(), rG1.get()) );

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

		cbr.addCommonBits( ret.get() ); 

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

#if GEOS_CHECK_COMMONBITS_VALIDITY
		// check that result is a valid geometry after the
		// reshift to orginal precision (see EnhancedPrecisionOp)
		using operation::valid::IsValidOp;
		using operation::valid::TopologyValidationError;
		IsValidOp ivo(ret.get());
		if ( ! ivo.isValid() )
		{
			TopologyValidationError* e = ivo.getValidationError();
			throw geos::util::TopologyException(
				"Result of overlay became invalid "
				"after re-addin common bits of operand "
				"coordinates: " + e->toString(),
			        e->getCoordinate());
		}
#endif // GEOS_CHECK_COMMONBITS_VALIDITY

		return ret;
	}
	catch (const geos::util::TopologyException& ex)
	{
        ::geos::ignore_unused_variable_warning(ex);
#if GEOS_DEBUG_BINARYOP
		std::cerr << "CBR: " << ex.what() << std::endl;
#endif
	}
#endif

	// Try with snapping
	//
	// TODO: possible optimization would be reusing the
	//       already common-bit-removed inputs and just
	//       apply geometry snapping, whereas the current
	//       SnapOp function does both.
// {
#if USE_SNAPPING_POLICY

#if GEOS_DEBUG_BINARYOP
	std::cerr << "Trying with snapping " << std::endl;
#endif

	try {
		ret = SnapOp(g0, g1, _Op);
#if GEOS_DEBUG_BINARYOP
	std::cerr << "SnapOp succeeded" << std::endl;
#endif
		return ret;
		
	}
	catch (const geos::util::TopologyException& ex)
	{
        ::geos::ignore_unused_variable_warning(ex);
#if GEOS_DEBUG_BINARYOP
		std::cerr << "SNAP: " << ex.what() << std::endl;
#endif
	}

#endif // USE_SNAPPING_POLICY }

// {
#if USE_PRECISION_REDUCTION_POLICY


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

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

		double maxScale = 1e16;

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


		for (double scale=maxScale; scale >= 1; scale /= 10)
		{
			PrecisionModel pm(scale);
			GeometryFactory::unique_ptr gf = GeometryFactory::create(&pm);
#if GEOS_DEBUG_BINARYOP
			std::cerr << "Trying with scale " << scale << std::endl;
#endif

			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 ( scale == 1 ) throw ex;
#if GEOS_DEBUG_BINARYOP
				std::cerr << "Reduced with scale (" << scale << "): "
				          << ex.what() << std::endl;
#endif
			}

		}

	}
	catch (const geos::util::TopologyException& ex)
	{
#if GEOS_DEBUG_BINARYOP
		std::cerr << "Reduced: " << ex.what() << std::endl;
#endif
        ::geos::ignore_unused_variable_warning(ex);
	}

#endif
// USE_PRECISION_REDUCTION_POLICY }





// {
#if USE_TP_SIMPLIFY_POLICY 

	// Try simplifying
	try
	{

		double maxTolerance = 0.04;
		double minTolerance = 0.01;
		double tolStep = 0.01;

		for (double tol = minTolerance; tol <= maxTolerance; tol += tolStep)
		{
#if GEOS_DEBUG_BINARYOP
			std::cerr << "Trying simplifying with tolerance " << tol << std::endl;
#endif

			GeomPtr rG0( simplify::TopologyPreservingSimplifier::simplify(g0, tol) );
			GeomPtr rG1( simplify::TopologyPreservingSimplifier::simplify(g1, tol) );

			try
			{
				ret.reset( _Op(rG0.get(), rG1.get()) );
				return ret;
			}
			catch (const geos::util::TopologyException& ex)
			{
				if ( tol >= maxTolerance ) throw ex;
#if GEOS_DEBUG_BINARYOP
				std::cerr << "Simplified with tolerance (" << tol << "): "
				          << ex.what() << std::endl;
#endif
			}

		}

		return ret;

	}
	catch (const geos::util::TopologyException& ex)
	{
#if GEOS_DEBUG_BINARYOP
		std::cerr << "Simplified: " << ex.what() << std::endl;
#endif
	}

#endif
// USE_TP_SIMPLIFY_POLICY }

	throw origException;
}


} // namespace geos::geom
} // namespace geos

#endif // GEOS_GEOM_BINARYOP_H