-- COGHomo2Impl.mesa: Basic geometrical operations on the two-sided plane in homogeneous coordinates (implem.)
-- last modified by Stolfi - September 12, 1983 4:45 pm
DIRECTORY
COGHomo2,
Real USING [SqRt, TrappingNaN, LargestNumber],
RealFns USING [AlmostZero],
COGCart2 USING [Point, PositiveAngle, Vector, ScaleFactors, Add];
COGHomo2Impl: CEDAR PROGRAM
IMPORTS Real, RealFns, COGHomo2, COGCart2
EXPORTS COGHomo2 =
BEGIN OPEN Real, RealFns, Homo: COGHomo2, Cart: COGCart2;
-- HOMOGENEOUS PREDICATES AND OPERATIONS - - - - - - - - - - - - - - - - -
IndPt:
PUBLIC
PROC [u: Point]
RETURNS [
BOOLEAN] =
{RETURN[u.w = 0.0 AND u.x = 0.0 AND u.y = 0.0]};
IndLn:
PUBLIC
PROC [r: Line]
RETURNS [
BOOLEAN] =
{RETURN[r.w = 0.0 AND r.x = 0.0 AND r.y = 0.0]};
Anti:
PUBLIC
PROC [a: Point]
RETURNS [Point] =
{RETURN [[-a.x, -a.y, -a.w]]};
Opp:
PUBLIC
PROC [r: Line]
RETURNS [
Line] =
INLINE
{RETURN [[-r.x, -r.y, -r.w]]};
Incident:
PUBLIC
PROC [p: Homo.Point, r: Homo.Line, tol: [-126..127] ← -8]
RETURNS [
BOOLEAN] =
{RETURN[ AlmostZero[r.x*p.x+r.y*p.y+r.w*p.w, tol]]};
LeftOf:
PUBLIC
PROC [p: Point, r: Line]
RETURNS [
R:
BOOLEAN] =
{RETURN[r.x*p.x+r.y*p.y+r.w*p.w > 0]};
RightOf:
PUBLIC
PROC [p: Point, r: Line]
RETURNS [
BOOLEAN] =
{RETURN[r.x*p.x+r.y*p.y+r.w*p.w < 0]};
SameSide:
PUBLIC
PROC [p, q: Homo.Point, r: Homo.Line]
RETURNS [
BOOLEAN] =
{RETURN[(r.x*p.x+r.y*p.y+r.w*p.w)*(r.x*q.x+r.y*q.y+r.w*q.w) >= 0]};
StrictlySameSide:
PUBLIC
PROC [p, q: Homo.Point, r: Homo.Line]
RETURNS [
BOOLEAN] =
{RETURN[(r.x*p.x+r.y*p.y+r.w*p.w)*(r.x*q.x+r.y*q.y+r.w*q.w) > 0]};
PPLine:
PUBLIC
PROC [p, q: Homo.Point]
RETURNS [Homo.Line] =
{RETURN [[p.y*q.w-p.w*q.y, -p.x*q.w+p.w*q.x, p.x*q.y-p.y*q.x]]};
Intercept:
PUBLIC
PROC [r, s: Homo.Line]
RETURNS [Homo.Point] =
{RETURN [[r.y*s.w-r.w*s.y, -r.x*s.w+r.w*s.x, r.x*s.y-r.y*s.x]]};
PosOr:
PUBLIC
PROC [p, q, r: Homo.Point]
RETURNS [
BOOLEAN] =
{RETURN [r.w*(p.x*q.y-q.x*p.y)-q.w*(p.x*r.y-r.x*p.y)+p.w*(q.x*r.y-r.x*q.y) > 0]};
NegOr:
PUBLIC
PROC [p, q, r: Homo.Point]
RETURNS [
BOOLEAN] =
{RETURN [r.w*(p.x*q.y-q.x*p.y)-q.w*(p.x*r.y-r.x*p.y)+p.w*(q.x*r.y-r.x*q.y) < 0]};
HWay:
PUBLIC
PROC [alpha:
REAL, p, q: Homo.Point]
RETURNS [Homo.Point] =
BEGIN
mu: REAL = (1.0-alpha)*SqRt[q.x*q.x+q.y*q.y+q.w*q.w];
nu: REAL = alpha*SqRt[p.x*p.x+p.y*p.y+p.w*p.w];
RETURN [[mu*p.x+nu*q.x, mu*p.y+nu*q.y, mu*p.w+nu*q.w]]
END;
HClos:
PUBLIC
PROC [p, q: Homo.Point]
RETURNS [
REAL] =
{RETURN [(p.x*q.x+p.y*q.y+p.w*q.w)/
SqRt[(p.x*p.x+p.y*p.y+p.w*p.w)*(q.x*q.x+q.y*q.y+q.w*q.w)]]};
HClosSq:
PUBLIC
PROC [p, q: Homo.Point]
RETURNS [
REAL] =
{prod: REAL =p.x*q.x+p.y*q.y+p.w*q.w;
RETURN [prod*prod/((p.x*p.x+p.y*p.y+p.w*p.w)*(q.x*q.x+q.y*q.y+q.w*q.w))]};
HLPDist:
PUBLIC
PROC [r: Homo.Line, p: Homo.Point]
RETURNS [
REAL] =
{RETURN [(r.x*p.x+r.y*p.y+r.w*p.w)/
SqRt[(r.x*r.x+r.y*r.y+r.w*r.w)*(p.x*p.x+p.y*p.y+p.w*p.w)]]};
HLPDistSq:
PUBLIC
PROC [r: Homo.Line, p: Homo.Point]
RETURNS [
REAL] =
{prod: REAL = r.x*p.x+r.y*p.y+r.w*p.w;
RETURN [prod*prod/((r.x*r.x+r.y*r.y+r.w*r.w)*(p.x*p.x+p.y*p.y+p.w*p.w))]};
-- OPERATIONS RELATED TO THE CARTESIAN PROJECTION - - - - - - - - - - - - -
CartPt:
PUBLIC
PROC [p: Point]
RETURNS [Cart.Point] =
{RETURN [[p.x/p.w, p.y/p.w]]};
Way:
PUBLIC
PROC [alpha:
REAL, p, q: Homo.Point]
RETURNS [Homo.Point] =
BEGIN
u: REAL = (1.0-alpha)*q.w;
v: REAL = alpha*p.w;
RETURN [[u*p.x+v*q.x, u*p.y+v*q.y, p.w*q.w]]
END;
Dist:
PUBLIC
PROC [p, q: Homo.Point]
RETURNS [
REAL] =
{RETURN [Length[Homo.Sub[p,q]]]};
DistSq:
PUBLIC
PROC [p, q: Homo.Point]
RETURNS [
REAL] =
{v: Homo.Vector=[q.x*p.w-p.x*q.w, q.y*p.w-p.y*q.w, q.w*p.w];
RETURN [(v.x*v.x+v.y*v.y)/(v.w*v.w)]};
Clos:
PUBLIC
PROC [p, q: Homo.Point]
RETURNS [
REAL] =
{RETURN [Shrt[Sub[p,q]]]};
ClosSq:
PUBLIC
PROC [p, q: Homo.Point]
RETURNS [
REAL] =
{v: Homo.Vector=[q.x*p.w-p.x*q.w, q.y*p.w-p.y*q.w, q.w*p.w];
RETURN [(v.w*v.w)/(v.x*v.x+v.y*v.y+v.w*v.w)]};
Add:
PUBLIC
PROC [a, b: Vector]
RETURNS [Vector] =
{RETURN [[a.x*b.w+b.x*a.w, a.y*b.w+b.y*a.w, a.w*b.w]]};
Sub:
PUBLIC
PROC [a, b: Vector]
RETURNS [Vector] =
{RETURN [[a.x*b.w-b.x*a.w, a.y*b.w-b.y*a.w, a.w*b.w]]};
Mul:
PUBLIC
PROC [a:
REAL, u: Vector]
RETURNS [w: Vector] =
{RETURN [[u.x*a, u.y*a, u.w]]};
Length:
PUBLIC
PROC [v: Homo.Vector]
RETURNS [
REAL] =
TRUSTED
{RETURN [SqRt[v.x*v.x+v.y*v.y]/v.w]};
LengthSq:
PUBLIC
PROC [v: Homo.Vector]
RETURNS [
REAL] =
{RETURN [(v.x*v.x+v.y*v.y)/(v.w*v.w)]};
Normalize:
PUBLIC
PROC [v: Homo.Vector]
RETURNS [Homo.Vector]
=
{RETURN [[r.x, r.y, SqRt[r.x*r.x+r.y*r.y]]]};
Shrt:
PUBLIC
PROC [v: Vector]
RETURNS [
REAL] =
{RETURN [v.w/SqRt[v.x*v.x+v.y*v.y+v.w*v.w]]};
ShrtSq:
PUBLIC
PROC [v: Vector]
RETURNS [
REAL] =
{RETURN [(v.w*v.w)/(v.x*v.x+v.y*v.y+v.w*v.w)]};
LPDist:
PUBLIC
PROC [r: Homo.Line, p: Homo.Point]
RETURNS [
REAL] =
BEGIN
prod: REAL = p.x*r.x+p.y*r.y+p.w*r.w;
mod: REAL = r.x*r.x+r.y*r.y;
RETURN [prod/(ABS[p.w]*SqRt[mod])]
END;
Bisector:
PUBLIC
PROC [p, q: Homo.Point]
RETURNS [Homo.Line] =
BEGIN
v: Homo.Vector = Homo.Sub [q, p];
m: Homo.Point = Homo.Add [p, q];
RETURN [[2.0*v.x*m.w, 2.0*v.y*m.w, -v.x*m.x-v.y*m.y]]
END;
PerpVector:
PUBLIC
PROC [r: Homo.Line]
RETURNS [Homo.Vector]
=
{RETURN [[r.x, r.y, 1.0]]};
ParVector:
PUBLIC
PROC [r: Homo.Line]
RETURNS [Homo.Vector]
=
{RETURN [[r.y, -r.x, 1.0]]};
PNLine:
PUBLIC PROC [p: Homo.Point, v: Homo.Vector]
RETURNS [Homo.Line] =
{RETURN [v.x*p.w, v.y*p.w, -v.x*p.x-v.y*p.y]};
{STOPPED HERE September 12, 1983 4:44 pm}
NorthOf:
PUBLIC
PROC [p, q: Homo.Point]
RETURNS [
R:
BOOL] =
{RETURN [p.y/p.w > q.y/q.w]};
SouthOf:
PUBLIC
PROC [p, q: Homo.Point]
RETURNS [
R:
BOOL] =
{RETURN [p.y/p.w > q.y/q.w]};
Below:
PUBLIC
PROC [p, q: Homo.Point]
RETURNS [
R:
BOOL] =
BEGIN
ydir: REAL;
ydir ← IF Homo.InfPt[p] AND Homo.InfPt[q] THEN q.y - p.y ELSE Homo.PtPtDir[p, q].y;
RETURN [ydir > 0.0]
END;
ToTheLeftOf:
PUBLIC
PROC [p, q: Homo.Point]
RETURNS [
R:
BOOL] =
BEGIN
xdir: REAL = IF Homo.InfPt[p] AND Homo.InfPt[q] THEN q.x - p.x ELSE Homo.PtPtDir[p, q].x;
RETURN [xdir > 0.0]
END;
Precedes:
PUBLIC
PROC [p, q: Homo.Point]
RETURNS [
R:
BOOL] =
BEGIN
dir: Cart.Vector =
IF Homo.InfPt[p] AND Homo.InfPt[q] THEN
[q.x - p.x, q.y - p.y]
ELSE
Homo.PtPtDir[p, q];
RETURN [dir.x > 0 OR dir.x = 0 AND dir.y > 0]
END;
CircumCenter:
PUBLIC
PROC [p, q, r: Homo.Point]
RETURNS [Homo.Vector] =
BEGIN
-- assumes p, q, r are in counterclockwise order
IF Homo.InfPt [p] THEN RETURN [InfiniteCC[p, q, r]]
ELSE IF Homo.InfPt [q] THEN RETURN [InfiniteCC[q, r, p]]
ELSE IF Homo.InfPt [r] THEN RETURN [InfiniteCC[r, p, q]]
ELSE RETURN [Homo.Intercept[Homo.Bisector[p, q], Homo.Bisector[q, r]]]
END;
InCircle:
PUBLIC
PROC [p, q, r, s: Homo.Point]
RETURNS [
R:
BOOL] =
BEGIN -- used to be ShouldConnect24
-- Checks whether the diagonal qs should be connected in the convex Delaunay
-- face pqrs, by intersecting the bisectors of pr and qs and checking which pair
-- is closest to the intersection.
a, b: Homo.Line;
o: Homo.Point;
IF q.w = 0.0 OR s.w = 0.0 THEN RETURN [FALSE];
IF p.w = 0.0 OR r.w = 0.0 THEN RETURN [TRUE];
a ← Homo.Bisector[p, r];
b ← Homo.Bisector[q, s];
o ← Homo.Intercept[a, b];
R ← Homo.DistSq[o, p] > Homo.DistSq [o, q]
END;
ScalePoint:
PUBLIC
PROC [p: Homo.Point, sf: Cart.ScaleFactors]
RETURNS [ps: Cart.Point] =
-- Transforms point p according to the scale factors sf. If p is infinite, will map it into a point in the same direction (as seen from the scaled origin) and at a very large distance. If p is indeterminate, will map it as if it were [0, 0].
{s: REAL = IF Homo.FinPt[p] THEN 1.0/p.w ELSE 1.0e5;
RETURN [[p.x*s*sf.sx + sf.tx, p.y*s*sf.sy + sf.ty]]};
ScaleVector:
PUBLIC
PROC [v: Homo.Vector, sf: Cart.ScaleFactors]
RETURNS [vs: Cart.Point] =
-- Transforms vector v according to the scale factors sf (will not displace origin). If p is infinite, will map it into a vector in the same direction (as seen from the scaled origin) and at a very large distance. If v is indeterminate, will map it as if it were [0, 0].
{s: REAL = IF Homo.FinPt[v] THEN 1.0/v.w ELSE 1.0e5;
RETURN [[v.x*s*sf.sx, v.y*s*sf.sy]]};
ScaleSeg:
PUBLIC
PROC [p, q: Homo.Point, sf: Cart.ScaleFactors]
RETURNS [ps, qs: Cart.Point] =
BEGIN
IF Homo.FinPt [p] THEN
{ps ← Homo.ScalePoint [p, sf];
IF Homo.FinPt [q] THEN
{qs ← Homo.ScalePoint [q, sf]}
ELSE
{qs ← Cart.Add[ps, Homo.ScaleVector [q, sf]]}
}
ELSE IF Homo.FinPt[q] THEN
{qs ← Homo.ScalePoint [q, sf];
ps ← Cart.Add [qs, Homo.ScaleVector [p, sf]]
}
ELSE
{ps ← qs ← [sf.tx, sf.ty]}
END;
END.