<<-- 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)]}; <<>> <<-- VECTORS >> <<>> 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)]}; -- LINES 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; -- LINES AND VECTORS 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.