/* private */
auto_ptr<Geometry>
GeometryPrecisionReducer::fixPolygonalTopology(const geom::Geometry& geom )
{
  /**
   * If precision model was *not* changed, need to flip
   * geometry to targetPM, buffer in that model, then flip back
   */
  auto_ptr<Geometry> tmp;
  auto_ptr<GeometryFactory> tmpFactory;

  const Geometry* geomToBuffer = &geom;

  if ( ! newFactory ) {
    tmpFactory = createFactory(*geom.getFactory(), targetPM);
    tmp.reset( tmpFactory->createGeometry(&geom) );
    geomToBuffer = tmp.get();
  }

  auto_ptr<Geometry> bufGeom ( geomToBuffer->buffer(0) );

  if ( ! newFactory ) {
    // a slick way to copy the geometry with the original precision factory
    bufGeom.reset( geom.getFactory()->createGeometry(bufGeom.get()) );
  }

  return bufGeom;
}

static int
checkSingleSidedBufferSuccess(geom::Geometry& gRes,
        geom::Geometry& gRealRes, double dist)
{
    int success = 1;
    do
    {

        if ( gRes.getGeometryTypeId() != gRealRes.getGeometryTypeId() )
        {
            std::cerr << "Expected result is of type "
                << gRes.getGeometryType()
                << "; obtained result is of type "
                << gRealRes.getGeometryType()
                << std::endl;
            success=0;
            break;
        }

        geos::xmltester::SingleSidedBufferResultMatcher matcher;
        if ( ! matcher.isBufferResultMatch(gRealRes,
                           gRes,
                           dist) )
        {
std::cerr << "SingleSidedBufferResultMatcher FAILED" << std::endl;
            success=0;
            break;
        }

    }
    while (0);

    return success;
}

bool
BufferResultMatcher::isBoundaryHausdorffDistanceInTolerance(
	const geom::Geometry& actualBuffer,
	const geom::Geometry& expectedBuffer, double distance)
{
	typedef std::auto_ptr<geom::Geometry> GeomPtr;

	using geos::algorithm::distance::DiscreteHausdorffDistance;

	GeomPtr actualBdy ( actualBuffer.getBoundary() );
	GeomPtr expectedBdy ( expectedBuffer.getBoundary() );

	DiscreteHausdorffDistance haus(*actualBdy, *expectedBdy);
	haus.setDensifyFraction(0.25);

	double maxDistanceFound = haus.orientedDistance();

	double expectedDistanceTol = fabs(distance) / MAX_HAUSDORFF_DISTANCE_FACTOR;
	if (expectedDistanceTol < MIN_DISTANCE_TOLERANCE)
	{
		expectedDistanceTol = MIN_DISTANCE_TOLERANCE;
	}

	if (maxDistanceFound > expectedDistanceTol)
	{
std::cerr << "maxDistanceFound: " << maxDistanceFound << " tolerated " << expectedDistanceTol << std::endl;
		return false;
	}

	return true;
}

bool
BufferResultMatcher::isBufferResultMatch(const geom::Geometry& actualBuffer,
	                         const geom::Geometry& expectedBuffer,
	                         double distance)
{
	if (actualBuffer.isEmpty() && expectedBuffer.isEmpty())
		return true;

	/**
	 * MD - need some more checks here - symDiffArea won't catch
	 * very small holes ("tears")
	 * near the edge of computed buffers (which can happen
	 * in current version of JTS (1.8)).
	 * This can probably be handled by testing
	 * that every point of the actual buffer is at least a certain
	 * distance away from the geometry boundary.
	 */
	if (! isSymDiffAreaInTolerance(actualBuffer, expectedBuffer))
	{
std::cerr << "isSymDiffAreaInTolerance failed" << std::endl;
		return false;
	}

	if (! isBoundaryHausdorffDistanceInTolerance(actualBuffer,
	           expectedBuffer, distance))
	{
std::cerr << "isBoundaryHasudorffDistanceInTolerance failed" << std::endl;
		return false;
	}

	return true;
}

bool
AbstractPreparedPolygonContains::isSingleShell( const geom::Geometry & geom)
{
    // handles single-element MultiPolygons, as well as Polygons
    if ( geom.getNumGeometries() != 1)
        return false;

    const geom::Polygon * poly = (Polygon *)geom.getGeometryN( 0);

    int numHoles = poly->getNumInteriorRing();
    return (numHoles == 0);
}

/* private */
void
SnapOverlayOp::removeCommonBits(const geom::Geometry& geom0,
                                const geom::Geometry& geom1,
				geom::GeomPtrPair& remGeom)
{
	cbr.reset(new precision::CommonBitsRemover());
	cbr->add(&geom0);
	cbr->add(&geom1);

	remGeom.first.reset( cbr->removeCommonBits(geom0.clone()) );
	remGeom.second.reset( cbr->removeCommonBits(geom1.clone()) );
}

bool 
AbstractPreparedPolygonContains::isSingleShell(const geom::Geometry& geom)
{
	// handles single-element MultiPolygons, as well as Polygons
	if (geom.getNumGeometries() != 1)
    {
		return false;
    }
	
    const geom::Geometry* g = geom.getGeometryN(0);
	const geom::Polygon* poly = static_cast<const Polygon*>(g);

    std::size_t numHoles = poly->getNumInteriorRing();
	return (0 == numHoles);
}

bool
RectangleIntersects::intersects(const geom::Geometry& geom)
{
	if (!rectEnv.intersects(geom.getEnvelopeInternal()))
		return false;

	// test envelope relationships
	EnvelopeIntersectsVisitor visitor(rectEnv);
	visitor.applyTo(geom);
	if (visitor.intersects()) return true;

	// test if any rectangle corner is contained in the target
	ContainsPointVisitor ecpVisitor(rectangle);
	ecpVisitor.applyTo(geom);
	if (ecpVisitor.containsPoint())
		return true;

	// test if any lines intersect
	LineIntersectsVisitor liVisitor(rectangle);
	liVisitor.applyTo(geom);
	if (liVisitor.intersects())
		return true;

	return false;
}

bool
BufferResultMatcher::isSymDiffAreaInTolerance(
	const geom::Geometry& actualBuffer,
	const geom::Geometry& expectedBuffer)
{
	typedef std::auto_ptr<geom::Geometry> GeomPtr;

	using namespace operation::overlay;
	using geos::geom::BinaryOp;

	double area = expectedBuffer.getArea();
	GeomPtr diff = BinaryOp(&actualBuffer, &expectedBuffer,
	                   overlayOp(OverlayOp::opSYMDIFFERENCE));

	double areaDiff = diff->getArea();
	
	// can't get closer than difference area = 0 ! 
	// This also handles case when symDiff is empty
	if (areaDiff <= 0.0) return true;

	if (area <= 0) return false; 
	double frac = areaDiff / area;

	bool ret = frac < MAX_RELATIVE_AREA_DIFFERENCE;
	if ( ! ret )
	{
std::cerr << "symDiffArea frac: " << frac << " tolerated " << MAX_RELATIVE_AREA_DIFFERENCE << std::endl;
	}
	return ret;
}

	void visit(const geom::Geometry &geom)
	{
		using GEOMETRY::algorithm::locate::SimplePointInAreaLocator;

		const geom::Polygon *poly;

		if ( !(poly=dynamic_cast<const geom::Polygon *>(&geom)) ) {
			return;
		}

		const geom::Envelope &elementEnv = *(geom.getEnvelopeInternal());

		if ( !rectEnv.intersects(elementEnv) ) {
			return;
		}

		// test each corner of rectangle for inclusion
		for (int i=0; i<4; i++)
		{

			const geom::Coordinate &rectPt=rectSeq.getAt(i);

			if ( !elementEnv.contains(rectPt) ) {
				continue;
			}

			// check rect point in poly (rect is known not to
			// touch polygon at this point)
			if ( SimplePointInAreaLocator::containsPointInPolygon(rectPt, poly) )
			{
				containsPointVar=true;
				return;
			}
		}
	}

/* private static */
void
GeometryNoder::extractSegmentStrings(const geom::Geometry& g,
                                     SegmentString::NonConstVect& to)
{
  SegmentStringExtractor ex(to);
	g.apply_ro(&ex);
}

/* public static */
void
DistanceToPoint::computeDistance(const geom::Geometry& geom,
                                          const geom::Coordinate& pt,
                                          PointPairDistance& ptDist)
{
	if ( const LineString* ls = dynamic_cast<const LineString*>(&geom) )
	{
		computeDistance(*ls, pt, ptDist);
	}
	else if ( const Polygon* pl = dynamic_cast<const Polygon*>(&geom) )
	{
		computeDistance(*pl, pt, ptDist);
	}
	else if ( const GeometryCollection* gc = dynamic_cast<const GeometryCollection*>(&geom) )
	{
		for (size_t i = 0; i < gc->getNumGeometries(); i++)
		{
			const Geometry* g = gc->getGeometryN(i);
			computeDistance(*g, pt, ptDist);
		}
	}
	else
	{
		// assume geom is Point
		ptDist.setMinimum(*(geom.getCoordinate()), pt);
	}
}

/*public static*/
double
GeometrySnapper::computeSizeBasedSnapTolerance(const geom::Geometry& g)
{
	const Envelope* env = g.getEnvelopeInternal();
	double minDimension = (std::min)(env->getHeight(), env->getWidth());
	double snapTol = minDimension * snapPrecisionFactor;
	return snapTol;
}

/* Could be an XMLTester class private but oh well.. */
static int
checkBufferSuccess(geom::Geometry const& gRes, geom::Geometry const& gRealRes, double dist)
{
    int success = 1;
    do
    {

        if ( gRes.getGeometryTypeId() != gRealRes.getGeometryTypeId() )
        {
            std::cerr << "Expected result is of type "
                << gRes.getGeometryType()
                << "; obtained result is of type "
                << gRealRes.getGeometryType()
                << std::endl;
            success=0;
            break;
        }

        // Is a buffer always an area ?
        if ( gRes.getDimension() != 2 )
        {
            std::cerr << "Don't know how to validate "
                << "result of buffer operation "
                << "when expected result is not an "
                << "areal type."
                << std::endl;
        }
        

        geos::xmltester::BufferResultMatcher matcher;
        if ( ! matcher.isBufferResultMatch(gRealRes,
                           gRes,
                           dist) )
        {
std::cerr << "BufferResultMatcher FAILED" << std::endl;
            success=0;
            break;
        }

    }
    while (0);

    return success;
}

/* public */
SharedPathsOp::SharedPathsOp(
    const geom::Geometry& g1, const geom::Geometry& g2)
  :
  _g1(g1),
  _g2(g2),
  _gf(*g1.getFactory())
{
  checkLinealInput(_g1);
  checkLinealInput(_g2);
}

bool
SingleSidedBufferResultMatcher::isBufferResultMatch(const geom::Geometry& actualBuffer,
	                         const geom::Geometry& expectedBuffer,
	                         double distance)
{
	bool aEmpty = actualBuffer.isEmpty();
	bool eEmpty = expectedBuffer.isEmpty();

	// Both empty succeeds
	if (aEmpty && eEmpty) return true;

	// One empty and not the other is a failure
	if (aEmpty || eEmpty)
	{
std::cerr << "isBufferResultMatch failed (one empty and one not)" << std::endl;
		return false;
	}


   // NOTE: we need to test hausdorff distance in both directions

	if (! isBoundaryHausdorffDistanceInTolerance(actualBuffer,
	           expectedBuffer, distance))
	{
std::cerr << "isBoundaryHasudorffDistanceInTolerance failed (actual,expected)" << std::endl;
		return false;
	}

	if (! isBoundaryHausdorffDistanceInTolerance(expectedBuffer,
	           actualBuffer, distance))
	{
std::cerr << "isBoundaryHasudorffDistanceInTolerance failed (expected,actual)" << std::endl;
		return false;
	}

	return true;
}

/*public static*/
double
GeometrySnapper::computeOverlaySnapTolerance(const geom::Geometry& g)
{
	double snapTolerance = computeSizeBasedSnapTolerance(g);

	/**
	 * Overlay is carried out in the precision model
	 * of the two inputs.
	 * If this precision model is of type FIXED, then the snap tolerance
	 * must reflect the precision grid size.
	 * Specifically, the snap tolerance should be at least
	 * the distance from a corner of a precision grid cell
	 * to the centre point of the cell.
	 */
	assert(g.getPrecisionModel());
	const PrecisionModel& pm = *(g.getPrecisionModel());
	if (pm.getType() == PrecisionModel::FIXED)
	{
		double fixedSnapTol = (1 / pm.getScale()) * 2 / 1.415;
		if ( fixedSnapTol > snapTolerance )
			snapTolerance = fixedSnapTol;
	}
	return snapTolerance;
}