G2dMatrixImpl.mesa
Copyright Ó 1984, 1992 by Xerox Corporation. All rights reserved.
Bloomenthal, July 2, 1992 4:55 pm PDT
DIRECTORY G2dBasic, G2dMatrix, Imager, RealFns;
G2dMatrixImpl: CEDAR MONITOR
IMPORTS RealFns
EXPORTS G2dMatrix
~ BEGIN
Pair:   TYPE ~ G2dBasic.Pair;
Triple:  TYPE ~ G2dBasic.Triple;
Matrix:  TYPE ~ G2dMatrix.Matrix;
Rectangle: TYPE ~ Imager.Rectangle;
singular:  PUBLIC ERROR ~ CODE;
Basic 3 by 3 Matrix Operations
Identity: PUBLIC PROC RETURNS [Matrix] ~ {
RETURN[[[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]];
};
Transpose: PUBLIC PROC [m: Matrix] RETURNS [t: Matrix] ~ {
t ¬ [[m.row1.x, m.row2.x, m.row3.x],
[m.row1.y, m.row2.y, m.row3.y],
[m.row1.z, m.row2.z, m.row3.z]];
};
Determinant: PUBLIC PROC [m: Matrix] RETURNS [d: REAL] ~ {
d ¬
m.row1.x*(m.row2.y*m.row3.z-m.row3.y*m.row2.z)-
m.row1.y*(m.row2.x*m.row3.z-m.row3.x*m.row2.z)+
m.row1.z*(m.row2.x*m.row3.y-m.row3.x*m.row2.y);
};
Adjoint: PUBLIC PROC [m: Matrix] RETURNS [r: Matrix] ~ {
r.row1.x ¬ m.row2.y*m.row3.z-m.row2.z*m.row3.y;
r.row2.x ¬ -(m.row2.x*m.row3.z-m.row2.z*m.row3.x);
r.row3.x ¬ m.row2.x*m.row3.y-m.row2.y*m.row3.x;
r.row1.y ¬ -(m.row1.y*m.row3.z-m.row1.z*m.row3.y);
r.row2.y ¬ m.row1.x*m.row3.z-m.row1.z*m.row3.x;
r.row3.y ¬ -(m.row1.x*m.row3.y-m.row1.y*m.row3.x);
r.row1.z ¬ m.row1.y*m.row2.z-m.row1.z*m.row2.y;
r.row2.z ¬ -(m.row1.x*m.row2.z-m.row1.z*m.row2.x);
r.row3.z ¬ m.row1.x*m.row2.y-m.row1.y*m.row2.x;
};
Invert: PUBLIC PROC [m: Matrix] RETURNS [r: Matrix] ~ {
det: REAL;
r ¬ Adjoint[m];
IF (det ¬ r.row1.x*m.row1.x+r.row1.y*m.row2.x+r.row1.z*m.row3.x) = 0.0
THEN ERROR singular
ELSE {
di: REAL ¬ 1.0/det;
r.row1.x ¬ di*r.row1.x;
r.row1.y ¬ di*r.row1.y;
r.row1.z ¬ di*r.row1.z;
r.row2.x ¬ di*r.row2.x;
r.row2.y ¬ di*r.row2.y;
r.row2.z ¬ di*r.row2.z;
r.row3.x ¬ di*r.row3.x;
r.row3.y ¬ di*r.row3.y;
r.row3.z ¬ di*r.row3.z;
};
};
Mul: PUBLIC PROC [left, rite: Matrix] RETURNS [ret: Matrix] ~ {
ret.row1.x ¬ left.row1.x*rite.row1.x+left.row1.y*rite.row2.x+left.row1.z*rite.row3.x;
ret.row1.y ¬ left.row1.x*rite.row1.y+left.row1.y*rite.row2.y+left.row1.z*rite.row3.y;
ret.row1.z ¬ left.row1.x*rite.row1.z+left.row1.y*rite.row2.z+left.row1.z*rite.row3.z;
ret.row2.x ¬ left.row2.x*rite.row1.x+left.row2.y*rite.row2.x+left.row2.z*rite.row3.x;
ret.row2.y ¬ left.row2.x*rite.row1.y+left.row2.y*rite.row2.y+left.row2.z*rite.row3.y;
ret.row2.z ¬ left.row2.x*rite.row1.z+left.row2.y*rite.row2.z+left.row2.z*rite.row3.z;
ret.row3.x ¬ left.row3.x*rite.row1.x+left.row3.y*rite.row2.x+left.row3.z*rite.row3.x;
ret.row3.y ¬ left.row3.x*rite.row1.y+left.row3.y*rite.row2.y+left.row3.z*rite.row3.y;
ret.row3.z ¬ left.row3.x*rite.row1.z+left.row3.y*rite.row2.z+left.row3.z*rite.row3.z;
};
Rotate: PUBLIC PROC [mat: Matrix, radians: REAL] RETURNS [Matrix] ~ {
cos: REAL ← RealFns.Cos[radians];
sin: REAL ← RealFns.Sin[radians];
RETURN[Mul[mat, [[cos, sin, 0.0], [-sin, cos, 0.0], [0.0, 0.0, 1.0]]]];
};
Scale: PUBLIC PROC [mat: Matrix, scale: Pair] RETURNS [Matrix] ~ {
RETURN[Mul[mat, [[scale.x, 0.0, 0.0], [0.0, scale.y, 0.0], [0.0, 0.0, 1.0]]]];
};
Translate: PUBLIC PROC [mat: Matrix, translate: Pair] RETURNS [Matrix] ~ {
RETURN[Mul[mat, [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [translate.x, translate.y, 1.0]]]];
};
Transform: PUBLIC PROC [p: Triple, mat: Matrix] RETURNS [Triple] ~ {
RETURN[[
p.x*mat.row1.x+p.y*mat.row2.x+p.z*mat.row3.x,
p.x*mat.row1.y+p.y*mat.row2.y+p.z*mat.row3.y,
p.x*mat.row1.z+p.y*mat.row2.z+p.z*mat.row3.z]];
};
Transformation Miscellany
UnitSquareToQuadrilateral: PUBLIC PROC [p1, p2, p3, p4: Pair] RETURNS [Matrix] ~ {
We wish to solve for N, the 3 by 3 perspective matrix that transforms a unit square,
{S0, S1, S2, S3}) to an arbitrary quadrilateral, {P0, P1, P2, P3}. Thus,
[Artwork node; type 'Artwork on' to command tool]
S is given by: {[0, 0, 1], [0, 1, 1], [1, 1, 1], [1, 0, 1]}; since the transform is homogeneous,
P is given by {[x1w, y1w, w], [x2w, y2w, w], [x3w, y3w, w], [x4w, y4w, w]}:
[Artwork node; type 'Artwork on' to command tool]
N is derived from the system of equations that result from the above.
Since N is homogeneous, it can be scaled arbitrarily; in particular, we can set N22 = 1.
For example, using P1, we obtain three equations:
(0, 0, 1) . (N00, N10, N20) = x1w, or N20 = x1w.
(0, 0, 1) . (N01, N11, N21) = y1w, or N21 = y1w.
(0, 0, 1) . (N02, N12, N22) = w, or N22 = w = 1.
These three simplify to:
N20 = x1 and N21 = y1.
The remaining six equations are obtained similarly. Cf. Heckbert, CG&A, Nov. 1986.
Det2: PROC [a, b, c, d: REAL] RETURNS [REAL] ~ INLINE {RETURN[a*d-b*c]};
x0: REAL ¬ p1.x;
y0: REAL ¬ p1.y;
x1: REAL ¬ p2.x;
y1: REAL ¬ p2.y;
x2: REAL ¬ p3.x;
y2: REAL ¬ p3.y;
x3: REAL ¬ p4.x;
y3: REAL ¬ p4.y;
px: REAL ¬ x0-x1+x2-x3;
py: REAL ¬ y0-y1+y2-y3;
IF RealFns.AlmostZero[px, -21] AND RealFns.AlmostZero[py, -21]
THEN RETURN[[[x2-x1, y2-y1, 0.0], [x1-x0, y1-y0, 0.0], [x0, y0, 1.0]]] -- affine
ELSE {                   -- perspective
dx1: REAL ¬ x3-x2;
dx2: REAL ¬ x1-x2;
dy1: REAL ¬ y3-y2;
dy2: REAL ¬ y1-y2;
del: REAL ¬ Det2[dx1, dx2, dy1, dy2];
a, b, c, d, e, f, g, h: REAL;
IF del = 0.0 THEN ERROR; -- should be Error["Degenerate quadrilateral"]
Compute matrix: |a d g|
     |b e h|
     |c f 1|
g ¬ Det2[px, dx2, py, dy2]/del;
g ← Det2[px, dx2, dy1, dy2]/del;
h ¬ Det2[dx1, px, dy1, py]/del;
a ¬ x3-x0+g*x3;
b ¬ x1-x0+h*x1;
c ¬ x0;
d ¬ y3-y0+g*y3;
e ¬ y1-y0+h*y1;
f ¬ y0;
RETURN[[[a, d, g], [b, e, h], [c, f, 1.0]]];
};
};
QuadrilateralToRectangle: PUBLIC PROC [p1, p2, p3, p4: Pair, rectangle: Rectangle]
RETURNS [transform: Matrix]
~ {
transform ¬ UnitSquareToQuadrilateral[p1, p2, p3, p4];
transform ¬ Invert[transform];
transform ¬ Scale[transform, [rectangle.w, rectangle.h]];
transform ¬ Translate[transform, [rectangle.x, rectangle.y]];
};
END.