home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 2 / Apprentice-Release2.iso / Source Code / Pascal / Libraries / GrafSys 2.0 / GrafSys 2.0 source / Matrix2.code < prev    next >
Encoding:
Text File  |  1993-08-23  |  11.5 KB  |  390 lines  |  [TEXT/PJMM]
  1. unit Matrix;
  2.  
  3. { fast, 4x4 Matrix multiplikation, fast Ax=b multiplikation for A ist 4x4 mtrx, x ist 4-vektor (x,y,z,1) }
  4.  
  5. { for use with the 3D Project © 1991 by Christian Franz }
  6.  
  7. {$IFC UNDEFINED UseFixedMath}
  8. {$SETC UseFixedMath := FALSE}
  9. {$ENDC}
  10.  
  11.  
  12. interface
  13.  
  14. {$IFC UseFixedMath = TRUE }
  15. uses
  16.     FixedMath;
  17.  
  18. type
  19.  
  20.     Matrix4 = record
  21.             M: array[1..4, 1..4] of fixed;
  22.             identity: boolean;
  23.         end;
  24.     Vector4 = array[1..4] of fixed;
  25.  
  26. var
  27.     one: fixed;
  28.  
  29. {$ELSEC}
  30.  
  31. type
  32.  
  33.     Matrix4 = record
  34.             M: array[1..4, 1..4] of real;
  35.             identityFlag: boolean;
  36.         end;
  37.     Vector4 = array[1..4] of real;
  38.  
  39. var
  40.     one: real;
  41.  
  42. {$ENDC}
  43.  
  44. (* InitMatrix : Init matrix module *)
  45. procedure InitMatrix;
  46.  
  47. (* identity : returns identity-matrix *)
  48. function Identity: Matrix4;
  49.  
  50. (* MMult : Multiply 4x4 by 4x4 Matrix. Returns C = AB *)
  51. function MMult (var A, B: Matrix4): Matrix4;
  52.  
  53. (* VMult : Multiply the Vector x with the Matrix A giving the Vector b (return value *)
  54. function VMult (x: Vector4; var A: Matrix4): Vector4;
  55.  
  56. function MatrixVectorMult (M: Matrix4; p: Vector4): Vector4;
  57. (* GeoBench's matrix * vector multiplication... To be optimized for speed *)
  58.  
  59. (* Transpose : Generate Transpose of Matrix A *)
  60. function Transpose (A: Matrix4): Matrix4;
  61.  
  62. (* VSub : Subtract two vectors: result := x - y *)
  63. function VSub (x, y: Vector4): Vector4;
  64.  
  65. (* VAdd : Add two vectors *)
  66. function VAdd (x, y: Vector4): Vector4;
  67.  
  68. (* Set Vector to coordinates *)
  69. procedure SetVector4 (var theVector: Vector4; x, y, z: real);
  70.  
  71. (* get vector coordinates *)
  72. procedure GetVector4 (theVector: Vector4; var x, y, z: real);
  73.  
  74. (* IsVisible determines if the plane defined through the three points KLM (as vector 4)     *)
  75. (* is visible form the eye. Visibility is determined as the value of the normal vector of     *)
  76. (* the plane defined by abc. it is visible if the z component of n is greater than zero         *)
  77. (* note that for this test you must label k,l,m clockwise !!!                                     *)
  78.  
  79. function IsVisible (k, l, m: Vector4): Boolean;
  80.  
  81. implementation
  82.  
  83. {$IFC UseFixedMath = FALSE }
  84.  
  85. (* InitMatrix : Init matrix module *)
  86. procedure InitMatrix;
  87.  
  88.     begin
  89.         one := 1;
  90.     end;
  91.  
  92. { this function will generate a 4x4 Identity Matrix }
  93.  
  94. function Identity: Matrix4;
  95.  
  96.     var
  97.         I: Matrix4;
  98.         j: INTEGER;
  99.  
  100.     begin
  101.         with I do begin
  102.             for j := 1 to 4 do begin
  103.                 M[j, 1] := 0;
  104.                 M[j, 2] := 0;
  105.                 M[j, 3] := 0;
  106.                 M[j, 4] := 0;
  107.                 M[j, j] := 1;
  108.             end;
  109.             IdentityFlag := TRUE;
  110.             Identity := I;
  111.         end;
  112.     end;
  113.  
  114. function MMult (var A, B: Matrix4): Matrix4;
  115.  
  116.     var
  117.         C: Matrix4;
  118.  
  119.     begin
  120.         if A.IdentityFlag then begin
  121.             MMult := B;
  122.             Exit(MMult);
  123.         end;
  124.         if B.IdentityFlag then begin
  125.             MMult := A;
  126.             Exit(MMult);
  127.         end;
  128.  
  129.         C.M[1, 1] := A.M[1, 1] * B.M[1, 1] + A.M[1, 2] * B.M[2, 1] + A.M[1, 3] * B.M[3, 1] + A.M[1, 4] * B.M[4, 1];
  130.         C.M[1, 2] := A.M[1, 1] * B.M[1, 2] + A.M[1, 2] * B.M[2, 2] + A.M[1, 3] * B.M[3, 2] + A.M[1, 4] * B.M[4, 2];
  131.         C.M[1, 3] := A.M[1, 1] * B.M[1, 3] + A.M[1, 2] * B.M[2, 3] + A.M[1, 3] * B.M[3, 3] + A.M[1, 4] * B.M[4, 3];
  132.         C.M[1, 4] := A.M[1, 1] * B.M[1, 4] + A.M[1, 2] * B.M[2, 4] + A.M[1, 3] * B.M[3, 4] + A.M[1, 4] * B.M[4, 4];
  133.  
  134.         C.M[2, 1] := A.M[2, 1] * B.M[1, 1] + A.M[2, 2] * B.M[2, 1] + A.M[2, 3] * B.M[3, 1] + A.M[2, 4] * B.M[4, 1];
  135.         C.M[2, 2] := A.M[2, 1] * B.M[1, 2] + A.M[2, 2] * B.M[2, 2] + A.M[2, 3] * B.M[3, 2] + A.M[2, 4] * B.M[4, 2];
  136.         C.M[2, 3] := A.M[2, 1] * B.M[1, 3] + A.M[2, 2] * B.M[2, 3] + A.M[2, 3] * B.M[3, 3] + A.M[2, 4] * B.M[4, 3];
  137.         C.M[2, 4] := A.M[2, 1] * B.M[1, 4] + A.M[2, 2] * B.M[2, 4] + A.M[2, 3] * B.M[3, 4] + A.M[2, 4] * B.M[4, 4];
  138.  
  139.         C.M[3, 1] := A.M[3, 1] * B.M[1, 1] + A.M[3, 2] * B.M[2, 1] + A.M[3, 3] * B.M[3, 1] + A.M[3, 4] * B.M[4, 1];
  140.         C.M[3, 2] := A.M[3, 1] * B.M[1, 2] + A.M[3, 2] * B.M[2, 2] + A.M[3, 3] * B.M[3, 2] + A.M[3, 4] * B.M[4, 2];
  141.         C.M[3, 3] := A.M[3, 1] * B.M[1, 3] + A.M[3, 2] * B.M[2, 3] + A.M[3, 3] * B.M[3, 3] + A.M[3, 4] * B.M[4, 3];
  142.         C.M[3, 4] := A.M[3, 1] * B.M[1, 4] + A.M[3, 2] * B.M[2, 4] + A.M[3, 3] * B.M[3, 4] + A.M[3, 4] * B.M[4, 4];
  143.  
  144.         C.M[4, 1] := A.M[4, 1] * B.M[1, 1] + A.M[4, 2] * B.M[2, 1] + A.M[4, 3] * B.M[3, 1] + A.M[4, 4] * B.M[4, 1];
  145.         C.M[4, 2] := A.M[4, 1] * B.M[1, 2] + A.M[4, 2] * B.M[2, 2] + A.M[4, 3] * B.M[3, 2] + A.M[4, 4] * B.M[4, 2];
  146.         C.M[4, 3] := A.M[4, 1] * B.M[1, 3] + A.M[4, 2] * B.M[2, 3] + A.M[4, 3] * B.M[3, 3] + A.M[4, 4] * B.M[4, 3];
  147.         C.M[4, 4] := A.M[4, 1] * B.M[1, 4] + A.M[4, 2] * B.M[2, 4] + A.M[4, 3] * B.M[3, 4] + A.M[4, 4] * B.M[4, 4];
  148.         C.IdentityFlag := FALSE;
  149.         MMult := C;
  150.     end;
  151.  
  152. function VMult (x: Vector4; var A: Matrix4): Vector4;
  153.  
  154.     var
  155.         b: Vector4;
  156.  
  157.     begin
  158.         if A.IdentityFlag then begin
  159.             VMult := x;
  160.             Exit(VMult);
  161.         end;
  162.         b[1] := x[1] * A.M[1, 1] + x[2] * A.M[2, 1] + x[3] * A.M[3, 1] + x[4] * A.M[4, 1];
  163.         b[2] := x[1] * A.M[1, 2] + x[2] * A.M[2, 2] + x[3] * A.M[3, 2] + x[4] * A.M[4, 2];
  164.         b[3] := x[1] * A.M[1, 3] + x[2] * A.M[2, 3] + x[3] * A.M[3, 3] + x[4] * A.M[4, 3];
  165.         b[4] := x[1] * A.M[1, 4] + x[2] * A.M[2, 4] + x[3] * A.M[3, 4] + x[4] * A.M[4, 4];
  166.         VMult := b;
  167.     end;
  168.  
  169. {SetVector4 : set coordinates of a point or vector to the passed params }
  170. procedure SetVector4 (var theVector: Vector4; x, y, z: real);
  171.     begin
  172.         theVector[1] := x;
  173.         theVector[2] := y;
  174.         theVector[3] := z;
  175.         theVector[4] := 1;
  176.     end;
  177.  
  178. {GetVector4 : return the settings of the inner definition of vector4 }
  179. procedure GetVector4 (theVector: Vector4; var x, y, z: real);
  180.  
  181.     begin
  182.         x := theVector[1];
  183.         y := theVector[2];
  184.         z := theVector[3];
  185.     end;
  186.  
  187. (* IsVisible determines if the plane defined through the three points XYZ (as vectors 4)     *)
  188. (* is visible form the eye. Visibility is determined as the value of the normal vector of     *)
  189. (* the plane defined by abc. it is visible if the z component of n is greater than zero         *)
  190. (* note that for this test you must label x,y,z clockwise !!!                                     *)
  191.  
  192. function IsVisible (k, l, m: Vector4): Boolean;
  193.     var
  194.         nz: real;
  195.         ax, bx, ay, by: real;
  196.  
  197.     begin
  198. {ax := (l[1] - k[1]);{}
  199. {bx := (m[1] - k[1]);{}
  200. {ay := (l[2] - k[2]);{}
  201. {by := (m[2] - k[2]);{}
  202. {nz := ax * by - ay * bx;{}
  203.         nz := (l[1] - k[1]) * (m[2] - k[2]) - (l[2] - k[2]) * (m[1] - k[1]);
  204.         IsVisible := nz >= 0;
  205.     end;
  206.  
  207.  
  208. {$ELSEC}
  209. { BEGIN FIXED-MATH routines }
  210.  
  211. (* InitMatrix : Init matrix module *)
  212. procedure InitMatrix;
  213.     begin
  214.         one := FixRatio(1, 1);
  215.     end;
  216.  
  217. { this function will generate a 4x4 Identity Matrix }
  218.  
  219. function Identity: Matrix4;
  220.  
  221.     var
  222.         I: Matrix4;
  223.         j: INTEGER;
  224.  
  225.     begin
  226.         with I do
  227.             for j := 1 to 4 do begin
  228.                 I[j, 1] := 0;
  229.                 I[j, 2] := 0;
  230.                 I[j, 3] := 0;
  231.                 I[j, 4] := 0;
  232.                 I[j, j] := FixRatio(1, 1);
  233.             end;
  234.         I.IdentityFlag := TRUE;
  235.         Identity := I;
  236.     end;
  237.  
  238. function MMult (var A, B: Matrix4): Matrix4;
  239.  
  240.     var
  241.         C: Matrix4;
  242.  
  243.     begin
  244.         if A.IdentityFlag then begin
  245.             MMult := B;
  246.             Exit(MMult);
  247.         end;
  248.         if B.IdentityFlag then begin
  249.             MMult := A;
  250.             Exit(MMult);
  251.         end;
  252.  
  253.         C.M[1, 1] := FixMul(A.M[1, 1], B.M[1, 1]) + FixMul(A.M[1, 2], B.M[2, 1]) + FixMul(A.M[1, 3], B.M[3, 1]) + FixMul(A.M[1, 4], B.M[4, 1]);
  254.         C.M[1, 2] := FixMul(A.M[1, 1], B.M[1, 2]) + FixMul(A.M[1, 2], B.M[2, 2]) + FixMul(A.M[1, 3], B.M[3, 2]) + FixMul(A.M[1, 4], B.M[4, 2]);
  255.         C.M[1, 3] := FixMul(A.M[1, 1], B.M[1, 3]) + FixMul(A.M[1, 2], B.M[2, 3]) + FixMul(A.M[1, 3], B.M[3, 3]) + FixMul(A.M[1, 4], B.M[4, 3]);
  256.         C.M[1, 4] := FixMul(A.M[1, 1], B.M[1, 4]) + FixMul(A.M[1, 2], B.M[2, 4]) + FixMul(A.M[1, 3], B.M[3, 4]) + FixMul(A.M[1, 4], B.M[4, 4]);
  257.  
  258.         C.M[2, 1] := FixMul(A.M[2, 1], B.M[1, 1]) + FixMul(A.M[2, 2], B.M[2, 1]) + FixMul(A.M[2, 3], B.M[3, 1]) + FixMul(A.M[2, 4], B.M[4, 1]);
  259.         C.M[2, 2] := FixMul(A.M[2, 1], B.M[1, 2]) + FixMul(A.M[2, 2], B.M[2, 2]) + FixMul(A.M[2, 3], B.M[3, 2]) + FixMul(A.M[2, 4], B.M[4, 2]);
  260.         C.M[2, 3] := FixMul(A.M[2, 1], B.M[1, 3]) + FixMul(A.M[2, 2], B.M[2, 3]) + FixMul(A.M[2, 3], B.M[3, 3]) + FixMul(A.M[2, 4], B.M[4, 3]);
  261.         C.M[2, 4] := FixMul(A.M[2, 1], B.M[1, 4]) + FixMul(A.M[2, 2], B.M[2, 4]) + FixMul(A.M[2, 3], B.M[3, 4]) + FixMul(A.M[2, 4], B.M[4, 4]);
  262.  
  263.         C.M[3, 1] := FixMul(A.M[3, 1], B.M[1, 1]) + FixMul(A.M[3, 2], B.M[2, 1]) + FixMul(A.M[3, 3], B.M[3, 1]) + FixMul(A.M[3, 4], B.M[4, 1]);
  264.         C.M[3, 2] := FixMul(A.M[3, 1], B.M[1, 2]) + FixMul(A.M[3, 2], B.M[2, 2]) + FixMul(A.M[3, 3], B.M[3, 2]) + FixMul(A.M[3, 4], B.M[4, 2]);
  265.         C.M[3, 3] := FixMul(A.M[3, 1], B.M[1, 3]) + FixMul(A.M[3, 2], B.M[2, 3]) + FixMul(A.M[3, 3], B.M[3, 3]) + FixMul(A.M[3, 4], B.M[4, 3]);
  266.         C.M[3, 4] := FixMul(A.M[3, 1], B.M[1, 4]) + FixMul(A.M[3, 2], B.M[2, 4]) + FixMul(A.M[3, 3], B.M[3, 4]) + FixMul(A.M[3, 4], B.M[4, 4]);
  267.  
  268.         C.M[4, 1] := FixMul(A.M[4, 1], B.M[1, 1]) + FixMul(A.M[4, 2], B.M[2, 1]) + FixMul(A.M[4, 3], B.M[3, 1]) + FixMul(A.M[4, 4], B.M[4, 1]);
  269.         C.M[4, 2] := FixMul(A.M[4, 1], B.M[1, 2]) + FixMul(A.M[4, 2], B.M[2, 2]) + FixMul(A.M[4, 3], B.M[3, 2]) + FixMul(A.M[4, 4], B.M[4, 2]);
  270.         C.M[4, 3] := FixMul(A.M[4, 1], B.M[1, 3]) + FixMul(A.M[4, 2], B.M[2, 3]) + FixMul(A.M[4, 3], B.M[3, 3]) + FixMul(A.M[4, 4], B.M[4, 3]);
  271.         C.M[4, 4] := FixMul(A.M[4, 1], B.M[1, 4]) + FixMul(A.M[4, 2], B.M[2, 4]) + FixMul(A.M[4, 3], B.M[3, 4]) + FixMul(A.M[4, 4], B.M[4, 4]);
  272.  
  273.         MMult := C;
  274.     end;
  275.  
  276. function VMult (x: Vector4; var A: Matrix4): Vector4;
  277.  
  278.     var
  279.         b: Vector4;
  280.  
  281.     begin
  282.         if A.IdentityFlag then begin
  283.             VMult := x;
  284.             Exit(VMult);
  285.         end;
  286.         b[1] := FixMul(x[1], A.M[1, 1]) + FixMul(x[2], A.M[2, 1]) + FixMul(x[3], A.M[3, 1]) + FixMul(x[4], A.M[4, 1]);
  287.         b[2] := FixMul(x[1], A.M[1, 2]) + FixMul(x[2], A.M[2, 2]) + FixMul(x[3], A.M[3, 2]) + FixMul(x[4], A.M[4, 2]);
  288.         b[3] := FixMul(x[1], A.M[1, 3]) + FixMul(x[2], A.M[2, 3]) + FixMul(x[3], A.M[3, 3]) + FixMul(x[4], A.M[4, 3]);
  289.         b[4] := FixMul(x[1], A.M[1, 4]) + FixMul(x[2], A.M[2, 4]) + FixMul(x[3], A.M[3, 4]) + FixMul(x[4], A.M[4, 4]);
  290.         VMult := b;
  291.     end;
  292.  
  293. {SetVector4 : set coordinates of a point or vector to the passed params }
  294. procedure SetVector4 (var theVector: Vector4; x, y, z: real);
  295.     begin
  296.         theVector[1] := X2Fix(x);
  297.         theVector[2] := X2Fix(y);
  298.         theVector[3] := X2Fix(z);
  299.         theVector[4] := FixRatio(1, 1);
  300.     end;
  301.  
  302. {GetVector4 : return the settings of the inner definition of vector4 }
  303. procedure GetVector4 (theVector: Vector4; var x, y, z: real);
  304.  
  305.     begin
  306.         x := Fix2X(theVector[1]);
  307.         y := Fix2X(theVector[2]);
  308.         z := Fix2X(theVector[3]);
  309.     end;
  310.  
  311. function IsVisible (k, l, m: Vector4): Boolean;
  312.     var
  313.         nz: Fixed;
  314. {ax, bx, ay, by: real;}
  315.  
  316.     begin
  317. {ax := (l[1] - k[1]);{}
  318. {bx := (m[1] - k[1]);{}
  319. {ay := (l[2] - k[2]);{}
  320. {by := (m[2] - k[2]);{}
  321. {nz := ax * by - ay * bx;{}
  322.         nz := FixMul((l[1] - k[1]), (m[2] - k[2])) - FixMul((l[2] - k[2]), (m[1] - k[1]));
  323.         IsVisible := nz >= 0;
  324.     end;
  325.  
  326. {$ENDC}
  327.  
  328. (* PROCEDURE THAT ARE THE SAME FOR BOTH VERSIONS *)
  329.  
  330. (* Transpose : Generate Transpose of Matrix A . This proc is identical to both fixed an float version*)
  331. function Transpose (A: Matrix4): Matrix4;
  332.  
  333.     var
  334.         T: Matrix4;
  335.         index: Integer;
  336.  
  337.     begin
  338.         with T do
  339.             for index := 1 to 4 do begin
  340.                 M[index, 1] := A.M[1, index];
  341.                 M[index, 2] := A.M[2, index];
  342.                 M[index, 3] := A.M[3, index];
  343.                 M[index, 4] := A.M[4, index];
  344.             end;
  345.         Transpose := T;
  346.     end;
  347.  
  348. (* VSub : Subtract two vectors: result := x - y *)
  349. function VSub (x, y: Vector4): Vector4;
  350.     var
  351.         theResult: Vector4;
  352.  
  353.     begin
  354.         theResult[1] := x[1] - y[1];
  355.         theResult[2] := x[2] - y[2];
  356.         theResult[3] := x[3] - y[3];
  357.         theResult[4] := one;
  358.         VSub := theResult;
  359.     end;
  360.  
  361.  
  362. (* VAdd : Add two vectors *)
  363. function VAdd (x, y: Vector4): Vector4;
  364.     var
  365.         theResult: Vector4;
  366.  
  367.     begin
  368.         theResult[1] := x[1] + y[1];
  369.         theResult[2] := x[2] + y[2];
  370.         theResult[3] := x[3] + y[3];
  371.         theResult[4] := one;
  372.         VAdd := theResult;
  373.     end;
  374.  
  375. (* mvm : from geobench the matrix times vector multiplication *)
  376. function MatrixVectorMult (M: Matrix4; p: Vector4): Vector4;
  377.     var
  378.         i, j: integer;
  379.         res: Vector4;
  380.  
  381.     begin
  382.         for i := 1 to 4 do begin
  383.             res[i] := 0;
  384.             for j := 1 to 4 do
  385.                 res[i] := res[i] + M.M[i, j] * p[j];
  386.         end;
  387.         MatrixVectorMult := res;
  388.  
  389.     end;
  390. end.