Creation Operations
CopyMatrix:
PUBLIC
PROC [in: Matrix, out: Matrix ←
NIL]
RETURNS [Matrix] ~ {
IF in = NIL THEN RETURN[NIL];
IF out = NIL THEN out ← NEW[MatrixRep];
out^ ← in^;
RETURN[out];
};
CopyMatrixSequence:
PUBLIC
PROC [in: MatrixSequence]
RETURNS [MatrixSequence] ~ {
copy: MatrixSequence ← NIL;
IF in #
NIL
THEN {
copy ← NEW[MatrixSequenceRep[in.maxLength]];
FOR n: NAT IN [0..in.maxLength) DO copy[n] ← CopyMatrix[in[n]]; ENDLOOP;
};
RETURN[copy];
};
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:
BOOL ←
TRUE, out: Matrix ←
NIL]
RETURNS [Matrix] ~ {
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 ← Vector3d.Normalize[axis];
[[sqx, sqy, sqz]] ← Vector3d.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];
};
MakeRotate:
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],
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:
BOOL ←
TRUE, 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:
BOOL ←
TRUE, 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];
};