overlayop.cpp

在Linux下做的QuadTree的程序

/**********************************************************************
 * $Id: OverlayOp.cpp 1932 2006-12-05 10:42:05Z strk $
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.refractions.net
 *
 * 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: operation/overlay/OverlayOp.java rev. 1.23
 *
 **********************************************************************/

#include <geos/operation/overlay/OverlayOp.h>
#include <geos/operation/overlay/OverlayResultValidator.h>
#include <geos/operation/overlay/ElevationMatrix.h>
#include <geos/operation/overlay/OverlayNodeFactory.h>
#include <geos/operation/overlay/PolygonBuilder.h>
#include <geos/operation/overlay/LineBuilder.h>
#include <geos/operation/overlay/PointBuilder.h>
#include <geos/geom/Geometry.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/geom/Polygon.h>
#include <geos/geom/LineString.h>
#include <geos/geom/Point.h>
#include <geos/geom/PrecisionModel.h>
#include <geos/geomgraph/Label.h>
#include <geos/geomgraph/Edge.h>
#include <geos/geomgraph/Node.h>
#include <geos/geomgraph/GeometryGraph.h>
#include <geos/geomgraph/EdgeEndStar.h>
#include <geos/geomgraph/DirectedEdgeStar.h>
#include <geos/geomgraph/DirectedEdge.h>
#include <geos/geomgraph/Position.h>
#include <geos/geomgraph/index/SegmentIntersector.h>
#include <geos/util/TopologyException.h>
#include <geos/precision/SimpleGeometryPrecisionReducer.h>
#include <geos/geomgraph/EdgeNodingValidator.h>

#include <cassert>
#include <functional>
#include <vector>
#include <sstream>
#include <memory> // for auto_ptr

#ifndef GEOS_DEBUG
#define GEOS_DEBUG 0
#endif

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

// A result validator using FuzzyPointLocator to
// check validity of result. Pretty expensive...
//#define ENABLE_OVERLAY_RESULT_VALIDATOR 1

// Edge noding validator, more lightweighted and
// catches robustness errors earlier
#define ENABLE_EDGE_NODING_VALIDATOR 1

// Define this to get some reports about validations
//#define GEOS_DEBUG_VALIDATION 1

// Other validators, not found in JTS
//#define ENABLE_OTHER_OVERLAY_RESULT_VALIDATORS 1

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

namespace geos {
namespace operation { // geos.operation
namespace overlay { // geos.operation.overlay

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

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


/* static public */
bool
OverlayOp::isResultOfOp(int loc0, int loc1, OverlayOp::OpCode opCode)
{
	if (loc0==Location::BOUNDARY) loc0=Location::INTERIOR;
	if (loc1==Location::BOUNDARY) loc1=Location::INTERIOR;
	switch (opCode) {
		case opINTERSECTION:
			return loc0==Location::INTERIOR && loc1==Location::INTERIOR;
		case opUNION:
			return loc0==Location::INTERIOR || loc1==Location::INTERIOR;
		case opDIFFERENCE:
			return loc0==Location::INTERIOR && loc1!=Location::INTERIOR;
		case opSYMDIFFERENCE:
			return (loc0==Location::INTERIOR && loc1!=Location::INTERIOR) 
				|| (loc0!=Location::INTERIOR && loc1==Location::INTERIOR);
	}
	return false;
}

OverlayOp::OverlayOp(const Geometry *g0, const Geometry *g1)

	:

	// this builds graphs in arg[0] and arg[1]
	GeometryGraphOperation(g0, g1),

	/*
	 * 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),
	graph(OverlayNodeFactory::instance()),
	resultPolyList(NULL),
	resultLineList(NULL),
	resultPointList(NULL)

{

#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 GEOS_DEBUG
	cerr<<elevationMatrix->print()<<endl;
#endif
#endif // USE_ELEVATION_MATRIX
#endif // COMPUTE_Z
}

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

/*public*/
Geometry*
OverlayOp::getResultGeometry(OverlayOp::OpCode funcCode)
	//throw(TopologyException *)
{
	computeOverlay(funcCode);
	return resultGeom;
}

/*private*/
void
OverlayOp::insertUniqueEdges(vector<Edge*> *edges)
{
	for_each(edges->begin(), edges->end(),
			bind1st(mem_fun(&OverlayOp::insertUniqueEdge), this));

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

}

/*private*/
void
OverlayOp::replaceCollapsedEdges()
{
	vector<Edge*>& edges=edgeList.getEdges();

	for(size_t i=0, nedges=edges.size(); i<nedges; ++i)
	{
		Edge *e=edges[i];
		assert(e);
		if (e->isCollapsed())
		{
#if GEOS_DEBUG
		cerr << " replacing collapsed edge " << i << endl;
#endif // GEOS_DEBUG
			//Debug.print(e);
			edges[i]=e->getCollapsedEdge();

			// should we keep this alive some more ?
			delete e;
		} 
	}
}

/*private*/
void
OverlayOp::copyPoints(int argIndex)
{
	NodeMap::container& nodeMap=arg[argIndex]->getNodeMap()->nodeMap;
	for ( NodeMap::const_iterator it=nodeMap.begin(), itEnd=nodeMap.end();
			it != itEnd; ++it )
	{
		Node* graphNode=it->second;
		assert(graphNode);

		Node* newNode=graph.addNode(graphNode->getCoordinate());
		assert(newNode);

		newNode->setLabel(argIndex,graphNode->getLabel()->getLocation(argIndex));
	}
}

/*private*/
void
OverlayOp::computeLabelling()
	//throw(TopologyException *) // and what else ?
{
	NodeMap::container& nodeMap=graph.getNodeMap()->nodeMap;

#if GEOS_DEBUG
	cerr<<"OverlayOp::computeLabelling(): at call time: "<<edgeList.print()<<endl;
#endif

#if GEOS_DEBUG
	cerr<<"OverlayOp::computeLabelling() scanning "<<nodeMap.size()<<" nodes from map:"<<endl;
#endif

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

/*private*/
void
OverlayOp::mergeSymLabels()
{
	NodeMap::container& nodeMap=graph.getNodeMap()->nodeMap;

#if GEOS_DEBUG
	cerr<<"OverlayOp::mergeSymLabels() scanning "<<nodeMap.size()<<" nodes from map:"<<endl;
#endif

	for ( NodeMap::const_iterator it=nodeMap.begin(), itEnd=nodeMap.end();
			it != itEnd; ++it )
	{
		Node *node=it->second;
		EdgeEndStar* ees=node->getEdges();
		assert(dynamic_cast<DirectedEdgeStar*>(ees));
		static_cast<DirectedEdgeStar*>(ees)->mergeSymLabels();
		//((DirectedEdgeStar*)node->getEdges())->mergeSymLabels();
#if GEOS_DEBUG
		cerr<<"     "<<node->print()<<endl;
#endif
		//node.print(System.out);
	}
}

/*private*/
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)

	NodeMap::container& nodeMap=graph.getNodeMap()->nodeMap;

#if GEOS_DEBUG
	cerr << "OverlayOp::updateNodeLabelling() scanning "
	     << nodeMap.size() << " nodes from map:" << endl;
#endif

	for ( NodeMap::const_iterator it=nodeMap.begin(), itEnd=nodeMap.end();
			it != itEnd; ++it )
	{
		Node *node=it->second;
		EdgeEndStar* ees = node->getEdges();
		assert( dynamic_cast<DirectedEdgeStar*>(ees) );
		DirectedEdgeStar* des = static_cast<DirectedEdgeStar*>(ees);
		Label &lbl = des->getLabel();
		node->getLabel()->merge(lbl);
#if GEOS_DEBUG
		cerr<<"     "<<node->print()<<endl;
#endif
	}
}

/*private*/
void
OverlayOp::labelIncompleteNodes()
{
	NodeMap::container& nodeMap=graph.getNodeMap()->nodeMap;

#if GEOS_DEBUG
	cerr<<"OverlayOp::labelIncompleteNodes() scanning "<<nodeMap.size()<<" nodes from map:"<<endl;
#endif

	for ( NodeMap::const_iterator it=nodeMap.begin(), itEnd=nodeMap.end();
			it != itEnd; ++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
		EdgeEndStar* ees = n->getEdges();
		assert( dynamic_cast<DirectedEdgeStar*>(ees) );
		DirectedEdgeStar* des = static_cast<DirectedEdgeStar*>(ees);

		des->updateLabelling(label);
		//((DirectedEdgeStar*)n->getEdges())->updateLabelling(label);
		//n.print(System.out);
	}
}

/*private*/
void
OverlayOp::labelIncompleteNode(Node *n, int targetIndex)
{
#if GEOS_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 GEOS_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
}

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

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

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

/*private*/
double
OverlayOp::getAverageZ(int targetIndex)
{
	if ( avgzcomputed[targetIndex] ) return avgz[targetIndex];

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

	// OverlayOp::getAverageZ(int) called with a ! polygon
	assert(targetGeom->getGeometryTypeId() == GEOS_POLYGON);

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

/*private*/
int
OverlayOp::mergeZ(Node *n, const Polygon *poly) const
{
	const LineString *ls;
	int found = 0;
	ls = poly->getExteriorRing();
	found = mergeZ(n, ls);
	if ( found ) return 1;
	for (size_t i=0, nr=poly->getNumInteriorRing(); i<nr; ++i)
	{
		ls = poly->getInteriorRingN(i);
		found = mergeZ(n, ls);
		if ( found ) return 1;
	}
	return 0;
}

/*private*/
int
OverlayOp::mergeZ(Node *n, const LineString *line) const
{
	const CoordinateSequence *pts = line->getCoordinatesRO();
	const Coordinate &p = n->getCoordinate();
	LineIntersector li;
	for(size_t i=1, size=pts->size(); i<size; ++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(LineIntersector::interpolateZ(p,
					p0, p1));
			}
			return 1;
		}
	}