bufferbuilder.cpp

在Linux下做的QuadTree的程序

/**********************************************************************
 * $Id: BufferBuilder.cpp 1998 2007-08-23 20:35:16Z pramsey $
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.refractions.net
 *
 * Copyright (C) 2005-2007 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/buffer/BufferBuilder.java rev. 1.23 (JTS-1.7)
 *
 **********************************************************************/

#include <geos/operation/overlay/PolygonBuilder.h> 
#include <geos/operation/overlay/OverlayNodeFactory.h> 
#include <geos/geom/GeometryFactory.h>
#include <geos/geom/Location.h>
#include <geos/geom/Geometry.h>
#include <geos/geom/Polygon.h>
#include <geos/geom/GeometryCollection.h>
#include <geos/operation/buffer/BufferBuilder.h>
#include <geos/operation/buffer/OffsetCurveBuilder.h>
#include <geos/operation/buffer/OffsetCurveSetBuilder.h>
#include <geos/operation/buffer/BufferSubgraph.h>
#include <geos/operation/buffer/SubgraphDepthLocater.h>
#include <geos/algorithm/LineIntersector.h>
#include <geos/noding/IntersectionAdder.h>
#include <geos/noding/SegmentString.h>
#include <geos/noding/MCIndexNoder.h>
#include <geos/geomgraph/Position.h>
#include <geos/geomgraph/PlanarGraph.h>
#include <geos/geomgraph/Label.h>
#include <geos/geomgraph/Node.h>
#include <geos/geomgraph/Edge.h>
#include <geos/util/GEOSException.h>
#include <geos/profiler.h>

#include <cassert>
#include <vector>
#include <iomanip>
#include <algorithm>
#include <iostream>

#ifndef GEOS_DEBUG
#define GEOS_DEBUG 0
#endif

#ifndef JTS_DEBUG 
#define JTS_DEBUG 0
#endif

//using namespace std;
using namespace geos::geom;
using namespace geos::geomgraph;
using namespace geos::noding;
using namespace geos::algorithm;
using namespace geos::operation::overlay;

namespace geos {
namespace operation { // geos.operation
namespace buffer { // geos.operation.buffer

#if PROFILE
static Profiler *profiler = Profiler::instance();
#endif

/**
* Compute the change in depth as an edge is crossed from R to L
*/
int
BufferBuilder::depthDelta(Label *label)
{
	int lLoc=label->getLocation(0, Position::LEFT);
	int rLoc=label->getLocation(0, Position::RIGHT);
	if (lLoc== Location::INTERIOR && rLoc== Location::EXTERIOR)
		return 1;
	else if (lLoc== Location::EXTERIOR && rLoc== Location::INTERIOR)
		return -1;
	return 0;
}

//static CGAlgorithms rCGA;
//CGAlgorithms *BufferBuilder::cga=&rCGA;

BufferBuilder::~BufferBuilder()
{
	delete li; // could be NULL
	delete intersectionAdder;
	//delete edgeList;
	for (size_t i=0; i<newLabels.size(); i++)
		delete newLabels[i];
}

/*public*/
Geometry*
BufferBuilder::buffer(const Geometry *g, double distance)
	// throw(GEOSException *)
{
	const PrecisionModel *precisionModel=workingPrecisionModel;
	if (precisionModel==NULL)
		precisionModel=g->getPrecisionModel();

	assert(precisionModel);
	assert(g);

	// factory must be the same as the one used by the input
	geomFact=g->getFactory();

	OffsetCurveBuilder curveBuilder(precisionModel, quadrantSegments);
	curveBuilder.setEndCapStyle(endCapStyle);
	OffsetCurveSetBuilder curveSetBuilder(*g, distance, curveBuilder);

	std::vector<SegmentString*>& bufferSegStrList=curveSetBuilder.getCurves();

#if GEOS_DEBUG
	std::cerr << "OffsetCurveSetBuilder got " << bufferSegStrList.size()
	          << " curves" << std::endl;
#endif
	// short-circuit test
	if (bufferSegStrList.size()<=0) {
		return createEmptyResultGeometry();
	}

#if GEOS_DEBUG
	std::cerr<<"BufferBuilder::buffer computing NodedEdges"<<std::endl;
#endif

#if JTS_DEBUG
std::cerr << "before noding: SegStr # " << bufferSegStrList.size() << std::endl;
for (size_t i = 0, n=bufferSegStrList.size(); i<n; i++)
{
 	SegmentString* segStr = bufferSegStrList[i];
	std::cerr << "SegStr " << i << ": pts # " << segStr->size()
		<< " nodes # " << segStr->getNodeList().size()
		<< std::endl;
}
#endif

	computeNodedEdges(bufferSegStrList, precisionModel);

#if JTS_DEBUG
std::cerr << "after noding: SegStr # " << bufferSegStrList.size() << std::endl;
for (size_t i = 0, n=bufferSegStrList.size(); i<n; i++)
{
 	SegmentString* segStr = bufferSegStrList[i];
	std::cerr << "SegStr " << i << ": pts # " << segStr->size()
		<< " nodes # " << segStr->getNodeList().size()
		<< std::endl;
}
#endif

#if GEOS_DEBUG > 1
	std::cerr << std::endl << edgeList << std::endl;
#endif

	Geometry* resultGeom=NULL;
	std::vector<Geometry*> *resultPolyList=NULL;
	std::vector<BufferSubgraph*> subgraphList;

	try {
		PlanarGraph graph(OverlayNodeFactory::instance());
		graph.addEdges(edgeList.getEdges());

		createSubgraphs(&graph, subgraphList);
#if GEOS_DEBUG
	std::cerr<<"Created "<<subgraphList.size()<<" subgraphs"<<std::endl;
#if GEOS_DEBUG > 1
	for (size_t i=0, n=subgraphList.size(); i<n; i++)
		std::cerr << std::setprecision(10) << *(subgraphList[i]) << std::endl;
#endif
#endif
		PolygonBuilder polyBuilder(geomFact);
		buildSubgraphs(subgraphList, polyBuilder);
		resultPolyList=polyBuilder.getPolygons();
#if GEOS_DEBUG
	std::cerr << "PolygonBuilder got " << resultPolyList->size()
	          << " polygons" << std::endl;
#if GEOS_DEBUG > 1
	for (size_t i=0, n=resultPolyList->size(); i<n; i++)
		std::cerr << (*resultPolyList)[i]->toString() << std::endl;
#endif
#endif
		// just in case ...
		if ( resultPolyList->empty() )
		{
			return createEmptyResultGeometry();
		}

		resultGeom=geomFact->buildGeometry(resultPolyList);
	} catch (const util::GEOSException& /* exc */) {
		for (size_t i=0, n=subgraphList.size(); i<n; i++)
			delete subgraphList[i];
		throw;
	} 

	for (size_t i=0, n=subgraphList.size(); i<n; i++)
		delete subgraphList[i];

	return resultGeom;
}

/*private*/
Noder*
BufferBuilder::getNoder(const PrecisionModel* pm)
{
	// this doesn't change workingNoder precisionModel!
	if (workingNoder != NULL) return workingNoder;

	// otherwise use a fast (but non-robust) noder

	if ( li ) // reuse existing IntersectionAdder and LineIntersector
	{
		li->setPrecisionModel(pm);
		assert(intersectionAdder!=NULL);
	}
	else
	{
		li = new LineIntersector(pm);
		intersectionAdder = new IntersectionAdder(*li);
	}

	MCIndexNoder* noder = new MCIndexNoder(intersectionAdder);

#if 0
	/* CoordinateArraySequence.cpp:84:
	 * virtual const geos::Coordinate& geos::CoordinateArraySequence::getAt(size_t) const:
	 * Assertion `pos<vect->size()' failed.
	 */
	//Noder* noder = new snapround::SimpleSnapRounder(*pm);

	Noder* noder = new IteratedNoder(pm);

	Noder noder = new SimpleSnapRounder(pm);
	Noder noder = new MCIndexSnapRounder(pm);
	Noder noder = new ScaledNoder(
		new MCIndexSnapRounder(new PrecisionModel(1.0)),
			pm.getScale());
#endif

	return noder;

}

/* private */
void
BufferBuilder::computeNodedEdges(SegmentString::NonConstVect& bufferSegStrList,
		const PrecisionModel *precisionModel) // throw(GEOSException)
{
	Noder *noder = getNoder(precisionModel);

	noder->computeNodes(&bufferSegStrList);

	SegmentString::NonConstVect* nodedSegStrings = \
			noder->getNodedSubstrings();


	for (SegmentString::NonConstVect::iterator
		i=nodedSegStrings->begin(), e=nodedSegStrings->end();
		i!=e;
		++i)
	{
		SegmentString* segStr = *i;
		const Label* oldLabel = static_cast<const Label*>(segStr->getData());

		CoordinateSequence* cs = CoordinateSequence::removeRepeatedPoints(segStr->getCoordinates());
		if ( cs->size() < 2 ) 
		{
			delete cs; // we need to take care of the memory here as cs is a new sequence
			return; // don't insert collapsed edges
		}
		// we need to clone SegmentString coordinates
		// as Edge will take ownership of them
		// TODO: find a way to transfer ownership instead
		// Who will own the edge ? FIXME: find out and handle that!
		Edge* edge = new Edge(cs, new Label(*oldLabel));

		// will take care of the Edge ownership
		insertEdge(edge);
	}

	if ( nodedSegStrings != &bufferSegStrList )
	{
		delete nodedSegStrings;
	}

	if ( noder != workingNoder ) delete noder;
}

/*private*/
void
BufferBuilder::insertEdge(Edge *e)
{
	//<FIX> MD 8 Oct 03  speed up identical edge lookup
	// fast lookup
	Edge *existingEdge=edgeList.findEqualEdge(e);
	// If an identical edge already exists, simply update its label
	if (existingEdge != NULL) {
		Label *existingLabel=existingEdge->getLabel();
		Label *labelToMerge=e->getLabel();

		// check if new edge is in reverse direction to existing edge
		// if so, must flip the label before merging it
		if (! existingEdge->isPointwiseEqual(e))