G3dMatrixImpl.mesa
Copyright Ó 1984, 1992 by Xerox Corporation. All rights reserved.
Bier, June 1982
Crow, October 31, 1985 12:31:44 pm PST
Bloomenthal, October 21, 1992 4:55 pm PDT
Ken Shoemake, August 30, 1989 8:17:57 pm PDT
See ``Principles of Interactive Computer Graphics,'' by Newman and Sproull.
Left-handed eyespace coordinate system (right, up, depth = into screen);
right-handed world coordinate system (longitude, latitude, altitude).
Ken Fishkin, November 6, 1991 1:30 pm PST
DIRECTORY G2dMatrix, G3dBasic, G3dMatrix, G3dPlane, G3dVector, Real, RealFns;
G3dMatrixImpl: CEDAR MONITOR
IMPORTS G2dMatrix, G3dPlane, G3dVector, Real, RealFns
EXPORTS G3dMatrix
~ BEGIN
singular: PUBLIC ERROR ~ CODE;
Type Declarations
Matrix:    TYPE ~ G3dMatrix.Matrix;
MatrixRep:   TYPE ~ G3dMatrix.MatrixRep;
MatrixSequence:  TYPE ~ G3dMatrix.MatrixSequence;
MatrixSequenceRep: TYPE ~ G3dMatrix.MatrixSequenceRep;
Viewport:    TYPE ~ G3dMatrix.Viewport;
Pair:     TYPE ~ G3dBasic.Pair;
Triple:    TYPE ~ G3dBasic.Triple;
Quad:     TYPE ~ G3dBasic.Quad;
Ray:     TYPE ~ G3dBasic.Ray;
Screen:    TYPE ~ G3dBasic.Screen;
origin:    Triple ~ G3dBasic.origin;
Creation Operations
CopyMatrix: PUBLIC PROC [matrix: Matrix, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
IF matrix = NIL THEN RETURN[NIL];
IF out = NIL THEN out ¬ NEW[MatrixRep];
out­ ¬ matrix­;
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,
unitize: BOOL ¬ FALSE,
out: Matrix ¬ NIL]
RETURNS [Matrix]
~ {
IF out = NIL THEN out ¬ NEW[MatrixRep];
IF unitize THEN {
v1 ¬ G3dVector.Unit[v1];
v2 ¬ G3dVector.Unit[v2];
v3 ¬ G3dVector.Unit[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];
};
MakeFromRows: PUBLIC PROC [r0, r1, r2, r3: Quad, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out ¬ NEW[MatrixRep];
out[0] ¬ [r0.x, r0.y, r0.z, r0.w];
out[1] ¬ [r1.x, r1.y, r1.z, r1.w];
out[2] ¬ [r2.x, r2.y, r2.z, r2.w];
out[3] ¬ [r3.x, r3.y, r3.z, r3.w];
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];
};
MakeRotate: PUBLIC PROC [
axis: Triple,
theta: REAL,
degrees: BOOL ¬ TRUE,
out: Matrix ¬ NIL]
RETURNS [Matrix]
~ {
This performs correctly for right-handed coordinate system and right-handed rotations.
sin: REAL ¬ IF degrees THEN RealFns.SinDeg[theta] ELSE RealFns.Sin[theta];
cos: REAL ¬ IF degrees THEN RealFns.CosDeg[theta] ELSE RealFns.Cos[theta];
cos1, sqx, sqy, sqz, xycos1, zxcos1, yzcos1, xsin, ysin, zsin: REAL;
axis ¬ G3dVector.Unit[axis];
[[sqx, sqy, sqz]] ¬ G3dVector.MulVectors[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];
};
MakeRotateAbout: PUBLIC PROC [
axis: Triple,
theta: REAL,
degrees: BOOL ¬ TRUE,
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],
   MakeRotate[axis, theta, degrees, rot], mul],
   base, out];
ReleaseMatrix[rot];
ReleaseMatrix[tran];
ReleaseMatrix[mul];
RETURN[out];
};
MakePerspective: PUBLIC PROC [nInv, fInv, fov: REAL, out: Matrix ¬ NIL]
RETURNS [Matrix]
~ {
IF fov = 0.0
THEN RETURN[Identity[out]]
ELSE {
d: REAL ¬ 1.0/nInv;
s: REAL ¬ d*RealFns.TanDeg[0.5*fov];
out ¬ Identity[out];
out[0][0] ¬ 1.0/s;
out[1][1] ¬ 1.0/s;
out[2][2] ¬ 1.0/(d*(1.0-d*fInv));
out[2][3] ¬ nInv;
out[3][2] ¬ -1.0/(1.0-d*fInv);
out[3][3] ¬ 0.0;
RETURN[out];
};
};
MakePerspective: PUBLIC PROC [nInv, fInv, fov: REAL, out: Matrix ← NIL]
RETURNS [Matrix]
~ {
tan: REAL ← RealFns.TanDeg[0.5*fov];
out ← Identity[out];
out[2][2] ← tan*(nInv+fInv)/(nInv-fInv);
out[2][3] ← tan;
out[3][2] ← -2.*tan/(nInv-fInv);
out[3][3] ← 0.;
RETURN[out];
};
HasPerspective: PUBLIC PROC [mat: Matrix] RETURNS [BOOL] ~ {
IF mat = NIL THEN RETURN[FALSE];
RETURN[mat[2][3] # 0.0]; used to be
RETURN[mat[0][3] # 0. OR mat[1][3] # 0. OR mat[2][3] # 0. OR mat[3][3] # 1.0];
};
MakeVectorTransform: PUBLIC PROC [pointTransform: Matrix, out: Matrix ¬ NIL]
RETURNS [n: Matrix]
~ {
IF out = NIL THEN out ¬ NEW[MatrixRep];
out ¬ Transpose[Adjoint[pointTransform], out];
FOR k: INT IN [0..4) DO out[3][k] ¬ out[k][3] ¬ 0.0; ENDLOOP;
out[3][3] ¬ 1.0;
RETURN[out];
};
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 = 0.0 THEN [0., 0.] ELSE 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]];
};
TransformQuad: PUBLIC PROC [q: Quad, mat: Matrix] RETURNS [Quad] ~ {
RETURN[[
q.x*mat[0][0]+q.y*mat[1][0]+q.z*mat[2][0]+q.w*mat[3][0],
q.x*mat[0][1]+q.y*mat[1][1]+q.z*mat[2][1]+q.w*mat[3][1],
q.x*mat[0][2]+q.y*mat[1][2]+q.z*mat[2][2]+q.w*mat[3][2],
q.x*mat[0][3]+q.y*mat[1][3]+q.z*mat[2][3]+q.w*mat[3][3]]];
};
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.
};
TransformPlane: PUBLIC PROC [plane: Quad, mat: Matrix, matInverse: Matrix ¬ NIL]
RETURNS [Quad]
~ {
m: Matrix ¬ IF matInverse # NIL THEN matInverse ELSE Adjoint[mat];
RETURN[PremultiplyPlaneByMatrix[plane, m]];
};
PremultiplyPlaneByMatrix: PUBLIC PROC [plane: Quad, m: Matrix] RETURNS [Quad] ~ {
RETURN[[
plane.x*m[0][0]+plane.y*m[0][1]+plane.z*m[0][2]+plane.w*m[0][3],
plane.x*m[1][0]+plane.y*m[1][1]+plane.z*m[1][2]+plane.w*m[1][3],
plane.x*m[2][0]+plane.y*m[2][1]+plane.z*m[2][2]+plane.w*m[2][3],
plane.x*m[3][0]+plane.y*m[3][1]+plane.z*m[3][2]+plane.w*m[3][3]]];
};
TransformByViewport: PUBLIC PROC [p: Pair, vp: Viewport] RETURNS [x: Pair] ~ {
x ¬ [vp.scale.y*p.x+vp.translate.x, vp.scale.y*p.y+vp.translate.y]; -- uniform scale
};
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: BOOL ¬ TRUE,
base: Triple ¬ origin,
out: Matrix ¬ NIL]
RETURNS [Matrix]
~ {
scratch: Matrix ¬ ObtainMatrix[];
out ¬ Mul[in, MakeRotateAbout[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: BOOL ¬ TRUE,
base: Triple ¬ origin,
out: Matrix ¬ NIL]
RETURNS [Matrix]
~ {
scratch: Matrix ¬ ObtainMatrix[];
out ¬ Mul[MakeRotateAbout[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: [0..4) IN [0..4) DO
FOR j: [0..4) IN [0..4) DO
out­[i][j] ¬ s*in[i][j];
ENDLOOP;
ENDLOOP;
RETURN[out];
};
Mul: PUBLIC PROC [left, rite: Matrix, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
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;
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 [ret: Matrix] ~ {
ret ¬ IF out = NIL THEN NEW[MatrixRep] ELSE out;
IF in[3][0] = 0.0 AND in[3][1] = 0.0 AND in[3][2] = 0.0 AND in[3][3] = 1.0
THEN { -- a little trick from Martin Newell (or was it Jim Blinn?)
inv3x3: G2dMatrix.Matrix ¬ G2dMatrix.Invert[[
[in[0][0], in[0][1], in[0][2]],
[in[1][0], in[1][1], in[1][2]],
[in[2][0], in[2][1], in[2][2]]]];
q: Triple ¬ G2dMatrix.Transform[[-in[3][0], -in[3][1], -in[3][2]], inv3x3];
ret[0] ¬ [inv3x3.row1.x, inv3x3.row1.y, inv3x3.row1.z, 0.0];
ret[1] ¬ [inv3x3.row2.x, inv3x3.row2.y, inv3x3.row2.z, 0.0];
ret[2] ¬ [inv3x3.row3.x, inv3x3.row3.y, inv3x3.row3.z, 0.0];
ret[3] ¬ [q.x, q.y, q.z, 1.0];
}
ELSE {
adj: Matrix ¬ Adjoint[in, ObtainMatrix[]];
det: REAL ¬ in[0][0]*adj[0][0]+in[0][1]*adj[1][0]+in[0][2]*adj[2][0]+in[0][3]*adj[3][0];
IF det = 0.0 THEN ERROR singular;
ret ¬ ScalarMul[adj, 1.0/det, ret];
ReleaseMatrix[adj];
};
};
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 [m: G2dMatrix.Matrix] RETURNS [REAL] ~ G2dMatrix.Determinant;
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: [0..4) IN [0..4) DO
out[i][i] ¬ in[i][i];
FOR j: [0..4) IN [0..i) DO
ij: REAL ¬ in[i][j];
ji: REAL ¬ in[j][i];
out[i][j] ¬ ji;
out[j][i] ¬ ij;
ENDLOOP;
ENDLOOP;
RETURN[out];
};
Atan: PROC [y, x: REAL] RETURNS [a: REAL] ~ INLINE {a ¬ RealFns.ArcTanDeg[y, x]};
Return the ArcTangent of y and x.
Asin: PROC [x: REAL] RETURNS [a: REAL] ~ INLINE {a ¬ Atan[x, RealFns.SqRt[1.0-x*x]]};
Return the ArcSin of x.
ExtractRotate: PUBLIC PROC [m: Matrix] RETURNS [Triple] ~ {
RETURN[IF ABS[m[0][2]] # 1.0
THEN [Atan[m[1][2], m[2][2]], -Asin[m[0][2]], Atan[m[0][1], m[0][0]]]
ELSE [Atan[m[1][0], m[2][0]], -Asin[m[0][2]], 0.]];
};
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];
};
Sequence Operations
CopyMatrixSequence: PUBLIC PROC [matrices: MatrixSequence]
RETURNS [MatrixSequence] ~ {
copy: MatrixSequence ¬ NIL;
IF matrices # NIL THEN {
copy ¬ NEW[MatrixSequenceRep[matrices.length]];
copy.length ¬ matrices.length;
FOR n: NAT IN [0..matrices.length) DO copy[n] ¬ CopyMatrix[matrices[n]]; ENDLOOP;
};
RETURN[copy];
};
AddToMatrixSequence: PUBLIC PROC [matrices: MatrixSequence, matrix: Matrix]
RETURNS [MatrixSequence]
~ {
IF matrices = NIL THEN matrices ¬ NEW[MatrixSequenceRep[1]];
IF matrices.length = matrices.maxLength THEN matrices ¬ LengthenMatrixSequence[matrices];
matrices[matrices.length] ¬ matrix;
matrices.length ¬ matrices.length+1;
RETURN[matrices];
};
LengthenMatrixSequence: PUBLIC PROC [matrices: MatrixSequence, amount: REAL ¬ 1.3]
RETURNS [new: MatrixSequence] ~ {
newLength: NAT ¬ MAX[Real.Ceiling[amount*matrices.maxLength], 3];
new ¬ NEW[MatrixSequenceRep[newLength]];
FOR i: NAT IN [0..matrices.length) DO new[i] ¬ matrices[i]; ENDLOOP;
new.length ¬ matrices.length;
};
Scratch Matrix 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;
};
Point Transformation and Visibility Testing
GetScreen: PUBLIC PROC [
point: Triple,
view: Matrix,
hasPerspective: BOOL,
viewport: Viewport ¬ []]
RETURNS [s: Screen]
~ {
IF hasPerspective
THEN {
q: Quad ¬ s.quad ¬ TransformH[point, view];
q.x ¬ s.quad.x ¬ q.x*viewport.aspectRecip; -- deform to square NDC clip-space
s.pos ¬ SELECT q.w FROM 0, 1 => [q.x, q.y], ENDCASE => [q.x/q.w, q.y/q.w];
s.l ¬ q.w+q.x < 0.0;  -- test left clipping
s.r ¬ q.w-q.x < 0.0;  -- test right clipping
s.b ¬ q.w+q.y < 0.0;  -- test bottom clipping
s.t ¬ q.w-q.y < 0.0;  -- test top clipping
s.n ¬ q.z < 0.0;   -- test near clipping
s.f ← q.w-q.z < 0.0;  -- test far clipping
IF (s.visible ¬ NOT (s.l OR s.r OR s.b OR s.t OR s.n)) AND viewport # [] THEN {
s.pos.x ¬ viewport.scale.x*s.pos.x+viewport.translate.x; -- aspectRecip*scale.x=scale.y,
s.pos.y ¬ viewport.scale.y*s.pos.y+viewport.translate.y; -- so uniform scale
};
}
ELSE {
z: REAL ¬ point.x*view[0][2]+point.y*view[1][2]+point.z*view[2][2]+view[3][2];
s.pos ¬ [
point.x*view[0][0]+point.y*view[1][0]+point.z*view[2][0]+view[3][0],
point.x*view[0][1]+point.y*view[1][1]+point.z*view[2][1]+view[3][1]];
IF (s.visible ¬ NOT (s.n ¬ z < 0.0)) AND viewport # [] THEN {
s.pos.x ¬ viewport.scale.y*s.pos.x+viewport.translate.x; -- uniform scale
s.pos.y ¬ viewport.scale.y*s.pos.y+viewport.translate.y;
};
};
s.intPos ¬ [Real.Round[s.pos.x], Real.Round[s.pos.y]];
};
Three-dimensional Information Derived from the Screen
TripleFromScreen: PUBLIC PROC [p: Pair, view: Matrix, viewInverse: Matrix ¬ NIL]
RETURNS [t: Triple]
~ {
x: Matrix ¬ IF viewInverse # NIL THEN viewInverse ELSE Invert[view, ObtainMatrix[]];
t ¬ Transform[[p.x, p.y, 0.0], x]; -- z is irrelevant
IF viewInverse = NIL THEN ReleaseMatrix[x];
};
TripleFromScreenAndPlane: PUBLIC PROC [p: Pair, plane: Quad, view: Matrix]
RETURNS [t: Triple]
~ {
vScreen: Quad ¬ [1.0, 0.0, 0.0, -p.x];  -- vertical plane through the given screen point
hScreen: Quad ¬ [0.0, 1.0, 0.0, -p.y];  -- horizontal plane through the given screen point
The screen planes must be transformed into world space; that is, transformed by the inverse of view. To transform plane F by matrix M is to premultiply F by the inverse of M. Thus, if M = view; transform plane F by M-1, yielding F'. Thus, F' = MF.
vWorld: Quad ¬ PremultiplyPlaneByMatrix[vScreen, view];
hWorld: Quad ¬ PremultiplyPlaneByMatrix[hScreen, view];
t ¬ G3dPlane.IntersectionOf3[[vWorld.x, vWorld.y, vWorld.z, vWorld.w, FALSE], [hWorld.x, hWorld.y, hWorld.z, hWorld.w, FALSE], [plane.x, plane.y, plane.z, plane.w, FALSE]];
};
LineFromScreen: PUBLIC PROC [p: Pair, view: Matrix, viewInverse: Matrix ¬ NIL]
RETURNS [r: Ray]
~ {
x: Matrix ¬ IF viewInverse # NIL THEN viewInverse ELSE Invert[view, ObtainMatrix[]];
p0: Triple ¬ TripleFromScreen[p, view, x];
p1: Triple ¬ TripleFromScreen[[0.0, 0.0], view, x];
r ¬ [p0, G3dVector.Sub[p1, p0]];
IF viewInverse = NIL THEN ReleaseMatrix[x];
};
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
The Wayside
dInverse = 0 for orthographic projection:
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]]);
Basically, with an infinite far clipping plane, d is irrelevant.
out[2][3] ¬ tan;
RETURN[out];
};