<<>> <> <> <> <> <> <> <<>> <> <> <> <> <<>> DIRECTORY G2dMatrix, G3dBasic, G3dMatrix, G3dPlane, G3dVector, Real, RealFns; G3dMatrixImpl: CEDAR MONITOR IMPORTS G2dMatrix, G3dPlane, G3dVector, Real, RealFns EXPORTS G3dMatrix ~ BEGIN singular: PUBLIC ERROR ~ CODE; <> 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; <> 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] ~ { <> 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]; }; }; <> <> <<~ {>> <> <> <> <> <> <> <> <<};>> <<>> HasPerspective: PUBLIC PROC [mat: Matrix] RETURNS [BOOL] ~ { IF mat = NIL THEN RETURN[FALSE]; <> 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]; }; <> 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]]]; }; <> <> <> <<>> 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]]]; }; <> <> 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 }; <> <> 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]; }; <<>> <> <<>> 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]; }; <> 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] ~ { <> <> <> <> <> 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]; <> }; 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]}; <> <<>> Asin: PROC [x: REAL] RETURNS [a: REAL] ~ INLINE {a ¬ Atan[x, RealFns.SqRt[1.0-x*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]; }; <> 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; }; <> 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; }; <> 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 <> 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]]; }; <> 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. .. <> 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 <> <> 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]; };