[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