[geos-commits] r2348 - in trunk/source: headers/geos/operation/buffer operation/buffer

svn_geos at osgeo.org svn_geos at osgeo.org
Tue Apr 14 04:52:30 EDT 2009


Author: strk
Date: 2009-04-14 04:52:30 -0400 (Tue, 14 Apr 2009)
New Revision: 2348

Modified:
   trunk/source/headers/geos/operation/buffer/OffsetCurveBuilder.h
   trunk/source/operation/buffer/OffsetCurveBuilder.cpp
Log:
Port OffsetCurveBuilder from JTS-1.9. Twenty time faster completion of fme.xml testcase !


Modified: trunk/source/headers/geos/operation/buffer/OffsetCurveBuilder.h
===================================================================
--- trunk/source/headers/geos/operation/buffer/OffsetCurveBuilder.h	2009-04-10 20:22:27 UTC (rev 2347)
+++ trunk/source/headers/geos/operation/buffer/OffsetCurveBuilder.h	2009-04-14 08:52:30 UTC (rev 2348)
@@ -4,6 +4,7 @@
  * GEOS - Geometry Engine Open Source
  * http://geos.refractions.net
  *
+ * Copyright (C) 2009  Sandro Santilli <strk at keybit.net>
  * Copyright (C) 2006-2007 Refractions Research Inc.
  *
  * This is free software; you can redistribute and/or modify it under
@@ -13,7 +14,7 @@
  *
  **********************************************************************
  *
- * Last port: operation/buffer/OffsetCurveBuilder.java rev 1.10
+ * Last port: operation/buffer/OffsetCurveBuilder.java rev. 1.30 (JTS-1.9)
  *
  **********************************************************************/
 
@@ -62,14 +63,6 @@
  */
 class OffsetCurveBuilder {
 public:
-	/** \brief
-	 * The default number of facets into which to divide a fillet
-	 * of 90 degrees.
-	 *
-	 * A value of 8 gives less than 2% max error in the buffer distance.
-	 * For a max error of < 1%, use QS = 12
-	 */
-	static const int DEFAULT_QUADRANT_SEGMENTS=8;
 
 	/*
 	 * @param nBufParams buffer parameters, this object will
@@ -91,8 +84,9 @@
 	 * @param lineList the std::vector to which CoordinateSequences will
 	 *                 be pushed_back
 	 */
-	void getLineCurve(const geom::CoordinateSequence* inputPts, double distance,
-		std::vector<geom::CoordinateSequence*>& lineList);
+	void getLineCurve(const geom::CoordinateSequence* inputPts,
+	                  double distance,
+	                  std::vector<geom::CoordinateSequence*>& lineList);
 
 	/**
 	 * This method handles the degenerate cases of single points and lines,
@@ -102,17 +96,76 @@
 	 *                 be pushed_back
 	 */
 	void getRingCurve(const geom::CoordinateSequence *inputPts, int side,
-		double distance, std::vector<geom::CoordinateSequence*>& lineList);
+	                  double distance,
+	                  std::vector<geom::CoordinateSequence*>& lineList);
 
 
 private:
 
-	static const double MIN_CURVE_VERTEX_FACTOR; //  1.0E-6;
+	/// The mitre will be beveled if it exceeds the mitre ratio limit.
+	//
+	/// @param offset0 the first offset segment
+	/// @param offset1 the second offset segment
+	/// @param distance the offset distance
+	///
+	void addMitreJoin(const geom::Coordinate& p,
+	                  const geom::LineSegment& offset0,
+	                  const geom::LineSegment& offset1,
+	                  double distance);
 
+	/// Adds a limited mitre join connecting the two reflex offset segments.
+	//
+	/// A limited mitre is a mitre which is beveled at the distance
+	/// determined by the mitre ratio limit.
+	///
+	/// @param offset0 the first offset segment
+	/// @param offset1 the second offset segment
+	/// @param distance the offset distance
+	/// @param mitreLimit the mitre limit ratio
+	///
+	void addLimitedMitreJoin(
+	                  const geom::LineSegment& offset0,
+	                  const geom::LineSegment& offset1,
+	                  double distance, double mitreLimit);
+
+	/// \brief
+	/// Adds a bevel join connecting the two offset segments
+	/// around a reflex corner.
+	// 
+	/// @param offset0 the first offset segment
+	/// @param offset1 the second offset segment
+	///
+	void addBevelJoin(const geom::LineSegment& offset0,
+	                  const geom::LineSegment& offset1);
+
+   
+	/**
+	 * Factor which controls how close curve vertices can be to be snapped
+	 */
+	static const double CURVE_VERTEX_SNAP_DISTANCE_FACTOR; //  1.0E-6;
+
 	static const double PI; //  3.14159265358979
 
 	static const double MAX_CLOSING_SEG_LEN; // 3.0
 
+	/**
+	 * Factor which controls how close offset segments can be to
+	 * skip adding a filler or mitre.
+	 */
+	static const double OFFSET_SEGMENT_SEPARATION_FACTOR; // 1.0E-3;
+
+	/**
+	 * Factor which controls how close curve vertices on inside turns
+	 * can be to be snapped
+	 */
+	static const double INSIDE_TURN_VERTEX_SNAP_DISTANCE_FACTOR; // 1.0E-3;
+
+	/** \brief
+	 * Factor which determines how short closing segs can be
+	 * for round buffers
+	 */
+	static const int MAX_CLOSING_SEG_FRACTION = 80;
+
 	algorithm::LineIntersector li;
 
 	/** \brief
@@ -143,6 +196,25 @@
 
 	const BufferParameters& bufParams; 
 
+	/// The Closing Segment Factor controls how long "closing
+	/// segments" are.  Closing segments are added at the middle of
+	/// inside corners to ensure a smoother boundary for the buffer
+	/// offset curve.  In some cases (particularly for round joins
+	/// with default-or-better quantization) the closing segments
+	/// can be made quite short.  This substantially improves
+	/// performance (due to fewer intersections being created).
+	/// 
+	/// A closingSegFactor of 0 results in lines to the corner vertex.
+	/// A closingSegFactor of 1 results in lines halfway
+	/// to the corner vertex.
+	/// A closingSegFactor of 80 results in lines 1/81 of the way
+	/// to the corner vertex (this option is reasonable for the very
+	/// common default situation of round joins and quadrantSegs >= 8).
+	///
+	/// The default is 1.
+	///
+	int closingSegFactor; // 1;
+
 	geom::Coordinate s0, s1, s2;
 
 	geom::LineSegment seg0;
@@ -155,60 +227,103 @@
 
 	int side;
 
-//	static geom::CoordinateSequence* copyCoordinates(geom::CoordinateSequence *pts);
-
 	void init(double newDistance);
 
-	//geom::CoordinateSequence* getCoordinates();
+	/**
+	 * Use a value which results in a potential distance error which is
+	 * significantly less than the error due to
+	 * the quadrant segment discretization.
+	 * For QS = 8 a value of 100 is reasonable.
+	 * This should produce a maximum of 1% distance error.
+	 */
+	static const double SIMPLIFY_FACTOR; // 100.0;
 
+	/** \brief
+	 * Computes the distance tolerance to use during input
+	 * line simplification.
+	 *
+	 * @param distance the buffer distance
+	 * @return the simplification tolerance
+	 */
+	double simplifyTolerance(double bufDistance);
+
 	void computeLineBufferCurve(const geom::CoordinateSequence& inputPts);
 
-	void computeRingBufferCurve(const geom::CoordinateSequence& inputPts, int side);
+	void computeRingBufferCurve(const geom::CoordinateSequence& inputPts,
+	                            int side);
 
-	//void addPt(const geom::Coordinate &pt);
+	void initSideSegments(const geom::Coordinate &nS1,
+	                      const geom::Coordinate &nS2, int nSide);
 
-	//void closePts();
+	void addNextSegment(const geom::Coordinate &p, bool addStartPoint);
 
-	void initSideSegments(const geom::Coordinate &nS1, const geom::Coordinate &nS2, int nSide);
+	void addCollinear(bool addStartPoint);
 
-	void addNextSegment(const geom::Coordinate &p, bool addStartPoint);
+	/// Adds the offset points for an outside (convex) turn
+	//
+	/// @param orientation 
+	/// @param addStartPoint  
+	///
+	void addOutsideTurn(int orientation, bool addStartPoint);
 
+	/// Adds the offset points for an inside (concave) turn
+	//
+	/// @param orientation 
+	/// @param addStartPoint  
+	///
+	void addInsideTurn(int orientation, bool addStartPoint);
+
 	/// Add last offset point
 	void addLastSegment();
 
-	/**
-	 * Compute an offset segment for an input segment on a given side and at a
-	 * given distance.
-	 * The offset points are computed in full double precision, for accuracy.
+	/** \brief
+	 * Compute an offset segment for an input segment on a given
+	 * side and at a given distance.
 	 *
+	 * The offset points are computed in full double precision,
+	 * for accuracy.
+	 *
 	 * @param seg the segment to offset
 	 * @param side the side of the segment the offset lies on
 	 * @param distance the offset distance
 	 * @param offset the points computed for the offset segment
 	 */
-	void computeOffsetSegment(const geom::LineSegment& seg, int side, double distance,
-			geom::LineSegment& offset);
+	void computeOffsetSegment(const geom::LineSegment& seg,
+	                          int side, double distance,
+	                          geom::LineSegment& offset);
 
-	/// Add an end cap around point p1, terminating a line segment coming from p0
-	void addLineEndCap(const geom::Coordinate &p0,const geom::Coordinate &p1);
+	/// \brief
+	/// Add an end cap around point p1, terminating a line segment
+	/// coming from p0
+	void addLineEndCap(const geom::Coordinate &p0,
+	                   const geom::Coordinate &p1);
 
 	/**
+	 * Adds points for a circular fillet around a reflex corner.
+	 * 
+	 * Adds the start and end points
+	 *
 	 * @param p base point of curve
 	 * @param p0 start point of fillet curve
 	 * @param p1 endpoint of fillet curve
+	 * @param direction the orientation of the fillet
+	 * @param radius the radius of the fillet
 	 */
 	void addFillet(const geom::Coordinate &p, const geom::Coordinate &p0,
-			const geom::Coordinate &p1, int direction, double distance);
+	               const geom::Coordinate &p1,
+	               int direction, double radius);
 
-	/** \brief
-	 * Adds points for a fillet. 
+	/** 
+	 * Adds points for a circular fillet arc between two specified angles.
+	 *
 	 * The start and end point for the fillet are not added -
 	 * the caller must add them if required.
 	 *
 	 * @param direction is -1 for a CW angle, 1 for a CCW angle
+	 * @param radius the radius of the fillet
 	 */
-	void addFillet(const geom::Coordinate &p, double startAngle, double endAngle,
-			int direction, double distance);
+	void addFillet(const geom::Coordinate &p, double startAngle,
+	               double endAngle, int direction, double radius);
 
 	/// Adds a CW circle around a point
 	void addCircle(const geom::Coordinate &p, double distance);

Modified: trunk/source/operation/buffer/OffsetCurveBuilder.cpp
===================================================================
--- trunk/source/operation/buffer/OffsetCurveBuilder.cpp	2009-04-10 20:22:27 UTC (rev 2347)
+++ trunk/source/operation/buffer/OffsetCurveBuilder.cpp	2009-04-14 08:52:30 UTC (rev 2348)
@@ -4,8 +4,9 @@
  * GEOS-Geometry Engine Open Source
  * http://geos.refractions.net
  *
+ * Copyright (C) 2009  Sandro Santilli <strk at keybit.net>
+ * Copyright (C) 2005 Refractions Research Inc.
  * Copyright (C) 2001-2002 Vivid Solutions Inc.
- * Copyright (C) 2005 Refractions Research Inc.
  *
  * This is free software; you can redistribute and/or modify it under
  * the terms of the GNU Lesser General Public Licence as published
@@ -14,7 +15,7 @@
  *
  **********************************************************************
  *
- * Last port: operation/buffer/OffsetCurveBuilder.java rev. 1.9
+ * Last port: operation/buffer/OffsetCurveBuilder.java rev. 1.30 (JTS-1.9)
  *
  **********************************************************************/
 
@@ -23,6 +24,7 @@
 #include <vector>
 
 #include <geos/algorithm/CGAlgorithms.h>
+#include <geos/algorithm/Angle.h>
 #include <geos/operation/buffer/OffsetCurveBuilder.h>
 #include <geos/operation/buffer/BufferOp.h>
 #include <geos/operation/buffer/BufferParameters.h>
@@ -31,13 +33,19 @@
 #include <geos/geom/CoordinateSequence.h>
 #include <geos/geom/Coordinate.h>
 #include <geos/geom/PrecisionModel.h>
+#include <geos/algorithm/NotRepresentableException.h>
+#include <geos/algorithm/HCoordinate.h>
 
 #include "OffsetCurveVertexList.h"
+#include "BufferInputLineSimplifier.h"
 
 #ifndef GEOS_DEBUG
 #define GEOS_DEBUG 0
 #endif
 
+// Define this to skip the input simplification step (for testing)
+//#define SKIP_INPUT_SIMPLIFICATION
+
 using namespace std;
 using namespace geos::geomgraph;
 using namespace geos::algorithm;
@@ -48,9 +56,12 @@
 namespace buffer { // geos.operation.buffer
 
 /*private data*/
-const double OffsetCurveBuilder::MIN_CURVE_VERTEX_FACTOR = 1.0E-6;
+const double OffsetCurveBuilder::CURVE_VERTEX_SNAP_DISTANCE_FACTOR = 1.0E-6;
 const double OffsetCurveBuilder::PI = 3.14159265358979;
 const double OffsetCurveBuilder::MAX_CLOSING_SEG_LEN = 3.0;
+const double OffsetCurveBuilder::OFFSET_SEGMENT_SEPARATION_FACTOR = 1.0E-3;
+const double OffsetCurveBuilder::INSIDE_TURN_VERTEX_SNAP_DISTANCE_FACTOR = 1.0E-3;
+const double OffsetCurveBuilder::SIMPLIFY_FACTOR = 100.0;
 
 /*public*/
 OffsetCurveBuilder::OffsetCurveBuilder(const PrecisionModel *newPrecisionModel,
@@ -62,6 +73,7 @@
 		distance(0.0),
 		precisionModel(newPrecisionModel),
 		bufParams(nBufParams),
+		closingSegFactor(1),
 		s0(),
 		s1(),
 		s2(),
@@ -72,7 +84,20 @@
 		side(0),
 		vertexLists()
 {
-	filletAngleQuantum=PI / 2.0 / bufParams.getQuadrantSegments();
+	// compute intersections in full precision, to provide accuracy
+	// the points are rounded as they are inserted into the curve line
+	filletAngleQuantum = PI / 2.0 / bufParams.getQuadrantSegments();
+
+	/**
+	 * Non-round joins cause issues with short closing segments,
+	 * so don't use them.  In any case, non-round joins
+	 * only really make sense for relatively small buffer distances.
+	 */
+	if (bufParams.getQuadrantSegments() >= 8
+	     && bufParams.getJoinStyle() == BufferParameters::JOIN_ROUND)
+	{
+		closingSegFactor = MAX_CLOSING_SEG_FRACTION;
+	}
 }
 
 /*public*/
@@ -93,7 +118,7 @@
 
 	init(distance);
 
-	if (inputPts->getSize() < 2) {
+	if (inputPts->getSize() <= 1) {
 		switch (bufParams.getEndCapStyle()) {
 			case BufferParameters::CAP_ROUND:
 				addCircle(inputPts->getAt(0), distance);
@@ -109,8 +134,11 @@
 	} else {
 		computeLineBufferCurve(*inputPts);
 	}
-	// NOTE: we take ownership of lineCoord here
+
+	// NOTE: we take ownership of lineCoord here ...
 	CoordinateSequence *lineCoord=vertexList->getCoordinates();
+
+	// ... and we give it away here
 	lineList.push_back(lineCoord);
 }
 
@@ -121,18 +149,20 @@
 		vector<CoordinateSequence*>& lineList)
 {
 	init(distance);
-	if (inputPts->getSize()<= 2)
+	if (inputPts->getSize() <= 2)
 	{
 		getLineCurve(inputPts, distance, lineList);
 		return;
 	}
-	// optimize creating ring for for zero distance
-	if (distance==0.0) {
+
+	// optimize creating ring for zero distance
+	if (distance == 0.0) {
 		vertexLists.push_back(vertexList);
 		vertexList = new OffsetCurveVertexList(); // is this needed ?
 		lineList.push_back(inputPts->clone());
 		return;
 	}
+
 	computeRingBufferCurve(*inputPts, side);
 
 	// this will be vertexList
@@ -144,41 +174,79 @@
 void
 OffsetCurveBuilder::init(double newDistance)
 {
-	distance=newDistance;
-	maxCurveSegmentError=distance*(1-cos(filletAngleQuantum/2.0));
+	distance = newDistance;
+	maxCurveSegmentError = distance * (1 - cos(filletAngleQuantum/2.0));
+
 	// Point list needs to be reset
 	// but if a previous point list exists
 	// we'd better back it up for final deletion
 	vertexLists.push_back(vertexList);
 	vertexList=new OffsetCurveVertexList();
 	vertexList->setPrecisionModel(precisionModel);
+
 	/**
-	 * Choose the min vertex separation as a small fraction of the offset distance.
+	 * Choose the min vertex separation as a small fraction of
+	 * the offset distance.
 	 */
-	vertexList->setMinimumVertexDistance(distance * MIN_CURVE_VERTEX_FACTOR);
+	vertexList->setMinimumVertexDistance(
+		distance * CURVE_VERTEX_SNAP_DISTANCE_FACTOR);
 }
 
+/* private */
+double
+OffsetCurveBuilder::simplifyTolerance(double bufDistance)
+{
+	return bufDistance / SIMPLIFY_FACTOR;
+}
+
 /*private*/
 void
 OffsetCurveBuilder::computeLineBufferCurve(const CoordinateSequence& inputPts)
 {
-	int n=inputPts.size()-1;
-	// compute points for left side of line
-	initSideSegments(inputPts[0], inputPts[1], Position::LEFT);
-	for (int i=2;i<= n;i++) {
-		addNextSegment(inputPts[i], true);
+	double distTol = simplifyTolerance(distance);
+
+	//--------- compute points for left side of line
+#ifndef SKIP_INPUT_SIMPLIFICATION
+	// Simplify the appropriate side of the line before generating
+	std::auto_ptr<CoordinateSequence> simp1_ = 
+		BufferInputLineSimplifier::simplify(inputPts, distTol);
+	const CoordinateSequence& simp1 = *simp1_;
+#else
+	// MD - used for testing only (to eliminate simplification)
+	const CoordinateSequence& simp1 = inputPts;
+#endif
+
+
+	int n1 = simp1.size() - 1;
+	initSideSegments(simp1[0], simp1[1], Position::LEFT);
+	for (int i = 2; i <= n1; ++i) {
+		addNextSegment(simp1[i], true);
 	}
 	addLastSegment();
 	// add line cap for end of line
-	addLineEndCap(inputPts[n-1], inputPts[n]);
-	// compute points for right side of line
-	initSideSegments(inputPts[n], inputPts[n-1], Position::LEFT);
-	for (int i=n-2; i>=0; i--) {
-		addNextSegment(inputPts[i], true);
+	addLineEndCap(simp1[n1-1], simp1[n1]);
+
+
+	//---------- compute points for right side of line
+#ifndef SKIP_INPUT_SIMPLIFICATION
+	// Simplify the appropriate side of the line before generating
+	std::auto_ptr<CoordinateSequence> simp2_ = 
+		BufferInputLineSimplifier::simplify(inputPts, -distTol);
+	const CoordinateSequence& simp2 = *simp2_;
+#else
+	// MD - used for testing only (to eliminate simplification)
+	const CoordinateSequence& simp2 = inputPts;
+#endif
+
+	int n2 = simp2.size() - 1;
+	initSideSegments(simp2[n2], simp2[n2-1], Position::LEFT);
+	for (int i = n2-2; i >= 0; --i) {
+		addNextSegment(simp2[i], true);
 	}
 	addLastSegment();
 	// add line cap for start of line
-	addLineEndCap(inputPts[1], inputPts[0]);
+	addLineEndCap(simp2[1], simp2[0]);
+
 	vertexList->closeRing();
 }
 
@@ -187,18 +255,32 @@
 OffsetCurveBuilder::computeRingBufferCurve(const CoordinateSequence& inputPts,
 	int side)
 {
-	int n=inputPts.size()-1;
-	initSideSegments(inputPts[n-1], inputPts[0], side);
-	for (int i=1; i<=n; i++) {
-		bool addStartPoint=i != 1;
-		addNextSegment(inputPts[i], addStartPoint);
+
+#ifndef SKIP_INPUT_SIMPLIFICATION
+	double distTol = simplifyTolerance(distance);
+	// ensure that correct side is simplified
+	if (side == Position::RIGHT)
+		distTol = -distTol;      
+	std::auto_ptr<CoordinateSequence> simp_ = 
+		BufferInputLineSimplifier::simplify(inputPts, distTol);
+	const CoordinateSequence& simp = *simp_;
+#else
+	const CoordinateSequence& simp = inputPts;
+#endif
+
+	int n = simp.size()-1;
+	initSideSegments(simp[n-1], simp[0], side);
+	for (int i = 1; i <= n; i++) {
+		bool addStartPoint = i != 1;
+		addNextSegment(simp[i], addStartPoint);
 	}
 	vertexList->closeRing();
 }
 
 /*private*/
 void
-OffsetCurveBuilder::initSideSegments(const Coordinate &nS1, const Coordinate &nS2, int nSide)
+OffsetCurveBuilder::initSideSegments(const Coordinate &nS1,
+		const Coordinate &nS2, int nSide)
 {
 	s1=nS1;
 	s2=nS2;
@@ -211,7 +293,8 @@
 void
 OffsetCurveBuilder::addNextSegment(const Coordinate &p, bool addStartPoint)
 {
-	// s0-s1-s2 are the coordinates of the previous segment and the current one
+	// s0-s1-s2 are the coordinates of the previous segment
+	// and the current one
 	s0=s1;
 	s1=s2;
 	s2=p;
@@ -222,72 +305,28 @@
 
 	// do nothing if points are equal
 	if (s1==s2) return;
+
 	int orientation=CGAlgorithms::computeOrientation(s0, s1, s2);
-	bool outsideTurn =(orientation==CGAlgorithms::CLOCKWISE
-						&& side==Position::LEFT)
-						||(orientation==CGAlgorithms::COUNTERCLOCKWISE 
-						&& side==Position::RIGHT);
+	bool outsideTurn =
+		(orientation==CGAlgorithms::CLOCKWISE
+		 && side==Position::LEFT)
+		||
+		(orientation==CGAlgorithms::COUNTERCLOCKWISE 
+		 && side==Position::RIGHT);
+
 	if (orientation==0)
-	{ // lines are collinear
-		li.computeIntersection(s0,s1,s1,s2);
-		int numInt=li.getIntersectionNum();
-
-		/**
-		 * if numInt is<2, the lines are parallel and in the same direction.
-		 * In this case the point can be ignored, since the offset lines will also be
-		 * parallel.
-		 */
-		if (numInt>= 2) {
-			/**
-			 * segments are collinear but reversing.  Have to add an "end-cap" fillet
-			 * all the way around to other direction
-			 * This case should ONLY happen for LineStrings, so the orientation is always CW.
-			 * Polygons can never have two consecutive segments which are parallel but
-			 * reversed, because that would be a self intersection.
-			 */
-			addFillet(s1, offset0.p1, offset1.p0, CGAlgorithms::CLOCKWISE, distance);
-		}
+	{
+		// lines are collinear
+		addCollinear(addStartPoint);
 	}
 	else if (outsideTurn)
 	{
-		// add a fillet to connect the endpoints of the offset segments
-		if (addStartPoint) vertexList->addPt(offset0.p1);
-		// TESTING-comment out to produce beveled joins
-		addFillet(s1, offset0.p1, offset1.p0, orientation, distance);
-		vertexList->addPt(offset1.p0);
+		addOutsideTurn(orientation, addStartPoint);
 	}
 	else
-	{ // inside turn
-
-		// add intersection point of offset segments (if any)
-		li.computeIntersection( offset0.p0, offset0.p1, offset1.p0, offset1.p1);
-		if (li.hasIntersection()) {
-			vertexList->addPt(li.getIntersection(0));
-		} else {
-			/**
-			 * If no intersection, it means the angle is so small and/or the offset so large
-			 * that the offsets segments don't intersect.
-			 * In this case we must add a offset joining curve to make sure the buffer line
-			 * is continuous and tracks the buffer correctly around the corner.
-			 * Note that the joining curve won't appear in the final buffer.
-			 *
-			 * The intersection test above is vulnerable to robustness errors;
-			 * i.e. it may be that the offsets should intersect very close to their
-			 * endpoints, but don't due to rounding.  To handle this situation
-			 * appropriately, we use the following test:
-			 * If the offset points are very close, don't add a joining curve
-			 * but simply use one of the offset points
-			 */
-			if (offset0.p1.distance(offset1.p0)<distance / 1000.0) {
-				vertexList->addPt(offset0.p1);
-			} else {
-				// add endpoint of this segment offset
-				vertexList->addPt(offset0.p1);
-				// <FIX> MD-add in centre point of corner, to make sure offset closer lines have correct topology
-				vertexList->addPt(s1);
-				vertexList->addPt(offset1.p0);
-			}
-		}
+	{
+		// inside turn
+		addInsideTurn(orientation, addStartPoint);
 	}
 }
 
@@ -307,7 +346,8 @@
 	double dx = seg.p1.x - seg.p0.x;
 	double dy = seg.p1.y - seg.p0.y;
 	double len = sqrt(dx * dx + dy * dy);
-	// u is the vector that is the length of the offset, in the direction of the segment
+	// u is the vector that is the length of the offset,
+	// in the direction of the segment
 	double ux = sideSign * distance * dx / len;
 	double uy = sideSign * distance * dy / len;
 	offset.p0.x = seg.p0.x - uy;
@@ -321,18 +361,22 @@
 OffsetCurveBuilder::addLineEndCap(const Coordinate &p0,const Coordinate &p1)
 {
 	LineSegment seg(p0, p1);
+
 	LineSegment offsetL;
 	computeOffsetSegment(seg, Position::LEFT, distance, offsetL);
 	LineSegment offsetR;
 	computeOffsetSegment(seg, Position::RIGHT, distance, offsetR);
+
 	double dx=p1.x-p0.x;
 	double dy=p1.y-p0.y;
 	double angle=atan2(dy, dx);
+
 	switch (bufParams.getEndCapStyle()) {
 		case BufferParameters::CAP_ROUND:
 			// add offset seg points with a fillet between them
 			vertexList->addPt(offsetL.p1);
-			addFillet(p1, angle+PI/2.0, angle-PI/2.0, CGAlgorithms::CLOCKWISE, distance);
+			addFillet(p1, angle+PI/2.0, angle-PI/2.0,
+			          CGAlgorithms::CLOCKWISE, distance);
 			vertexList->addPt(offsetR.p1);
 			break;
 		case BufferParameters::CAP_FLAT:
@@ -341,14 +385,18 @@
 			vertexList->addPt(offsetR.p1);
 			break;
 		case BufferParameters::CAP_SQUARE:
-			// add a square defined by extensions of the offset segment endpoints
+			// add a square defined by extensions of the offset
+			// segment endpoints
 			Coordinate squareCapSideOffset;
 			squareCapSideOffset.x=fabs(distance)*cos(angle);
 			squareCapSideOffset.y=fabs(distance)*sin(angle);
-			Coordinate squareCapLOffset(offsetL.p1.x+squareCapSideOffset.x,
-					offsetL.p1.y+squareCapSideOffset.y);
-			Coordinate squareCapROffset(offsetR.p1.x+squareCapSideOffset.x,
-					offsetR.p1.y+squareCapSideOffset.y);
+
+			Coordinate squareCapLOffset(
+				offsetL.p1.x+squareCapSideOffset.x,
+				offsetL.p1.y+squareCapSideOffset.y);
+			Coordinate squareCapROffset(
+				offsetR.p1.x+squareCapSideOffset.x,
+				offsetR.p1.y+squareCapSideOffset.y);
 			vertexList->addPt(squareCapLOffset);
 			vertexList->addPt(squareCapROffset);
 			break;
@@ -358,43 +406,52 @@
 /*private*/
 void
 OffsetCurveBuilder::addFillet(const Coordinate &p, const Coordinate &p0,
-	const Coordinate &p1, int direction, double distance)
+	const Coordinate &p1, int direction, double radius)
 {
-	double dx0=p0.x-p.x;
-	double dy0=p0.y-p.y;
-	double startAngle=atan2(dy0, dx0);
-	double dx1=p1.x-p.x;
-	double dy1=p1.y-p.y;
-	double endAngle=atan2(dy1, dx1);
-	if (direction==CGAlgorithms::CLOCKWISE) {
-		if (startAngle<= endAngle) startAngle += 2.0*PI;
-	} else {    // direction==COUNTERCLOCKWISE
-		if (startAngle>=endAngle) startAngle-=2.0*PI;
+	double dx0 = p0.x - p.x;
+	double dy0 = p0.y - p.y;
+	double startAngle = atan2(dy0, dx0);
+	double dx1 = p1.x - p.x;
+	double dy1 = p1.y - p.y;
+	double endAngle = atan2(dy1, dx1);
+
+	if (direction == CGAlgorithms::CLOCKWISE) {
+		if (startAngle <= endAngle) startAngle += 2.0 * PI;
 	}
+	else {    // direction==COUNTERCLOCKWISE
+		if (startAngle >= endAngle) startAngle -= 2.0 * PI;
+	}
+
 	vertexList->addPt(p0);
-	addFillet(p, startAngle, endAngle, direction, distance);
+	addFillet(p, startAngle, endAngle, direction, radius);
 	vertexList->addPt(p1);
 }
 
 /*private*/
 void
 OffsetCurveBuilder::addFillet(const Coordinate &p, double startAngle,
-	double endAngle, int direction, double distance)
+	double endAngle, int direction, double radius)
 {
-	int directionFactor=direction==CGAlgorithms::CLOCKWISE ? -1 : 1;
-	double totalAngle=fabs(startAngle-endAngle);
-	int nSegs=(int) (totalAngle / filletAngleQuantum+0.5);
-	if (nSegs<1) return;   // no segments because angle is less than increment-nothing to do!
+	int directionFactor = direction == CGAlgorithms::CLOCKWISE ? -1 : 1;
+
+	double totalAngle = fabs(startAngle - endAngle);
+	int nSegs = (int) (totalAngle / filletAngleQuantum + 0.5);
+
+	// no segments because angle is less than increment-nothing to do!
+	if (nSegs<1) return;
+
 	double initAngle, currAngleInc;
+
 	// choose angle increment so that each segment has equal length
-	initAngle=0.0;
-	currAngleInc=totalAngle / nSegs;
-	double currAngle=initAngle;
+	initAngle = 0.0;
+	currAngleInc = totalAngle / nSegs;
+
+	double currAngle = initAngle;
 	Coordinate pt;
-	while (currAngle<totalAngle) {
-		double angle=startAngle+directionFactor*currAngle;
-		pt.x=p.x+distance*cos(angle);
-		pt.y=p.y+distance*sin(angle);
+	while (currAngle < totalAngle) {
+		double angle = startAngle + directionFactor * currAngle;
+		pt.x = p.x + radius * cos(angle);
+		pt.y = p.y + radius * sin(angle);
 		vertexList->addPt(pt);
 		currAngle += currAngleInc;
 	}
@@ -406,8 +463,7 @@
 OffsetCurveBuilder::addCircle(const Coordinate &p, double distance)
 {
 	// add start point
-	Coordinate pt(p);
-	pt.x+=distance;
+	Coordinate pt(p.x + distance, p.y);
 	vertexList->addPt(pt);
 	addFillet(p, 0.0, 2.0*PI, -1, distance);
 }
@@ -424,6 +480,279 @@
 	vertexList->addPt(Coordinate(p.x+distance, p.y+distance));
 }
 
+/* private */
+void
+OffsetCurveBuilder::addCollinear(bool addStartPoint)
+{
+	/**
+	 * This test could probably be done more efficiently,
+	 * but the situation of exact collinearity should be fairly rare.
+	 */
+
+	li.computeIntersection(s0,s1,s1,s2);
+	int numInt=li.getIntersectionNum();
+
+	/**
+	 * if numInt is<2, the lines are parallel and in the same direction.
+	 * In this case the point can be ignored, since the offset lines
+	 * will also be parallel.
+	 */
+	if (numInt>= 2)
+	{
+		/**
+		 * Segments are collinear but reversing. 
+		 * Add an "end-cap" fillet
+		 * all the way around to other direction
+		 *
+		 * This case should ONLY happen for LineStrings,
+		 * so the orientation is always CW (Polygons can never
+		 * have two consecutive segments which are parallel but
+		 * reversed, because that would be a self intersection).
+		 */
+		if (  bufParams.getJoinStyle() == BufferParameters::JOIN_BEVEL
+		   || bufParams.getJoinStyle() == BufferParameters::JOIN_MITRE)
+		{
+			if (addStartPoint) vertexList->addPt(offset0.p1);
+			vertexList->addPt(offset1.p0);
+		}
+		else
+		{
+			addFillet(s1, offset0.p1, offset1.p0,
+				  CGAlgorithms::CLOCKWISE, distance);
+		}
+	}
+}
+
+/* private */
+void
+OffsetCurveBuilder::addOutsideTurn(int orientation, bool addStartPoint)
+{
+	/**
+	 * Heuristic: If offset endpoints are very close together,
+	 * just use one of them as the corner vertex.
+	 * This avoids problems with computing mitre corners in the case
+	 * where the two segments are almost parallel
+	 * (which is hard to compute a robust intersection for).
+	 */
+
+	if (offset0.p1.distance(offset1.p0) <
+		distance*OFFSET_SEGMENT_SEPARATION_FACTOR)
+	{
+		vertexList->addPt(offset0.p1);
+		return;
+	}
+
+	if (bufParams.getJoinStyle() == BufferParameters::JOIN_MITRE)
+	{
+		addMitreJoin(s1, offset0, offset1, distance);
+	}
+	else if (bufParams.getJoinStyle() == BufferParameters::JOIN_BEVEL)
+	{
+		addBevelJoin(offset0, offset1);
+	}
+	else
+	{
+		// add a circular fillet connecting the endpoints 
+		// of the offset segments
+		if (addStartPoint) vertexList->addPt(offset0.p1);
+
+		// TESTING - comment out to produce beveled joins
+		addFillet(s1, offset0.p1, offset1.p0, orientation, distance);
+		vertexList->addPt(offset1.p0);
+	}
+}
+
+/* private */
+void
+OffsetCurveBuilder::addInsideTurn(int orientation, bool addStartPoint)
+{
+	// add intersection point of offset segments (if any)
+	li.computeIntersection(offset0.p0, offset0.p1, offset1.p0, offset1.p1);
+	if (li.hasIntersection())
+	{
+		vertexList->addPt(li.getIntersection(0));
+		return;
+	}
+
+	// If no intersection is detected, it means the angle is so small
+	// and/or the offset so large that the offsets segments don't
+	// intersect. In this case we must add a "closing segment" to make
+	// sure the buffer curve is continuous,
+	// fairly smooth (e.g. no sharp reversals in direction)
+	// and tracks the buffer correctly around the corner.
+	// The curve connects the endpoints of the segment offsets to points
+	// which lie toward the centre point of the corner.
+	// The joining curve will not appear in the final buffer outline,
+	// since it is completely internal to the buffer polygon.
+	//
+	// In complex buffer cases the closing segment may cut across many
+	// other segments in the generated offset curve. 
+	// In order to improve the performance of the noding, the closing
+	// segment should be kept as short as possible.
+	// (But not too short, since that would defeat it's purpose).
+	// This is the purpose of the closingSegFactor heuristic value.
+
+       /**
+        * The intersection test above is vulnerable to robustness errors;
+	* i.e. it may be that the offsets should intersect very close to
+	* their endpoints, but aren't reported as such due to rounding.
+	* To handle this situation appropriately, we use the following test:
+	* If the offset points are very close, don't add closing segments
+	* but simply use one of the offset points
+        */
+
+	if (offset0.p1.distance(offset1.p0) <
+		distance * INSIDE_TURN_VERTEX_SNAP_DISTANCE_FACTOR)
+	{
+		vertexList->addPt(offset0.p1);
+	}
+	else
+	{
+		// add endpoint of this segment offset
+		vertexList->addPt(offset0.p1);
+
+		// Add "closing segment" of required length.
+		if ( closingSegFactor > 0 )
+		{
+			Coordinate mid0(
+ (closingSegFactor*offset0.p1.x + s1.x)/(closingSegFactor + 1),
+ (closingSegFactor*offset0.p1.y + s1.y)/(closingSegFactor + 1)
+			);
+			vertexList->addPt(mid0);
+
+			Coordinate mid1(
+ (closingSegFactor*offset1.p0.x + s1.x)/(closingSegFactor + 1),
+ (closingSegFactor*offset1.p0.y + s1.y)/(closingSegFactor + 1)
+			);
+			vertexList->addPt(mid1);
+		}
+		else
+		{
+			// This branch is not expected to be used
+			// except for testing purposes.
+			// It is equivalent to the JTS 1.9 logic for
+			// closing segments (which results in very poor
+			// performance for large buffer distances)
+			vertexList->addPt(s1);
+		}
+
+		// add start point of next segment offset
+		vertexList->addPt(offset1.p0);
+	}
+}
+
+/* private */
+void
+OffsetCurveBuilder::addMitreJoin(const geom::Coordinate& p,
+	                  const geom::LineSegment& offset0,
+	                  const geom::LineSegment& offset1,
+	                  double distance)
+{
+	bool isMitreWithinLimit = true;
+	Coordinate intPt;
+
+        /**
+         * This computation is unstable if the offset segments
+	 * are nearly collinear.
+         * Howver, this situation should have been eliminated earlier
+	 * by the check for whether the offset segment endpoints are
+	 * almost coincident
+         */
+	try
+	{
+		HCoordinate::intersection(offset0.p0, offset0.p1,
+		                          offset1.p0, offset1.p1,
+		                          intPt);
+
+		double mitreRatio = distance <= 0.0 ? 1.0
+			: intPt.distance(p) / fabs(distance);
+
+		if (mitreRatio > bufParams.getMitreLimit())
+			isMitreWithinLimit = false;
+	}
+	catch (const NotRepresentableException& e)
+	{
+                intPt = Coordinate(0,0);
+                isMitreWithinLimit = false;
+        }
+
+	if (isMitreWithinLimit)
+	{
+                vertexList->addPt(intPt);
+        }
+        else
+	{
+		addLimitedMitreJoin(offset0, offset1, distance,
+		                    bufParams.getMitreLimit());
+		//addBevelJoin(offset0, offset1);
+	}
+}
+
+/* private */
+void
+OffsetCurveBuilder::addLimitedMitreJoin(
+	                  const geom::LineSegment& offset0,
+	                  const geom::LineSegment& offset1,
+	                  double distance, double mitreLimit)
+{
+	const Coordinate& basePt = seg0.p1;
+
+	double ang0 = Angle::angle(basePt, seg0.p0);
+	//double ang1 = Angle::angle(basePt, seg1.p1); // unused in JTS, bug ?
+
+	// oriented angle between segments
+	double angDiff = Angle::angleBetweenOriented(seg0.p0, basePt, seg1.p1);
+	// half of the interior angle
+	double angDiffHalf = angDiff / 2;
+
+	// angle for bisector of the interior angle between the segments
+	double midAng = Angle::normalize(ang0 + angDiffHalf);
+	// rotating this by PI gives the bisector of the reflex angle
+	double mitreMidAng = Angle::normalize(midAng + PI);
+
+	// the miterLimit determines the distance to the mitre bevel
+	double mitreDist = mitreLimit * distance;
+	// the bevel delta is the difference between the buffer distance
+	// and half of the length of the bevel segment
+	double bevelDelta = mitreDist * fabs(sin(angDiffHalf));
+	double bevelHalfLen = distance - bevelDelta;
+
+	// compute the midpoint of the bevel segment
+	double bevelMidX = basePt.x + mitreDist * cos(mitreMidAng);
+	double bevelMidY = basePt.y + mitreDist * sin(mitreMidAng);
+	Coordinate bevelMidPt(bevelMidX, bevelMidY);
+
+	// compute the mitre midline segment from the corner point to
+	// the bevel segment midpoint
+	LineSegment mitreMidLine(basePt, bevelMidPt);
+
+	// finally the bevel segment endpoints are computed as offsets from
+	// the mitre midline
+	Coordinate bevelEndLeft;
+	mitreMidLine.pointAlongOffset(1.0, bevelHalfLen, bevelEndLeft);
+	Coordinate bevelEndRight;
+	mitreMidLine.pointAlongOffset(1.0, -bevelHalfLen, bevelEndRight);
+
+	if (side == Position::LEFT) {
+		vertexList->addPt(bevelEndLeft);
+		vertexList->addPt(bevelEndRight);
+	}
+	else {
+		vertexList->addPt(bevelEndRight);
+		vertexList->addPt(bevelEndLeft);
+	}
+
+}
+
+/* private */
+void
+OffsetCurveBuilder::addBevelJoin( const geom::LineSegment& offset0,
+	                  const geom::LineSegment& offset1)
+{
+	vertexList->addPt(offset0.p1);
+	vertexList->addPt(offset1.p0);
+}
+
 } // namespace geos.operation.buffer
 } // namespace geos.operation
 } // namespace geos



More information about the geos-commits mailing list