-- CGMatrixImpl.mesa
-- Last changed by Doug Wyatt, August 23, 1982 4:14 pm

DIRECTORY
  CGMatrix USING [ErrorType, Id, Ref, Rep],
  CGStorage USING [qZone],
  GraphicsBasic USING [Transformation, Vec],
  Real USING [RealException];

CGMatrixImpl: CEDAR MONITOR
IMPORTS CGStorage, Real
EXPORTS CGMatrix = { OPEN CGMatrix, GraphicsBasic;

repZone: ZONE = CGStorage.qZone;

lastid: Id _ 1; -- Monitored global!

identityMatrix: Rep = [m: [a: 1, b: 0, c: 0, d: 1, e: 0, f: 0],
  id: 0, trivial: TRUE, rectangular: TRUE, determinant: 1];

Error: PUBLIC ERROR[type: ErrorType] = CODE;

NewID: ENTRY PROC RETURNS[Id] = INLINE {
  IF lastid<LAST[Id] THEN lastid _ lastid + 1
  ELSE lastid _ 1; -- just wrap around
  RETURN[lastid];
  };

SetHints: PROC[self: Ref] = {
  trivial: BOOLEAN _ self.m.a=1 AND self.m.b=0 AND self.m.c=0 AND self.m.d=1;
  self.id _ IF trivial THEN 1 ELSE NewID[];
  self.trivial _ trivial;
  self.rectangular _ trivial OR (self.m.b=0 AND self.m.c=0) OR (self.m.a=0 AND self.m.d=0);
  self.determinant _ self.m.a*self.m.d - self.m.b*self.m.c;
  };

New: PUBLIC PROC RETURNS[Ref] = {
  self: Ref _ repZone.NEW[Rep _ identityMatrix];
  RETURN[self];
  };

Make: PUBLIC PROC[m: Transformation] RETURNS[Ref] = {
  self: Ref _ repZone.NEW[Rep _ [m: m,
    id: 0, trivial: FALSE, rectangular: FALSE, determinant: 0]];
  SetHints[self];
  RETURN[self];
  };

MapH: PUBLIC PROC[self: Ref, v: Vec, h: [0..1]] RETURNS[Vec] = {
  x,y: REAL;
  x _ v.x*self.m.a + v.y*self.m.c;
  y _ v.x*self.m.b + v.y*self.m.d;
  RETURN[IF h=1 THEN [x + self.m.e, y + self.m.f] ELSE [x,y]];
  };

InvH: PUBLIC PROC[self: Ref, v: Vec, h: [0..1]] RETURNS[Vec] = {
  ENABLE Real.RealException => GOTO Singular;
  u: Vec _ IF h=1 THEN Inv0[self,v.x-self.m.e,v.y-self.m.f] ELSE Inv0[self,v.x,v.y];
  RETURN[u]; EXITS Singular => ERROR Error[singularMatrix];
  };

Inv0: PROC[self: Ref, x,y: REAL] RETURNS[Vec] = {
  det: REAL _ self.determinant;
  RETURN[[(x*self.m.d-y*self.m.c)/det,(y*self.m.a-x*self.m.b)/det]];
  };

Concat: PUBLIC PROC[self: Ref, a,b,c,d: REAL] = {
  IF self.trivial THEN {
    self.m.a _ a; self.m.b _ b;
    self.m.c _ c; self.m.d _ d;
    }
  ELSE {
    aa: REAL _ self.m.a; bb: REAL _ self.m.b;
    cc: REAL _ self.m.c; dd: REAL _ self.m.d;
    self.m.a _ a*aa + b*cc; self.m.b _ a*bb + b*dd;
    self.m.c _ c*aa + d*cc; self.m.d _ c*bb + d*dd;
    };
  SetHints[self];
  };

ConcatTransformation: PUBLIC PROC[self: Ref, m: Transformation] = {
  IF self.trivial THEN {
    self.m.a _ m.a; self.m.b _ m.b;
    self.m.c _ m.c; self.m.d _ m.d;
    self.m.e _ self.m.e + m.e;
    self.m.f _ self.m.f + m.f;
    }
  ELSE {
    aa: REAL _ self.m.a; bb: REAL _ self.m.b;
    cc: REAL _ self.m.c; dd: REAL _ self.m.d;
    self.m.a _ m.a*aa + m.b*cc; self.m.b _ m.a*bb + m.b*dd;
    self.m.c _ m.c*aa + m.d*cc; self.m.d _ m.c*bb + m.d*dd;
    self.m.e _ self.m.e + m.e*aa + m.f*cc;
    self.m.f _ self.m.f + m.e*bb + m.f*dd;
    };
  SetHints[self];
  };

Invert: PUBLIC PROC[self: Ref] = {
  aa: REAL _ self.m.a; bb: REAL _ self.m.b;
  cc: REAL _ self.m.c; dd: REAL _ self.m.d;
  ee: REAL _ self.m.e; ff: REAL _ self.m.f;
  det: REAL _ aa*dd - bb*cc;
  self.m.a _ dd/det;
  self.m.b _ -bb/det;
  self.m.c _ -cc/det;
  self.m.d _ aa/det;
  self.m.e _ (cc*ff-dd*ee)/det;
  self.m.f _ (bb*ee-aa*ff)/det;
  SetHints[self];
  };

Copy: PUBLIC PROC[self: Ref] RETURNS[Ref] = {
  RETURN[repZone.NEW[Rep _ self^]];
  };

}.