Matrix3dImpl.mesa
Copyright © 1984 by Xerox Corporation. All rights reserved.
Bier, June 1982
Crow, October 31, 1985 12:31:44 pm PST
Bloomenthal, May 21, 1986 6:34:06 pm PDT
Cf ``Principles of Interactive Computer Graphics,'' by Newman and Sproull.
Left-handed eyespace coordinate system (right, up, depth);
right-handed world coordinate system (longitude, latitude, altitude).
DIRECTORY Matrix3d, RealFns, Rope, Vector3d;
Matrix3dImpl: CEDAR MONITOR
IMPORTS RealFns, Vector3d
EXPORTS Matrix3d
~ BEGIN
Pair:   TYPE ~ Vector3d.Pair;
Triple:  TYPE ~ Vector3d.Triple;
Quad:   TYPE ~ Vector3d.Quad;
Row:   TYPE ~ Matrix3d.Row;
Col:   TYPE ~ Matrix3d.Col;
Matrix:  TYPE ~ Matrix3d.Matrix;
MatrixRep: TYPE ~ Matrix3d.MatrixRep;
origin:  Triple ~ Matrix3d.origin;
singular:  PUBLIC ERROR ~ CODE;
Creation Operations
Copy: PUBLIC PROC [in: Matrix, out: Matrix ← NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out NEW[MatrixRep];
out^ ← in^;
RETURN[out];
};
Identity: PUBLIC PROC [out: Matrix ← NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out NEW[MatrixRep];
out^ ← [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]];
RETURN[out];
};
MakeFromTriad: PUBLIC PROC [v1, v2, v3: Triple, p: Triple ← origin, out: Matrix ← NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out NEW[MatrixRep];
v1 ← Vector3d.Normalize[v1];
v2 ← Vector3d.Normalize[v2];
v3 ← Vector3d.Normalize[v3];
out^ ← [[v1.x,v1.y,v1.z,0.0], [v2.x,v2.y,v2.z,0.0], [v3.x,v3.y,v3.z,0.0], [p.x,p.y,p.z,1.0]];
RETURN[out];
};
MakeScale: PUBLIC PROC [s: Triple, out: Matrix ← NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out NEW[MatrixRep];
out^ ← [[s.x, 0.0, 0.0, 0.0], [0.0, s.y, 0.0, 0.0], [0.0, 0.0, s.z, 0.0], [0.0, 0.0, 0.0, 1.0]];
RETURN[out];
};
MakeTranslate: PUBLIC PROC [p: Triple, out: Matrix ← NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out NEW[MatrixRep];
out^ ← [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [p.x, p.y, p.z, 1.0]];
RETURN[out];
};
MakePureRotate: PUBLIC PROC [axis: Triple, theta: REAL, degrees: BOOLTRUE, out: Matrix ← NIL] RETURNS [Matrix] ~ {
sin: REALIF degrees THEN RealFns.SinDeg[theta] ELSE RealFns.Sin[theta];
cos: REALIF degrees THEN RealFns.CosDeg[theta] ELSE RealFns.Cos[theta];
cos1, sqx, sqy, sqz, xycos1, zxcos1, yzcos1, xsin, ysin, zsin: REAL;
axis ← Vector3d.Normalize[axis];
[[sqx, sqy, sqz]] ← Vector3d.MulC[axis, axis];
cos1 ← 1.0-cos;
xycos1 ← axis.x*axis.y*cos1; xsin ← axis.x*sin;
yzcos1 ← axis.y*axis.z*cos1; ysin ← axis.y*sin;
zxcos1 ← axis.x*axis.z*cos1; zsin ← axis.z*sin;
IF out = NIL THEN out NEW[MatrixRep];
out^ ← [
[sqx+(1.0-sqx)*cos, xycos1-zsin,   zxcos1+ysin,   0.0],
[xycos1+zsin,  sqy+(1.0-sqy)*cos, yzcos1-xsin,   0.0],
[zxcos1-ysin,  yzcos1+xsin,   sqz+(1.0-sqz)*cos, 0.0],
[0.0,     0.0,     0.0,     1.0]];
RETURN[out];
};
MakeRotate: PUBLIC PROC [axis: Triple, theta: REAL, degrees: BOOLTRUE, base: Triple ← origin, out: Matrix ← NIL] RETURNS [Matrix] ~ {
rot: Matrix ← ObtainMatrix[];
tran: Matrix ← ObtainMatrix[];
mul: Matrix ← ObtainMatrix[];
out ← Translate[
   Mul[
   MakeTranslate[[-base.x, -base.y, -base.z], tran],
   MakePureRotate[axis, theta, degrees, rot], mul],
   base, out];
ReleaseMatrix[rot];
ReleaseMatrix[tran];
ReleaseMatrix[mul];
RETURN[out];
};
MakePerspective: PUBLIC PROC [dInv, fInv, fov: REAL, out: Matrix ← NIL] RETURNS [Matrix]
~ {
tan: REAL ← RealFns.TanDeg[0.5*fov];
out ← Identity[out];
out[2][2] ← tan*(1.0+fInv/MAX[0.0001, ABS[dInv]]);
out[2][3] ← tan;
RETURN[out];
};
HasPerspective: PUBLIC PROC [mat: Matrix] RETURNS [BOOL] ~ {
IF mat = NIL THEN RETURN[FALSE];
RETURN[mat[2][3] # 0.0];
};
Transformation Operations
TransformH: PUBLIC PROC [p: Triple, mat: Matrix] RETURNS [Quad] ~ {
RETURN[[
p.x*mat[0][0]+p.y*mat[1][0]+p.z*mat[2][0]+mat[3][0],
p.x*mat[0][1]+p.y*mat[1][1]+p.z*mat[2][1]+mat[3][1],
p.x*mat[0][2]+p.y*mat[1][2]+p.z*mat[2][2]+mat[3][2],
p.x*mat[0][3]+p.y*mat[1][3]+p.z*mat[2][3]+mat[3][3]]];
};
The following four transformation procedures test the transformation matrix for perspective (i.e.,
the matrix has perspective iff the last column is not [0, 0, 0, 1]. If a perspective element exists,
the return coordinates are first divided by w.
Transform: PUBLIC PROC [p: Triple, mat: Matrix] RETURNS [Triple] ~ {
q: Quad;
IF NOT HasPerspective[mat] THEN RETURN[[
p.x*mat[0][0]+p.y*mat[1][0]+p.z*mat[2][0]+mat[3][0],
p.x*mat[0][1]+p.y*mat[1][1]+p.z*mat[2][1]+mat[3][1],
p.x*mat[0][2]+p.y*mat[1][2]+p.z*mat[2][2]+mat[3][2]]];
q ← TransformH[p, mat];
RETURN[IF q.w = 1.0 THEN [q.x, q.y, q.z] ELSE [q.x/q.w, q.y/q.w, q.z/q.w]];
};
TransformD: PUBLIC PROC [p: Triple, mat: Matrix] RETURNS [Pair] ~ {
w: REAL;
pp: Pair ← [
p.x*mat[0][0]+p.y*mat[1][0]+p.z*mat[2][0]+mat[3][0],
p.x*mat[0][1]+p.y*mat[1][1]+p.z*mat[2][1]+mat[3][1]];
IF NOT HasPerspective[mat] THEN RETURN[pp];
w ← p.x*mat[0][3]+p.y*mat[1][3]+p.z*mat[2][3]+mat[3][3];
RETURN[IF w = 1.0 THEN pp ELSE [pp.x/w, pp.y/w]];
};
TransformPair: PUBLIC PROC [p: Pair, mat: Matrix] RETURNS [Triple] ~ {
w: REAL;
t: Triple ← [
p.x*mat[0][0]+p.y*mat[1][0]+mat[3][0],
p.x*mat[0][1]+p.y*mat[1][1]+mat[3][1],
p.x*mat[0][2]+p.y*mat[1][2]+mat[3][2]];
IF NOT HasPerspective[mat] THEN RETURN[t];
w ← p.x*mat[0][3]+p.y*mat[1][3]+mat[3][3];
RETURN[IF w = 1.0 THEN t ELSE [t.x/w, t.y/w, t.z/w]];
};
TransformPairD: PUBLIC PROC [p: Pair, mat: Matrix] RETURNS [Pair] ~ {
w: REAL;
pp: Pair ← [p.x*mat[0][0]+p.y*mat[1][0]+mat[3][0], p.x*mat[0][1]+p.y*mat[1][1]+mat[3][1]];
IF NOT HasPerspective[mat] THEN RETURN[pp];
w ← p.x*mat[0][3]+p.y*mat[1][3]+mat[3][3];
RETURN[IF w = 1.0 THEN pp ELSE [pp.x/w, pp.y/w]];
};
If no differential scaling exists within the transformation matrix, TransformVec may be used.
Otherwise, use TransformVecDiffS.
TransformVec: PUBLIC PROC [vec: Triple, mat: Matrix] RETURNS [Triple] ~ {
RETURN[[
vec.x*mat[0][0]+vec.y*mat[1][0]+vec.z*mat[2][0],
vec.x*mat[0][1]+vec.y*mat[1][1]+vec.z*mat[2][1],
vec.x*mat[0][2]+vec.y*mat[1][2]+vec.z*mat[2][2]]];
};
TransformVecDiffS: PUBLIC PROC [vec: Triple, mat: Matrix] RETURNS [Triple] ~ {
RETURN[Transform[vec, Cofactors[mat]]]; -- or check if sum of squares of rows are all equal.
};
Concatenation Operations
The next 3 procs post-multiply the matrix, thus transforming in the reference coordinate space.
Scale: PUBLIC PROC [in: Matrix, s: REAL, out: Matrix ← NIL] RETURNS [Matrix] ~ {
scratch: Matrix ← ObtainMatrix[];
out ← Mul[in, MakeScale[[s, s, s], scratch], out];
ReleaseMatrix[scratch];
RETURN[out];
};
DiffScale: PUBLIC PROC [in: Matrix, s: Triple, out: Matrix ← NIL] RETURNS [Matrix] ~ {
scratch: Matrix ← ObtainMatrix[];
out ← Mul[in, MakeScale[s, scratch], out];
ReleaseMatrix[scratch];
RETURN[out];
};
Rotate: PUBLIC PROC [in: Matrix, axis: Triple, theta: REAL, degrees: BOOLTRUE, base: Triple ← origin, out: Matrix ← NIL] RETURNS [Matrix] ~ {
scratch: Matrix ← ObtainMatrix[];
out ← Mul[in, MakeRotate[axis, theta, degrees, base, scratch], out];
ReleaseMatrix[scratch];
RETURN[out];
};
Translate: PUBLIC PROC [in: Matrix, p: Triple, out: Matrix ← NIL] RETURNS [Matrix] ~ {
scratch: Matrix ← ObtainMatrix[];
out ← Mul[in, MakeTranslate[p, scratch], out];
ReleaseMatrix[scratch];
RETURN[out];
};
The next 3 procs pre-multiply the matrix, thus transforming in the local coordinate space.
LocalScale: PUBLIC PROC [in: Matrix, s: REAL, out: Matrix ← NIL] RETURNS [Matrix] ~ {
scratch: Matrix ← ObtainMatrix[];
out ← Mul[MakeScale[[s, s, s], scratch], in, out];
ReleaseMatrix[scratch];
RETURN[out];
};
LocalDiffScale: PUBLIC PROC [in: Matrix, s: Triple, out: Matrix ← NIL] RETURNS [Matrix] ~ {
scratch: Matrix ← ObtainMatrix[];
out ← Mul[MakeScale[s, scratch], in, out];
ReleaseMatrix[scratch];
RETURN[out];
};
LocalRotate: PUBLIC PROC [in: Matrix, axis: Triple, theta: REAL, degrees: BOOLTRUE, base: Triple ← origin, out: Matrix ← NIL] RETURNS [Matrix] ~ {
scratch: Matrix ← ObtainMatrix[];
out ← Mul[MakeRotate[axis, theta, degrees, base, scratch], in, out];
ReleaseMatrix[scratch];
RETURN[out];
};
LocalTranslate: PUBLIC PROC [in: Matrix, p: Triple, out: Matrix ← NIL] RETURNS [Matrix] ~ {
scratch: Matrix ← ObtainMatrix[];
out ← Mul[MakeTranslate[p, scratch], in];
ReleaseMatrix[scratch];
RETURN[out];
};
Mathematical Operations
ScalarMul: PROC [in: Matrix, s: REAL, out: Matrix ← NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out NEW[MatrixRep];
FOR i: Row IN Row DO
FOR j: Col IN Col DO
out^[i][j] ← s*in[i][j];
ENDLOOP;
ENDLOOP;
RETURN[out];
};
Mul: PUBLIC PROC [left, rite: Matrix, out: Matrix ← NIL] RETURNS [Matrix] ~ {
ret: Matrix ← IF out # NIL THEN out ELSE NEW[MatrixRep];
FOR i: INT IN [0..3] DO
FOR j: INT IN [0..3] DO
ret[i][j] ← left[i][0]*rite[0][j]+left[i][1]*rite[1][j]+left[i][2]*rite[2][j]+left[i][3]*rite[3][j];
ENDLOOP;
ENDLOOP;
RETURN[ret];
};
Mul: PUBLIC PROC [left, rite: Matrix, out: Matrix ← NIL] RETURNS [Matrix] ~ {
ret: MatrixRep;
IF left = NIL OR rite = NIL THEN RETURN[out];
ret[0][0] ← left[0][0]*rite[0][0]+left[0][1]*rite[1][0]+left[0][2]*rite[2][0]+left[0][3]*rite[3][0];
ret[0][1] ← left[0][0]*rite[0][1]+left[0][1]*rite[1][1]+left[0][2]*rite[2][1]+left[0][3]*rite[3][1];
ret[0][2] ← left[0][0]*rite[0][2]+left[0][1]*rite[1][2]+left[0][2]*rite[2][2]+left[0][3]*rite[3][2];
ret[0][3] ← left[0][0]*rite[0][3]+left[0][1]*rite[1][3]+left[0][2]*rite[2][3]+left[0][3]*rite[3][3];
ret[1][0] ← left[1][0]*rite[0][0]+left[1][1]*rite[1][0]+left[1][2]*rite[2][0]+left[1][3]*rite[3][0];
ret[1][1] ← left[1][0]*rite[0][1]+left[1][1]*rite[1][1]+left[1][2]*rite[2][1]+left[1][3]*rite[3][1];
ret[1][2] ← left[1][0]*rite[0][2]+left[1][1]*rite[1][2]+left[1][2]*rite[2][2]+left[1][3]*rite[3][2];
ret[1][3] ← left[1][0]*rite[0][3]+left[1][1]*rite[1][3]+left[1][2]*rite[2][3]+left[1][3]*rite[3][3];
ret[2][0] ← left[2][0]*rite[0][0]+left[2][1]*rite[1][0]+left[2][2]*rite[2][0]+left[2][3]*rite[3][0];
ret[2][1] ← left[2][0]*rite[0][1]+left[2][1]*rite[1][1]+left[2][2]*rite[2][1]+left[2][3]*rite[3][1];
ret[2][2] ← left[2][0]*rite[0][2]+left[2][1]*rite[1][2]+left[2][2]*rite[2][2]+left[2][3]*rite[3][2];
ret[2][3] ← left[2][0]*rite[0][3]+left[2][1]*rite[1][3]+left[2][2]*rite[2][3]+left[2][3]*rite[3][3];
ret[3][0] ← left[3][0]*rite[0][0]+left[3][1]*rite[1][0]+left[3][2]*rite[2][0]+left[3][3]*rite[3][0];
ret[3][1] ← left[3][0]*rite[0][1]+left[3][1]*rite[1][1]+left[3][2]*rite[2][1]+left[3][3]*rite[3][1];
ret[3][2] ← left[3][0]*rite[0][2]+left[3][1]*rite[1][2]+left[3][2]*rite[2][2]+left[3][3]*rite[3][2];
ret[3][3] ← left[3][0]*rite[0][3]+left[3][1]*rite[1][3]+left[3][2]*rite[2][3]+left[3][3]*rite[3][3];
IF out = NIL THEN out NEW[MatrixRep];
out^ ← ret;
RETURN[out];
};
Cat: PUBLIC PROC [m1, m2: Matrix, m3, m4, m5, m6, out: Matrix ← NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out NEW[MatrixRep];
out ← Mul[m1, m2, out];
IF m3 # NIL THEN out ← Mul[out, m3, out];
IF m4 # NIL THEN out ← Mul[out, m4, out];
IF m5 # NIL THEN out ← Mul[out, m5, out];
IF m6 # NIL THEN out ← Mul[out, m6, out];
RETURN[out];
};
Invert: PUBLIC PROC [in: Matrix, out: Matrix ← NIL] RETURNS [Matrix] ~ {
det: REAL ← Determinant[in];
scratch: Matrix;
IF det = 0.0 THEN ERROR singular;
scratch ← ObtainMatrix[];
out ← ScalarMul[Adjoint[in, scratch], 1.0/det, out];
ReleaseMatrix[scratch];
RETURN[out];
};
Adjoint: PUBLIC PROC [in: Matrix, out: Matrix ← NIL] RETURNS [Matrix] ~ {
scratch: Matrix ← ObtainMatrix[];
out ← Transpose[Cofactors[in, scratch], out];
ReleaseMatrix[scratch];
RETURN[out];
or: transpose; chk for no perspec; check constant sum of squares for rows
};
Det3x3: PROC [a0, b0, c0, a1, b1, c1, a2, b2, c2: REAL] RETURNS [REAL] ~ {
RETURN[a0*(b1*c2-c1*b2)-b0*(a1*c2-c1*a2)+c0*(a1*b2-b1*a2)];
};
Cofactors: PUBLIC PROC [in: Matrix, out: Matrix ← NIL] RETURNS [Matrix] ~ {
m: Matrix ← in;
IF out = NIL THEN out ← NEW[MatrixRep];
out[0][0] ← Det3x3[m[1][1],m[1][2],m[1][3],m[2][1],m[2][2],m[2][3],m[3][1],m[3][2],m[3][3]];
out[0][1] ← -Det3x3[m[1][0],m[1][2],m[1][3],m[2][0],m[2][2],m[2][3],m[3][0],m[3][2],m[3][3]];
out[0][2] ← Det3x3[m[1][0],m[1][1],m[1][3],m[2][0],m[2][1],m[2][3],m[3][0],m[3][1],m[3][3]];
out[0][3] ← -Det3x3[m[1][0],m[1][1],m[1][2],m[2][0],m[2][1],m[2][2],m[3][0],m[3][1],m[3][2]];
out[1][0] ← -Det3x3[m[0][1],m[0][2],m[0][3],m[2][1],m[2][2],m[2][3],m[3][1],m[3][2],m[3][3]];
out[1][1] ← Det3x3[m[0][0],m[0][2],m[0][3],m[2][0],m[2][2],m[2][3],m[3][0],m[3][2],m[3][3]];
out[1][2] ← -Det3x3[m[0][0],m[0][1],m[0][3],m[2][0],m[2][1],m[2][3],m[3][0],m[3][1],m[3][3]];
out[1][3] ← Det3x3[m[0][0],m[0][1],m[0][2],m[2][0],m[2][1],m[2][2],m[3][0],m[3][1],m[3][2]];
out[2][0] ← Det3x3[m[0][1],m[0][2],m[0][3],m[1][1],m[1][2],m[1][3],m[3][1],m[3][2],m[3][3]];
out[2][1] ← -Det3x3[m[0][0],m[0][2],m[0][3],m[1][0],m[1][2],m[1][3],m[3][0],m[3][2],m[3][3]];
out[2][2] ← Det3x3[m[0][0],m[0][1],m[0][3],m[1][0],m[1][1],m[1][3],m[3][0],m[3][1],m[3][3]];
out[2][3] ← -Det3x3[m[0][0],m[0][1],m[0][2],m[1][0],m[1][1],m[1][2],m[3][0],m[3][1],m[3][2]];
out[3][0] ← -Det3x3[m[0][1],m[0][2],m[0][3],m[1][1],m[1][2],m[1][3],m[2][1],m[2][2],m[2][3]];
out[3][1] ← Det3x3[m[0][0],m[0][2],m[0][3],m[1][0],m[1][2],m[1][3],m[2][0],m[2][2],m[2][3]];
out[3][2] ← -Det3x3[m[0][0],m[0][1],m[0][3],m[1][0],m[1][1],m[1][3],m[2][0],m[2][1],m[2][3]];
out[3][3] ← Det3x3[m[0][0],m[0][1],m[0][2],m[1][0],m[1][1],m[1][2],m[2][0],m[2][1],m[2][2]];
RETURN[out];
};
Determinant: PUBLIC PROC [m: Matrix] RETURNS [REAL] ~ {
RETURN[
m[0][0]*Det3x3[m[1][1],m[1][2],m[1][3],m[2][1],m[2][2],m[2][3],m[3][1],m[3][2],m[3][3]] -
m[0][1]*Det3x3[m[1][0],m[1][2],m[1][3],m[2][0],m[2][2],m[2][3],m[3][0],m[3][2],m[3][3]] +
m[0][2]*Det3x3[m[1][0],m[1][1],m[1][3],m[2][0],m[2][1],m[2][3],m[3][0],m[3][1],m[3][3]] -
m[0][3]*Det3x3[m[1][0],m[1][1],m[1][2],m[2][0],m[2][1],m[2][2],m[3][0],m[3][1],m[3][2]]];
};
Transpose: PUBLIC PROC [in: Matrix, out: Matrix ← NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out NEW[MatrixRep];
FOR i: Row IN Row DO
FOR j: Col IN Col DO
out[i][j] ← in[j][i];
ENDLOOP;
ENDLOOP;
RETURN[out];
};
InTermsOf: PUBLIC PROC [A, B: Matrix, out: Matrix ← NIL] RETURNS [Matrix] ~ {
scratch: Matrix ← ObtainMatrix[];
out ← Invert[A, scratch];
out ← Mul[scratch, B, out];
ReleaseMatrix[scratch];
RETURN[out];
};
Invert3x3: PUBLIC PROC [m: Matrix, out: Matrix ← NIL] RETURNS [Matrix] ~ {
d: REAL ← Det3x3[m[0][0],m[0][1],m[0][2],m[1][0],m[1][1],m[1][2],m[2][0],m[2][1],m[2][2]];
IF out = NIL THEN out ← NEW[MatrixRep];
out[0][0] ← m[1][1]*m[2][2]-m[1][2]*m[2][1];
out[0][1] ← m[0][2]*m[2][1]-m[0][1]*m[2][2];
out[0][2] ← m[0][1]*m[1][2]-m[0][2]*m[1][1];
out[1][0] ← m[1][2]*m[2][0]-m[1][0]*m[2][2];
out[1][1] ← m[0][0]*m[2][2]-m[0][2]*m[2][0];
out[1][2] ← m[0][2]*m[1][0]-m[0][0]*m[1][2];
out[2][0] ← m[1][0]*m[2][1]-m[1][1]*m[2][0];
out[2][1] ← m[0][1]*m[2][0]-m[0][0]*m[2][1];
out[2][2] ← m[0][0]*m[1][1]-m[0][1]*m[1][0];
IF d # 0.0 THEN {
di: REAL ← 1.0/d;
FOR i: NAT IN[0..2] DO FOR j: NAT IN[0..2] DO out[i][j] ← di*out[i][j]; ENDLOOP; ENDLOOP;
};
RETURN[out];
};
Miscellaneous Operations
nScratchMatrices: NAT ~ 6;
scratchMatrices: ARRAY [0..nScratchMatrices) OF Matrix ← ALL[NIL];
ObtainMatrix: PUBLIC ENTRY PROC RETURNS [Matrix] ~ {
FOR i: NAT IN [0..nScratchMatrices) DO
matrix: Matrix ~ scratchMatrices[i];
IF matrix # NIL THEN {
scratchMatrices[i] ← NIL;
RETURN [matrix];
};
ENDLOOP;
RETURN [NEW[MatrixRep]];
};
ReleaseMatrix: PUBLIC ENTRY PROC [matrix: Matrix] ~ {
FOR i: NAT IN [0..nScratchMatrices) DO
IF scratchMatrices[i] = NIL THEN {
scratchMatrices[i] ← matrix;
RETURN;
};
ENDLOOP;
};
END.
..
References:
by Crow:  /Cyan/Imaging/6.0/ThreeDWorld/Linear3DImpl.mesa
by Bier:  /Cedar/CedarChest6.0/SolidModeler/Matrix3dImpl.mesa
by Wyatt: /Cedar/Cedar6.0/Imager/ImagerTransformationImpl.mesa