PEAlgebraImpl.mesa
Written by Darlene Plebon on August 16, 1983 3:32 pm
Useful linear algebra routines.
DIRECTORY
PEAlgebra,
PETypes USING [Point];
PEAlgebraImpl:
CEDAR
PROGRAM
IMPORTS PEAlgebra
EXPORTS PEAlgebra =
BEGIN OPEN PEAlgebra, PETypes;
ClosestPointOnLineSegment:
PUBLIC PROCEDURE [point, p0, p1: Point]
RETURNS [closestPoint: Point] = {
This routine returns the point on the line segment [p0,p1] which is closest to the specified point.
distanceP0P1, distanceP0ClosestPoint, distanceP1ClosestPoint: REAL;
closestPoint ← Projection[point, p0, p1];
distanceP0P1 ← DistanceSquared[p0, p1];
distanceP0ClosestPoint ← DistanceSquared[p0, closestPoint];
distanceP1ClosestPoint ← DistanceSquared[p1, closestPoint];
IF distanceP0ClosestPoint > distanceP0P1
OR distanceP1ClosestPoint > distanceP0P1
THEN
IF distanceP0ClosestPoint < distanceP1ClosestPoint THEN closestPoint ← p0
ELSE closestPoint ← p1;
};
Projection:
PUBLIC PROCEDURE [point, p0, p1: Point]
RETURNS [colinearPoint: Point] = {
This procedure projects the specified point onto the line defined by p0 and p1. The returned point is colinear with p0 and p1.
diffX: REAL ← p1.x - p0.x;
diffY: REAL ← p1.y - p0.y;
diffX2: REAL ← Sqr[diffX];
diffY2: REAL ← Sqr[diffY];
diffX2PlusDiffY2: REAL ← diffX2 + diffY2;
diffXY: REAL ← p0.y*p1.x - p0.x*p1.y;
diffXDiffY: REAL ← diffX * diffY;
IF diffX2PlusDiffY2 = 0.0 THEN colinearPoint ← p0
ELSE {
colinearPoint.x ← (diffXDiffY*point.y + diffX2*point.x - diffXY*diffY) / diffX2PlusDiffY2;
colinearPoint.y ← (diffXDiffY*point.x + diffY2*point.y + diffXY*diffX) / diffX2PlusDiffY2;
};
};
Intersection:
PUBLIC
PROCEDURE [line1A, line1B, line2A, line2B: Point]
RETURNS [intersecting:
BOOLEAN, p: Point] = {
This routine determines the intersection of two lines, each of which is defined by two points on the line. This routine returns an indication of whether the two lines intersect, along with the point of intersection.
line1Diff, line2Diff: Point;
t1, t2, denominator: REAL;
line1Diff.x ← line1B.x - line1A.x;
line1Diff.y ← line1B.y - line1A.y;
line2Diff.x ← line2B.x - line2A.x;
line2Diff.y ← line2B.y - line2A.y;
t1 ← line1A.y*line1B.x - line1B.y*line1A.x;
t2 ← line2A.y*line2B.x - line2B.y*line2A.x;
denominator ← line1Diff.y*line2Diff.x - line2Diff.y*line1Diff.x;
intersecting ← denominator # 0.0;
IF intersecting
THEN {
p.x ← (t2*line1Diff.x - t1*line2Diff.x) / denominator;
p.y ← (t2*line1Diff.y - t1*line2Diff.y) / denominator;
};
};
END.