DIRECTORY
Angles2d, Imager, ImagerTransformation, Lines2d, Lines2dTypes, RealFns, Vectors2d;

Lines2dImpl: CEDAR PROGRAM
IMPORTS Angles2d, Imager, ImagerTransformation, RealFns, Vectors2d
EXPORTS Lines2d =

BEGIN
Point: TYPE = Lines2dTypes.Point;
Edge: TYPE = REF EdgeObj;
EdgeObj: TYPE = Lines2dTypes.EdgeObj;
Line: TYPE = REF LineObj;
LineObj: TYPE = Lines2dTypes.LineObj;
Ray: TYPE = REF RayObj;
RayObj: TYPE = Lines2dTypes.RayObj;
Vector: TYPE = Lines2dTypes.Vector;

CreateEmptyLine: PUBLIC PROC RETURNS [line: Line] = {
line _ NEW[LineObj];
};

CopyLine: PUBLIC PROC [from: Line, to: Line] = {
to.c _ from.c;
to.s _ from.s;
to.theta _ from.theta;
to.d _ from.d;
};


EqualLine: PUBLIC PROC [a: Line, b: Line] RETURNS [BOOL] = {
RETURN[a.d = b.d AND a.theta = b.theta];
};

AlmostEqualLine: PUBLIC PROC [a: Line, b: Line, errorDegrees: REAL, errorDistance: REAL] RETURNS [BOOL] = {
RETURN[ABS[a.theta - b.theta] < errorDegrees AND ABS[a.d-b.d] < errorDistance];
};

FillLineFromPoints: PUBLIC PROC [v1, v2: Point, line: Line] = {
epsilon: REAL = 1.0e-3; -- changed from 1.0e-5 in August 1986, by Bier
-- changed from 1.0e-6 in August 1985 by Bier
x2Minusx1: REAL _ v2.x - v1.x;
y2Minusy1: REAL _ v2.y - v1.y;
IF ABS[x2Minusx1] < epsilon THEN {-- vertical line
IF v2.y > v1.y THEN {-- line goes up
line.theta _ 90.0;
line.s _ 1;
line.d _ -v1.x}
ELSE { -- line goes down
line.theta _ -90;
line.s _ -1;
line.d _ v1.x};
line.c _ 0;
}
ELSE {
line.theta _ RealFns.ArcTanDeg[y2Minusy1, x2Minusx1];
line.c _ RealFns.CosDeg[line.theta];
line.s _ RealFns.SinDeg[line.theta];
line.d _ v1.y*line.c - v1.x*line.s;
};
}; -- end of FillLineFromPoints

VectorTooSmall: PUBLIC SIGNAL = CODE;
FillLineFromPointAndVector: PUBLIC PROC [pt: Point, vec: Vector, line: Line] = {
epsilon: REAL = 1.0e-3;
IF ABS[vec.x] < epsilon AND ABS[vec.y] < epsilon THEN {
SIGNAL VectorTooSmall;
line.theta _ 90.0;
}
ELSE line.theta _ RealFns.ArcTanDeg[vec.y, vec.x];
line.c _ RealFns.CosDeg[line.theta];
line.s _ RealFns.SinDeg[line.theta];
line.d _ pt.y*line.c - pt.x*line.s;
};

FillLineFromCoefficients: PUBLIC PROC [sineOfTheta, cosineOfTheta, distance: REAL, line: Line] = {
line.s _ sineOfTheta;
line.c _ cosineOfTheta;
line.d _ distance;
line.theta _ RealFns.ArcTanDeg[sineOfTheta, cosineOfTheta];
}; -- end of FillLineFromCoefficients

FillLineFromPointAndAngle: PUBLIC PROC [pt: Point, degrees: REAL, line: Line] = {
line.theta _ Angles2d.Normalize[degrees];
line.c _ RealFns.CosDeg[line.theta];
line.s _ RealFns.SinDeg[line.theta];
line.d _ pt.y*line.c - pt.x*line.s;
}; -- end of FillLineFromPointAndAngle

FillLineNormalToLineThruPoint: PUBLIC PROC [line: Line, pt: Point, normalLine: Line] = {
normalLine.s _ line.c;
normalLine.c _ -line.s;
normalLine.d _ -pt.y*line.s - pt.x*line.c;
normalLine.theta _ Angles2d.Add[line.theta, 90];
}; -- end of FillLineAsNormal

FillLineLeftOfLine: PUBLIC PROC [line: Line, dist: REAL, parallelLine: Line] = {
parallelLine.s _ line.s;
parallelLine.c _ line.c;
parallelLine.d _ line.d + dist;
parallelLine.theta _ line.theta;
};

FillLineRightOfLine: PUBLIC PROC [line: Line, dist: REAL, parallelLine: Line] = {
parallelLine.s _ line.s;
parallelLine.c _ line.c;
parallelLine.d _ line.d - dist;
parallelLine.theta _ line.theta;
};

FillLineTransform: PUBLIC PROC [fixed: Line, transform: ImagerTransformation.Transformation, line: Line] = {
point, newPoint: Point;
direction, newDirection: Vector;
point _ PointOnLine[fixed];
direction _ DirectionOfLine[fixed];
newPoint _ ImagerTransformation.Transform[transform, point];
newDirection _ ImagerTransformation.TransformVec[transform, direction];
FillLineFromPointAndVector[newPoint, newDirection, line];
};

LineFromPoints: PUBLIC PROC [v1, v2: Point] RETURNS [line: Line] = {
line _ CreateEmptyLine[];
FillLineFromPoints[v1, v2, line];
};
LineFromPointAndVector: PUBLIC PROC [pt: Point, vec: Vector] RETURNS [line: Line] = {
pt2: Point;
pt2 _ Vectors2d.Add[pt, vec];
line _ LineFromPoints[pt, pt2];
};
LineFromCoefficients: PUBLIC PROC [sineOfTheta, cosineOfTheta, distance: REAL] RETURNS [line: Line] = {
line _ CreateEmptyLine[];
FillLineFromCoefficients[sineOfTheta, cosineOfTheta, distance, line];
};
LineFromPointAndAngle: PUBLIC PROC  [pt: Point, degrees: REAL] RETURNS [line: Line] = {
line _ CreateEmptyLine[];
FillLineFromPointAndAngle[pt, degrees, line];
};

LineNormalToLineThruPoint: PUBLIC PROC [line: Line, pt: Point] RETURNS [normalLine: Line] = {
normalLine _ CreateEmptyLine[];
FillLineNormalToLineThruPoint[line, pt, normalLine];
};

LineLeftOfLine: PUBLIC PROC [line: Line, dist: REAL] RETURNS [parallelLine: Line] = {
parallelLine _ CreateEmptyLine[];
FillLineLeftOfLine[line, dist, parallelLine];
};

LineRightOfLine: PUBLIC PROC [line: Line, dist: REAL] RETURNS [parallelLine: Line] = {
parallelLine _ CreateEmptyLine[];
FillLineRightOfLine[line, dist, parallelLine];
};

LineTransform: PUBLIC PROC [fixed: Line, transform: ImagerTransformation.Transformation] RETURNS [rotatedLine: Line] = {
rotatedLine _ CreateEmptyLine[];
FillLineTransform[fixed, transform, rotatedLine];
};

DrawLine: PUBLIC PROC [dc: Imager.Context, line: Line, clippedBy: Imager.Rectangle, strokeWidth: REAL _ 1.0] = {
count: NAT;
ray: Ray;
params: ARRAY[1..2] OF REAL;
p1, p2, basePoint: Point;
direction: Vector;
DoDrawLine: PROC = {
Imager.SetXY[dc, [p1.x, p1.y]];
Imager.Move[dc];
Imager.SetStrokeEnd[dc, round];
Imager.SetStrokeWidth[dc, strokeWidth];
Imager.MaskVector[dc, [0.0, 0.0], [p2.x - p1.x, p2.y - p1.y]];
};
p1 _ [clippedBy.x, clippedBy.y];
p2 _ [clippedBy.x + clippedBy.w, clippedBy.y + clippedBy.h];
basePoint _ PointOnLine[line];
direction _ DirectionOfLine[line];
ray _ CreateRay[basePoint, direction];
[count, params] _ LineRayMeetsBox[ray, p1.x, p1.y, p2.x, p2.y];
IF count = 2 THEN {
p1 _ EvalRay[ray, params[1]];
p2 _ EvalRay[ray, params[2]];
Imager.DoSave[dc, DoDrawLine];
};
};

CreateEmptyEdge: PUBLIC PROC RETURNS [edge: Edge] = {
edge _ NEW[EdgeObj];
edge.line _ CreateEmptyLine[];
};
CopyEdge: PUBLIC PROC [from: Edge, to: Edge] = {
CopyLine[from.line, to.line];
to.start _ from.start;
to.end _ from.end;
to.startIsFirst _ from.startIsFirst;
}; -- end of CopyEdge
FillEdge: PUBLIC PROC [v1, v2: Point, edge: Edge] = {
y2Minusy1: REAL;
FillLineFromPoints[v1, v2, edge.line];
y2Minusy1 _ v2.y - v1.y;
IF y2Minusy1 = 0 THEN 
IF v2.x > v1.x THEN {edge.end _ v2; edge.start _ v1; edge.startIsFirst _ TRUE}
ELSE {edge.end _ v1; edge.start _ v2; edge.startIsFirst _ FALSE}
ELSE
IF v2.y > v1.y THEN {edge.end _ v2; edge.start _ v1; edge.startIsFirst _ TRUE}
ELSE {edge.end _ v1; edge.start _ v2; edge.startIsFirst _ FALSE};
}; -- end of FillEdge
FillEdgeTransform: PUBLIC PROC [fixed: Edge, transform: ImagerTransformation.Transformation, edge: Edge] = {
start, end: Point;
start _ ImagerTransformation.Transform[m: transform, v: fixed.start];
end _ ImagerTransformation.Transform[m: transform, v: fixed.end];
FillEdge[start, end, edge];
};
CreateEdge: PUBLIC PROC [v1, v2: Point] RETURNS [edge: Edge] = {
edge _ CreateEmptyEdge[];
FillEdge[v1, v2, edge];
}; -- end of CreateEdge
EdgeTransform: PUBLIC PROC [fixed: Edge, transform: ImagerTransformation.Transformation] RETURNS [edge: Edge] = {
edge _ CreateEmptyEdge[];
FillEdgeTransform[fixed, transform, edge];
};
CreateRay: PUBLIC PROC [basePoint: Point, direction: Vector] RETURNS [ray: Ray] = {
ray _ NEW[RayObj _ [basePoint, direction]];
};

CreateRayFromPoints: PUBLIC PROC [p1, p2: Point] RETURNS [ray: Ray] = {
ray _ NEW[RayObj _ [p1, Vectors2d.Sub[p2, p1]]];
};

AlmostEqual: PROC [r1, r2, almostZero: REAL] RETURNS [BOOL] = {
RETURN[ABS[r1 - r2] < almostZero];
};

LineRayMeetsBox: PUBLIC PROC [ray: Ray, xmin, ymin, xmax, ymax: REAL] RETURNS [count: NAT, params: ARRAY[1..2] OF REAL] = {
almostZero: REAL _ 1.0e-3;
x, y, t: REAL;
count _ 0;
IF ABS[ray.d.y] > almostZero THEN { -- intersection occurs
t _ (ymax-ray.p.y)/ray.d.y;
x _ ray.p.x + t*ray.d.x;
IF x >=xmin-almostZero AND x<= xmax+almostZero THEN { -- hits box
count _ count + 1;
params[count] _ t;
};
t _ (ymin-ray.p.y)/ray.d.y;
x _ ray.p.x + t*ray.d.x;
IF x >=xmin-almostZero AND x<= xmax+almostZero THEN { -- hits box
count _ count + 1;
params[count] _ t;
};
};
IF ABS[ray.d.x] > almostZero THEN { -- intersection occurs
IF count < 2 THEN {
t _ (xmax-ray.p.x)/ray.d.x;
IF count = 0 OR (count = 1 AND NOT AlmostEqual[t, params[1], almostZero]) THEN {
y _ ray.p.y + t*ray.d.y;
IF y >=ymin-almostZero AND y<= ymax+almostZero THEN { -- hits box
count _ count + 1;
params[count] _ t;
};
};
};
IF count < 2 THEN {
t _ (xmin-ray.p.x)/ray.d.x;
IF count = 0 OR (count = 1 AND NOT AlmostEqual[t, params[1], almostZero]) THEN {
y _ ray.p.y + t*ray.d.y;
IF y >=ymin-almostZero AND y<= ymax+almostZero THEN { -- hits box
count _ count + 1;
params[count] _ t;
};
};
};
};
IF count = 2 THEN {
IF params[2] <  params[1] THEN {
temp: REAL _ params[1];
params[1] _ params[2];
params[2] _ temp;
};
}; -- make sure hits are sorted
}; -- end of LineRayMeetsBox

EvalRay: PUBLIC PROC [ray: Ray, param: REAL] RETURNS [point: Point] = {
point.x _ ray.p.x + param*ray.d.x;
point.y _ ray.p.y + param*ray.d.y;
};
AlmostZero: PROC [r: REAL] RETURNS [BOOL] = {
epsilon: REAL = 1.0e-5;
RETURN[ABS[r] < epsilon];
};


LineMeetsLine: PUBLIC PROC [line1, line2: Line] RETURNS [intersection: Point, parallel: BOOL] = {
determinant: REAL;
epsilon: REAL = 4E-4;
parallel _ FALSE;
IF Angles2d.AlmostParallel[line1.theta, line2.theta, epsilon] THEN {parallel _ TRUE; RETURN};
determinant _ line2.s*line1.c - line1.s*line2.c;
intersection.x _ (line2.c*line1.d - line1.c*line2.d)/determinant;
intersection.y _ (line2.s*line1.d - line1.s*line2.d)/determinant;
};

LineMeetsYAxis: PUBLIC PROC [line: Line] RETURNS [yInt: REAL, parallel: BOOL] = {
IF line.theta = 90 OR line.theta = -90 THEN parallel _ TRUE
ELSE {-- we just want the y Intercept which is calculated at line creation time for now.
parallel _ FALSE;
yInt _ line.yInt;}
};

LineMeetsEdge: PUBLIC PROC [line: Line, edge: Edge] RETURNS [intersection: Point, noHit: BOOL] = {
edgeLine: Line _ edge.line;
parallel: BOOL;
[intersection, parallel] _ LineMeetsLine[edgeLine, line];
IF parallel THEN {noHit _ TRUE; RETURN};
noHit _ NOT OnEdge[intersection, edge];
};

EdgeMeetsEdge: PUBLIC PROC  [e1, e2: Edge] RETURNS [intersection: Point, noHit: BOOL] = {
e1Line: Line _ e1.line;
[intersection, noHit] _ LineMeetsEdge[e1Line, e2];
IF noHit THEN RETURN;
noHit _ NOT OnEdge[intersection, e1];
};
SignedLineDistance: PUBLIC PROC [pt: Point, line: Line] RETURNS [d: REAL] = {
d _ pt.y*line.c - pt.x*line.s - line.d;
}; -- SignedLineDistance

LineDistance: PUBLIC PROC [pt: Point, line: Line] RETURNS [d: REAL] = {
d _ ABS[pt.y*line.c - pt.x*line.s - line.d];
}; -- LineDistance


DropPerpendicular: PUBLIC PROC [pt: Point, line: Line] RETURNS [projectedPt: Point] = {
D: REAL _ pt.y*line.c - pt.x*line.s - line.d;
projectedPt.x _ pt.x + D*line.s;
projectedPt.y _ pt.y - D*line.c;
};

PointOnLine: PUBLIC PROC [line: Line] RETURNS [pt: Point] = {
IF ABS[line.c] > ABS[line.s] THEN {
pt.x _ 0.0;
pt.y _ line.d/line.c;
}
ELSE {
pt.y _ 0.0;
pt.x _ -line.d/line.s;
};
};
DirectionOfLine: PUBLIC PROC [line: Line] RETURNS [direction: Vector] = {
direction.x _ line.c;
direction.y _ line.s;
};
NearestEndpoint: PUBLIC PROC [pt: Point, edge: Edge] RETURNS [endpoint: Point] = {
IF ABS[pt.x-edge.start.x] <= ABS[pt.x-edge.end.x] THEN
IF ABS[pt.y-edge.start.y] <= ABS[pt.y-edge.end.y] THEN RETURN[edge.start]
ELSE GOTO DoMath
ELSE
IF ABS[pt.y-edge.start.y] > ABS[pt.y-edge.end.y] THEN RETURN[edge.end]
ELSE GOTO DoMath;
EXITS
DoMath =>
IF DistanceSquaredPointToPoint[pt, edge.start] < DistanceSquaredPointToPoint[pt, edge.end]
THEN endpoint _ edge.start
ELSE endpoint _ edge.end;
};

DistanceSquaredToNearestEndpoint: PUBLIC PROC [pt: Point, edge: Edge] RETURNS [distanceSquared: REAL] = {
distance2ToPLo, distance2ToPHi: REAL;
distance2ToPLo _ DistanceSquaredPointToPoint[pt, edge.start];
distance2ToPHi _ DistanceSquaredPointToPoint[pt, edge.end];
RETURN[MIN[distance2ToPLo, distance2ToPHi]];
};

NearestPointOnEdge: PUBLIC PROC [pt: Point, edge: Edge] RETURNS [onEdge: Point] = {
projectedPt: Point _ DropPerpendicular[pt, edge.line];
IF LinePointOnEdge[projectedPt, edge] THEN onEdge _ projectedPt
ELSE onEdge _ NearestEndpoint[pt, edge];
};

DistancePointToEdge: PUBLIC PROC [pt: Point, edge: Edge] RETURNS [distance: REAL] = {
projectedPt: Point _ DropPerpendicular[pt, edge.line];
nearEndpoint: Point;
IF LinePointOnEdge[projectedPt, edge] THEN distance _ ABS[LineDistance[pt, edge.line]]
ELSE {
nearEndpoint _ NearestEndpoint[pt, edge];
distance _ DistancePointToPoint[pt, nearEndpoint];
};
};

DistanceSquaredPointToEdge: PUBLIC PROC [pt: Point, edge: Edge] RETURNS [distanceSquared: REAL] = {
projectedPt: Point _ DropPerpendicular[pt, edge.line];
IF LinePointOnEdge[projectedPt, edge]
THEN {distanceSquared _ LineDistance[pt, edge.line];
 distanceSquared _ distanceSquared*distanceSquared}
ELSE distanceSquared _ DistanceSquaredToNearestEndpoint[pt, edge];
};

OnEdge: PUBLIC PROC [pt: Point, edge: Edge] RETURNS [BOOL] = {
d2: REAL;
d2 _ DistanceSquaredPointToEdge[pt, edge];
RETURN[AlmostZero[d2]];
};
LinePointOnEdge: PUBLIC PROC [pt: Point, edge: Edge] RETURNS [BOOL] = {
IF ABS[edge.end.x - edge.start.x] <= ABS[edge.end.y - edge.start.y] THEN  -- line is more vertical or has zero length
RETURN[Between[pt.y, edge.start.y, edge.end.y]]
ELSE -- line is more horizontal
RETURN[Between[pt.x, edge.start.x, edge.end.x]];
}; -- end of LinePointOnEdge
DistancePointToPoint: PUBLIC PROC [p1, p2: Point] RETURNS [distance: REAL] = {
distance _ RealFns.SqRt[(p2.y-p1.y)*(p2.y-p1.y) + (p2.x-p1.x)*(p2.x-p1.x)];
};
DistanceSquaredPointToPoint: PUBLIC PROC [p1, p2: Point] RETURNS [distance: REAL] = {
distance _ (p2.y-p1.y)*(p2.y-p1.y) + (p2.x-p1.x)*(p2.x-p1.x);
};
PointLeftOfLine: PUBLIC PROC [distance: REAL, pOnLine: Point, line: Line] RETURNS [point: Point] = {

lineParallel, linePerp: Line;
parallel: BOOL;
lineParallel _ CreateEmptyLine[];
linePerp _ CreateEmptyLine[];
FillLineFromCoefficients[line.s, line.c, line.d + distance, lineParallel];
FillLineNormalToLineThruPoint[line, pOnLine, linePerp];
[point, parallel] _ LineMeetsLine[lineParallel, linePerp];
IF parallel THEN ERROR; -- perpendicular lines are not parallel
};


Between: PRIVATE PROC [test, a, b: REAL] RETURNS [BOOL] = {
SELECT a FROM
< b => RETURN [a <= test AND test <= b];
= b => RETURN [test = b];
> b => RETURN [b <= test AND test <= a];
ENDCASE => ERROR;
};

END.

:File: Lines2dImpl.mesa
Author: Eric Bier on June 4, 1985 5:04:38 pm PDT
Last edited by Bier on June 24, 1987 11:34:21 am PDT
Contents:  Routines for finding the intersections of various types of lines and line segments in Gargoyle.
Pier, August 8, 1986 12:14:59 pm PDT
Bier, March 21, 1988 12:02:31 pm PST
Making Lines

to.slope _ from.slope;
to.yInt _ from.yInt;
Are these lines mathematically identical?
Returns TRUE if the lines differ in slope by no more than errorDegrees and differ in distance from the origin by no more than errorDistance.
Recall y*c - x*s - d = 0;
Calculates the different parts of a line given an ordered pair of points (the tail and the head).  Trig lines are directed in sense since 0 <= line.theta <= 180 implies that v1 is lower than or to the right of) v2.
Deal with very short and vertical segments
Notice that zero length lines are considered vertical.
we have -x*s = d. where s = 1. Plug in v1.  -v1.x*s = d
we have -x*s = d. where s = -1. Plug in v1.  -v1.x*s = d
Otherwise, use trig functions.
d _ y1c - x1s.  Subsitute in a point to find d.
recall y*c - x*s - d = 0;
Calculates the different parts of a line given c, s and d.
IF cosineOfTheta # 0 THEN { -- find its slope and y intercept.
line.slope _ sineOfTheta/cosineOfTheta;
y intercept occurs when x = 0, ie when y*c = d.  y = d/c;
line.yInt _ line.d/line.c};
Find a line which is perpendicular to "line" and passes thru "pt".  Useful for dropping perpendiculars.
If line has the form: y*cos(theta) - x*sin(theta) - d = 0, then normalLine will have the form:
y*cos(theta+90) - x*sin(theta+90) - D = 0;
or -y*sin(theta) - (x*cos(theta)) - D = 0;
to find D, we substitute in pt:
D = -pt.y*sin(theta) - pt.x*cos(theta);
IF normalLine.c #0 THEN { -- compute slope and yInt
normalLine.slope _ normalLine.s/normalLine.c;
line.yInt _ normalLine.d/normalLine.c
};
parallelLine.slope _ line.slope;
parallelLine.slope _ line.slope;
Makes a new line that results by transforming line by transform.
Imager.Trans[dc];
Draw2d.Line[dc, [0.0, 0.0], [p2.x - p1.x, p2.y - p1.y], solid];
Making Edges
Making Rays
We can take advantage of the horizontal and vertical lines of the box to do an easy intersection test.  Note that we are really testing for line intersections rather than ray intersections.
The top line has equation y = ymax.  If ray.d.y = 0, we don't hit this line.  Otherwise, we use y(t) = ymax = ray.p.y+t*ray.d.y;  Solve for t:  t = (ymax-ray.p.y)/ray.d.y.
Top Line
Bottom Line
Right Line
Left Line
Intersections
To ensure no errors of more than 0.072 screen dots in a picture of size 14 inches by 14 inches, our angles in radians must be accurate to (theta*1008 < 0.072) theta <  7.142857e-5).  In degrees, this is 4.092559e-3.  I will use 4e-4 for extra accuracy.  e-5 results in determinant = 0.0 for Window.script.  (Bier, January 7, 1987)
If line1 is of the form: c1*y - s1*x - d1 = 0;
and line2 of the form: c2*y - s2*x - d2 = 0;
then we solve simultaneously.
x = (c2d1 - c1d2)/(s2c1 -s1c2);
y = (s2d1 - s1d2)/(s2c1 - s1c2);
determinant should not be zero since the lines are not parallel.
Find the intersection of line with the line of seg.  See if this point is within the bounds of seg.
Find the intersection of e1, with e2.  See if this point is within the bounds of e1 and e2.
Direction and Distance for Lines
Because of our choice or representation for a Line, plugging the point into the line equation gives us the signed distance.
ie.  distance = y*cos(theta) - x*sin(theta) - d;
Because of our choice or representation for a Line, plugging the point into the line equation gives us the signed distance.
ie.  distance = y*cos(theta) - x*sin(theta) - d;
We drop a normal from the point onto the line and find where it hits.
The line equation of the normal we drop can be found using FillLineAsNormal above.
We will have line equations:
c*y - s*x - d = 0.  The vector v = [c, s] is the unit vector parallel to line.  The vector
l = [-s, c] is the unit vector 90 degrees counter-clockwise of v.  If pt is distance D from line (along l), then the new point we want is pt-D*l.
This routine takes 4 mults, 4 adds.
Finds any old point on line and returns it.

Returns the direction vector of line.
Distance for Edges
Look for an obvious winner first.  If that fails, do math.
perpendicular distance if possible, else distance to nearest endpoint.
Perpendicular distance if possible, else distance to nearest endpoint.
Assumes pt is on edge.line.  Is it on edge?
Distance for Points
point is a point to the left of the directed line, on the normal to the line which intersects the line at pOnLine.  If distance is negative, the point will be to the right of the directed line.
Method:  The point we want will be at the intersection these two lines
1) The line parallel to "line" but distance to its left
2) The line perpendicular to "line" at pOnLine.
We can generate both of these easily as follows:
UTILITY FUNCTIONS
Κ§˜Icodeš
œ™Kš
œ0™0Kš
œ4™4š
œj™jK™$K™$—K˜š
Οk	˜	Kš
œR˜R—K˜šœ
œ
˜Kšœ;˜BKšœ
˜—K˜š
˜Kšœœ˜!Kšœœœ	˜Kšœ	œ˜%Kšœœœ	˜Kšœ	œ˜%Kšœœœ˜Kšœœ˜#Kšœœ˜#K˜—Ihead™K™šΟnœœ
œ
œ˜5Kšœœ
˜Kš
œ˜K˜—šžœœ
œ˜0Kš
œ˜Kš
œ˜Kš
œ˜Kš
œ˜Kš
œ™Kš
œ™Kš
œ˜K˜—K˜š
ž	œœ
œœœ˜<K™)Kšœœ˜(K˜K˜—šžœœ
œ"œœœœ˜kKš
œŒ
™Œ
Kšœ
œ#œ
œ˜OK˜—K˜šžœœ
œ ˜?Kš
œ™Kš
œΦ
™Φ
šœ	œΟc˜FKš
Ÿ-˜-—Kš
Πbc*™*Kšœœ˜Kšœœ˜Kš
Οb6™6šœ
œœŸ˜2šœ
œŸ˜$Kš
œ˜Kš
œ˜Kš
œ7™7Kš
œ˜—šœŸ˜Kš
œ˜Kš
œ˜Kš
œ8™8Kš
œ˜—Kš
œ˜Kš
œ
˜
—Kš
‘™šœ˜Kš
œ5˜5Kš
œ$˜$Kš
œ$˜$Kš
œ/™/Kš
œ#˜#Kš
œ˜—KšœŸ˜K˜—Kšžœœ
œœ
˜%šžœœ
œ)˜PKšœ	œ
˜š
œ
œœ
œœ˜7Kšœ˜Kš
œ˜Kš
œ
˜
—Kšœ.˜2Kš
œ$˜$Kš
œ$˜$Kš
œ#˜#K˜K˜—šžœœ
œ(œ˜bKš
œ™Kš
œ:™:Kš
œ˜Kš
œ˜Kš
œ˜Kš
œ;˜;šœœ%™>Kš
œ'™'Kš
œ9™9Kš
œ™—KšœŸ"˜%K˜—šžœœ
œœ˜QKš
œ)˜)Kš
œ$˜$Kš
œ$˜$Kš
œ#˜#KšœŸ#˜&K˜—šžœœ
œ.˜XK™gKš
œ^™^Kš
œ*™*Kš
œ*™*Kš
œ™Kš
œ'™'Kš
œ˜Kš
œ˜Kš
œ*˜*Kš
œ0˜0šœœŸ™3Kš
œ-™-Kš
œ%™%Kš
œ™—KšœŸ˜—K˜šžœœ
œœ˜PKš
œ˜Kš
œ˜Kš
œ˜Kš
œ ˜ Kš
œ ™ K˜K˜—šžœœ
œœ˜QKš
œ˜Kš
œ˜Kš
œ˜Kš
œ ˜ Kš
œ ™ K˜K˜—šžœœ
œN˜lKš
œ˜Kš
œ ˜ Kš
œ˜Kš
œ#˜#Kš
œ<˜<Kš
œG˜GKš
œ9˜9K˜—K˜šžœœ
œœ˜DKš
œ˜Kš
œ!˜!Kš
œ˜—šžœœ
œœ˜UKš
œ˜Kš
œ˜Kš
œ˜K˜—š
žœœ
œ(œœ˜gKš
œ˜Kš
œE˜EKš
œ˜—š
žœœ
œœœ˜WKš
œ˜Kš
œ-˜-K˜K˜—šžœœ
œœ˜]Kš
œ˜Kš
œ4˜4Kš
œ˜K˜—š
žœœ
œœœ˜UKš
œ!˜!Kš
œ-˜-Kš
œ˜K˜—š
žœœ
œœœ˜VKš
œ!˜!Kš
œ.˜.Kš
œ˜K˜—šž
œœ
œ?œ˜xKš
œ@™@Kš
œ ˜ Kš
œ1˜1K˜—K˜šžœœ
œLœ˜pKšœœ
˜K˜	Kšœœœ
œ
˜Kš
œ˜K˜šž
œœ˜Kš
œ˜Kš
œ™Kš
œ˜Kš
œ˜Kš
œ'˜'Kš
œ>˜>Kš
œ?™?K˜—Kš
œ ˜ Kš
œ<˜<Kš
œ˜Kš
œ"˜"Kš
œ&˜&Kš
œ?˜?šœœ˜Kš
œ˜Kš
œ˜K˜K˜—K˜K˜—L™šžœœ
œ
œ˜5Kšœœ
˜Kš
œ˜Kš
œ˜—šžœœ
œ˜0Kš
œ˜Kš
œ˜Kš
œ˜Kš
œ$˜$KšœŸ˜—šžœœ
œ ˜5Kšœœ
˜Kš
œ&˜&Kš
œ˜šœœ
˜Kšœ
œ6œ
˜NKšœ6œ
˜@—š
˜Kšœ
œ6œ
˜NKšœ6œ˜A—KšœŸ˜—šžœœ
œN˜lK˜K•StartOfExpansion4[m: 
ImagerTransformation.Transformation, v: 
VEC]š
œE˜EK–4[m: 
ImagerTransformation.Transformation, v: 
VEC]š
œA˜AKš
œ˜K˜—šž
œœ
œœ˜@Kš
œ˜Kš
œ˜KšœŸ˜—šž
œœ
œ?œ˜qKš
œ˜Kš
œ*˜*K˜—L™šž	œœ
œ'œ˜SKšœœ"˜+K˜—šžœœ
œœ˜GKšœœ'˜0K˜K˜—š
žœœœœœ˜?Kšœ
œ˜"K˜K˜—šžœœ
œ$œœ	œ
œœ
œ˜{K™½
Kšœœ
˜Kšœ	œ
˜˜
Kš
œ«
™«
Kšœ
œœŸ˜:™Kš
œ˜Kš
œ˜šœœœŸ˜AK˜Kš
œ˜K˜——™Kš
œ˜Kš
œ˜šœœœŸ˜AK˜Kš
œ˜K˜—K˜—Kšœ
œœŸ˜:™
šœœ˜Kš
œ˜š
œœœ
œ(œ˜PKš
œ˜šœœœŸ˜AK˜Kš
œ˜K˜—K˜—K˜——™	šœœ˜Kš
œ˜š
œœœ
œ(œ˜PKš
œ˜šœœœŸ˜AK˜Kš
œ˜K˜—K˜—K˜—K˜—šœœ˜šœœ˜ Kšœœ
˜Kš
œ˜Kš
œ˜K˜—KšœŸ˜——KšœŸ˜—š
žœœ
œœœ˜GKš
œ"˜"Kš
œ"˜"K˜—L™
š
ž
œœœœœ˜-Kšœ	œ
˜Kšœ
œ˜K˜K˜—K˜š
ž
œœ
œœ!œ˜aKš
œΚ™ΚKš
œ.™.Kš
œ,™,Kš
œ™Kš
œ™Kš
œ ™ Kšœ
œ
˜Kšœ	œ˜Kšœœ
˜Kšœ<œ
œœ˜]š
œ0˜0Kš
œ@™@—Kš
œA˜AKš
œA˜AKš
œ˜K˜—šžœœ
œœœœ˜QKšœœœ˜;šœŸR˜XKšœœ
˜Kš
œ˜—Kš
œ˜K˜—š
ž
œœ
œœœ˜bKš
œc™cKš
œ˜Kšœ
œ
˜Kš
œ9˜9Kšœ
œ
œœ˜(Kšœœ˜'Kš
œ˜K˜—š
ž
œœ
œœœ˜YKš
œ[™[Kš
œ˜Kš
œ2˜2Kšœœ
œ
˜Kšœœ˜%Kš
œ˜—L™ š
žœœ
œœœ˜MKš
œ{™{Kš
œ0™0Kš
œ'˜'KšœŸ˜K˜—š
žœœ
œœœ˜GKš
œ{™{Kš
œ0™0Kšœœ%˜,KšœŸ˜K˜K˜—šžœœ
œœ˜WKš
œE™EKš
œR™RKš
œ™Kš
œμ
™μ
Kšœ‘™#Kšœœ&˜-Kš
œ ˜ Kš
œ ˜ Kš
œ˜K˜—šžœœ
œœ˜=K™+šœ
œœ	œ˜#Kš
œ˜Kš
œ˜K˜
—šœ˜K˜Kš
œ˜K˜—K˜K™—šžœœ
œœ˜IK™%Kš
œ˜Kš
œ˜K˜—L™šžœœ
œœ˜RKš
œ:™:šœ
œœ˜6Kš
œ
œœœ
œ˜IKšœ
œ˜—š
˜Kš
œ
œœœ
œ
˜FKšœ
œ˜—š
˜š
œ	˜	šœX˜ZKšœ˜Kšœ˜———Kš
œ˜K˜—š
ž œœ
œœœ˜iKšœ œ
˜%Kš
œ=˜=Kš
œ;˜;Kšœ
œ"˜,Kš
œ˜—K˜šžœœ
œœ˜SKš
œ6˜6Kšœ$œ˜?Kšœ$˜(K˜K˜—š
žœœ
œœœ˜UKš
œF™FKš
œ6˜6Kš
œ˜Kšœ$œœ˜Všœ˜Kš
œ)˜)Kš
œ2˜2Kš
œ˜—Kš
œ˜K˜—š
žœœ
œœœ˜cKš
œF™FKš
œ6˜6Kšœ#˜%Kšœ0˜4Kš
œ3˜3Kšœ>˜BKš
œ˜K˜—šžœœœœ˜>Kšœ
Οu
œœ
˜	Kšœ
’
œ(˜*Kšœ
’
œ˜K˜—š
žœœ
œœœ˜GKš
œ+™+Kš	œ
œœœŸ+˜uKšœ)˜/Kšœ
Ÿ˜Kšœ*˜0KšœŸ˜—L™š
žœœ
œœœ˜NKš
œK˜KKš
œ˜—š
žœœ
œœœ˜UKš
œ=˜=Kš
œ˜—š
žœœ
œœœ˜dKš
œΑ
™Α
Kš
œF™FKš
œ7™7Kš
œ/™/Kš
œ0™0K˜Kš
œ˜Kšœ
œ
˜Kš
œ!˜!Kš
œ˜Kš
œJ˜JKš
œ7˜7Kš
œ:˜:Kšœ
œ
œŸ'˜?Kš
œ˜—K˜Kš
œ™K˜šžœœ
œœœœ˜;šœ˜
Kšœœœ˜(Kšœœ˜Kšœœœ˜(Kšœœ
˜Kš
œ˜——K˜Kšœ
˜K˜—…—8~__