INTERFACE HLR3;

(* 
  Oriented projective geometry in three dimensions.
  Created 93-04-18 by Marcos C. Carrard.
  Based on H3.pas by J. Stolfi.
  Last edited by stolfi *)
   
IMPORT LR3, LR4, LR4x4, LR6;

TYPE
  Point = RECORD c: LR4.T END;  (* c[0..3] are coords [w,x,y,z]. *)
  Plane = RECORD f: LR4.T END;  (* f[0..3] are coeffs <W,X,Y,Z>. *)
  Line  = RECORD k: LR6.T END;  (* k[0..5] are Plücker coords [wx,wy,xy,wz,xz,yz] *)
  Sign =  [-1..1];
  
PROCEDURE FromCartesian(READONLY c: LR3.T): Point;
  (*
    Point on "hither" half of space (i.e, with positive weight)
    whose Cartesian coordinates are "c". *)
    
PROCEDURE ToCartesian(READONLY p: Point): LR3.T;
  (*
    Cartesian coordinates of point "p" (which must be finite). *)

PROCEDURE Side(READONLY p: Point; READONLY Q: Plane): Sign; 
  (*
    Returns the position of point "p" relative to plane "Q": 
    0 on the plane, +1 in positive halfspace, -1 in negative halfspace.
    May give inconsistent results for points very close to the plane. *)

PROCEDURE Orient(READONLY p, q, r, s: Point): Sign;
  (*
    Returns the orientation (handedness) of the tetrahedron "p q r s".
    Specifically, the tetrahedron "[1 0 0 0] [0 1 0 0] [0 0 1 0] [0 0 0 1]"
    has positive orientation. Returns 0 iff the points are coplanar.
    May give inconsistent results for points very close to coplanar. *)
    
PROCEDURE LineFromTwoPoints(READONLY p, q: Point): Line;
  (*
    The line from "p" to "q". *)

PROCEDURE PlaneFromThreePoints(READONLY p, q, r: Point): Plane;
  (*
    The plane through "p", "q", and "r". *)
    
PROCEDURE PlaneFromLineAndPoint(READONLY n: Line; READONLY r: Point): Plane;
  (*
    The plane through "n" and "r". *)
    
PROCEDURE LineFromTwoPlanes(READONLY P, Q: Plane): Line;
  (*
    The line where "P" meets "Q". *)

PROCEDURE PointFromThreePlanes(READONLY P, Q, R: Plane): Point;
  (*
    The point where "P", "Q", and "R" meet. *)
    
PROCEDURE PointFromLineAndPlane(READONLY n: Line; R: Plane): Point;
  (*
    The point where "n" meets "R". *)

PROCEDURE Dir(READONLY frm, tto: Point): LR3.T;
  (*
    Direction (a unit-length vector) of point "tto" seen from point "frm".
    Works even if one of them is at infinity.  Does not work if both
    are at infinity, or coincident, or antipodal. *)
    
PROCEDURE Dist(READONLY a, b: Point): LONGREAL;
  (*
    Distance between "a" and "b", which must lie in the front half-plane. *)
    
PROCEDURE DistSqr(READONLY a, b: Point): LONGREAL;
  (*
    Distance squared between "a" and "b", which must lie in the front half-plane. *)
    
PROCEDURE Normal(READONLY P: Plane): LR3.T;
  (*
    The normal direction of plane "P", on the hither side, pointing into 
    "P"'s positive halfspace.  Assumes "P" is not at infinity. *)
    
(* PROJECTIVE MAPS *)

TYPE
  PMap =  RECORD dir: LR4x4.T; inv: LR4x4.T END;  
    (* 
      "dir" is the map's matrix, "inv" is its inverse. *)

PROCEDURE MapPoint(READONLY p: Point; READONLY m: PMap): Point;
  (*
    Applies projective map "m" to point "p". *)

PROCEDURE InvMapPoint(READONLY p: Point; READONLY m: PMap): Point;
  (*
    Applies the inverse of projective map "m" to point "p". *)

PROCEDURE MapPlane(READONLY P: Plane; READONLY m: PMap): Plane;
  (*
    Applies projective map "m" to plane "P" *)

PROCEDURE InvMapPlane(READONLY P: Plane; READONLY m: PMap): Plane;
  (*
    Applies the inverse of projective map "m" to plane "P" *)

PROCEDURE CompMap(READONLY m, n: PMap): PMap;
  (* 
    Returns the composition of "m" and "n", applied in that order *)
    
PROCEDURE InvMap(READONLY m: PMap): PMap;
  (*
    Returns the inverse of map "m". *)

PROCEDURE PerspMap(READONLY obs, foc, upp: Point): PMap;
  (* 
    Computes a perspective transformation with given viewing parameters:
      "obs" = the scenesys coords of imagesys (0,0,d) (the observer);
      "foc" = the scenesys coords of imagesys (0,0,0) (the image focus);
      "up"  = the scenesys coords of some point 
               with imagesys X = 0, Y > 0 (up reference). *)

END HLR3.