DIRECTORY JunoStorage USING [Point, Coords], JunoMatrix; JunoMatrixImpl: PROGRAM EXPORTS JunoMatrix = BEGIN OPEN Stor: JunoStorage, JunoMatrix; Identity: PUBLIC PROC [m: REF Matrix _ NIL] RETURNS [mid: REF Matrix] = BEGIN IF m = NIL THEN m _ NEW [Matrix]; FOR i: INTEGER IN [1..3] DO FOR j: INTEGER IN [1..3] DO m^[i][j] _ IF i = j THEN 1 ELSE 0 ENDLOOP ENDLOOP; RETURN[m] END; m1: REF Matrix _ NEW[Matrix]; -- Work matrices mInv: REF Matrix _ NEW[Matrix]; GetFrameMatrix: PUBLIC PROC [frame: Frame, m: REF Matrix _ NIL] RETURNS [mr: REF Matrix] = BEGIN IF m = NIL THEN m _ NEW [Matrix]; m^[1][3] _ frame.org.coords.x; m^[2][3] _ frame.org.coords.y; m^[3][3] _ 1; IF frame.hor=NIL THEN {m^[1][1] _ 300; m^[2][1] _ 0} ELSE {m^[1][1] _ frame.hor.coords.x-frame.org.coords.x; m^[2][1] _ frame.hor.coords.y-frame.org.coords.y}; m^[3][1] _ 0; IF frame.ver=NIL THEN {m^[1][2] _ -m^[2][1]; m^[2][2] _ m^[1][1]} ELSE {m^[1][2] _ frame.ver.coords.x-frame.org.coords.x; m^[2][2] _ frame.ver.coords.y-frame.org.coords.y}; m^[3][2] _ 0; RETURN[m] END; InvertMatrix: PUBLIC PROC [m: REF Matrix, work: REF Matrix _ NIL] RETURNS [mInv: REF Matrix, singular: BOOL] = BEGIN -- Inverts m^ into mInv^ by pivoting three times; or sets "singular" flag. -- Uses work if not NIL, otherwise allocates new matrix for result. i, j, k, l: INTEGER; c: ARRAY [1..3] OF INTEGER; pivoted: ARRAY [1..3] OF BOOLEAN _ [FALSE, FALSE, FALSE]; p: REAL; -- k is the row in which we are pivoting. -- l is the column in which we are pivoting. -- i and j are miscellaneous row and column indices respectively -- c[i] is the column of the pivot in the ith row. -- p is the reciprocal of the pivot element; also used as temp for swapping. IF work = NIL THEN work _ NEW[Matrix]; FOR k IN [1..3] DO -- set l so m^[k,l] is largest of m^[k,1], m^[k, 2], m^[k, 3], excluding -- columns in which we have already pivoted. p _ 0; FOR j IN [1 .. 3] DO IF ABS[m^[k][j]] >= p AND NOT pivoted[j] THEN {l _ j; p _ ABS[m^[k][l]]} ENDLOOP; -- We will pivot at m^[k,l], if it is not too small: IF ABS[m^[k][l]] < .0001 THEN RETURN[NIL, TRUE]; c[k] _ l; pivoted[l] _ TRUE; p _ 1.0 / m^[k][l]; m^[k][l] _ 1.0; -- divide everything in pivot row by the pivot element: FOR j IN [1..3] DO m^[k][j] _ m^[k][j] * p ENDLOOP; FOR i IN [1..3] DO IF i # k THEN FOR j IN [1..3] DO IF j # l THEN -- for each m^[i,j] outside the pivot row and column m^[i][j] _ m^[i][j] - m^[i][l] * m^[k][j]; -- note that m^[k,j] was already * p. ENDLOOP ENDLOOP; -- Finally process pivot column: FOR i IN [1..3] DO IF i # k THEN m^[i][l] _ -m^[i][l] * p ENDLOOP; ENDLOOP; -- Now we permute rows and columns: FOR i IN [1..3] DO FOR j IN [1..3] DO work^[c[i]][j] _ m^[i][c[j]] ENDLOOP ENDLOOP; RETURN[work, FALSE] END; MultiplyMatrix: PUBLIC PROC [ma, mb: REF Matrix, mc: REF Matrix _ NIL] RETURNS [mr: REF Matrix] = BEGIN i, j, k: INTEGER; sum: REAL; IF mc = NIL THEN mc _ NEW [Matrix]; FOR i IN [1..3] DO FOR j IN [1..3] DO sum _ 0.0; FOR k IN [1..3] DO sum _ sum + ma^[i][k] * mb^[k][j] ENDLOOP; mc^[i][j] _ sum ENDLOOP ENDLOOP; RETURN[mc] END; ComputeTransform: PUBLIC PROC [src, dest: Frame, mm: REF Matrix _ NIL] RETURNS [mat: REF Matrix, singular: BOOL] = BEGIN -- We want -- mat * [ src.org, src.hor, src.ver ] = [ dest.org, dest.hor, dest.ver ], -- where the points are viewed as -- column vectors with third component 1. Hence we compute the inverse -- of src and multipy on the left by dest. But the pairs -- (dest.hor, src.hor), (dest.ver, src.ver) may be missing, in which case they -- are filled in by default to make the transformation a -- translation (if both are missing) or a Euclidean motion (if just (c, sc) is missing). -- The matrix mm is optional, and is used to store the result: if NIL, a new one is allocated. -- Uses m1, mInv as work areas. sing: BOOL; mm _ GetFrameMatrix[src, mm]; [mInv, sing] _ InvertMatrix[mm, mInv]; IF sing THEN RETURN [NIL, sing]; m1 _ GetFrameMatrix[dest, m1]; mm _ MultiplyMatrix[m1, mInv, mm]; RETURN [mm, sing] END; MapCoords: PUBLIC PROC [coords: Coords, mat: REF Matrix] RETURNS [cMap: Coords] = BEGIN -- Transforms coords.x, y by the given matrix cMap.x _ coords.x * mat^[1][1] + coords.y * mat^[1][2] + mat^[1][3]; cMap.y _ coords.x * mat^[2][1] + coords.y * mat^[2][2] + mat^[2][3] END; END. Ž JunoMatrixImpl.mesa Last Edited by: Stolfi June 13, 1984 8:59:59 am PDT Tools for coordinate transformations and matrix algebra. ÊN˜šœ™™Jšœ7™7—Jšœ9™9šœÏk ˜ Jšœœ˜3——šœÏbœ˜šœ˜Jšœ ˜ ——J˜šœ˜Jšœœ#˜(šÏn œœœœ œœœ ˜HJš'œœœœœœ œœœœœœœœœœœœœœ˜®—JšŸœœ œ Ïc˜/JšŸœœ œ ˜ šŸœœœœ œœœ ˜\JšŸœœœœœbœ œœ,œƒœ œœ9œƒœœ˜­—šŸ œœœœœ œœœœ˜pJš¬Ÿœœ Kœ Dœœ œœœœœœœœœ œ *œ -œ Aœ 3œ Mœœœœœœœœ Iœ -œœœœœœœœ œ œ œ 5œœœœœœœ œ3 8œœœœœœœœ œœœœœœœ 5œ2 &œœœ !œœœœœœœœ $œœœœœœœœœœœœ˜â —šŸœœœ œ œ œœœ ˜dJš+Ÿœ œ œœœœœœœœœœœœœœ#œœœœœ˜—šŸœœœœ œœœœ˜vJš)Ÿœ œ kœ Iœ Ðcr ¡ œ Rœ ;œ kœ ZÏr ¢ œœMœœœœQœ œ˜ß—š Ÿ œœœœ œ˜SJšŸœ .œœ˜È——Jšœœ˜J˜—…—: