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; 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]]; }; UnitSquareToQuadrilateral: PUBLIC PROC [p1, p2, p3, p4: Pair] RETURNS [Matrix] ~ { 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"] g Det2[px, dx2, py, 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.  G2dMatrixImpl.mesa Copyright 1984, 1992 by Xerox Corporation. All rights reserved. Bloomenthal, July 2, 1992 4:55 pm PDT Basic 3 by 3 Matrix Operations Transformation Miscellany 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. Compute matrix: |a d g| |b e h| |c f 1| g _ Det2[px, dx2, dy1, dy2]/del; "cedarcode" styleNewlineDelimiter J e6B%JJk &/Jbl JJ JJJJ J!J #JJ headln  +J6J