GeometryGraph.cpp

/**********************************************************************
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.osgeo.org
 *
 * Copyright (C) 2011 Sandro Santilli <strk@kbt.io>
 * Copyright (C) 2005-2006 Refractions Research Inc.
 * Copyright (C) 2001-2002 Vivid Solutions 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: geomgraph/GeometryGraph.java r428 (JTS-1.12+)
 *
 **********************************************************************/

#include <geos/algorithm/CGAlgorithms.h>
#include <geos/algorithm/BoundaryNodeRule.h>

#include <geos/util/UnsupportedOperationException.h>

#include <geos/geomgraph/GeometryGraph.h>
#include <geos/geomgraph/Node.h>
#include <geos/geomgraph/Edge.h>
#include <geos/geomgraph/Label.h>
#include <geos/geomgraph/Position.h>

#include <geos/geomgraph/index/SimpleMCSweepLineIntersector.h>
#include <geos/geomgraph/index/SegmentIntersector.h>
#include <geos/geomgraph/index/EdgeSetIntersector.h>

#include <geos/geom/CoordinateArraySequence.h>
#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/Location.h>
#include <geos/geom/Point.h>
#include <geos/geom/Envelope.h>
#include <geos/geom/LinearRing.h>
#include <geos/geom/LineString.h>
#include <geos/geom/Polygon.h>
#include <geos/geom/MultiPoint.h>
#include <geos/geom/MultiLineString.h>
#include <geos/geom/MultiPolygon.h>
#include <geos/geom/GeometryCollection.h>
#include <geos/util/Interrupt.h>

#include <geos/inline.h>

#include <vector>
#include <memory> // unique_ptr
#include <cassert>
#include <typeinfo>

#ifndef GEOS_DEBUG
#define GEOS_DEBUG 0
#endif

#ifndef GEOS_INLINE
# include "geos/geomgraph/GeometryGraph.inl"
#endif

using namespace std;
using namespace geos::geomgraph::index;
using namespace geos::algorithm;
using namespace geos::geom;

namespace geos {
namespace geomgraph { // geos.geomgraph

/*
 * This method implements the Boundary Determination Rule
 * for determining whether
 * a component (node or edge) that appears multiple times in elements
 * of a MultiGeometry is in the boundary or the interior of the Geometry
 *
 * The SFS uses the "Mod-2 Rule", which this function implements
 *
 * An alternative (and possibly more intuitive) rule would be
 * the "At Most One Rule":
 *    isInBoundary = (componentCount == 1)
 */
bool
GeometryGraph::isInBoundary(int boundaryCount)
{
	// the "Mod-2 Rule"
	return boundaryCount%2==1;
}

int
GeometryGraph::determineBoundary(int boundaryCount)
{
	return isInBoundary(boundaryCount)?Location::BOUNDARY : Location::INTERIOR;
}


EdgeSetIntersector*
GeometryGraph::createEdgeSetIntersector()
{
	// various options for computing intersections, from slowest to fastest

	//private EdgeSetIntersector esi = new SimpleEdgeSetIntersector();
	//private EdgeSetIntersector esi = new MonotoneChainIntersector();
	//private EdgeSetIntersector esi = new NonReversingChainIntersector();
	//private EdgeSetIntersector esi = new SimpleSweepLineIntersector();
	//private EdgeSetIntersector esi = new MCSweepLineIntersector();

	//return new SimpleEdgeSetIntersector();
	return new SimpleMCSweepLineIntersector();
}

/*public*/
vector<Node*>*
GeometryGraph::getBoundaryNodes()
{
	if ( ! boundaryNodes.get() )
	{
		boundaryNodes.reset(new vector<Node*>());
		getBoundaryNodes(*(boundaryNodes.get()));
	}
	return boundaryNodes.get();
}

/*public*/
CoordinateSequence*
GeometryGraph::getBoundaryPoints()
{

	if ( ! boundaryPoints.get() )
	{
		// Collection will be destroied by GeometryGraph dtor
		vector<Node*>* coll = getBoundaryNodes();
		boundaryPoints.reset(new CoordinateArraySequence(coll->size()));
		size_t i=0;
		for (vector<Node*>::iterator it=coll->begin(), endIt=coll->end();
			it!=endIt; ++it)
		{
			Node *node=*it;
			boundaryPoints->setAt(node->getCoordinate(), i++);
		}
	}

	// We keep ownership of this, will be destroyed by destructor
	return boundaryPoints.get();
}

Edge*
GeometryGraph::findEdge(const LineString *line)
{
	return lineEdgeMap.find(line)->second;
}

void
GeometryGraph::computeSplitEdges(vector<Edge*> *edgelist)
{
#if GEOS_DEBUG
	cerr<<"["<<this<<"] GeometryGraph::computeSplitEdges() scanning "<<edges->size()<<" local and "<<edgelist->size()<<" provided edges"<<endl;
#endif
	for (vector<Edge*>::iterator i=edges->begin(), endIt=edges->end();
		i!=endIt; ++i)
	{
		Edge *e=*i;
#if GEOS_DEBUG
		cerr<<"   "<<e->print()<<" adding split edges from arg"<<endl;
#endif
		e->eiList.addSplitEdges(edgelist);
	}
#if GEOS_DEBUG
	cerr<<"["<<this<<"] GeometryGraph::computeSplitEdges() completed "<<endl;
#endif
}

void
GeometryGraph::add(const Geometry *g)
	//throw (UnsupportedOperationException *)
{
	if (g->isEmpty()) return;

	// check if this Geometry should obey the Boundary Determination Rule
	// all collections except MultiPolygons obey the rule
	if ( dynamic_cast<const MultiPolygon*>(g) )
		useBoundaryDeterminationRule = false;


	if ( const Polygon* x1 = dynamic_cast<const Polygon*>(g) )
		addPolygon(x1);

	// LineString also handles LinearRings
	else if ( const LineString* x2 = dynamic_cast<const LineString*>(g) )
		addLineString(x2);

	else if ( const Point* x3 = dynamic_cast<const Point*>(g) )
		addPoint(x3);

	else if ( const GeometryCollection* x4 =
            dynamic_cast<const GeometryCollection*>(g) )
  {
		addCollection(x4);
  }

	else {
		string out=typeid(*g).name();
		throw util::UnsupportedOperationException("GeometryGraph::add(Geometry *): unknown geometry type: "+out);
	}
}

void
GeometryGraph::addCollection(const GeometryCollection *gc)
{
	for (size_t i=0, n=gc->getNumGeometries(); i<n; ++i)
	{
		const Geometry *g=gc->getGeometryN(i);
		add(g);
	}
}

/*
 * Add a Point to the graph.
 */
void
GeometryGraph::addPoint(const Point *p)
{
	const Coordinate& coord=*(p->getCoordinate());
	insertPoint(argIndex, coord, Location::INTERIOR);
}

/*
 * The left and right topological location arguments assume that the ring
 * is oriented CW.
 * If the ring is in the opposite orientation,
 * the left and right locations must be interchanged.
 */
void
GeometryGraph::addPolygonRing(const LinearRing *lr, int cwLeft, int cwRight)
	// throw IllegalArgumentException (see below)
{
	// skip empty component (see bug #234)
	if ( lr->isEmpty() ) return;

	const CoordinateSequence *lrcl = lr->getCoordinatesRO();

	CoordinateSequence* coord=CoordinateSequence::removeRepeatedPoints(lrcl);
	if (coord->getSize()<4) {
		hasTooFewPointsVar=true;
		invalidPoint=coord->getAt(0); // its now a Coordinate
		delete coord;
		return;
	}
	int left=cwLeft;
	int right=cwRight;

	/*
	 * the isCCW call might throw an
	 * IllegalArgumentException if degenerate ring does
	 * not contain 3 distinct points.
	 */
	try
	{
		if (CGAlgorithms::isCCW(coord)) {
			left=cwRight;
			right=cwLeft;
		}
	}
	catch(...)
	{
		delete coord;
		throw;
	}

	Edge *e=new Edge(coord, Label(argIndex, Location::BOUNDARY, left, right));
	lineEdgeMap[lr]=e;
	insertEdge(e);
	insertPoint(argIndex,coord->getAt(0), Location::BOUNDARY);
}

void
GeometryGraph::addPolygon(const Polygon *p)
{
	const LineString* ls;
	const LinearRing* lr;

	ls = p->getExteriorRing();
	assert(dynamic_cast<const LinearRing*>(ls));
	lr = static_cast<const LinearRing*>(ls);
	addPolygonRing(lr, Location::EXTERIOR, Location::INTERIOR);
	for (size_t i=0, n=p->getNumInteriorRing(); i<n; ++i)
	{
		// Holes are topologically labelled opposite to the shell, since
		// the interior of the polygon lies on their opposite side
		// (on the left, if the hole is oriented CW)
		ls = p->getInteriorRingN(i);
		assert(dynamic_cast<const LinearRing*>(ls));
		lr = static_cast<const LinearRing*>(ls);
		addPolygonRing(lr, Location::INTERIOR, Location::EXTERIOR);
	}
}

void
GeometryGraph::addLineString(const LineString *line)
{
	CoordinateSequence* coord=CoordinateSequence::removeRepeatedPoints(line->getCoordinatesRO());
	if(coord->getSize()<2) {
		hasTooFewPointsVar=true;
		invalidPoint=coord->getAt(0);
		delete coord;
		return;
	}

	Edge *e=new Edge(coord, Label(argIndex, Location::INTERIOR));
	lineEdgeMap[line]=e;
	insertEdge(e);

	/*
	 * Add the boundary points of the LineString, if any.
	 * Even if the LineString is closed, add both points as if they
	 * were endpoints.
	 * This allows for the case that the node already exists and is
	 * a boundary point.
	 */
	assert(coord->size() >= 2); // found LineString with single point
	insertBoundaryPoint(argIndex, coord->getAt(0));
	insertBoundaryPoint(argIndex, coord->getAt(coord->getSize()-1));
}

/*
 * Add an Edge computed externally.  The label on the Edge is assumed
 * to be correct.
 */
void
GeometryGraph::addEdge(Edge *e)
{
	insertEdge(e);
	const CoordinateSequence* coord=e->getCoordinates();
	// insert the endpoint as a node, to mark that it is on the boundary
	insertPoint(argIndex,coord->getAt(0),Location::BOUNDARY);
	insertPoint(argIndex,coord->getAt(coord->getSize()-1),Location::BOUNDARY);
}

/*
 * Add a point computed externally.  The point is assumed to be a
 * Point Geometry part, which has a location of INTERIOR.
 */
void
GeometryGraph::addPoint(Coordinate& pt)
{
	insertPoint(argIndex,pt,Location::INTERIOR);
}

template <class T, class C>
void collect_intersecting_edges(const Envelope *env, T start, T end, C &to)
{
  for (T i=start; i != end; ++i)
  {
    Edge *e = *i;
    if ( e->getEnvelope()->intersects(env) ) to.push_back(e);
  }
}

/*public*/
SegmentIntersector*
GeometryGraph::computeSelfNodes(LineIntersector &li,
  bool computeRingSelfNodes, const Envelope *env)
{
	return computeSelfNodes(li, computeRingSelfNodes, false, env);
}

SegmentIntersector*
GeometryGraph::computeSelfNodes(LineIntersector &li,
  bool computeRingSelfNodes, bool isDoneIfProperInt, const Envelope *env)
{
	SegmentIntersector *si = new SegmentIntersector(&li, true, false);
	si->setIsDoneIfProperInt(isDoneIfProperInt);
	unique_ptr<EdgeSetIntersector> esi(createEdgeSetIntersector());

	typedef vector<Edge*> EC;
	EC *se = edges;
	EC self_edges_copy;

	if ( env && ! env->covers(parentGeom->getEnvelopeInternal()) ) {
		collect_intersecting_edges(env, se->begin(), se->end(), self_edges_copy);
    //cerr << "(computeSelfNodes) Self edges reduced from " << se->size() << " to " << self_edges_copy.size() << endl;
		se = &self_edges_copy;
	}

	bool isRings = dynamic_cast<const LinearRing*>(parentGeom)
	    || dynamic_cast<const Polygon*>(parentGeom)
	    || dynamic_cast<const MultiPolygon*>(parentGeom);

	bool computeAllSegments = computeRingSelfNodes || ! isRings;

	esi->computeIntersections(se, si, computeAllSegments);

#if GEOS_DEBUG
	cerr << "SegmentIntersector # tests = " << si->numTests << endl;
#endif // GEOS_DEBUG

	addSelfIntersectionNodes(argIndex);
	return si;
}

SegmentIntersector*
GeometryGraph::computeEdgeIntersections(GeometryGraph *g,
	LineIntersector *li, bool includeProper, const Envelope *env)
{
#if GEOS_DEBUG
	cerr<<"GeometryGraph::computeEdgeIntersections call"<<endl;
#endif
	SegmentIntersector *si=new SegmentIntersector(li, includeProper, true);

	si->setBoundaryNodes(getBoundaryNodes(), g->getBoundaryNodes());
	unique_ptr<EdgeSetIntersector> esi(createEdgeSetIntersector());

	typedef vector<Edge*> EC;

	EC self_edges_copy;
	EC other_edges_copy;

	EC *se = edges;
	EC *oe = g->edges;
	if ( env && ! env->covers(parentGeom->getEnvelopeInternal()) ) {
		collect_intersecting_edges(env, se->begin(), se->end(), self_edges_copy);
    //cerr << "Self edges reduced from " << se->size() << " to " << self_edges_copy.size() << endl;
		se = &self_edges_copy;
	}
	if ( env && ! env->covers(g->parentGeom->getEnvelopeInternal()) ) {
		collect_intersecting_edges(env, oe->begin(), oe->end(), other_edges_copy);
    //cerr << "Other edges reduced from " << oe->size() << " to " << other_edges_copy.size() << endl;
		oe = &other_edges_copy;
	}
	esi->computeIntersections(se, oe, si);
#if GEOS_DEBUG
	cerr<<"GeometryGraph::computeEdgeIntersections returns"<<endl;
#endif
	return si;
}

void
GeometryGraph::insertPoint(int p_argIndex, const Coordinate& coord,
	int onLocation)
{
#if GEOS_DEBUG > 1
	cerr<<"GeometryGraph::insertPoint("<<coord.toString()<<" called"<<endl;
#endif
	Node *n=nodes->addNode(coord);
	Label& lbl = n->getLabel();
	if ( lbl.isNull() )
	{
		n->setLabel(p_argIndex, onLocation);
	}
	else
	{
		lbl.setLocation(p_argIndex, onLocation);
	}
}

/*
 * Adds points using the mod-2 rule of SFS.  This is used to add the boundary
 * points of dim-1 geometries (Curves/MultiCurves).  According to the SFS,
 * an endpoint of a Curve is on the boundary
 * iff if it is in the boundaries of an odd number of Geometries
 */
void
GeometryGraph::insertBoundaryPoint(int p_argIndex, const Coordinate& coord)
{
	Node *n=nodes->addNode(coord);
	// nodes always have labels
	Label& lbl = n->getLabel();

	// the new point to insert is on a boundary
	int boundaryCount=1;

	// determine the current location for the point (if any)
	int loc = lbl.getLocation(p_argIndex,Position::ON);
	if (loc==Location::BOUNDARY) boundaryCount++;

	// determine the boundary status of the point according to the
	// Boundary Determination Rule
	int newLoc = determineBoundary(boundaryNodeRule, boundaryCount);
	lbl.setLocation(p_argIndex,newLoc);
}

/*private*/
void
GeometryGraph::addSelfIntersectionNodes(int p_argIndex)
{
	for (vector<Edge*>::iterator i=edges->begin(), endIt=edges->end();
		i!=endIt; ++i)
	{
		Edge *e=*i;
		int eLoc = e->getLabel().getLocation(p_argIndex);
		EdgeIntersectionList &eiL = e->eiList;
		for (EdgeIntersectionList::iterator
			eiIt=eiL.begin(), eiEnd=eiL.end();
			eiIt!=eiEnd; ++eiIt)
		{
			EdgeIntersection *ei=*eiIt;
			addSelfIntersectionNode(p_argIndex, ei->coord, eLoc);
			GEOS_CHECK_FOR_INTERRUPTS();
		}
	}
}

/*private*/
void
GeometryGraph::addSelfIntersectionNode(int p_argIndex,
	const Coordinate& coord, int loc)
{
	// if this node is already a boundary node, don't change it
	if (isBoundaryNode(p_argIndex,coord)) return;
	if (loc==Location::BOUNDARY && useBoundaryDeterminationRule)
	{
		insertBoundaryPoint(p_argIndex,coord);
	}
	else
	{
		insertPoint(p_argIndex,coord,loc);
	}
}

vector<Edge*> *
GeometryGraph::getEdges()
{
	return edges;
}

bool
GeometryGraph::hasTooFewPoints()
{
	return hasTooFewPointsVar;
}

const Coordinate&
GeometryGraph::getInvalidPoint()
{
	return invalidPoint;
}

GeometryGraph::GeometryGraph(int newArgIndex,
		const geom::Geometry *newParentGeom)
	:
	PlanarGraph(),
	parentGeom(newParentGeom),
	useBoundaryDeterminationRule(true),
    boundaryNodeRule(algorithm::BoundaryNodeRule::getBoundaryOGCSFS()),
	argIndex(newArgIndex),
	hasTooFewPointsVar(false)
{
	if (parentGeom!=nullptr) add(parentGeom);
}

GeometryGraph::GeometryGraph(int newArgIndex,
		const geom::Geometry *newParentGeom,
		const algorithm::BoundaryNodeRule& bnr)
	:
	PlanarGraph(),
	parentGeom(newParentGeom),
	useBoundaryDeterminationRule(true),
	boundaryNodeRule(bnr),
	argIndex(newArgIndex),
	hasTooFewPointsVar(false)
{
	if (parentGeom!=nullptr) add(parentGeom);
}

GeometryGraph::GeometryGraph()
	:
	PlanarGraph(),
	parentGeom(nullptr),
	useBoundaryDeterminationRule(true),
    boundaryNodeRule(algorithm::BoundaryNodeRule::getBoundaryOGCSFS()),
	argIndex(-1),
	hasTooFewPointsVar(false)
{
}


/* public static */
int
GeometryGraph::determineBoundary(
	             const algorithm::BoundaryNodeRule& boundaryNodeRule,
	                                            int boundaryCount)
{
	return boundaryNodeRule.isInBoundary(boundaryCount)
		? Location::BOUNDARY : Location::INTERIOR;
}

} // namespace geos.geomgraph
} // namespace geos