overlayop.cpp

来自「一个很好的vc底层代码」· C++ 代码 · 共 941 行 · 第 1/2 页

/**********************************************************************
 * $Id: OverlayOp.cpp,v 1.37.2.1 2005/05/23 16:41:45 strk Exp $
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.refractions.net
 *
 * 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.
 *
 **********************************************************************/

#include <geos/opOverlay.h>
#include <stdio.h>
#include <geos/util.h>

#define DEBUG 0
#define COMPUTE_Z 1
#define USE_ELEVATION_MATRIX 1
#define USE_INPUT_AVGZ 0

namespace geos {

Geometry*
OverlayOp::overlayOp(const Geometry *geom0,const Geometry *geom1,int opCode)
	// throw(TopologyException *)
{
	OverlayOp gov(geom0, geom1);
	return gov.getResultGeometry(opCode);
}

bool
OverlayOp::isResultOfOp(Label *label,int opCode)
{
	int loc0=label->getLocation(0);
	int loc1=label->getLocation(1);
	return isResultOfOp(loc0,loc1,opCode);
}


/*
 * This method will handle arguments of Location.NULL correctly
 *
 * @return true if the locations correspond to the opCode
 */
bool
OverlayOp::isResultOfOp(int loc0,int loc1,int opCode)
{
	if (loc0==Location::BOUNDARY) loc0=Location::INTERIOR;
	if (loc1==Location::BOUNDARY) loc1=Location::INTERIOR;
	switch (opCode) {
		case INTERSECTION:
			return loc0==Location::INTERIOR && loc1==Location::INTERIOR;
	case UNION:
		return loc0==Location::INTERIOR || loc1==Location::INTERIOR;
	case DIFFERENCE:
		return loc0==Location::INTERIOR && loc1!=Location::INTERIOR;
	case SYMDIFFERENCE:
		return (loc0==Location::INTERIOR && loc1!=Location::INTERIOR) 
			|| (loc0!=Location::INTERIOR && loc1==Location::INTERIOR);
	}
	return false;
}

OverlayOp::OverlayOp(const Geometry *g0, const Geometry *g1): GeometryGraphOperation(g0,g1)
{
	graph=new PlanarGraph(new OverlayNodeFactory());

	/*
	 * Use factory of primary geometry.
	 * Note that this does NOT handle mixed-precision arguments
	 * where the second arg has greater precision than the first.
	 */
	geomFact=g0->getFactory();
	resultGeom=NULL;
	edgeList=new EdgeList();
	resultPolyList=NULL;
	resultLineList=NULL;
	resultPointList=NULL;
	ptLocator=new PointLocator();
#if COMPUTE_Z
#if USE_INPUT_AVGZ
	avgz[0] = DoubleNotANumber;
	avgz[1] = DoubleNotANumber;
	avgzcomputed[0] = false;
	avgzcomputed[1] = false;
#endif // USE_INPUT_AVGZ

	Envelope env(*(g0->getEnvelopeInternal()));
	env.expandToInclude(g1->getEnvelopeInternal());
#if USE_ELEVATION_MATRIX
	elevationMatrix = new ElevationMatrix(env, 3, 3);
	elevationMatrix->add(g0);
	elevationMatrix->add(g1);
#if DEBUG
	cerr<<elevationMatrix->print()<<endl;
#endif
#endif // USE_ELEVATION_MATRIX
#endif // COMPUTE_Z
}

OverlayOp::~OverlayOp()
{
	delete graph;
	delete edgeList;
	delete resultPolyList;
	delete resultLineList;
	delete resultPointList;
	delete ptLocator;
	for (unsigned int i=0; i<dupEdges.size(); i++)
		delete dupEdges[i];
#if USE_ELEVATION_MATRIX
	delete elevationMatrix;
#endif
}

Geometry*
OverlayOp::getResultGeometry(int funcCode)
	//throw(TopologyException *)
{
	computeOverlay(funcCode); // this can throw TopologyException *
	return resultGeom;
}

PlanarGraph*
OverlayOp::getGraph()
{
	return graph;
}

void
OverlayOp::insertUniqueEdges(vector<Edge*> *edges)
{
	for(int i=0;i<(int)edges->size();i++) {
		Edge *e=(*edges)[i];
		insertUniqueEdge(e);
	}

#if DEBUG
	cerr<<"OverlayOp::insertUniqueEdges("<<edges->size()<<"): "<<endl;
	for(int i=0;i<(int)edges->size();i++) {
		Edge *e=(*edges)[i];
		if ( ! e ) cerr <<" NULL"<<endl;
		cerr <<" "<< e->print() << endl;
	}
#endif // DEBUG

}

/*
 * If edges which have undergone dimensional collapse are found,
 * replace them with a new edge which is a L edge
 */
void
OverlayOp::replaceCollapsedEdges()
{
	vector<Edge*> *newEdges=new vector<Edge*>();
	vector<Edge*> *oldEdges=new vector<Edge*>();
	for(int i=0;i<(int)edgeList->getEdges()->size();i++) {
		Edge *e=edgeList->get(i);
		if (e->isCollapsed()) {
			//Debug.print(e);
			newEdges->push_back(e->getCollapsedEdge());
			delete e;
		} else {
			//instead of removing from edgeList
			oldEdges->push_back(e);
		}
	}
	oldEdges->insert(oldEdges->end(),newEdges->begin(),newEdges->end());
	edgeList->getEdges()->assign(oldEdges->begin(),oldEdges->end());
	delete oldEdges;
	delete newEdges;
}

/*
 * Copy all nodes from an arg geometry into this graph.
 * The node label in the arg geometry overrides any previously computed
 * label for that argIndex.
 * (E.g. a node may be an intersection node with
 * a previously computed label of BOUNDARY,
 * but in the original arg Geometry it is actually
 * in the interior due to the Boundary Determination Rule)
 */
void
OverlayOp::copyPoints(int argIndex)
{
	map<Coordinate,Node*,CoordLT> *nodeMap=(*arg)[argIndex]->getNodeMap()->nodeMap;
	map<Coordinate,Node*,CoordLT>::iterator	it=nodeMap->begin();
	for (;it!=nodeMap->end();it++) {
		Node *graphNode=it->second;
		Node *newNode=graph->addNode(graphNode->getCoordinate());
		newNode->setLabel(argIndex,graphNode->getLabel()->getLocation(argIndex));
	}
}

/*
 * Compute initial labelling for all DirectedEdges at each node.
 * In this step, DirectedEdges will acquire a complete labelling
 * (i.e. one with labels for both Geometries)
 * only if they
 * are incident on a node which has edges for both Geometries
 */
void
OverlayOp::computeLabelling()
	//throw(TopologyException *) // and what else ?
{
	map<Coordinate,Node*,CoordLT> *nodeMap=graph->getNodeMap()->nodeMap;

#if DEBUG
	cerr<<"OverlayOp::computeLabelling(): at call time: "<<edgeList->print()<<endl;
#endif

#if DEBUG
	cerr<<"OverlayOp::computeLabelling() scanning "<<nodeMap->size()<<" nodes from map:"<<endl;
#endif

	map<Coordinate,Node*,CoordLT>::iterator	it=nodeMap->begin();
	for (;it!=nodeMap->end();it++) {
		Node *node=it->second;
#if DEBUG
		cerr<<"     "<<node->print()<<" has "<<node->getEdges()->getEdges()->size()<<" edgeEnds"<<endl;
#endif
		node->getEdges()->computeLabelling(arg);
	}
#if DEBUG
	cerr<<"OverlayOp::computeLabelling(): after edge labelling: "<<edgeList->print()<<endl;
#endif
	mergeSymLabels();
#if DEBUG
	cerr<<"OverlayOp::computeLabelling(): after labels sym merging: "<<edgeList->print()<<endl;
#endif
	updateNodeLabelling();
#if DEBUG
	cerr<<"OverlayOp::computeLabelling(): after node labeling update: "<<edgeList->print()<<endl;
#endif
}

/*
 * For nodes which have edges from only one Geometry incident on them,
 * the previous step will have left their dirEdges with no labelling
 * for the other Geometry. 
 * However, the sym dirEdge may have a labelling for the other
 * Geometry, so merge the two labels.
 */
void
OverlayOp::mergeSymLabels()
{
	map<Coordinate,Node*,CoordLT> *nodeMap=graph->getNodeMap()->nodeMap;

#if DEBUG
	cerr<<"OverlayOp::mergeSymLabels() scanning "<<nodeMap->size()<<" nodes from map:"<<endl;
#endif

	map<Coordinate,Node*,CoordLT>::iterator	it=nodeMap->begin();
	for (;it!=nodeMap->end();it++) {
		Node *node=it->second;
		((DirectedEdgeStar*)node->getEdges())->mergeSymLabels();
#if DEBUG
		cerr<<"     "<<node->print()<<endl;
#endif
		//node.print(System.out);
	}
}

void
OverlayOp::updateNodeLabelling()
{
	// update the labels for nodes
	// The label for a node is updated from the edges incident on it
	// (Note that a node may have already been labelled
	// because it is a point in one of the input geometries)
	map<Coordinate,Node*,CoordLT> *nodeMap=graph->getNodeMap()->nodeMap;
#if DEBUG
	cerr<<"OverlayOp::updateNodeLabelling() scanning "<<nodeMap->size()<<" nodes from map:"<<endl;
#endif
	map<Coordinate,Node*,CoordLT>::iterator	it=nodeMap->begin();
	for (;it!=nodeMap->end();it++) {
		Node *node=it->second;
		Label *lbl=((DirectedEdgeStar*)node->getEdges())->getLabel();
		node->getLabel()->merge(lbl);
#if DEBUG
		cerr<<"     "<<node->print()<<endl;
#endif
	}
}

/*
 * Incomplete nodes are nodes whose labels are incomplete.
 * (e.g. the location for one Geometry is NULL).
 * These are either isolated nodes,
 * or nodes which have edges from only a single Geometry incident on them.
 *
 * Isolated nodes are found because nodes in one graph which don't intersect
 * nodes in the other are not completely labelled by the initial process
 * of adding nodes to the nodeList.
 * To complete the labelling we need to check for nodes that lie in the
 * interior of edges, and in the interior of areas.
 *
 * When each node labelling is completed, the labelling of the incident
 * edges is updated, to complete their labelling as well.
 */
void
OverlayOp::labelIncompleteNodes()
{
	map<Coordinate,Node*,CoordLT> *nodeMap=graph->getNodeMap()->nodeMap;
#if DEBUG
	cerr<<"OverlayOp::labelIncompleteNodes() scanning "<<nodeMap->size()<<" nodes from map:"<<endl;
#endif
	map<Coordinate,Node*,CoordLT>::iterator	it=nodeMap->begin();
	for (;it!=nodeMap->end();it++) {
		Node *n=it->second;
		Label *label=n->getLabel();
		if (n->isIsolated()) {
			if (label->isNull(0))
				labelIncompleteNode(n,0);
			else
				labelIncompleteNode(n, 1);
		}
		// now update the labelling for the DirectedEdges incident on this node
		((DirectedEdgeStar*)n->getEdges())->updateLabelling(label);
		//n.print(System.out);
	}
}

/*
 * Label an isolated node with its relationship to the target geometry.
 */
void
OverlayOp::labelIncompleteNode(Node *n, int targetIndex)
{
#if DEBUG
	cerr<<"OverlayOp::labelIncompleteNode("<<n->print()<<", "<<targetIndex<<")"<<endl;
#endif
	const Geometry *targetGeom = (*arg)[targetIndex]->getGeometry();
	int loc=ptLocator->locate(n->getCoordinate(), targetGeom);
	n->getLabel()->setLocation(targetIndex,loc);

#if DEBUG
	cerr<<"   after location set: "<<n->print()<<endl;
#endif

#if COMPUTE_Z
	/*
	 * If this node has been labeled INTERIOR of a line
	 * or BOUNDARY of a polygon we must merge
	 * Z values of the intersected segment.
	 * The intersection point has been already computed
	 * by LineIntersector invoked by CGAlgorithms::isOnLine
	 * invoked by PointLocator.
	 */
	const LineString *line = dynamic_cast<const LineString *>(targetGeom);
	if ( loc == Location::INTERIOR && line )
	{
		mergeZ(n, line);
	}
	const Polygon *poly = dynamic_cast<const Polygon *>(targetGeom);
	if ( loc == Location::BOUNDARY && poly )
	{
		mergeZ(n, poly);
	}
#if USE_INPUT_AVGZ
	if ( loc == Location::INTERIOR && poly )
	{
		n->addZ(getAverageZ(targetIndex));
	}
#endif // USE_INPUT_AVGZ
#endif // COMPUTE_Z
}

double
OverlayOp::getAverageZ(const Polygon *poly)
{
	double totz = 0.0;
	int zcount = 0;

	const CoordinateSequence *pts =
		poly->getExteriorRing()->getCoordinatesRO();
	for (int i=0; i<pts->getSize(); i++)
	{
		const Coordinate &c = pts->getAt(i);
		if ( !ISNAN(c.z) )
		{
			totz += c.z;
			zcount++;
		}
	}

	if ( zcount ) return totz/zcount;
	else return DoubleNotANumber;
}

/*
 * This caches result to avoid multiple scans
 */
double
OverlayOp::getAverageZ(int targetIndex)
{
	if ( avgzcomputed[targetIndex] ) return avgz[targetIndex];

	const Geometry *targetGeom = (*arg)[targetIndex]->getGeometry();

	Assert::isTrue(targetGeom->getGeometryTypeId() == GEOS_POLYGON, 
		"OverlayOp::getAverageZ(int) called with a ! polygon");

	avgz[targetIndex] = getAverageZ((const Polygon *)targetGeom);
	avgzcomputed[targetIndex] = true;
	return avgz[targetIndex];
}

/*
 * Merge Z values of node with those of the segment or vertex in
 * the given Polygon it is on.
 */
int
OverlayOp::mergeZ(Node *n, const Polygon *poly) const
{
	const LineString *ls;
	int found = 0;
	ls = (const LineString *)poly->getExteriorRing();
	found = mergeZ(n, ls);
	if ( found ) return 1;
	for (int i=0; i<poly->getNumInteriorRing(); i++)
	{
		ls = (const LineString *)poly->getInteriorRingN(i);
		found = mergeZ(n, ls);
		if ( found ) return 1;
	}
	return 0;
}

/*
 * Merge Z values of node with those of the segment or vertex in
 * the given LineString it is on.
 * @returns 1 if an intersection is found, 0 otherwise.
 */
int
OverlayOp::mergeZ(Node *n, const LineString *line) const
{
	const CoordinateSequence *pts = line->getCoordinatesRO();
	const Coordinate &p = n->getCoordinate();
	RobustLineIntersector li;
	for(int i=1;i<pts->getSize();i++) {
		const Coordinate &p0=pts->getAt(i-1);
		const Coordinate &p1=pts->getAt(i);	
		li.computeIntersection(p, p0, p1);
		if (li.hasIntersection()) {
			if ( p == p0 ) n->addZ(p0.z);
			else if ( p == p1 ) n->addZ(p1.z);
			else {
				//n->addZ(p0.z);
				//n->addZ(p1.z);
				n->addZ(LineIntersector::interpolateZ(p,
					p0, p1));
			}
			return 1;
		}
	}
	return 0;
}

/*
 * Find all edges whose label indicates that they are in the result area(s),
 * according to the operation being performed.  Since we want polygon shells
 * to be oriented CW, choose dirEdges with the interior of the result on the
 * RHS.
 * Mark them as being in the result.
 * Interior Area edges are the result of dimensional collapses.