PESimplifyImpl.mesa
Copyright (C) 1984 by Xerox Corporation. All rights reserved.
Written by Pauline Ts'o on July 12, 1984 11:05:26 am PDT
Last Edited by: Tso, August 14, 1984 10:42:22 am PDT
Beta spline control point simplification routines for the Path Editor.
DIRECTORY
PEBetaSpline USING [BetaSplinePoints],
PERefresh USING [DisableTrajectoryRefresh, EnableTrajectoryRefresh],
PESimplify,
PETrajectoryOps USING [InsertSegment, InsertVertex, LinkSegments, UnlinkSegments],
PETypes,
RealFns USING [ArcTanDeg, CosDeg, Ln, SinDeg, SqRt];
PESimplifyImpl:
CEDAR
PROGRAM
IMPORTS PEBetaSpline, PERefresh, PETrajectoryOps, RealFns
EXPORTS PESimplify =
BEGIN
WeightedMidpoint:
PUBLIC
PROCEDURE [activeTrajectory: PETypes.TrajectoryNode, delta:
REAL]
RETURNS [filtered:
BOOLEAN ←
FALSE, collinearSegNode: PETypes.SegmentNode, newControlPolygon: PETypes.Trajectory, newSpline: PETypes.Trajectory, oldSegList: PETypes.SegmentList ←
NIL] = {
This routine computes new control points of a Beta spline given a list of control points with corresponding beta values (bias and tension).
bias: REAL;
center: PETypes.Point;
equiCurveSegNode: PETypes.SegmentNode;
factor1, factor2, factor3: REAL;
firstSegment: PETypes.Segment;
newSegmentList: PETypes.SegmentList;
newSplineSegmentList: PETypes.SegmentList;
numEquiCurveSegs: CARDINAL;
point: PETypes.Point;
prevSegment: PETypes.SegmentNode;
radius: REAL;
segment: PETypes.SegmentNode;
tempSegList: PETypes.SegmentList ← NIL;
tension: REAL;
vertex1, vertex2, vertex3, vertex4: PETypes.Vertex;
Test to see if any control points are collinear with its two neighbors. If so, remove one of these.
collinearSegNode ← ComputeCollinear[activeTrajectory, delta];
IF collinearSegNode # NIL THEN RETURN;
Test to see if any control points have the same curvature within a neighborhood of at least four. If so, remove one of these and reditribute the rest around the osculating circle.
[numEquiCurveSegs, equiCurveSegNode] ← EquiCurvatureSegments[activeTrajectory, delta];
IF equiCurveSegNode #
NIL
AND numEquiCurveSegs > 4
THEN{
PERefresh.DisableTrajectoryRefresh[activeTrajectory.first];
[center, radius] ← OsculatingCircle[equiCurveSegNode, numEquiCurveSegs];
newSegmentList ← NewCurveSegs[equiCurveSegNode, numEquiCurveSegs, center, radius];
prevSegment ← equiCurveSegNode.preceding;
[tempSegList, oldSegList] ← PETrajectoryOps.UnlinkSegments[activeTrajectory.first.segments, equiCurveSegNode, numEquiCurveSegs];
activeTrajectory.first.segments ← PETrajectoryOps.LinkSegments[tempSegList, prevSegment, newSegmentList];
PERefresh.EnableTrajectoryRefresh[activeTrajectory.first];
RETURN;
};
filtered ← TRUE;
Take pairs of points, find the weighted average between the two, attach to head of new control polygon list.
firstSegment ← activeTrajectory.first.segments.first;
segment ← activeTrajectory.first.segments.preceding;
vertex1 ← segment.first.vertices.first;
vertex2 ← segment.preceding.first.vertices.first;
vertex3 ← segment.preceding.preceding.first.vertices.first;
vertex4 ← segment.preceding.preceding.preceding.first.vertices.first;
newSegmentList ← ComputeEndPoint[vertex1, vertex2, vertex3, vertex4, tempSegList];
tempSegList ← newSegmentList;
-- The currently implemented code recomputes n new positions for n control points by using a weighted neighborhood averaging with a neighborhood of two.
WHILE
segment.preceding.first # firstSegment
DO
vertex1 ← segment.first.vertices.first;
vertex2 ← segment.preceding.first.vertices.first;
vertex3 ← segment.preceding.preceding.first.vertices.first;
vertex4 ← segment.preceding.preceding.preceding.first.vertices.first;
[factor1, factor2, factor3] ← ComputeWeightingFactors[vertex1, vertex2, vertex3];
point.x ← factor1*vertex1.point.x + factor2*vertex2.point.x + factor3*vertex3.point.x;
point.y ← factor1*vertex1.point.y + factor2*vertex2.point.y + factor3*vertex3.point.y;
bias ← factor1*vertex1.bias + factor2*vertex2.bias + factor3*vertex3.bias;
tension ← factor1*vertex1.tension + factor2*vertex2.tension + factor3*vertex3.tension;
newSegmentList ← LinkSegments[tempSegList, point, bias, tension];
tempSegList ← newSegmentList;
segment ← segment.preceding;
ENDLOOP;
vertex1 ← segment.preceding.first.vertices.first;
vertex2 ← segment.first.vertices.first;
vertex3 ← segment.rest.first.vertices.first;
vertex4 ← segment.rest.rest.first.vertices.first;
newSegmentList ← ComputeEndPoint[vertex1, vertex2, vertex3, vertex4, newSegmentList];
newSegmentList.first.type ← moveTo;
newSegmentList.first.fp ← NIL;
IF newSegmentList # NIL THEN newControlPolygon ← NEW[PETypes.TrajectoryRec ← [segments: newSegmentList]];
Compute the Beta-spline that goes with the new control polygon and make that into a trajectory node also.
[newSplineSegmentList] ← PEBetaSpline.BetaSplinePoints[newSegmentList, newSegmentList, FALSE];
IF newSplineSegmentList #
NIL
THEN{
newSplineSegmentList.first.type ← moveTo;
newSplineSegmentList.first.fp ← NIL;
newSpline ← NEW[PETypes.TrajectoryRec ← [segments: newSplineSegmentList]];
};
};
ComputeWeightingFactors:
PRIVATE
PROCEDURE [vertex1: PETypes.Vertex, vertex2: PETypes.Vertex, vertex3: PETypes.Vertex]
RETURNS [factor1:
REAL, factor2:
REAL, factor3:
REAL] = {
The returned parameter, factor, should be between 0.0 and 1.0
weight1, weight2, weight3: REAL;
denom: REAL;
-- Using the natural log because the effects of tension and bias are non-linear.
weight1 ← 0.1*(1.0 + RealFns.Ln[vertex1.tension + vertex1.bias]);
weight2 ← 0.8*(1.0 + RealFns.Ln[vertex2.tension + vertex2.bias]);
weight3 ← 0.1*(1.0 + RealFns.Ln[vertex3.tension + vertex3.bias]);
denom ← weight1 + weight2 + weight3;
IF denom # 0.0
THEN {
factor1 ← weight1/denom;
factor2 ← weight2/denom;
factor3 ← weight3/denom;
}
ELSE {factor1 ← 0.1; factor2 ← 0.8; factor3 ← 0.1};
};
LinkSegments:
PRIVATE
PROCEDURE [changedSegList: PETypes.SegmentList, point: PETypes.Point, bias:
REAL, tension:
REAL]
RETURNS [newSegmentList: PETypes.SegmentList] = {
This routine inserts a vertex at the head of a vertex list, inserts the vertex list into a segment node, and inserts the segment at the head of a segment list which it returns.
changedVertexList: PETypes.VertexList;
newSegment: PETypes.Segment;
newVertex: PETypes.Vertex;
newVertexList: PETypes.VertexList;
newVertex ← NEW[PETypes.VertexRec ← [point: point, bias: bias, tension: tension]];
newVertexList ← NEW[PETypes.VertexNodeRec];
newVertexList ← NIL;
[changedVertexList, ] ← PETrajectoryOps.InsertVertex[newVertexList, newVertex, NIL];
newSegment ← NEW[PETypes.SegmentRec ← [vertices: changedVertexList]];
[newSegmentList, ] ← PETrajectoryOps.InsertSegment[changedSegList, newSegment, NIL];
};
ComputeCollinear:
PRIVATE
PROCEDURE [activeTrajectory: PETypes.TrajectoryNode, delta:
REAL]
RETURNS [collinearSegNode: PETypes.SegmentNode] = {
This routine returns a pointer to the control point that is the most closely collinear (up to some error metric, delta) with its neighbors on either side.
firstSegment: PETypes.Segment;
lastAcceptedDist: REAL ← delta;
perpendicularDist: REAL;
segment: PETypes.SegmentNode;
segment0, segment1, segment2: PETypes.Vertex;
tempSegList: PETypes.SegmentList ← NIL;
u: PETypes.Point;
vector1, vector2: PETypes.Point;
vector1Length: REAL;
xComponent, yComponent: REAL;
firstSegment ← activeTrajectory.first.segments.rest.first;
segment ← activeTrajectory.first.segments.preceding.preceding.preceding;
WHILE
segment.first # firstSegment
DO
segment0 ← segment.rest.first.vertices.first;
segment1 ← segment.first.vertices.first;
segment2 ← segment.preceding.first.vertices.first;
vector1.x ← segment2.point.x - segment0.point.x;
vector1.y ← segment2.point.y - segment0.point.y;
vector2.x ← segment1.point.x - segment0.point.x;
vector2.y ← segment1.point.y - segment0.point.y;
vector1Length ← RealFns.SqRt[vector1.x*vector1.x + vector1.y*vector1.y];
u.x ← vector1.x/vector1Length;
u.y ← vector1.y/vector1Length;
xComponent ← vector2.x*u.x + vector2.y*u.y;
yComponent ← vector2.y*u.x - vector2.x*u.y;
IF (xComponent/vector1Length) >= 0.0
AND (xComponent/vector1Length) <= 1.0
THEN{
perpendicularDist ← ABS[yComponent];
IF perpendicularDist < lastAcceptedDist
THEN{
collinearSegNode ← segment;
lastAcceptedDist ← perpendicularDist;
};
};
segment ← segment.preceding;
ENDLOOP;
};
EquiCurvatureSegments:
PRIVATE
PROCEDURE [activeTrajectory: PETypes.TrajectoryNode, delta:
REAL]
RETURNS [saveCount:
CARDINAL, savePointer: PETypes.SegmentNode] = {
This routine returns a list of pointers to control points that have approximately the same curvature (to some error metric, delta) with its neighbors on either side.
count: INT;
currentCos, nextCos: REAL;
currentSine, nextSine: REAL;
currentTheta, nextTheta: REAL;
epsilon: REAL ← delta;
firstSegment: PETypes.Segment;
length1, length2, length3: REAL;
point1, point2: PETypes.Point;
pointer: PETypes.SegmentNode ← NIL;
sameBias, sameTension: BOOLEAN;
segment: PETypes.SegmentNode;
vertex1, vertex2, vertex3, vertex4: PETypes.Vertex;
firstSegment ← activeTrajectory.first.segments.first;
segment ← activeTrajectory.first.segments.preceding.preceding.preceding;
count ← 1;
saveCount ← 1;
point1.x ← segment.first.vertices.first.point.x - segment.rest.first.vertices.first.point.x;
point1.y ← segment.first.vertices.first.point.y - segment.rest.first.vertices.first.point.y;
point2.x ← segment.preceding.first.vertices.first.point.x - segment.first.vertices.first.point.x;
point2.y ← segment.preceding.first.vertices.first.point.y - segment.first.vertices.first.point.y;
[currentCos, currentSine, length1, length2] ← ComputeCurvatureElements[point1, point2];
WHILE
segment.first # firstSegment
DO
vertex1 ← segment.rest.first.vertices.first;
vertex2 ← segment.first.vertices.first;
vertex3 ← segment.preceding.first.vertices.first;
vertex4 ← segment.preceding.preceding.first.vertices.first;
IF vertex1.bias = vertex2.bias AND vertex2.bias = vertex3.bias AND vertex3.bias = vertex4.bias THEN sameBias ← TRUE;
IF vertex1.tension = vertex2.tension AND vertex2.tension = vertex3.tension AND vertex3.tension = vertex4.tension THEN sameTension ← TRUE;
Calculating the angles between segments and the length of each segment
point1.x ← vertex3.point.x - vertex2.point.x;
point1.y ← vertex3.point.y - vertex2.point.y;
point2.x ← vertex4.point.x - vertex3.point.x;
point2.y ← vertex4.point.y - vertex3.point.y;
[nextCos, nextSine, , length3] ← ComputeCurvatureElements[point1, point2];
currentTheta ← RealFns.ArcTanDeg[currentSine, currentCos];
IF currentSine < 0 THEN currentTheta ← pi + currentTheta;
IF currentTheta = 0.0 THEN IF currentCos < 0.0 THEN currentTheta ← pi;
nextTheta ← RealFns.ArcTanDeg[nextSine, nextCos];
IF nextSine < 0 THEN nextTheta ← pi + nextTheta;
IF nextTheta = 0.0 THEN IF nextCos < 0.0 THEN nextTheta ← pi;
The check for equicurvature segments involves testing for equal length segments, equal angles, bias, and tension between segments. Yet to be implemented, checking for segments that have more than a 360 degree sweep.
IF
ABS[nextTheta - currentTheta] < epsilon
THEN{
IF (
ABS[length1-length2] < epsilon)
AND (
ABS[length2-length3] < epsilon)
AND (
ABS[length1-length3] < epsilon)
THEN{
IF sameBias
AND sameTension
THEN
{
count ← count + 1;
pointer ← segment.preceding;
};
};
}
ELSE IF count > saveCount
THEN {
-- At the end of a run of equicurvature segments
saveCount ← count;
savePointer ← pointer;
count ← 1;
pointer ← segment.preceding;
};
segment ← segment.preceding;
currentCos ← nextCos;
currentSine ← nextSine;
length1 ← length2;
length2 ← length3;
sameBias ← FALSE;
sameTension ← FALSE;
ENDLOOP;
IF count > saveCount
THEN {
-- Pick up last equicurvature segments, if any
saveCount ← count;
savePointer ← pointer;
};
IF saveCount > 1
THEN{
saveCount ← saveCount + 2;
savePointer ← savePointer.preceding;
};
};
ComputeCurvatureElements:
PRIVATE
PROCEDURE [a: PETypes.Point, b: PETypes.Point]
RETURNS [cosTheta:
REAL, sineTheta:
REAL, aLength:
REAL, bLength:
REAL] = {
This routine approximates the curvature of a trajectory by using dq/ds. Cosq is calculated by the formula: cosq = A.B/(||A|| ||B||). Sin q is calculated by the formula: sinq = (AxB)/(||A|| ||B||), where A, B are vectors with tails at the origin. Assumes that vectors a and b have tails at the origin.
aCrossB: REAL;
aDotB: REAL;
aLength ← RealFns.SqRt[a.x*a.x + a.y*a.y];
bLength ← RealFns.SqRt[b.x*b.x + b.y*b.y];
aDotB ← a.x*b.x + a.y*b.y;
cosTheta ← aDotB/(aLength*bLength);
aCrossB ← a.x*b.y - a.y*b.x;
sineTheta ← aCrossB/(aLength*bLength);
};
ComputeEndPoint:
PRIVATE
PROCEDURE [vertex1: PETypes.Vertex, vertex2: PETypes.Vertex, vertex3: PETypes.Vertex, vertex4: PETypes.Vertex, oldSegmentList: PETypes.SegmentList]
RETURNS [newSegmentList: PETypes.SegmentNode] = {
This routine computes the endpoints of the simplified control polygon. Presently, the new head and tail segments bisect the old head and tail segments, respectively, and have the same length as their respective new neighboring segments.
bias: REAL;
factor, factor1, factor2, factor3: REAL;
midPoint: PETypes.Point;
midSegLength: REAL;
newPoint: PETypes.Point;
nextPoint, afterNextPoint: PETypes.Point;
oldSegLength: REAL;
point: PETypes.Point;
tension: REAL;
midPoint.x ← vertex2.point.x + .5*(vertex1.point.x - vertex2.point.x);
midPoint.y ← vertex2.point.y + .5*(vertex1.point.y - vertex2.point.y);
bias ← vertex2.bias + .5*(vertex1.bias - vertex2.bias);
tension ← vertex2.tension + .5*(vertex1.tension - vertex2.tension);
[factor1, factor2, factor3] ← ComputeWeightingFactors[vertex1, vertex2, vertex3];
nextPoint.x ← factor1*vertex1.point.x + factor2*vertex2.point.x + factor3*vertex3.point.x;
nextPoint.y ← factor1*vertex1.point.y + factor2*vertex2.point.y + factor3*vertex3.point.y;
[factor1, factor2, factor3] ← ComputeWeightingFactors[vertex2, vertex3, vertex4];
afterNextPoint.x ← factor1*vertex2.point.x + factor2*vertex3.point.x + factor3*vertex4.point.x;
afterNextPoint.y ← factor1*vertex2.point.y + factor2*vertex3.point.y + factor3*vertex4.point.y;
midSegLength ← RealFns.SqRt[(nextPoint.x-midPoint.x)*(nextPoint.x-midPoint.x) + (nextPoint.y-midPoint.y)*(nextPoint.y-midPoint.y)];
oldSegLength ← RealFns.SqRt[(nextPoint.x-afterNextPoint.x)*(nextPoint.x-afterNextPoint.x) + (nextPoint.y-afterNextPoint.y)*(nextPoint.y-afterNextPoint.y)];
factor ← MAX[oldSegLength/midSegLength, 1.0];
newPoint.x ← nextPoint.x + factor*(midPoint.x - nextPoint.x);
newPoint.y ← nextPoint.y + factor*(midPoint.y - nextPoint.y);
newSegmentList ← LinkSegments[oldSegmentList, newPoint, bias, tension];
};
NewCurveSegs:
PRIVATE
PROCEDURE [firstVertex: PETypes.SegmentNode, numVertices:
INT, center: PETypes.Point, radius:
REAL]
RETURNS [newSegmentList: PETypes.SegmentNode] = {
This routine removes a control point from a portion of the control polygon that is judged to have n adjacent points of equal curvature. The routine evenly spreads n-1 points around the osculating circle, with the first and last points remaining unchanged from the original first and last points.
angle: REAL;
cosStep: REAL;
i: INT ← 0;
lastVertex: PETypes.SegmentNode ← firstVertex;
newPoint: PETypes.Point;
numPts: REAL;
point: PETypes.Point;
sineStep: REAL;
stepAngle: REAL;
tempSegList: PETypes.SegmentList ← NIL;
theta: REAL;
First find the angle between the vector from the origin to the first vertex of the equicurvature area and the vector from the origin to the last vertex of the equicurvature area.
WHILE i < numVertices-1
DO
lastVertex ← lastVertex.rest;
i ← i+1;
ENDLOOP;
point.x ← lastVertex.first.vertices.first.point.x - center.x;
point.y ← lastVertex.first.vertices.first.point.y - center.y;
theta ← PolygonAngle[center, firstVertex, numVertices];
numPts ← ABS[(2*pi)/theta];
stepAngle ← 2*pi/(numPts-2);
IF theta > 0 THEN stepAngle ← -stepAngle;
The angle between new points is the angle found above divided by one less the original number of equicurvature points.
angle ← 0.0;
i ← 1;
WHILE i <= numVertices-2
DO
cosStep ← RealFns.CosDeg[angle];
sineStep ← RealFns.SinDeg[angle];
newPoint.x ← point.x*cosStep - point.y*sineStep + center.x;
newPoint.y ← point.y*cosStep + point.x*sineStep + center.y;
newSegmentList ← LinkSegments[tempSegList, newPoint, 1.0, 0.0];
tempSegList ← newSegmentList;
angle ← angle + stepAngle;
i ← i+1;
ENDLOOP;
newSegmentList ← LinkSegments[tempSegList, firstVertex.first.vertices.first.point, 1.0, 0.0];
};
OsculatingCircle:
PRIVATE
PROCEDURE [firstVertex: PETypes.SegmentNode, numVertices:
INT]
RETURNS [center: PETypes.Point, radius:
REAL] = {
This routine approximates the center and radius of the osculating circle at two adjacent points in the control polygon.
denom: REAL;
i: INT ← 1;
numer: REAL;
point1, point2, point3, point4: PETypes.Point;
segment: PETypes.SegmentNode;
sumCenter: PETypes.Point;
sumRadius: REAL ← 0.0;
t: REAL;
tempPoint: PETypes.Point;
IF numVertices < 4 THEN RETURN;
sumCenter.x ← 0.0;
sumCenter.y ← 0.0;
segment ← firstVertex;
point1 ← segment.first.vertices.first.point;
point2 ← segment.rest.first.vertices.first.point;
point3 ← segment.rest.rest.first.vertices.first.point;
point4 ← segment.rest.rest.rest.first.vertices.first.point;
WHILE i <= numVertices-3
DO
tempPoint.x ← point1.x + point3.x - point2.x;
tempPoint.y ← point1.y + point3.y - point2.y;
denom ← -point1.x + 3*point2.x - 3*point3.x + point4.x;
numer ← point2.x - point3.x;
IF denom = 0.0
THEN{
denom ← -point1.y + 3*point2.y - 3*point3.y + point4.y;
numer ← point2.y - point3.y;
};
t ← numer/denom;
sumCenter.x ← sumCenter.x + (point2.x + t*(tempPoint.x - point2.x));
sumCenter.y ← sumCenter.y + (point2.y + t*(tempPoint.y - point2.y));
center.x ← point2.x + t*(tempPoint.x - point2.x);
center.y ← point2.y + t*(tempPoint.y - point2.y);
sumRadius ← sumRadius + RealFns.SqRt[(center.x-point2.x)*(center.x-point2.x) + (center.y-point2.y)*(center.y-point2.y)];
segment ← segment.rest;
point1 ← segment.first.vertices.first.point;
point2 ← segment.rest.first.vertices.first.point;
point3 ← segment.rest.rest.first.vertices.first.point;
point4 ← segment.rest.rest.rest.first.vertices.first.point;
i ← i+1;
ENDLOOP;
center.x ← sumCenter.x/(numVertices-3);
center.y ← sumCenter.y/(numVertices-3);
radius ← sumRadius/numVertices;
};
PolygonAngle:
PRIVATE
PROCEDURE [center: PETypes.Point, firstVertex: PETypes.SegmentNode, numVertices:
INT]
RETURNS [theta:
REAL] = {
This routine finds the angle between two vectors whose heads are two adjacent vertices of a regular polygon and whose tails are both the center of that polygon. Used in the simplification algorithm when re-distributing control points around a circular arc.
cosTheta: REAL;
i: INT ← 1;
point1, point2: PETypes.Point;
segment: PETypes.SegmentNode;
sineTheta: REAL;
sumAngles: REAL ← 0.0;
segment ← firstVertex;
point1.x ← segment.first.vertices.first.point.x - center.x;
point1.y ← segment.first.vertices.first.point.y - center.y;
point2.x ← segment.rest.first.vertices.first.point.x - center.x;
point2.y ← segment.rest.first.vertices.first.point.y - center.y;
WHILE i < numVertices
DO
[cosTheta, sineTheta, ] ← ComputeCurvatureElements[point1, point2];
sumAngles ← sumAngles + RealFns.ArcTanDeg[sineTheta, cosTheta];
segment ← segment.rest;
point1 ← point2;
point2.x ← segment.rest.first.vertices.first.point.x - center.x;
point2.y ← segment.rest.first.vertices.first.point.y - center.y;
i ← i+1;
ENDLOOP;
theta ← sumAngles/(numVertices-1);
};
END.