File: Matrix3d.mesa
Last edited by Bier on January 8, 1985 12:49:48 pm PST
Author: Eric Allan Bier in June 1982
Contents: Definitions module for the 4 by 4 homogeneous 3d transform package. This package assumes a matrix of the form:
[nx ox ax px]  where n = [nx ny nz]T is the "normal" column vector,
[ny oy ay py]  o = [ox oy oz]T is the "orthogonal" column vector
[nz oz az pz]  a = [ax ay az]T is the "axial" column vector,
[0 0 0 1]   and p = [px, py, pz]T is the "position" column
see Richard Paul's "Robot Manipulators" Chapter 1 for more information
DIRECTORY
SV2d,
SV3d,
SVVector3d;
Matrix3d: DEFINITIONS =
BEGIN
Matrix4by4: TYPE = SV3d.Matrix4by4;
Plane: TYPE = SV3d.Plane;
Point3d: TYPE = SV3d.Point3d;
Point2d: TYPE = SV2d.Point2d;
Vector: TYPE = SVVector3d.Vector;
Identity: PROC [] RETURNS [identityMat: Matrix4by4];
Update: PROC [mat: Matrix4by4, point: Point3d] RETURNS [newPoint: Point3d];
UpdateVectorWithInverse: PROC [inverse: Matrix4by4, vec: Vector] RETURNS [newVec: Vector];
Because of the wonders of mathematics, if a shape is transformed by a matrix transform A, then its vector should be transformed by A transpose inverse. If A involves only even scaling, rotation and translation then using A (minus translational components) will produce a vector with the proper direction (though not scaled as it would be with A transpose inverse). However, if A has differential scaling, then A transpose inverse must be used. If you already have A inverse handy, then don't bother to recompute it. Use UpdateVectorWithInverse. Otherwise use UpdateVectorComputeInverse. If there is no differential scaling in A and the magnitude of the result is not important, use UpdateVectorEvenScaling.
UpdateVectorComputeInverse: PROC [mat: Matrix4by4, vec: Vector] RETURNS [newVec: Vector];
UpdateVectorEvenScaling: PROC [mat: Matrix4by4, vec: Vector] RETURNS [newVec: Vector];
UpdatePlaneWithInverse: PROC [inverse: Matrix4by4, plane: Plane] RETURNS [newPlane: Plane];
UpdatePlaneComputeInverse: PROC [mat: Matrix4by4, plane: Plane] RETURNS [newPlane: Plane];
PerspectiveTrans: PROC [oldPoint: Point3d, focalDistance: REAL ← 1800]
RETURNS [newPoint: Point2d];
IllegalDepth: SIGNAL;
The next five procedures perform an efficient matrix multiplication as though mat were being multiplied to its left (a transform in reference coordinates) by a matrix which translates or rotates respectively.
Scale: PROC [mat: Matrix4by4, sx, sy, sz: REAL] RETURNS [transMat: Matrix4by4];
Translate: PROC [mat: Matrix4by4, dx, dy, dz: REAL] RETURNS [transMat: Matrix4by4];
RotateAboutXAxis: PROC [mat: Matrix4by4, degrees: REAL] RETURNS [transMat: Matrix4by4];
RotateAboutYAxis: PROC [mat: Matrix4by4, degrees: REAL] RETURNS [transMat: Matrix4by4];
RotateAboutZAxis: PROC [mat: Matrix4by4, degrees: REAL] RETURNS [transMat: Matrix4by4];
The next five procedures perform an efficient matrix multiplication as though mat were being multiplied on its right (a transform in local coordinates) by a matrix which translates or rotates respectively.
LocalScale: PROC [mat: Matrix4by4, sx, sy, sz: REAL] RETURNS [transMat: Matrix4by4];
LocalTranslate: PROC [mat: Matrix4by4, dx, dy, dz: REAL] RETURNS [transMat: Matrix4by4];
LocalRotateAboutXAxis: PROC [mat: Matrix4by4, degrees: REAL] RETURNS [transMat: Matrix4by4];
LocalRotateAboutYAxis: PROC [mat: Matrix4by4, degrees: REAL] RETURNS [transMat: Matrix4by4];
LocalRotateAboutZAxis: PROC [mat: Matrix4by4, degrees: REAL] RETURNS [transMat: Matrix4by4];

The next six procedures allow the user direct access to the transform matrices and their manipulation.
MatMult: PROC [left, right: Matrix4by4] RETURNS [transMat: Matrix4by4];
MatMult is almost a general 4 by 4 matrix operation. However, it assumes that the last row of each matrix is [0 0 0 1].
MakeScaleMat: PROC [sx, sy, sz: REAL] RETURNS [scale: Matrix4by4];
MakeTranslateMat: PROC [dx, dy, dz: REAL] RETURNS [trans: Matrix4by4];
MakeRotateXMat: PROC [degrees: REAL] RETURNS [rotx: Matrix4by4];
MakeRotateYMat: PROC [degrees: REAL] RETURNS [roty: Matrix4by4];
MakeRotateZMat: PROC [degrees: REAL] RETURNS [rotz: Matrix4by4];
MakeMatFromZandXAxis: PROC [zAxis, xAxis: Vector, origin: Point3d] RETURNS [mat: Matrix4by4];
MakeMatFromYandXAxis: PROC [yAxis, xAxis: Vector, origin: Point3d] RETURNS [mat: Matrix4by4];
Uses the normalized version of the zAxis given as it is. Projects xAxis onto the plane of zAxis and normalizes. Finds y by cross product.
MakeHorizontalMatFromZAxis: PROC [zAxis: Vector, origin: Point3d] RETURNS [mat: Matrix4by4];
Uses zAxis as it is. Finds a horizontal x axis orthogonal to zAxis. If zAxis is vertical, then there are infinitely many. Chooses [1,0,0] in this case.
MakeMatFromAxes: PROC [xAxis, yAxis, zAxis: Vector, origin: Point3d] RETURNS [mat: Matrix4by4];
Uses the axes just as they are given. Doesn't even check to see if they are orthogonal.
ScaleFromMatrix: PROC [mat: Matrix4by4] RETURNS [sx, sy, sz: REAL];
Finds out how much this matrix will scale an object it is applied to.
OriginOfMatrix: PROC [mat: Matrix4by4] RETURNS [origin: Point3d];
XAxisOfMatrix: PROC [mat: Matrix4by4] RETURNS [xAxis: Vector];
YAxisOfMatrix: PROC [mat: Matrix4by4] RETURNS [yAxis: Vector];
ZAxisOfMatrix: PROC [mat: Matrix4by4] RETURNS [zAxis: Vector];

Because of the orthogonality of the transform matrices, matrix inversion is relatively easy.
InverseNoScaling: PROC [mat: Matrix4by4] RETURNS [inverse: Matrix4by4];
DegenerateInverse: SIGNAL;
WARNING: This procedure inverts the translation and rotation components of a matrix, but does not invert the scaling component. Hence, if you take a matrix A which scales by 10, translates, and rotates. Applying the matrix and its inverse to an object will leave in its original position, scaled by 100! Of course, if the scaling factor is 1, this doesn't matter, and this procedure is very fast.
Cofactor: PROC [mat: Matrix4by4, row, col: NAT] RETURNS [cof: REAL];
Returns the determinant of the 3 by 3 matrix which results by crossing off row row and column col in mat.
Determinant: PROC [mat: Matrix4by4] RETURNS [det: REAL];
Transpose: PROC [mat: Matrix4by4] RETURNS [MatT: Matrix4by4];
Inverse: PROC [mat: Matrix4by4] RETURNS [inverse: Matrix4by4];
Inverts any matrix of the given form (even with scaling)

Once you can create and manipulate objects in the world (reference) coordinate system, it is useful to be able to express one coordinate system in terms of another. Consider a solar system example. If you have the center of the sun as world coordinate system, and earth and moon, as two other coordinate system expressed in terms of world coordinates (recall that a 4 by 4 matrix represents the relationship between coordinate systems), then if you want the moon to rotate around the earth, you might try:

MoonInTermsOfEarth ← WorldToLocal[earth,moon];
FOREVER DO:
 MoonInTermsOfEarth ← RotateAboutY[MoonInTermsOfEarth, degrees];
 moon ← MatMult[earth, MoonInTermsOfEarth];
 DrawSphereCenteredAt[moon];
ENDLOOP;
WorldToLocal: PROC [AinWorld,BinWorld: Matrix4by4] RETURNS [BinTermsOfA: Matrix4by4];
Normalizing
UnitNormalize: PROC [mat: Matrix4by4] RETURNS [unitMat: Matrix4by4];
Finds a matrix transform which provides the same translation and rotation of mat, but with a scaling of 1. While we're at it, make sure the axes of unitMat are perpendicular.
Normalize: PROC [mat: Matrix4by4] RETURNS [normMat: Matrix4by4];
Finds a matrix with the same translation, rotation and scaling as mat, but with the axes guaranteed to be orthogonal.
NormalizeUniformScale: PROC [mat: Matrix4by4] RETURNS [normMat: Matrix4by4];
Finds a matrix with the same translation, and rotation as mat, with the axes guaranteed to be orthogonal and with the scaling along all three axes set equal to the maximum scaling of the three axes of mat.
NormalizeNoRotation: PROC [mat: Matrix4by4] RETURNS [normMat: Matrix4by4];
Finds a matrix with the same translation and scaling as mat, with the axes guaranteed to be orthogonal, removing rotational components.
NormalizeNoTranslation: PROC [mat: Matrix4by4] RETURNS [normMat: Matrix4by4];
Finds a matrix with the same scaling and rotation as mat, with the axes guaranteed to be orthogonal, removing translational components.
TidyUpVector: PROC [vec: Vector] RETURNS [fixedVec: Vector];
If a component is almost a multiple of .5, make it exactly so.
TidyUpMatrix: PROC [mat: Matrix4by4] RETURNS [fixedMat: Matrix4by4];
If a component is almost a multiple of .5, make it exactly so.

END.