overlayop.cpp

在Linux下做的QuadTree的程序

	return 0;
}

/*private*/
void
OverlayOp::findResultAreaEdges(OverlayOp::OpCode opCode)
{
	vector<EdgeEnd*> *ee=graph.getEdgeEnds();
	for(size_t i=0, e=ee->size(); i<e; ++i)
	{
		DirectedEdge *de=(DirectedEdge*) (*ee)[i];
		// mark all dirEdges with the appropriate label
		Label *label=de->getLabel();
		if (label->isArea()
			&& !de->isInteriorAreaEdge()
			&& isResultOfOp(label->getLocation(0,Position::RIGHT),
					label->getLocation(1,Position::RIGHT),
					opCode)
			)
		{
				de->setInResult(true);
				//Debug.print("in result "); Debug.println(de);
		}
	}
}

/*private*/
void
OverlayOp::cancelDuplicateResultEdges()
{
	// remove any dirEdges whose sym is also included
	// (they "cancel each other out")
	vector<EdgeEnd*> *ee=graph.getEdgeEnds();
	for(size_t i=0, eesize=ee->size(); i<eesize; ++i)
	{
		DirectedEdge *de=static_cast<DirectedEdge*>( (*ee)[i] );
		DirectedEdge *sym=de->getSym();
		if (de->isInResult() && sym->isInResult()) {
			de->setInResult(false);
			sym->setInResult(false);
			//Debug.print("cancelled "); Debug.println(de); Debug.println(sym);
		}
	}
}

/*public*/
bool
OverlayOp::isCoveredByLA(const Coordinate& coord)
{
	if (isCovered(coord,resultLineList)) return true;
	if (isCovered(coord,resultPolyList)) return true;
	return false;
}

/*public*/
bool
OverlayOp::isCoveredByA(const Coordinate& coord)
{
	if (isCovered(coord,resultPolyList)) return true;
	return false;
}

/*private*/
bool
OverlayOp::isCovered(const Coordinate& coord,vector<Geometry*> *geomList)
{
	for(size_t i=0, n=geomList->size(); i<n; ++i)
	{
		Geometry *geom=(*geomList)[i];
		int loc=ptLocator.locate(coord,geom);
		if (loc!=Location::EXTERIOR) return true;
	}
	return false;
}

/*private*/
bool
OverlayOp::isCovered(const Coordinate& coord,vector<LineString*> *geomList)
{
	for(size_t i=0, n=geomList->size(); i<n; ++i)
	{
		Geometry *geom=(Geometry*)(*geomList)[i];
		int loc=ptLocator.locate(coord,geom);
		if (loc!=Location::EXTERIOR) return true;
	}
	return false;
}

/*private*/
bool
OverlayOp::isCovered(const Coordinate& coord,vector<Polygon*> *geomList)
{
	for(size_t i=0, n=geomList->size(); i<n; ++i)
	{
		Geometry *geom=(Geometry*)(*geomList)[i];
		int loc=ptLocator.locate(coord,geom);
		if (loc!=Location::EXTERIOR) return true;
	}
	return false;
}

/*private*/
Geometry*
OverlayOp::computeGeometry(vector<Point*> *nResultPointList,
                              vector<LineString*> *nResultLineList,
                              vector<Polygon*> *nResultPolyList)
{
	size_t nPoints=nResultPointList->size();
	size_t nLines=nResultLineList->size();
	size_t nPolys=nResultPolyList->size();

	vector<Geometry*> *geomList=new vector<Geometry*>();
	geomList->reserve(nPoints+nLines+nPolys);

	// element geometries of the result are always in the order P,L,A
	geomList->insert(geomList->end(),
			nResultPointList->begin(),
			nResultPointList->end());

	geomList->insert(geomList->end(),
			nResultLineList->begin(),
			nResultLineList->end());

	geomList->insert(geomList->end(),
			nResultPolyList->begin(),
			nResultPolyList->end());
			
	// build the most specific geometry possible
	Geometry *g=geomFact->buildGeometry(geomList);
	return g;
}

/*private*/
void
OverlayOp::computeOverlay(OverlayOp::OpCode opCode)
	//throw(TopologyException *)
{

	// copy points from input Geometries.
	// This ensures that any Point geometries
	// in the input are considered for inclusion in the result set
	copyPoints(0);
	copyPoints(1);

	// node the input Geometries
	delete arg[0]->computeSelfNodes(li,false);
	delete arg[1]->computeSelfNodes(li,false);

#if GEOS_DEBUG
	cerr<<"OverlayOp::computeOverlay: computed SelfNodes"<<endl;
#endif

	// compute intersections between edges of the two input geometries
	delete arg[0]->computeEdgeIntersections(arg[1], &li,true);

#if GEOS_DEBUG
	cerr<<"OverlayOp::computeOverlay: computed EdgeIntersections"<<endl;
	cerr<<"OverlayOp::computeOverlay: li: "<<li.toString()<<endl;
#endif


	vector<Edge*> baseSplitEdges;
	arg[0]->computeSplitEdges(&baseSplitEdges);
	arg[1]->computeSplitEdges(&baseSplitEdges);

	// add the noded edges to this result graph
	insertUniqueEdges(&baseSplitEdges);
	computeLabelsFromDepths();
	replaceCollapsedEdges();
	//Debug.println(edgeList);

#ifdef ENABLE_EDGE_NODING_VALIDATOR // {
	if ( resultPrecisionModel->isFloating() )
	{
		// Will throw TopologyException if noding is found to be invalid
		EdgeNodingValidator nv(edgeList.getEdges());

#ifdef GEOS_DEBUG_VALIDATION // {
		try {
			nv.checkValid();
		} catch (const util::TopologyException& ex) {
			cout << "EdgeNodingValidator found noding invalid: " << ex.what() << endl;
			throw ex;
		}
		cout << "EdgeNodingValidator accepted the noding" << endl;

#else // }{
		nv.checkValid();
#endif // }
	}
#ifdef GEOS_DEBUG_VALIDATION // {
	else
	{
		cout << "Did not run EdgeNodingValidator as the precision model is not floating" << endl;
	}
#endif // GEOS_DEBUG_VALIDATION }
#endif // ENABLE_EDGE_NODING_VALIDATOR }

	graph.addEdges(edgeList.getEdges());

	// this can throw TopologyException *
	computeLabelling();

	//Debug.printWatch();
	labelIncompleteNodes();
	//Debug.printWatch();
	//nodeMap.print(System.out);


	/*
	 * The ordering of building the result Geometries is important.
	 * Areas must be built before lines, which must be built
	 * before points.
	 * This is so that lines which are covered by areas are not
	 * included explicitly, and similarly for points.
	 */
	findResultAreaEdges(opCode);
	cancelDuplicateResultEdges();

	PolygonBuilder polyBuilder(geomFact);
	
	// might throw a TopologyException *
	polyBuilder.add(&graph);

	vector<Geometry*> *gv=polyBuilder.getPolygons();
	size_t gvsize=gv->size();
	resultPolyList=new vector<Polygon*>(gvsize);
	for(size_t i=0; i<gvsize; ++i) {
		(*resultPolyList)[i]=(Polygon*)(*gv)[i];
	}
	delete gv;

	LineBuilder lineBuilder(this,geomFact,&ptLocator);
	resultLineList=lineBuilder.build(opCode);

	PointBuilder pointBuilder(this,geomFact,&ptLocator);
	resultPointList=pointBuilder.build(opCode);

	// gather the results from all calculations into a single
	// Geometry for the result set
	resultGeom=computeGeometry(resultPointList,resultLineList,resultPolyList);

	checkObviouslyWrongResult(opCode);


#if USE_ELEVATION_MATRIX
	elevationMatrix->elevate(resultGeom);
#endif // USE_ELEVATION_MATRIX
	
}

/*protected*/
void
OverlayOp::insertUniqueEdge(Edge *e)
{
	//Debug.println(e);
#if GEOS_DEBUG
	cerr<<"OverlayOp::insertUniqueEdge("<<e->print()<<")"<<endl;
#endif

	//<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)
	{

#if GEOS_DEBUG
		cerr<<"  found identical edge, should merge Z"<<endl;
#endif
		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))
		{
//			labelToMerge=new Label(e->getLabel());
			labelToMerge->flip();
		}
		Depth &depth=existingEdge->getDepth();

		// if this is the first duplicate found for this edge,
		// initialize the depths
		if (depth.isNull())
		{
			depth.add(*existingLabel);
		}

		depth.add(*labelToMerge);

		existingLabel->merge(*labelToMerge);
		//Debug.print("inserted edge: "); Debug.println(e);
		//Debug.print("existing edge: "); Debug.println(existingEdge);
		dupEdges.push_back(e);
	}
	else
	{  // no matching existing edge was found
#if GEOS_DEBUG
		cerr<<"  no matching existing edge"<<endl;
#endif
		// add this new edge to the list of edges in this graph
		//e.setName(name+edges.size());
		//e.getDepth().add(e.getLabel());
		edgeList.add(e);
	}
}

/*private*/
void
OverlayOp::computeLabelsFromDepths()
{
	for(size_t j=0, s=edgeList.getEdges().size(); j<s; ++j)
	{
		Edge *e=edgeList.get(j);
		Label *lbl=e->getLabel();
		Depth &depth=e->getDepth();

		/*
		 * Only check edges for which there were duplicates,
		 * since these are the only ones which might
		 * be the result of dimensional collapses.
		 */
		if (depth.isNull()) continue;

		depth.normalize();
		for (int i=0;i<2;i++)
		{
			if (!lbl->isNull(i) && lbl->isArea() && !depth.isNull(i))
			{
				/*
				 * if the depths are equal, this edge is the result of
				 * the dimensional collapse of two or more edges.
				 * It has the same location on both sides of the edge,
				 * so it has collapsed to a line.
				 */
				if (depth.getDelta(i)==0) {
					lbl->toLine(i);
				} else {
					/*
					 * This edge may be the result of a dimensional collapse,
					 * but it still has different locations on both sides.  The
					 * label of the edge must be updated to reflect the resultant
					 * side locations indicated by the depth values.
					 */
					assert(!depth.isNull(i,Position::LEFT)); // depth of LEFT side has not been initialized
					lbl->setLocation(i,Position::LEFT,depth.getLocation(i,Position::LEFT));
					assert(!depth.isNull(i,Position::RIGHT)); // depth of RIGHT side has not been initialized
					lbl->setLocation(i,Position::RIGHT,depth.getLocation(i,Position::RIGHT));
				}
			}
		}
	}
}

struct PointCoveredByAny: public geom::CoordinateFilter
{
	const vector<const Geometry*>& geoms;
	PointLocator locator;
	PointCoveredByAny(const vector<const Geometry*>& nGeoms)
		: geoms(nGeoms)
	{}

	void filter_ro(const Coordinate* coord)
	{
		for (size_t i=0, n=geoms.size(); i<n; ++i)
		{
			int loc = locator.locate(*coord, geoms[i]);
			if ( loc == Location::INTERIOR ||
				loc == Location::BOUNDARY )
			{
				return;
			}
		}
		throw util::TopologyException("Obviously wrong result: "
			"A point on first geom boundary isn't covered "
			"by either result or second geom");
	}
};

void
OverlayOp::checkObviouslyWrongResult(OverlayOp::OpCode opCode)
{
	using std::cerr;
	using std::endl;

	assert(resultGeom);


#ifdef ENABLE_OTHER_OVERLAY_RESULT_VALIDATORS

	if ( opCode == opINTERSECTION
		&& arg[0]->getGeometry()->getDimension() == Dimension::A
		&& arg[1]->getGeometry()->getDimension() == Dimension::A )
	{
		// Area of result must be less or equal area of smaller geom
		double area0 = arg[0]->getGeometry()->getArea();
		double area1 = arg[1]->getGeometry()->getArea();
		double minarea = min(area0, area1);
		double resultArea = resultGeom->getArea();

		if ( resultArea > minarea )
		{
			throw util::TopologyException("Obviously wrong result: "
				"area of intersection "
				"result is bigger then minimum area between "
				"input geometries");
		}
	}

	else if ( opCode == opDIFFERENCE 
		&& arg[0]->getGeometry()->getDimension() == Dimension::A
		&& arg[1]->getGeometry()->getDimension() == Dimension::A )
	{
		// Area of result must be less or equal area of first geom
		double area0 = arg[0]->getGeometry()->getArea();
		double resultArea = resultGeom->getArea();

		if ( resultArea > area0 )
		{
			throw util::TopologyException("Obviously wrong result: "
				"area of difference "
				"result is bigger then area of first "
				"input geometry");
		}

		// less obvious check:
		// each vertex in first geom must be either covered by
		// result or second geom
		vector<const Geometry*> testGeoms; testGeoms.reserve(2);
		testGeoms.push_back(resultGeom);
		testGeoms.push_back(arg[1]->getGeometry());
		PointCoveredByAny tester(testGeoms);
		arg[0]->getGeometry()->apply_ro(&tester);
	}

	// Add your tests here
#endif

#ifdef ENABLE_OVERLAY_RESULT_VALIDATOR
	// This only work for FLOATING precision
	if ( resultPrecisionModel->isFloating() )
	{
		OverlayResultValidator validator( *(arg[0]->getGeometry()),
			*(arg[1]->getGeometry()), *(resultGeom));
		bool isvalid = validator.isValid(opCode);
		if ( ! isvalid )
		{
#ifdef GEOS_DEBUG_VALIDATION
			cout << "OverlayResultValidator considered result INVALID" << endl;
#endif
			throw util::TopologyException(
				"Obviously wrong result: "
				"OverlayResultValidator didn't like "
				"the result: \n"
				"Invalid point: " +
				validator.getInvalidLocation().toString() +
				string("\nInvalid result: ") +
				resultGeom->toString());
		}
#ifdef GEOS_DEBUG_VALIDATION
		else
		{
			cout << "OverlayResultValidator considered result valid" << endl;
		}
#endif
	}
#ifdef GEOS_DEBUG_VALIDATION
	else
	{
		cout << "Did not run OverlayResultValidator as the precision model is not floating" << endl;
	}
#endif // ndef GEOS_DEBUG
#endif
}

} // namespace geos.operation.overlay
} // namespace geos.operation
} // namespace geos

/**********************************************************************
 * $Log$
 * Revision 1.77  2006/07/05 20:19:29  strk
 * added checks for obviously wrong result of difference and intersection ops
 *
 * Revision 1.76  2006/06/14 15:38:16  strk
 * Fixed just-introduced bug
 *
 * Revision 1.75  2006/06/14 15:03:32  strk
 * * source/operation/overlay/OverlayOp.cpp: use NodeMap::container and related typedefs, removed (int) casts, optimized loops.
 *
 * Revision 1.74  2006/06/13 21:42:55  strk
 * trimmed cvs log, cleanups
 *
 * Revision 1.73  2006/06/12 11:29:24  strk
 * unsigned int => size_t
 *
 * Revision 1.72  2006/06/09 07:42:13  strk
 * * source/geomgraph/GeometryGraph.cpp, source/operation/buffer/OffsetCurveSetBuilder.cpp, source/operation/overlay/OverlayOp.cpp, source/operation/valid/RepeatedPointTester.cpp: Fixed warning after Polygon ring accessor methods changed to work with size_t. Small optimizations in loops.
 *
 * Revision 1.71  2006/06/05 15:36:34  strk
 * Given OverlayOp funx code enum a name and renamed values to have a lowercase prefix. Drop all of noding headers from installed header set.
 *
 *
 **********************************************************************/