Last Edited by: Stolfi June 13, 1984 8:59:59 am PDT

JunoStorage USING [Point, Coords],
JunoMatrix;

JunoMatrix

BEGIN
IF m = NIL THEN m ← NEW [Matrix];
FOR i: INTEGER IN [1..3] DO FOR j: INTEGER IN [1..3] DO
m^[i][j] ← IF i = j THEN 1 ELSE 0
ENDLOOP ENDLOOP;
RETURN[m]
END;

BEGIN
IF m = NIL THEN m ← NEW [Matrix];
m^[1][3] ← frame.org.coords.x;
m^[2][3] ← frame.org.coords.y;
m^[3][3] ← 1;

IF frame.hor=NIL THEN
{m^[1][1] ← 300;
m^[2][1] ← 0}
ELSE
{m^[1][1] ← frame.hor.coords.x-frame.org.coords.x;
m^[2][1] ← frame.hor.coords.y-frame.org.coords.y};
m^[3][1] ← 0;
IF frame.ver=NIL THEN
{m^[1][2] ← -m^[2][1];
m^[2][2] ← m^[1][1]}
ELSE
{m^[1][2] ← frame.ver.coords.x-frame.org.coords.x;
m^[2][2] ← frame.ver.coords.y-frame.org.coords.y};
m^[3][2] ← 0;
RETURN[m]
END;

BEGIN
-- Inverts m^ into mInv^ by pivoting three times; or sets "singular" flag.
-- Uses work if not NIL, otherwise allocates new matrix for result.
i, j, k, l: INTEGER;
c: ARRAY [1..3] OF INTEGER;
pivoted: ARRAY [1..3] OF BOOLEAN ← [FALSE, FALSE, FALSE];
p: REAL;
-- k is the row in which we are pivoting.
-- l is the column in which we are pivoting.
-- i and j are miscellaneous row and column indices respectively
-- c[i] is the column of the pivot in the ith row.
-- p is the reciprocal of the pivot element; also used as temp for swapping.
IF work = NIL THEN work ← NEW[Matrix];
FOR k IN [1..3] DO
-- set l so m^[k,l] is largest of m^[k,1], m^[k, 2], m^[k, 3], excluding
-- columns in which we have already pivoted.
p ← 0;
FOR j IN [1 .. 3]
DO IF ABS[m^[k][j]] >= p AND NOT pivoted[j] THEN {l ← j; p ← ABS[m^[k][l]]} ENDLOOP;
-- We will pivot at m^[k,l], if it is not too small:
IF ABS[m^[k][l]] < .0001 THEN RETURN[NIL, TRUE];
c[k] ← l; pivoted[l] ← TRUE;
p ← 1.0 / m^[k][l]; m^[k][l] ← 1.0;
-- divide everything in pivot row by the pivot element:
FOR j IN [1..3] DO m^[k][j] ← m^[k][j] * p ENDLOOP;
FOR i IN [1..3] DO
IF i # k THEN
  FOR j IN [1..3] DO
  IF j # l THEN -- for each m^[i,j] outside the pivot row and column
  m^[i][j] ← m^[i][j] - m^[i][l] * m^[k][j]; -- note that m^[k,j] was already * p.
ENDLOOP ENDLOOP;
-- Finally process pivot column:
FOR i IN [1..3] DO IF i # k THEN m^[i][l] ← -m^[i][l] * p ENDLOOP;
ENDLOOP;
-- Now we permute rows and columns:
FOR i IN [1..3] DO FOR j IN [1..3] DO work^[c[i]][j] ← m^[i][c[j]] ENDLOOP ENDLOOP;
RETURN[work, FALSE]
END;

BEGIN
i, j, k: INTEGER;
sum: REAL;
IF mc = NIL THEN mc ← NEW [Matrix];
FOR i IN [1..3] DO FOR j IN [1..3] DO
sum ← 0.0;
FOR k IN [1..3] DO sum ← sum + ma^[i][k] * mb^[k][j] ENDLOOP;
mc^[i][j] ← sum
ENDLOOP ENDLOOP;
RETURN[mc]
END;

BEGIN
-- We want
-- mat * [ src.org, src.hor, src.ver ] = [ dest.org, dest.hor, dest.ver ],
-- where the points are viewed as
-- column vectors with third component 1. Hence we compute the inverse
-- of src and multipy on the left by dest. But the pairs
-- (dest.hor, src.hor), (dest.ver, src.ver) may be missing, in which case they
-- are filled in by default to make the transformation a
-- translation (if both are missing) or a Euclidean motion (if just (c, sc) is missing).
-- The matrix mm is optional, and is used to store the result: if NIL, a new one is allocated.
-- Uses m1, mInv as work areas.
sing: BOOL;
mm ← GetFrameMatrix[src, mm];
[mInv, sing] ← InvertMatrix[mm, mInv];
IF sing THEN RETURN [NIL, sing];
m1 ← GetFrameMatrix[dest, m1];
mm ← MultiplyMatrix[m1, mInv, mm];
RETURN [mm, sing]
END;

BEGIN
-- Transforms coords.x, y by the given matrix
cMap.x ← coords.x * mat^[1][1] + coords.y * mat^[1][2] + mat^[1][3];
cMap.y ← coords.x * mat^[2][1] + coords.y * mat^[2][2] + mat^[2][3]
END;