home *** CD-ROM | disk | FTP | other *** search
/ Delphi 5 for Professionals / DELPHI5.iso / AddOns / Components / TEECHART / Src Code / GEOMETRY.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1998-10-24  |  52.1 KB  |  1,593 lines
  1. unit Geometry;
  2.  
  3. // This unit contains many needed types, functions and procedures for
  4. // quaternion, vector and matrix arithmetics. It is specifically designed
  5. // for geometric calculations within R3 (affine vector space)
  6. // and R4 (homogeneous vector space).
  7. //
  8. // Identifiers containing no dimensionality (as affine or homogeneous)
  9. // and no datatype (integer..extended) are supposed as R4 representation
  10. // with 'single' floating point type (examples are TVector, TMatrix,
  11. // and TQuaternion). The default data type is 'single' ('GLFloat' for OpenGL)
  12. // and used in all routines (except conversions).
  13. //
  14. // Routines with an open array as argument can either take
  15. // Func([1,2,3,4,..]) or Func(Vect). The latter is prefered, since
  16. // no extra stack operations is required.
  17. // NOTE: Be careful while passing open array elements! If you pass more elements
  18. // than there's room in the result, then the behaviour will be unpredictable.
  19. //
  20. // If not otherwise stated, all angles are given in radians
  21. // (instead of degrees). Use RadToDeg or DegToRad from Math.pas
  22. // to convert between them.
  23. //
  24. // Geometry.pas was assembled from different sources (like GraphicGems)
  25. // and relevant books or based on self written code, respectivly.
  26. //
  27. // NOTE: Some aspects need to be considered when using Delphi and pure
  28. //       assembler code. Delphi esnures that the direction flag is always
  29. //       cleared while entering a function and expects it cleared on return.
  30. //       The registeres EDI, ESI and EBX (as well as the stack management
  31. //       registers EBP and ESP) must not be changed! EAX, ECX and EDX are
  32. //       freely available.
  33. //
  34. // last change : 03. January 1998
  35. // Done by: Dipl.-Ing. Mike Lischke (Lischke@hotmail.com)
  36.  
  37. interface
  38.  
  39. uses Math;
  40.  
  41. type // data types needed for 3D graphics calculation,
  42.      // included are 'C like' aliases for each type (to be
  43.      // conformal with OpenGL types)
  44.  
  45.      PByte            = ^Byte;
  46.      PWord            = ^Word;
  47.      PInteger         = ^Integer;
  48.      PFloat           = ^Single;
  49.      PDouble          = ^Double;
  50.      PExtended        = ^Extended;
  51.      PPointer         = ^Pointer;
  52.  
  53.      // types to specify continous streams of a specific type
  54.      // switch off range checking to access values beyond the limits
  55.      PByteVector = ^TByteVector;
  56.      TByteVector = array[Word] of Byte;
  57.  
  58.      PWordVector = ^TWordVector;
  59.      TWordVector = array[Word] of Word;
  60.  
  61.      PIntVector = ^TIntVector;
  62.      TIntVector = array[Word] of Integer;
  63.  
  64.      PFloatVector = ^TFloatVector;
  65.      TFloatVector = array[Word] of Single;
  66.  
  67.      PDblVector = ^TDblVector;
  68.      TDblVector = array[Word] of Double;
  69.  
  70.      // common vector and matrix types
  71.      // indices correspond like: x -> 0
  72.      //                          y -> 1
  73.      //                          z -> 2
  74.      //                          w -> 3
  75.  
  76.      PHomogeneousByteVector = ^THomogeneousByteVector;
  77.      THomogeneousByteVector = array[0..3] of Byte;
  78.      TVector4b              = THomogeneousByteVector;
  79.  
  80.      PHomogeneousWordVector = ^THomogeneousWordVector;
  81.      THomogeneousWordVector = array[0..3] of Word;
  82.      TVector4w              = THomogeneousWordVector;
  83.  
  84.      PHomogeneousIntVector  = ^THomogeneousIntVector;
  85.      THomogeneousIntVector  = array[0..3] of Integer;
  86.      TVector4i              = THomogeneousIntVector;
  87.  
  88.      PHomogeneousFltVector  = ^THomogeneousFltVector;
  89.      THomogeneousFltVector  = array[0..3] of Single;
  90.      TVector4f              = THomogeneousFltVector;
  91.  
  92.      PHomogeneousDblVector  = ^THomogeneousDblVector;
  93.      THomogeneousDblVector  = array[0..3] of Double;
  94.      TVector4d              = THomogeneousDblVector;
  95.  
  96.      PHomogeneousExtVector  = ^THomogeneousExtVector;
  97.      THomogeneousExtVector  = array[0..3] of Extended;
  98.      TVector4e              = THomogeneousExtVector;
  99.  
  100.      PHomogeneousPtrVector  = ^THomogeneousPtrVector;
  101.      THomogeneousPtrVector  = array[0..3] of Pointer;
  102.      TVector4p              = THomogeneousPtrVector;
  103.  
  104.      PAffineByteVector      = ^TAffineByteVector;
  105.      TAffineByteVector      = array[0..2] of Byte;
  106.      TVector3b              = TAffineByteVector;
  107.  
  108.      PAffineWordVector      = ^TAffineWordVector;
  109.      TAffineWordVector      = array[0..2] of Word;
  110.      TVector3w              = TAffineWordVector;
  111.  
  112.      PAffineIntVector       = ^TAffineIntVector;
  113.      TAffineIntVector       = array[0..2] of Integer;
  114.      TVector3i              = TAffineIntVector;
  115.  
  116.      PAffineFltVector       = ^TAffineFltVector;
  117.      TAffineFltVector       = array[0..2] of Single;
  118.      TVector3f              = TAffineFltVector;
  119.  
  120.      PAffineDblVector       = ^TAffineDblVector;
  121.      TAffineDblVector       = array[0..2] of Double;
  122.      TVector3d              = TAffineDblVector;
  123.  
  124.      PAffineExtVector       = ^TAffineExtVector;
  125.      TAffineExtVector       = array[0..2] of Extended;
  126.      TVector3e              = TAffineExtVector;
  127.  
  128.      PAffinePtrVector       = ^TAffinePtrVector;
  129.      TAffinePtrVector       = array[0..2] of Pointer;
  130.      TVector3p              = TAffinePtrVector;
  131.  
  132.      // some simplified names
  133.      PVector                = ^TVector;
  134.      TVector                = THomogeneousFltVector;
  135.  
  136.      PHomogeneousVector     = ^THomogeneousVector;
  137.      THomogeneousVector     = THomogeneousFltVector;
  138.  
  139.      PAffineVector          = ^TAffineVector;
  140.      TAffineVector          = TAffineFltVector;
  141.  
  142.      PVectorArray           = ^TVectorArray;
  143.      TVectorArray           = array[Word] of TAffineVector;
  144.  
  145.      // matrices
  146.      THomogeneousByteMatrix = array[0..3] of THomogeneousByteVector;
  147.      TMatrix4b              = THomogeneousByteMatrix;
  148.  
  149.      THomogeneousWordMatrix = array[0..3] of THomogeneousWordVector;
  150.      TMatrix4w              = THomogeneousWordMatrix;
  151.  
  152.      THomogeneousIntMatrix  = array[0..3] of THomogeneousIntVector;
  153.      TMatrix4i              = THomogeneousIntMatrix;
  154.  
  155.      THomogeneousFltMatrix  = array[0..3] of THomogeneousFltVector;
  156.      TMatrix4f              = THomogeneousFltMatrix;
  157.  
  158.      THomogeneousDblMatrix  = array[0..3] of THomogeneousDblVector;
  159.      TMatrix4d              = THomogeneousDblMatrix;
  160.  
  161.      THomogeneousExtMatrix  = array[0..3] of THomogeneousExtVector;
  162.      TMatrix4e              = THomogeneousExtMatrix;
  163.  
  164.      TAffineByteMatrix      = array[0..2] of TAffineByteVector;
  165.      TMatrix3b              = TAffineByteMatrix;
  166.  
  167.      TAffineWordMatrix      = array[0..2] of TAffineWordVector;
  168.      TMatrix3w              = TAffineWordMatrix;
  169.  
  170.      TAffineIntMatrix       = array[0..2] of TAffineIntVector;
  171.      TMatrix3i              = TAffineIntMatrix;
  172.  
  173.      TAffineFltMatrix       = array[0..2] of TAffineFltVector;
  174.      TMatrix3f              = TAffineFltMatrix;
  175.  
  176.      TAffineDblMatrix       = array[0..2] of TAffineDblVector;
  177.      TMatrix3d              = TAffineDblMatrix;
  178.  
  179.      TAffineExtMatrix       = array[0..2] of TAffineExtVector;
  180.      TMatrix3e              = TAffineExtMatrix;
  181.  
  182.      // some simplified names
  183.      PMatrix                = ^TMatrix;
  184.      TMatrix                = THomogeneousFltMatrix;
  185.  
  186.      PHomogeneousMatrix     = ^THomogeneousMatrix;
  187.      THomogeneousMatrix     = THomogeneousFltMatrix;
  188.  
  189.      PAffineMatrix          = ^TAffineMatrix;
  190.      TAffineMatrix          = TAffineFltMatrix;
  191.  
  192.      // q=([x,y,z],w)
  193.      TQuaternion        = record
  194.                             case Integer of
  195.                               0 : (Axis   : TAffineVector;
  196.                                    Angle  : Single);
  197.                               1 : (Vector : TVector);
  198.                           end;
  199.  
  200.      TRectangle         = record
  201.                             Left, Top, Width, Height: Integer;
  202.                           end;
  203.  
  204.      TTransType = (ttScaleX,ttScaleY,ttScaleZ,
  205.                    ttShearXY,ttShearXZ,ttShearYZ,
  206.                    ttRotateX,ttRotateY,ttRotateZ,
  207.                    ttTranslateX,ttTranslateY,ttTranslateZ,
  208.                    ttPerspectiveX,ttPerspectiveY,ttPerspectiveZ,ttPerspectiveW);
  209.  
  210.      // used to describe a sequence of transformations in following order:
  211.      // [Sx][Sy][Sz][ShearXY][ShearXZ][ShearZY][Rx][Ry][Rz][Tx][Ty][Tz][P(x,y,z,w)]
  212.      // constants are declared for easier access (see MatrixDecompose below)
  213.      TTransformations  = array[TTransType] of Extended;
  214.  
  215.     
  216. const // useful constants
  217.  
  218.       // standard vectors
  219.       XVector    : TAffineVector = (1,0,0);
  220.       YVector    : TAffineVector = (0,1,0);
  221.       ZVector    : TAffineVector = (0,0,1);
  222.       NullVector : TAffineVector = (0,0,0);
  223.  
  224.       IdentityMatrix : TMatrix = ((1,0,0,0),
  225.                                   (0,1,0,0),
  226.                                   (0,0,1,0),
  227.                                   (0,0,0,1));
  228.       EmptyMatrix    : TMatrix = ((0,0,0,0),
  229.                                   (0,0,0,0),
  230.                                   (0,0,0,0),
  231.                                   (0,0,0,0));
  232.       // some very small numbers
  233.       EPSILON  = 1E-100;
  234.       EPSILON2 = 1E-50;
  235.  
  236. //------------------------------------------------------------------------------
  237.  
  238. // vector functions
  239. function  VectorAdd(V1,V2: TVector): TVector;
  240. function  VectorAffineAdd(V1,V2: TAffineVector): TAffineVector;
  241. function  VectorAffineCombine(V1,V2: TAffineVector; F1,F2: Single): TAffineVector;
  242. function  VectorAffineDotProduct(V1,V2: TAffineVector): Single;
  243. function  VectorAffineLerp(V1,V2: TAffineVector; t: Single): TAffineVector;
  244. function  VectorAffineSubtract(V1,V2: TAffineVector): TAffineVector;
  245. function  VectorAngle(V1,V2: TAffineVector): Single;
  246. function  VectorCombine(V1,V2: TVector; F1,F2: Single): TVector;
  247. function  VectorCrossProduct(V1,V2: TAffineVector): TAffineVector;
  248. function  VectorDotProduct(V1,V2: TVector): Single;
  249. function  VectorLength(V: array of Single): Single;
  250. function  VectorLerp(V1,V2: TVector; t: Single): TVector;
  251. procedure VectorNegate(V: array of Single);
  252. function  VectorNorm(V: array of Single): Single; 
  253. function  VectorNormalize(V: array of Single): Single;
  254. function  VectorPerpendicular(V,N: TAffineVector): TAffineVector;
  255. function  VectorReflect(V, N: TAffineVector): TAffineVector;
  256. procedure VectorScale(V: array of Single; Factor: Single);
  257. function  VectorSubtract(V1,V2: TVector): TVector;
  258.  
  259. // matrix functions
  260. function  CreateRotationMatrixX(Sine, Cosine: Single): TMatrix;
  261. function  CreateRotationMatrixY(Sine, Cosine: Single): TMatrix;
  262. function  CreateRotationMatrixZ(Sine, Cosine: Single): TMatrix;
  263. function  CreateScaleMatrix(V: TAffineVector): TMatrix;
  264. function  CreateTranslationMatrix(V: TAffineVector): TMatrix;
  265. procedure MatrixAdjoint(var M: TMatrix);
  266. function  MatrixAffineDeterminant(M: TAffineMatrix): Single;
  267. procedure MatrixAffineTranspose(var M: TAffineMatrix);
  268. function  MatrixDeterminant(M: TMatrix): Single;
  269. procedure MatrixInvert(var M: TMatrix);
  270. function  MatrixMultiply(M1, M2: TMatrix): TMatrix;
  271. procedure MatrixScale(var M: TMatrix; Factor: Single);
  272. procedure MatrixTranspose(var M: TMatrix);
  273.  
  274. // quaternion functions
  275. function  QuaternionConjugate(Q: TQuaternion): TQuaternion;
  276. function  QuaternionFromPoints(V1,V2: TAffineVector): TQuaternion;
  277. function  QuaternionMultiply(qL, qR: TQuaternion): TQuaternion;
  278. function  QuaternionSlerp(QStart,QEnd: TQuaternion; Spin: Integer; t: Single): TQuaternion;
  279. function  QuaternionToMatrix(Q: TQuaternion): TMatrix;
  280. procedure QuaternionToPoints(Q: TQuaternion; var ArcFrom, ArcTo: TVector);
  281.  
  282. // mixed functions
  283. function  ConvertRotation(Angles: TAffineVector): TVector;
  284. function  CreateRotationMatrix(Axis: TAffineVector; Angle: Single): TMatrix;
  285. function  MatrixDecompose(M: TMatrix; var Tran: TTransformations): Boolean;
  286. function  VectorAffineTransform(V: TAffineVector; M: TAffineMatrix): TAffineVector;
  287. function  VectorTransform(V: TVector; M: TMatrix): TVector;
  288.  
  289. // miscellaneous functions
  290. function  MakeAffineDblVector(V: array of Double): TAffineDblVector;
  291. function  MakeDblVector(V: array of Double): THomogeneousDblVector; 
  292. function  MakeAffineVector(V: array of Single): TAffineVector;
  293. function  MakeVector(V: array of Single): TVector;
  294. function  PointInPolygon(xp, yp : array of Single; x,y: Single): Boolean;
  295. function  VectorAffineDblToFlt(V: TAffineDblVector): TAffineVector;
  296. function  VectorDblToFlt(V: THomogeneousDblVector): THomogeneousVector;
  297. function  VectorAffineFltToDbl(V: TAffineVector): TAffineDblVector;
  298. function  VectorFltToDbl(V: TVector): THomogeneousDblVector;
  299.  
  300. //------------------------------------------------------------------------------
  301.  
  302. implementation
  303.  
  304. const // FPU status flags (high order byte)
  305.       C0 =   1;
  306.       C1 =   2;
  307.       C2 =   4;
  308.       C3 = $40;
  309.  
  310.       // to be used as descriptive indices
  311.       X  =   0;
  312.       Y  =   1;
  313.       Z  =   2;
  314.       W  =   3;
  315.  
  316. //------------------------------------------------------------------------------
  317.  
  318. function MakeAffineDblVector(V: array of Double): TAffineDblVector; register; assembler;
  319.  
  320. // create a vector from given values
  321. // EAX contains address of V
  322. // ECX contains address to result vector
  323. // EDX contains highest index of V
  324.  
  325. asm
  326.               PUSH EDI
  327.               PUSH ESI
  328.               MOV EDI,ECX
  329.               MOV ESI,EAX
  330.               MOV ECX,6
  331.               REP MOVSD
  332.               POP ESI
  333.               POP EDI
  334. end;
  335.  
  336. //------------------------------------------------------------------------------
  337.  
  338. function MakeDblVector(V: array of Double): THomogeneousDblVector; register; assembler;
  339.  
  340. // create a vector from given values
  341. // EAX contains address of V
  342. // ECX contains address to result vector
  343. // EDX contains highest index of V
  344.  
  345. asm
  346.               PUSH EDI
  347.               PUSH ESI
  348.               MOV EDI,ECX
  349.               MOV ESI,EAX
  350.               MOV ECX,8
  351.               REP MOVSD
  352.               POP ESI
  353.               POP EDI
  354. end;
  355.  
  356. //------------------------------------------------------------------------------
  357.  
  358. function MakeAffineVector(V: array of Single): TAffineVector; register; assembler;
  359.  
  360. // create a vector from given values
  361. // EAX contains address of V
  362. // ECX contains address to result vector
  363. // EDX contains highest index of V
  364.  
  365. asm
  366.               PUSH EDI
  367.               PUSH ESI
  368.               MOV EDI,ECX
  369.               MOV ESI,EAX
  370.               MOV ECX,3
  371.               REP MOVSD
  372.               POP ESI
  373.               POP EDI
  374. end;
  375.  
  376. //------------------------------------------------------------------------------
  377.  
  378. function MakeVector(V: array of Single): TVector; register; assembler;
  379.  
  380. // create a vector from given values
  381. // EAX contains address of V
  382. // ECX contains address to result vector
  383. // EDX contains highest index of V
  384.  
  385. asm
  386.               PUSH EDI
  387.               PUSH ESI
  388.               MOV EDI,ECX
  389.               MOV ESI,EAX
  390.               MOV ECX,4
  391.               REP MOVSD
  392.               POP ESI
  393.               POP EDI
  394. end;
  395.  
  396. //------------------------------------------------------------------------------
  397.  
  398. function VectorLength(V: array of Single): Single; register; assembler;
  399.  
  400. // calculates the length of a vector following the equation: sqrt(x*x+y*y+...)
  401. // Note: The parameter of this function is declared as open array. Thus
  402. // there's no restriction about the number of the components of the vector.
  403. //
  404. // EAX contains the pointer to the data and EDX the number of array members.
  405. // The result is returned in ST(0) and will be automatically converted, if a
  406. // non-Single value is needed.
  407.  
  408. asm
  409.               FLDZ                          // initialize sum
  410. @@Loop:       FLD  DWORD PTR [EAX+4*EDX]    // load a component
  411.               FMUL ST,ST
  412.               FADDP
  413.               SUB  EDX,1
  414.               JNL  @@Loop
  415.               FSQRT
  416. end;
  417.  
  418. //------------------------------------------------------------------------------
  419.  
  420. function VectorAngle(V1,V2: TAffineVector): Single; register; assembler;
  421.  
  422. // calculates the cosine of the angle between Vector1 and Vector2
  423. // Result = DotProduct(V1,V2) / (Length(V1)*Length(V2))
  424. //
  425. // EAX contains Address of Vector1
  426. // EDX contains Address of Vector2
  427.  
  428. asm
  429.               FLD DWORD PTR [EAX]           // V1[0]
  430.               FLD ST                        // Single V1[0]
  431.               FMUL ST,ST                    // V1[0]**2 (prep. for divisor)
  432.               FLD DWORD PTR [EDX]           // V2[0]
  433.               FMUL ST(2),ST                 // ST(2):=V1[0]*V2[0]
  434.               FMUL ST,ST                    // V2[0]**2 (prep. for divisor)
  435.               FLD DWORD PTR [EAX+4]         // V1[1]
  436.               FLD ST                        // Single V1[1]
  437.               FMUL ST,ST                    // ST(0):=V1[1]**2
  438.               FADDP ST(3),ST                // ST(2):=V1[0]**2+V1[1]**2
  439.               FLD DWORD PTR [EDX+4]         // V2[1]
  440.               FMUL ST(1),ST                 // ST(1):=V1[1]*V2[1]
  441.               FMUL ST,ST                    // ST(0):=V2[1]**2
  442.               FADDP ST(2),ST                // ST(1):=V2[0]**2+V2[1]**2
  443.               FADDP ST(3),ST                // ST(2):=V1[0]*V2[0]+V1[1]*V2[1]
  444.               FLD DWORD PTR [EAX+8]         // load V2[1]
  445.               FLD ST                        // same calcs go here
  446.               FMUL ST,ST                    // (compare above)
  447.               FADDP ST(3),ST
  448.               FLD DWORD PTR [EDX+8]
  449.               FMUL ST(1),ST
  450.               FMUL ST,ST
  451.               FADDP ST(2),ST
  452.               FADDP ST(3),ST
  453.               FMULP                         // ST(0):=(V1[0]**2+V1[1]**2+V1[2])*
  454.                                             //        (V2[0]**2+V2[1]**2+V2[2])
  455.               FSQRT                         // sqrt(ST(0))
  456.               FDIVP                         // ST(0)=Result:=ST(1)/ST(0)
  457.   // the result is expected in ST(0), if it's invalid, an error is raised
  458. end;
  459.  
  460. //------------------------------------------------------------------------------
  461.  
  462. function PointInPolygon(xp, yp : array of Single; x,y: Single): Boolean;
  463.  
  464. // The code below is from Wm. Randolph Franklin <wrf@ecse.rpi.edu>
  465. // with some minor modifications for speed.  It returns 1 for strictly
  466. // interior points, 0 for strictly exterior, and 0 or 1 for points on
  467. // the boundary.
  468.  
  469. var I, J : Integer;
  470.  
  471. begin
  472.   Result:=False;
  473.   if High(XP) <> High(YP) then Exit;
  474.   J:=High(XP);
  475.   for I:=0 to High(XP) do
  476.   begin
  477.     if ((((yp[I]<=y) and (y<yp[J])) OR ((yp[J]<=y) and (y<yp[I]))) and
  478.         (x < (xp[J]-xp[I])*(y-yp[I])/(yp[J]-yp[I])+xp[I]))
  479.     then Result:=not Result;
  480.     J:=I+1;
  481.   end;
  482. end;
  483.  
  484. //------------------------------------------------------------------------------
  485.  
  486. function QuaternionConjugate(Q: TQuaternion): TQuaternion; register; assembler;
  487.  
  488. // return the conjugate of a quaternion
  489. // EAX contains address of Q
  490. // EDX contains address of result
  491.  
  492. asm
  493.               FLD DWORD PTR [EAX]
  494.               FCHS
  495.               WAIT
  496.               FSTP DWORD PTR [EDX]
  497.               FLD DWORD PTR [EAX+4]
  498.               FCHS
  499.               WAIT
  500.               FSTP DWORD PTR [EDX+4]
  501.               FLD DWORD PTR [EAX+8]
  502.               FCHS
  503.               WAIT
  504.               FSTP DWORD PTR [EDX+8]
  505.               MOV EAX,[EAX+12]
  506.               MOV [EDX+12],EAX
  507. end;
  508.  
  509. //------------------------------------------------------------------------------
  510.  
  511. function QuaternionFromPoints(V1,V2: TAffineVector): TQuaternion; register; assembler;
  512.  
  513. // construct a unit quaternion from two points on unit sphere
  514. // EAX contains address of V1
  515. // ECX contains address to result
  516. // EDX contains address of V2
  517.  
  518. asm
  519.   {Result.Vector[X]:= V1[Y]*V2[Z] - V1[Z]*V2[Y];
  520.   Result.Vector[Y]:= V1[Z]*V2[X] - V1[X]*V2[Z];
  521.   Result.Vector[Z]:= V1[X]*V2[Y] - V1[Y]*V2[X];
  522.   Result.Angle:= V1[X]*V2[X] + V1[Y]*V2[Y] + V1[Z]*V2[Z];}
  523.               FLD DWORD PTR [EDX+8]    // first load both vectors onto FPU register stack
  524.               FLD DWORD PTR [EDX+4]
  525.               FLD DWORD PTR [EDX+0]
  526.               FLD DWORD PTR [EAX+8]
  527.               FLD DWORD PTR [EAX+4]
  528.               FLD DWORD PTR [EAX+0]
  529.  
  530.               FLD ST(1)                // ST(0):=V1[Y]
  531.               FMUL ST,ST(6)            // ST(0):=V1[Y]*V2[Z]
  532.               FLD ST(3)                // ST(0):=V1[Z]
  533.               FMUL ST,ST(6)            // ST(0):=V1[Z]*V2[Y]
  534.               FSUBP ST(1),ST           // ST(0):=ST(1)-ST(0)
  535.               WAIT
  536.               FSTP DWORD [ECX]         // Result.Vector[X]:=ST(0)
  537.               FLD ST(2)                // ST(0):=V1[Z]
  538.               FMUL ST,ST(4)            // ST(0):=V1[Z]*V2[X]
  539.               FLD ST(1)                // ST(0):=V1[X]
  540.               FMUL ST,ST(7)            // ST(0):=V1[X]*V2[Z]
  541.               FSUBP ST(1),ST           // ST(0):=ST(1)-ST(0)
  542.               WAIT
  543.               FSTP DWORD [ECX+4]       // Result.Vector[Y]:=ST(0)
  544.               FLD ST                   // ST(0):=V1[X]
  545.               FMUL ST,ST(5)            // ST(0):=V1[X]*V2[Y]
  546.               FLD ST(2)                // ST(0):=V1[Y]
  547.               FMUL ST,ST(5)            // ST(0):=V1[Y]*V2[X]
  548.               FSUBP ST(1),ST           // ST(0):=ST(1)-ST(0)
  549.               WAIT
  550.               FSTP DWORD [ECX+8]       // Result.Vector[Z]:=ST(0)
  551.               FMUL ST,ST(3)            // ST(0):=V1[X]*V2[X]
  552.               FLD ST(1)                // ST(0):=V1[Y]
  553.               FMUL ST,ST(5)            // ST(0):=V1[Y]*V2[Y]
  554.               FADDP ST(1),ST           // ST(0):=V1[X]*V2[X]+V1[Y]*V2[Y]
  555.               FLD ST(2)                // etc...
  556.               FMUL ST,ST(6)
  557.               FADDP ST(1),ST
  558.               WAIT
  559.               FSTP DWORD PTR [ECX+12]  // Result.Angle:=ST(0)
  560.               FFREE ST                 // clear FPU register stack
  561.               FFREE ST(1)
  562.               FFREE ST(2)
  563.               FFREE ST(3)
  564.               FFREE ST(4)
  565. end;
  566.  
  567. //------------------------------------------------------------------------------
  568.  
  569. function QuaternionMultiply(qL, qR: TQuaternion): TQuaternion; register; 
  570.  
  571. // Return quaternion product qL * qR.  Note: order is important!
  572. // To combine rotations, use the product QuaternionMuliply(qSecond, qFirst),
  573. // which gives the effect of rotating by qFirst then qSecond.
  574.  
  575. begin
  576.   Result.Angle:=qL.Angle*qR.Angle-qL.Vector[X]*qR.Vector[X]-
  577.                 qL.Vector[Y]*qR.Vector[Y]-qL.Vector[Z]*qR.Vector[Z];
  578.   Result.Vector[X]:=qL.Angle*qR.Vector[X]+qL.Vector[X]*qR.Angle+
  579.                     qL.Vector[Y]*qR.Vector[Z]-qL.Vector[Z]*qR.Vector[Y];
  580.   Result.Vector[Y]:=qL.Angle*qR.Vector[Y]+qL.Vector[Y]*qR.Angle+
  581.                     qL.Vector[Z]*qR.Vector[X]-qL.Vector[X]*qR.Vector[Z];
  582.   Result.Vector[Z]:=qL.Angle*qR.Vector[Z]+qL.Vector[Z]*qR.Angle+
  583.                     qL.Vector[X]*qR.Vector[Y]-qL.Vector[Y]*qR.Vector[X];
  584. end;
  585.  
  586. //------------------------------------------------------------------------------
  587.  
  588. function QuaternionToMatrix(Q: TQuaternion): TMatrix; register;
  589.  
  590. // Construct rotation matrix from (possibly non-unit) quaternion.
  591. // Assumes matrix is used to multiply column vector on the left:
  592. // vnew = mat vold.  Works correctly for right-handed coordinate system
  593. // and right-handed rotations.
  594.  
  595. var Norm,S,
  596.     XS,YS,ZS,
  597.     WX,WY,WZ,
  598.     XX,XY,XZ,
  599.     YY,YZ,ZZ   : Single;
  600.  
  601. begin
  602.   Norm:=Q.Vector[X]*Q.Vector[X]+Q.Vector[Y]*Q.Vector[Y]+Q.Vector[Z]*Q.Vector[Z]+Q.Angle*Q.Angle;
  603.   if Norm > 0 then S:=2/Norm
  604.               else S:=0;
  605.               
  606.   XS:=Q.Vector[X]*S;  YS:=Q.Vector[Y]*S;  ZS:=Q.Vector[Z]*S;
  607.   WX:=Q.Angle*XS;     WY:=Q.Angle*YS;     WZ:=Q.Angle*ZS;
  608.   XX:=Q.Vector[X]*XS; XY:=Q.Vector[X]*YS; XZ:=Q.Vector[X]*ZS;
  609.   YY:=Q.Vector[Y]*YS; YZ:=Q.Vector[Y]*ZS; ZZ:=Q.Vector[Z]*ZS;
  610.  
  611.   Result[X,X]:=1-(YY+ZZ); Result[Y,X]:=XY+WZ;     Result[Z,X]:=XZ-WY;     Result[W,X]:=0;
  612.   Result[X,Y]:=XY-WZ;     Result[Y,Y]:=1-(XX+ZZ); Result[Z,Y]:=YZ+WX;     Result[W,Y]:=0;
  613.   Result[X,Z]:=XZ+WY;     Result[Y,Z]:=YZ-WX;     Result[Z,Z]:=1-(XX+YY); Result[W,Z]:=0;
  614.   Result[X,W]:=0;         Result[Y,W]:=0;         Result[Z,W]:=0;         Result[W,W]:=1;
  615. end;
  616.  
  617. //------------------------------------------------------------------------------
  618.  
  619. procedure QuaternionToPoints(Q: TQuaternion; var ArcFrom, ArcTo: TVector); register;
  620.  
  621. // convert a unit quaternion into two points on a unit sphere
  622.  
  623. var S: Single;
  624.  
  625. begin
  626.     S:=Sqrt(Q.Vector[X]*Q.Vector[X]+Q.Vector[Y]*Q.Vector[Y]);
  627.     if s = 0 then ArcFrom:=MakeVector([0,1,0,0])
  628.              else ArcFrom:=MakeVector([-Q.Vector[Y]/S,Q.Vector[X]/S,0,0]);
  629.     ArcTo[X]:=Q.Angle*ArcFrom[X]-Q.Vector[Z]*ArcFrom[Y];
  630.     ArcTo[Y]:=Q.Angle*ArcFrom[Y]+Q.Vector[Z]*ArcFrom[X];
  631.     ArcTo[Z]:=Q.Vector[X]*ArcFrom[Y]-Q.Vector[Y]*ArcFrom[X];
  632.     if Q.Angle < 0 then ArcFrom:=MakeVector([-ArcFrom[X],-ArcFrom[Y],0,0]);
  633. end;
  634.  
  635. //------------------------------------------------------------------------------
  636.  
  637. function VectorNorm(V: array of Single): Single; assembler; register;
  638.  
  639. // calculate norm of a vector which is defined as norm=x*x+y*y+...
  640. // EAX contains address of V
  641. // EDX contains highest index in V
  642. // result is passed in ST(0)
  643.  
  644. asm
  645.               FLDZ                          // initialize sum
  646. @@Loop:       FLD  DWORD PTR [EAX+4*EDX]    // load a component
  647.               FMUL ST,ST                    // make square
  648.               FADDP                         // add previous calculated sum
  649.               SUB  EDX,1
  650.               JNL  @@Loop
  651. end;
  652.  
  653. //------------------------------------------------------------------------------
  654.  
  655. function VectorNormalize(V: array of Single): Single; assembler; register;
  656.  
  657. // transform a vector to unit length and return length
  658. // EAX contains address of V
  659. // EDX contains the highest index in V
  660.  
  661. asm
  662.               PUSH EBX
  663.               MOV ECX,EDX                   // save size of V
  664.               CALL VectorLength             // calculate length of vector
  665.               FTST                          // test if length = 0
  666.               MOV EBX,EAX                   // save parameter address
  667.               FSTSW AX                      // get test result
  668.               TEST AH,C3                    // check the test result
  669.               JNZ @@Finish
  670.               SUB EBX,4                     // simplyfied address calculation
  671.               INC ECX
  672.               FLD1                          // calculate reciprocal of length
  673.               FDIV ST,ST(1)
  674. @@1:          FLD ST                        // double reciprocal
  675.               FMUL DWORD PTR [EBX+4*ECX]    // scale component
  676.               WAIT
  677.               FSTP DWORD PTR [EBX+4*ECX]    // store result
  678.               LOOP @@1
  679.               FSTP ST                       // remove reciprocal from FPU stack
  680. @@Finish:     POP EBX
  681. end;
  682.  
  683. //------------------------------------------------------------------------------
  684.  
  685. function VectorAffineSubtract(V1,V2: TAffineVector): TAffineVector; register;
  686.  
  687. // return the difference of v1 minus v2
  688.  
  689. begin
  690.   Result[X]:=V1[X]-V2[X];
  691.   Result[Y]:=V1[Y]-V2[Y];
  692.   Result[Z]:=V1[Z]-V2[Z];
  693. end;
  694.  
  695. //------------------------------------------------------------------------------
  696.  
  697. function VectorReflect(V, N: TAffineVector): TAffineVector; register;
  698.  
  699. // reflect vector V against N (assumes N is normalized)
  700.  
  701. var Dot : Single;
  702.  
  703. begin
  704.    Dot:=VectorAffineDotProduct(V,N);
  705.    Result[X]:=V[X]-2*Dot*N[X];
  706.    Result[Y]:=V[Y]-2*Dot*N[Y];
  707.    Result[Z]:=V[Z]-2*Dot*N[Z];
  708. end;
  709.  
  710. //------------------------------------------------------------------------------
  711.  
  712. procedure VectorScale(V: array of Single; Factor: Single); assembler; register;
  713.  
  714. // return a vector scaled by a factor
  715. // EAX contains address of V
  716. // EDX contains highest index in V
  717. // Factor is located on the stack (don't ask me why not in ECX)
  718.  
  719. asm
  720.   {for I:=Low(V) to High(V) do V[I]:=V[I]*Factor;}
  721.               FLD DWORD PTR [Factor]        // load factor
  722. @@Loop:       FLD DWORD PTR [EAX+4*EDX]     // load a component
  723.               FMUL ST,ST(1)                 // multiply it with the factor
  724.               WAIT
  725.               FSTP DWORD PTR [EAX+4*EDX]    // store the result
  726.               DEC EDX                       // do the entire array
  727.               JNS @@Loop
  728.               FSTP ST(0)                    // clean the FPU stack
  729. end;
  730.  
  731. //------------------------------------------------------------------------------
  732.  
  733. procedure VectorNegate(V: array of Single); assembler; register;
  734.  
  735. // return a negated vector
  736. // EAX contains address of V
  737. // EDX contains highest index in V
  738.  
  739. asm
  740.   {V[X]:=-V[X];
  741.   V[Y]:=-V[Y];
  742.   V[Z]:=-V[Z];}
  743. @@Loop:       FLD DWORD PTR [EAX+4*EDX]
  744.               FCHS
  745.               WAIT
  746.               FSTP DWORD PTR [EAX+4*EDX]
  747.               DEC EDX
  748.               JNS @@Loop
  749. end;
  750.  
  751. //------------------------------------------------------------------------------
  752.  
  753. function VectorAdd(V1,V2: TVector): TVector; register;
  754.  
  755. // return the sum of two vectors
  756.  
  757. begin
  758.   Result[X]:=V1[X]+V2[X];
  759.   Result[Y]:=V1[Y]+V2[Y];
  760.   Result[Z]:=V1[Z]+V2[Z];
  761.   Result[W]:=V1[W]+V2[W];
  762. end;
  763.  
  764. //------------------------------------------------------------------------------
  765.  
  766. function VectorAffineAdd(V1,V2: TAffineVector): TAffineVector; register;
  767.  
  768. // return the sum of two vectors
  769.  
  770. begin
  771.   Result[X]:=V1[X]+V2[X];
  772.   Result[Y]:=V1[Y]+V2[Y];
  773.   Result[Z]:=V1[Z]+V2[Z];
  774. end;
  775.  
  776. //------------------------------------------------------------------------------
  777.  
  778. function VectorSubtract(V1,V2: TVector): TVector; register;
  779.  
  780. // return the difference of two vectors
  781.  
  782. begin
  783.   Result[X]:=V1[X]-V2[X];
  784.   Result[Y]:=V1[Y]-V2[Y];
  785.   Result[Z]:=V1[Z]-V2[Z];
  786.   Result[W]:=V1[W]-V2[W];
  787. end;
  788.  
  789. //------------------------------------------------------------------------------
  790.  
  791. function VectorDotProduct(V1,V2: TVector): Single; register;
  792.  
  793. begin
  794.   Result:=V1[X]*V2[X]+V1[Y]*V2[Y]+V1[Z]*V2[Z]+V1[W]*V2[W];
  795. end;
  796.  
  797. //------------------------------------------------------------------------------
  798.  
  799. function VectorAffineDotProduct(V1,V2: TAffineVector): Single; register;
  800.  
  801. begin
  802.   Result:=V1[X]*V2[X]+V1[Y]*V2[Y]+V1[Z]*V2[Z];
  803. end;
  804.  
  805. //------------------------------------------------------------------------------
  806.  
  807. function VectorCrossProduct(V1,V2: TAffineVector): TAffineVector; register;
  808.  
  809. begin
  810.   Result[X]:=V1[Y]*V2[Z]-V1[Z]*V2[Y];
  811.   Result[Y]:=V1[Z]*V2[X]-V1[X]*V2[Z];
  812.   Result[Z]:=V1[X]*V2[Y]-V1[Y]*V2[X];
  813. end;
  814.  
  815. //------------------------------------------------------------------------------
  816.  
  817. function VectorPerpendicular(V, N: TAffineVector): TAffineVector; register; 
  818.  
  819. // calculate a vector perpendicular to N (N is assumed to be of unit length)
  820. // subtract out any component parallel to N
  821.  
  822. var Dot : Single;
  823.  
  824. begin
  825.    Dot:=VectorAffineDotProduct(V,N);
  826.    Result[X]:=V[X]-Dot*N[X];
  827.    Result[Y]:=V[Y]-Dot*N[Y];
  828.    Result[Z]:=V[Z]-Dot*N[Z];
  829. end;
  830.  
  831. //------------------------------------------------------------------------------
  832.  
  833. function VectorTransform(V: TVector; M: TMatrix): TVector; register;
  834.  
  835. // transform a homogeneous vector by multiplying it with a matrix
  836.  
  837. var TV : TVector;
  838.  
  839. begin
  840.   TV[X]:=V[X]*M[X,X]+V[Y]*M[Y,X]+V[Z]*M[Z,X]+V[W]*M[W,X];
  841.   TV[Y]:=V[X]*M[X,Y]+V[Y]*M[Y,Y]+V[Z]*M[Z,Y]+V[W]*M[W,Y];
  842.   TV[Z]:=V[X]*M[X,Z]+V[Y]*M[Y,Z]+V[Z]*M[Z,Z]+V[W]*M[W,Z];
  843.   TV[W]:=V[X]*M[X,W]+V[Y]*M[Y,W]+V[Z]*M[Z,W]+V[W]*M[W,W];
  844.   Result:=TV
  845. end;
  846.  
  847. //------------------------------------------------------------------------------
  848.  
  849. function VectorAffineTransform(V: TAffineVector; M: TAffineMatrix): TAffineVector; register;
  850.  
  851. // transform an affine vector by multiplying it with a matrix
  852.  
  853. var TV : TAffineVector;
  854.  
  855. begin
  856.   TV[X]:=V[X]*M[X,X]+V[Y]*M[Y,X]+V[Z]*M[Z,X];
  857.   TV[Y]:=V[X]*M[X,Y]+V[Y]*M[Y,Y]+V[Z]*M[Z,Y];
  858.   TV[Z]:=V[X]*M[X,Z]+V[Y]*M[Y,Z]+V[Z]*M[Z,Z];
  859.   Result:=TV;
  860. end;
  861.  
  862. //------------------------------------------------------------------------------
  863.  
  864. function MatrixAffineDeterminant(M: TAffineMatrix): Single; register;
  865.  
  866. // determinant of a 3x3 matrix
  867.  
  868. begin
  869.   Result:=M[X,X]*(M[Y,Y]*M[Z,Z]-M[Z,Y]*M[Y,Z])-
  870.           M[X,Y]*(M[Y,X]*M[Z,Z]-M[Z,X]*M[Y,Z])+
  871.       M[X,Z]*(M[Y,X]*M[Z,Y]-M[Z,X]*M[Y,Y]);
  872. end;
  873.  
  874. //------------------------------------------------------------------------------
  875.  
  876. function MatrixDetInternal(a1,a2,a3,b1,b2,b3,c1,c2,c3: Single): Single;
  877.  
  878. // internal version for the determinant of a 3x3 matrix
  879.  
  880. begin
  881.   Result:= a1 * (b2 * c3 - b3 * c2) -
  882.            b1 * (a2 * c3 - a3 * c2) +
  883.            c1 * (a2 * b3 - a3 * b2);
  884. end;
  885.  
  886. //------------------------------------------------------------------------------
  887.  
  888. procedure MatrixAdjoint(var M: TMatrix); register;
  889.  
  890. // Adjoint of a 4x4 matrix - used in the computation of the inverse
  891. // of a 4x4 matrix
  892.  
  893. var a1,a2,a3,a4,
  894.     b1,b2,b3,b4,
  895.     c1,c2,c3,c4,
  896.     d1,d2,d3,d4  : Single;
  897.  
  898.  
  899. begin
  900.     a1:= M[0,0]; b1:= M[0,1];
  901.     c1:= M[0,2]; d1:= M[0,3];
  902.     a2:= M[1,0]; b2:= M[1,1];
  903.     c2:= M[1,2]; d2:= M[1,3];
  904.     a3:= M[2,0]; b3:= M[2,1];
  905.     c3:= M[2,2]; d3:= M[2,3];
  906.     a4:= M[3,0]; b4:= M[3,1];
  907.     c4:= M[3,2]; d4:= M[3,3];
  908.  
  909.     // row column labeling reversed since we transpose rows & columns
  910.     M[X,X]:= MatrixDetInternal(b2,b3,b4,c2,c3,c4,d2,d3,d4);
  911.     M[Y,X]:=-MatrixDetInternal(a2,a3,a4,c2,c3,c4,d2,d3,d4);
  912.     M[Z,X]:= MatrixDetInternal(a2,a3,a4,b2,b3,b4,d2,d3,d4);
  913.     M[W,X]:=-MatrixDetInternal(a2,a3,a4,b2,b3,b4,c2,c3,c4);
  914.  
  915.     M[X,Y]:=-MatrixDetInternal(b1,b3,b4,c1,c3,c4,d1,d3,d4);
  916.     M[Y,Y]:= MatrixDetInternal(a1,a3,a4,c1,c3,c4,d1,d3,d4);
  917.     M[Z,Y]:=-MatrixDetInternal(a1,a3,a4,b1,b3,b4,d1,d3,d4);
  918.     M[W,Y]:= MatrixDetInternal(a1,a3,a4,b1,b3,b4,c1,c3,c4);
  919.  
  920.     M[X,Z]:= MatrixDetInternal(b1,b2,b4,c1,c2,c4,d1,d2,d4);
  921.     M[Y,Z]:=-MatrixDetInternal(a1,a2,a4,c1,c2,c4,d1,d2,d4);
  922.     M[Z,Z]:= MatrixDetInternal(a1,a2,a4,b1,b2,b4,d1,d2,d4);
  923.     M[W,Z]:=-MatrixDetInternal(a1,a2,a4,b1,b2,b4,c1,c2,c4);
  924.  
  925.     M[X,W]:=-MatrixDetInternal(b1,b2,b3,c1,c2,c3,d1,d2,d3);
  926.     M[Y,W]:= MatrixDetInternal(a1,a2,a3,c1,c2,c3,d1,d2,d3);
  927.     M[Z,W]:=-MatrixDetInternal(a1,a2,a3,b1,b2,b3,d1,d2,d3);
  928.     M[W,W]:= MatrixDetInternal(a1,a2,a3,b1,b2,b3,c1,c2,c3);
  929. end;
  930.  
  931. //------------------------------------------------------------------------------
  932.  
  933. function MatrixDeterminant(M: TMatrix): Single; register;
  934.  
  935. // Determinant of a 4x4 matrix
  936.  
  937. var a1, a2, a3, a4,
  938.     b1, b2, b3, b4,
  939.     c1, c2, c3, c4,
  940.     d1, d2, d3, d4  : Single;
  941.  
  942. begin
  943.   a1:=M[X,X]; b1:=M[X,Y]; c1:=M[X,Z]; d1:=M[X,W];
  944.   a2:=M[Y,X]; b2:=M[Y,Y]; c2:=M[Y,Z]; d2:=M[Y,W];
  945.   a3:=M[Z,X]; b3:=M[Z,Y]; c3:=M[Z,Z]; d3:=M[Z,W];
  946.   a4:=M[W,X]; b4:=M[W,Y]; c4:=M[W,Z]; d4:=M[W,W];
  947.  
  948.   Result:=a1*MatrixDetInternal(b2,b3,b4,c2,c3,c4,d2,d3,d4)-
  949.           b1*MatrixDetInternal(a2,a3,a4,c2,c3,c4,d2,d3,d4)+
  950.           c1*MatrixDetInternal(a2,a3,a4,b2,b3,b4,d2,d3,d4)-
  951.           d1*MatrixDetInternal(a2,a3,a4,b2,b3,b4,c2,c3,c4);
  952. end;
  953.  
  954. //------------------------------------------------------------------------------
  955.  
  956. procedure MatrixScale(var M: TMatrix; Factor: Single); register;
  957.  
  958. // multiply all elements of a 4x4 matrix with a factor
  959.  
  960. var I,J : Integer;
  961.  
  962. begin
  963.   for I:=0 to 3 do
  964.     for J:=0 to 3 do M[I,J]:=M[I,J]*Factor;
  965. end;
  966.  
  967. //------------------------------------------------------------------------------
  968.  
  969. procedure MatrixInvert(var M: TMatrix); register;
  970.  
  971. // find the inverse of a 4x4 matrix
  972.  
  973. var Det : Single;
  974.  
  975. begin
  976.   Det:=MatrixDeterminant(M);
  977.   if Abs(Det) < EPSILON then M:=IdentityMatrix
  978.                         else
  979.   begin
  980.     MatrixAdjoint(M);
  981.     MatrixScale(M,1/Det);
  982.   end;
  983. end;
  984.  
  985. //------------------------------------------------------------------------------
  986.  
  987. procedure MatrixTranspose(var M: TMatrix); register;
  988.  
  989. // compute transpose of 4x4 matrix
  990.  
  991. var I,J : Integer;
  992.     TM  : TMatrix;
  993.  
  994. begin
  995.   for I:=0 to 3 do
  996.     for J:=0 to 3 do TM[J,I]:=M[I,J];
  997.   M:=TM;
  998. end;
  999.  
  1000. //------------------------------------------------------------------------------
  1001.  
  1002. procedure MatrixAffineTranspose(var M: TAffineMatrix); register;
  1003.  
  1004. // compute transpose of 3x3 matrix
  1005.  
  1006. var I,J : Integer;
  1007.     TM  : TAffineMatrix;
  1008.  
  1009. begin
  1010.   for I:=0 to 2 do
  1011.     for J:=0 to 2 do TM[J,I]:=M[I,J];
  1012.   M:=TM;
  1013. end;
  1014.  
  1015. //------------------------------------------------------------------------------
  1016.  
  1017. function MatrixMultiply(M1, M2: TMatrix): TMatrix; register;
  1018.  
  1019. // multiply two 4x4 matrices
  1020.  
  1021. var I,J : Integer;
  1022.     TM  : TMatrix;
  1023.  
  1024. begin
  1025.   for I:=0 to 3 do
  1026.     for J:=0 to 3 do
  1027.       TM[I,J]:=M1[I,X]*M2[X,J]+
  1028.                M1[I,Y]*M2[Y,J]+
  1029.                M1[I,Z]*M2[Z,J]+
  1030.                M1[I,W]*M2[W,J];
  1031.   Result:=TM;
  1032. end;
  1033.  
  1034. //------------------------------------------------------------------------------
  1035.  
  1036. function CreateRotationMatrix(Axis: TAffineVector; Angle: Single): TMatrix; register; 
  1037.  
  1038. // Create a rotation matrix along the given axis by the given angle in radians.
  1039. // The axis is a set of direction cosines.
  1040.  
  1041. var cosine, sine,
  1042.     one_minus_cosine : Extended;
  1043.  
  1044. begin
  1045.     SinCos(Angle,Sine,Cosine);
  1046.     one_minus_cosine:=1-cosine;
  1047.  
  1048.     Result[X,X]:=SQR(Axis[X])+(1-SQR(axis[X]))*cosine;
  1049.     Result[X,Y]:=Axis[X]*Axis[Y]*one_minus_cosine+Axis[Z]*sine;
  1050.     Result[X,Z]:=Axis[X]*Axis[Z]*one_minus_cosine-Axis[Y]*sine;
  1051.     Result[X,W]:=0;
  1052.  
  1053.     Result[Y,X]:=Axis[X]*Axis[Y]*one_minus_cosine-Axis[Z]*sine;
  1054.     Result[Y,Y]:=SQR(Axis[Y])+(1-SQR(axis[Y]))*cosine;
  1055.     Result[Y,Z]:=Axis[Y]*Axis[Z]*one_minus_cosine+Axis[X]*sine;
  1056.     Result[Y,W]:=0;
  1057.  
  1058.     Result[Z,X]:=Axis[X]*Axis[Z]*one_minus_cosine+Axis[Y]*sine;
  1059.     Result[Z,Y]:=Axis[Y]*Axis[Z]*one_minus_cosine-Axis[X]*sine;
  1060.     Result[Z,Z]:=SQR(Axis[Z])+(1-SQR(axis[Z]))*cosine;
  1061.     Result[Z,W]:=0;
  1062.  
  1063.     Result[W,X]:=0;
  1064.     Result[W,Y]:=0;
  1065.     Result[W,Z]:=0;
  1066.     Result[W,W]:=1;
  1067. end;
  1068.  
  1069. //------------------------------------------------------------------------------
  1070.  
  1071. function ConvertRotation(Angles: TAffineVector): TVector; register;
  1072.  
  1073. { Turn a triplet of rotations about x, y, and z (in that order) into an
  1074.    equivalent rotation around a single axis (all in radians).
  1075.  
  1076.    Rotation of the angle t about the axis (X,Y,Z) is given by:
  1077.  
  1078.      | X^2+(1-X^2) Cos(t),    XY(1-Cos(t)) + Z Sin(t), XZ(1-Cos(t))-Y Sin(t) |
  1079.  M = | XY(1-Cos(t))-Z Sin(t), Y^2+(1-Y^2) Cos(t),      YZ(1-Cos(t))+X Sin(t) |
  1080.      | XZ(1-Cos(t))+Y Sin(t), YZ(1-Cos(t))-X Sin(t),   Z^2+(1-Z^2) Cos(t)    |
  1081.  
  1082.    Rotation about the three axes (angles a1, a2, a3) can be represented as
  1083.    the product of the individual rotation matrices:
  1084.  
  1085.       | 1  0       0       | | Cos(a2) 0 -Sin(a2) | |  Cos(a3) Sin(a3) 0 |
  1086.       | 0  Cos(a1) Sin(a1) |*| 0       1  0       |*| -Sin(a3) Cos(a3) 0 |
  1087.       | 0 -Sin(a1) Cos(a1) | | Sin(a2) 0  Cos(a2) | |  0       0       1 |
  1088.          Mx                       My                     Mz
  1089.  
  1090.    We now want to solve for X, Y, Z, and t given 9 equations in 4 unknowns.
  1091.    Using the diagonal elements of the two matrices, we get:
  1092.  
  1093.       X^2+(1-X^2) Cos(t) = M[0][0]
  1094.       Y^2+(1-Y^2) Cos(t) = M[1][1]
  1095.       Z^2+(1-Z^2) Cos(t) = M[2][2]
  1096.  
  1097.    Adding the three equations, we get:
  1098.  
  1099.       X^2 + Y^2 + Z^2 - (M[0][0] + M[1][1] + M[2][2]) =
  1100.      - (3 - X^2 - Y^2 - Z^2) Cos(t)
  1101.  
  1102.    Since (X^2 + Y^2 + Z^2) = 1, we can rewrite as:
  1103.  
  1104.       Cos(t) = (1 - (M[0][0] + M[1][1] + M[2][2])) / 2
  1105.  
  1106.    Solving for t, we get:
  1107.  
  1108.       t = Acos(((M[0][0] + M[1][1] + M[2][2]) - 1) / 2)
  1109.  
  1110.     We can substitute t into the equations for X^2, Y^2, and Z^2 above
  1111.     to get the values for X, Y, and Z.  To find the proper signs we note
  1112.     that:
  1113.  
  1114.     2 X Sin(t) = M[1][2] - M[2][1]
  1115.     2 Y Sin(t) = M[2][0] - M[0][2]
  1116.     2 Z Sin(t) = M[0][1] - M[1][0]
  1117. }
  1118.  
  1119. var Axis1, Axis2 : TAffineVector;
  1120.     M, M1, M2    : TMatrix;
  1121.     cost, cost1,
  1122.     sint,
  1123.     s1, s2, s3   : Single;
  1124.     I            : Integer;
  1125.  
  1126.  
  1127. begin
  1128.   // see if we are only rotating about a single axis 
  1129.   if Abs(Angles[X]) < EPSILON then
  1130.   begin
  1131.     if Abs(Angles[Y]) < EPSILON then
  1132.     begin
  1133.       Result:=MakeVector([0,0,1,Angles[Z]]);
  1134.       Exit;
  1135.     end
  1136.     else
  1137.       if Abs(Angles[Z]) < EPSILON then
  1138.       begin
  1139.         Result:=MakeVector([0,1,0,Angles[Y]]);
  1140.         Exit;
  1141.       end
  1142.    end
  1143.    else
  1144.      if (Abs(Angles[Y]) < EPSILON) and
  1145.         (Abs(Angles[Z]) < EPSILON) then
  1146.      begin
  1147.        Result:=MakeVector([1,0,0,Angles[X]]);
  1148.        Exit;
  1149.      end;
  1150.  
  1151.   // make the rotation matrix
  1152.   Axis1:=MakeAffineVector([1,0,0]);
  1153.   M:=CreateRotationMatrix(Axis1,Angles[X]);
  1154.  
  1155.   Axis2:=MakeAffineVector([0,1,0]);
  1156.   M2:=CreateRotationMatrix(Axis2,Angles[Y]);
  1157.   M1:=MatrixMultiply(M,M2);
  1158.  
  1159.   Axis2:=MakeAffineVector([0,0,1]);
  1160.   M2:=CreateRotationMatrix(Axis2,Angles[Z]);
  1161.   M:=MatrixMultiply(M1,M2);
  1162.  
  1163.   cost:=((M[X,X]+M[Y,Y]+M[Z,Z])-1)/2;
  1164.   if cost < -1 then cost:=-1
  1165.                else
  1166.     if cost > 1-EPSILON then
  1167.     begin
  1168.       // Bad angle - this would cause a crash
  1169.       Result:=MakeVector([1,0,0,0]);
  1170.       Exit;
  1171.     end;
  1172.  
  1173.   cost1:=1-cost;
  1174.   Result:=Makevector([Sqrt((M[X,X]-cost)/cost1),
  1175.                       Sqrt((M[Y,Y]-cost)/cost1),
  1176.                       sqrt((M[Z,Z]-cost)/cost1),
  1177.                       arccos(cost)]);
  1178.  
  1179.   sint:=2*Sqrt(1-cost*cost); // This is actually 2 Sin(t)
  1180.  
  1181.   // Determine the proper signs
  1182.   for I:=0 to 7 do
  1183.   begin
  1184.     if (I and 1) > 1 then s1:=-1 else s1:=1;
  1185.     if (I and 2) > 1 then s2:=-1 else s2:=1;
  1186.     if (I and 4) > 1 then s3:=-1 else s3:=1;
  1187.     if (Abs(s1*Result[X]*sint-M[Y,Z]+M[Z,Y]) < EPSILON2) and
  1188.        (Abs(s2*Result[Y]*sint-M[Z,X]+M[X,Z]) < EPSILON2) and
  1189.        (Abs(s3*Result[Z]*sint-M[X,Y]+M[Y,X]) < EPSILON2) then
  1190.         begin
  1191.           // We found the right combination of signs
  1192.           Result[X]:=Result[X]*s1;
  1193.           Result[Y]:=Result[Y]*s2;
  1194.           Result[Z]:=Result[Z]*s3;
  1195.           Exit;
  1196.         end;
  1197.   end;
  1198. end;
  1199.  
  1200. //------------------------------------------------------------------------------
  1201.  
  1202. function CreateRotationMatrixX(Sine, Cosine: Single): TMatrix; register;
  1203.  
  1204. // create matrix for rotation about x-axis
  1205.  
  1206. begin
  1207.   Result:=EmptyMatrix;
  1208.   Result[0,0]:=1;
  1209.   Result[1,1]:=Cosine;
  1210.   Result[1,2]:=Sine;
  1211.   Result[2,1]:=-Sine;
  1212.   Result[2,2]:=Cosine;
  1213.   Result[3,3]:=1;
  1214. end;
  1215.  
  1216. //------------------------------------------------------------------------------
  1217.  
  1218. function CreateRotationMatrixY(Sine, Cosine: Single): TMatrix; register;
  1219.  
  1220. // create matrix for rotation about y-axis
  1221.  
  1222. begin
  1223.   Result:=EmptyMatrix;
  1224.   Result[0,0]:=Cosine;
  1225.   Result[0,2]:=-Sine;
  1226.   Result[1,1]:=1;
  1227.   Result[2,0]:=Sine;
  1228.   Result[2,2]:=Cosine;
  1229.   Result[3,3]:=1;
  1230. end;
  1231.  
  1232. //------------------------------------------------------------------------------
  1233.  
  1234. function CreateRotationMatrixZ(Sine, Cosine: Single): TMatrix; register;
  1235.  
  1236. // create matrix for rotation about z-axis
  1237.  
  1238. begin
  1239.   Result:=EmptyMatrix;
  1240.   Result[0,0]:=Cosine;
  1241.   Result[0,1]:=Sine;
  1242.   Result[1,0]:=-Sine;
  1243.   Result[1,1]:=Cosine;
  1244.   Result[2,2]:=1;
  1245.   Result[3,3]:=1;
  1246. end;
  1247.  
  1248. //------------------------------------------------------------------------------
  1249.  
  1250. function CreateScaleMatrix(V: TAffineVector): TMatrix; register;
  1251.  
  1252. // create scaling matrix
  1253.  
  1254. begin
  1255.   Result:=EmptyMatrix;
  1256.   Result[X,X]:=V[X];
  1257.   Result[Y,Y]:=V[Y];
  1258.   Result[Z,Z]:=V[Z];
  1259.   Result[W,W]:=1;
  1260. end;
  1261.  
  1262. //------------------------------------------------------------------------------
  1263.  
  1264. function CreateTranslationMatrix(V: TAffineVector): TMatrix; register;
  1265.  
  1266. // create translation matrix
  1267.  
  1268. begin
  1269.   Result:=IdentityMatrix;
  1270.   Result[W,X]:=V[X];
  1271.   Result[W,Y]:=V[Y];
  1272.   Result[W,Z]:=V[Z];
  1273. end;
  1274.  
  1275. //------------------------------------------------------------------------------
  1276.  
  1277. function Lerp(Start,Stop,t: Single): Single;
  1278.  
  1279. // calculate linear interpolation between start and stop at point t
  1280.  
  1281. begin
  1282.   Result:=Start+(Stop-Start)*t;
  1283. end;
  1284.  
  1285. //------------------------------------------------------------------------------
  1286.  
  1287. function VectorAffineLerp(V1,V2: TAffineVector; t: Single): TAffineVector;
  1288.  
  1289. // calculate linear interpolation between vector1 and vector2 at point t
  1290.  
  1291. begin
  1292.   Result[X]:=Lerp(V1[X],V2[X],t);
  1293.   Result[Y]:=Lerp(V1[Y],V2[Y],t);
  1294.   Result[Z]:=Lerp(V1[Z],V2[Z],t);
  1295. end;
  1296.  
  1297. //------------------------------------------------------------------------------
  1298.  
  1299. function VectorLerp(V1,V2: TVector; t: Single): TVector;
  1300.  
  1301. // calculate linear interpolation between vector1 and vector2 at point t
  1302.  
  1303. begin
  1304.   Result[X]:=Lerp(V1[X],V2[X],t);
  1305.   Result[Y]:=Lerp(V1[Y],V2[Y],t);
  1306.   Result[Z]:=Lerp(V1[Z],V2[Z],t);
  1307.   Result[W]:=Lerp(V1[W],V2[W],t);
  1308. end;
  1309.  
  1310. //------------------------------------------------------------------------------
  1311.  
  1312. function QuaternionSlerp(QStart,QEnd: TQuaternion; Spin: Integer; t: Single): TQuaternion;
  1313.  
  1314. // spherical linear interpolation of unit quaternions with spins
  1315. // QStart, QEnd - start and end unit quaternions
  1316. // t            - interpolation parameter (0 to 1)
  1317. // Spin         - number of extra spin rotations to involve
  1318.  
  1319. var beta,                   // complementary interp parameter
  1320.     theta,                  // angle between A and B
  1321.     sint, cost,             // sine, cosine of theta
  1322.     phi         : Single;   // theta plus spins
  1323.     bflip       : Boolean;  // use negation of B?
  1324.  
  1325.  
  1326. begin
  1327.   // cosine theta
  1328.   cost:=VectorAngle(QStart.Axis,QEnd.Axis);
  1329.  
  1330.   // if QEnd is on opposite hemisphere from QStart, use -QEnd instead
  1331.   if cost < 0 then
  1332.   begin
  1333.     cost:=-cost;
  1334.     bflip:=True;
  1335.   end
  1336.   else bflip:=False;
  1337.  
  1338.   // if QEnd is (within precision limits) the same as QStart,
  1339.   // just linear interpolate between QStart and QEnd.
  1340.   // Can't do spins, since we don't know what direction to spin.
  1341.  
  1342.   if (1-cost) < EPSILON then beta:=1-t
  1343.                     else
  1344.   begin
  1345.     // normal case
  1346.     theta:=arccos(cost);
  1347.     phi:=theta+Spin*Pi;
  1348.     sint:=sin(theta);
  1349.     beta:=sin(theta-t*phi)/sint;
  1350.     t:=sin(t*phi)/sint;
  1351.   end;
  1352.  
  1353.   if bflip then t:=-t;
  1354.  
  1355.   // interpolate
  1356.   Result.Vector[X]:=beta*QStart.Vector[X]+t*QEnd.Vector[X];
  1357.   Result.Vector[Y]:=beta*QStart.Vector[Y]+t*QEnd.Vector[Y];
  1358.   Result.Vector[Z]:=beta*QStart.Vector[Z]+t*QEnd.Vector[Z];
  1359.   Result.Angle:=beta*QStart.Angle+t*QEnd.Angle;
  1360. end;
  1361.  
  1362. //------------------------------------------------------------------------------
  1363.  
  1364. function VectorAffineCombine(V1,V2: TAffineVector; F1,F2: Single): TAffineVector;
  1365.  
  1366. // make a linear combination of two vectors and return the result
  1367.  
  1368. begin
  1369.   Result[X]:=(F1*V1[X])+(F2*V2[X]);
  1370.   Result[Y]:=(F1*V1[Y])+(F2*V2[Y]);
  1371.   Result[Z]:=(F1*V1[Z])+(F2*V2[Z]);
  1372. end;
  1373.  
  1374. //------------------------------------------------------------------------------
  1375.  
  1376. function VectorCombine(V1,V2: TVector; F1,F2: Single): TVector;
  1377.  
  1378. // make a linear combination of two vectors and return the result
  1379.  
  1380. begin
  1381.   Result[X]:=(F1*V1[X])+(F2*V2[X]);
  1382.   Result[Y]:=(F1*V1[Y])+(F2*V2[Y]);
  1383.   Result[Z]:=(F1*V1[Z])+(F2*V2[Z]);
  1384.   Result[W]:=(F1*V1[W])+(F2*V2[W]);
  1385. end;
  1386.  
  1387. //------------------------------------------------------------------------------
  1388.  
  1389. function MatrixDecompose(M: TMatrix; var Tran: TTransformations): Boolean; register;
  1390.  
  1391. // Author: Spencer W. Thomas, University of Michigan
  1392. //
  1393. // MatrixDecompose - Decompose a non-degenerate 4x4 transformation matrix into
  1394. // the sequence of transformations that produced it.
  1395. //
  1396. // The coefficient of each transformation is returned in the corresponding
  1397. // element of the vector Tran.
  1398. //
  1399. // Returns true upon success, false if the matrix is singular.
  1400.  
  1401. var I,J      : Integer;
  1402.     LocMat,
  1403.     pmat,
  1404.     invpmat,
  1405.     tinvpmat : TMatrix;
  1406.     prhs,
  1407.     psol     : TVector;
  1408.     Row      : ARRAY[0..2] OF TAffineVector;
  1409.  
  1410. begin
  1411.   Result:=False;
  1412.   locmat:=M;
  1413.   // normalize the matrix
  1414.   if locmat[W,W] = 0 then Exit;
  1415.   for I:=0 to 3 do
  1416.     for J:=0 to 3 do
  1417.       locmat[I,J]:=locmat[I,J]/locmat[W,W];
  1418.  
  1419.   // pmat is used to solve for perspective, but it also provides
  1420.   // an easy way to test for singularity of the upper 3x3 component.
  1421.  
  1422.   pmat:=locmat;
  1423.   for I:=0 to 2 do pmat[I,W]:=0;
  1424.   pmat[W,W]:=1;
  1425.  
  1426.   if MatrixDeterminant(pmat) = 0 then Exit;
  1427.  
  1428.   // First, isolate perspective.  This is the messiest.
  1429.   if (locmat[X,W] <> 0) or
  1430.      (locmat[Y,W] <> 0) or
  1431.      (locmat[Z,W] <> 0) then
  1432.   begin
  1433.     // prhs is the right hand side of the equation.
  1434.     prhs[X]:=locmat[X,W];
  1435.     prhs[Y]:=locmat[Y,W];
  1436.     prhs[Z]:=locmat[Z,W];
  1437.     prhs[W]:=locmat[W,W];
  1438.  
  1439.     // Solve the equation by inverting pmat and multiplying
  1440.     // prhs by the inverse.  (This is the easiest way, not
  1441.     // necessarily the best.)
  1442.  
  1443.     invpmat:=pmat;
  1444.     MatrixInvert(invpmat);
  1445.     MatrixTranspose(invpmat);
  1446.     psol:=VectorTransform(prhs,tinvpmat);
  1447.  
  1448.     // stuff the answer away
  1449.     Tran[ttPerspectiveX]:=psol[X];
  1450.     Tran[ttPerspectiveY]:=psol[Y];
  1451.     Tran[ttPerspectiveZ]:=psol[Z];
  1452.     Tran[ttPerspectiveW]:=psol[W];
  1453.  
  1454.     // clear the perspective partition
  1455.     locmat[X,W]:=0;
  1456.     locmat[Y,W]:=0;
  1457.     locmat[Z,W]:=0;
  1458.     locmat[W,W]:=1;
  1459.   end
  1460.   else
  1461.   begin
  1462.     // no perspective
  1463.     Tran[ttPerspectiveX]:=0;
  1464.     Tran[ttPerspectiveY]:=0;
  1465.     Tran[ttPerspectiveZ]:=0;
  1466.     Tran[ttPerspectiveW]:=0;
  1467.   end;
  1468.  
  1469.   // next take care of translation (easy)
  1470.   for I:=0 to 2 do
  1471.   begin
  1472.     Tran[TTransType(Ord(ttTranslateX)+I)]:=locmat[W,I];
  1473.     locmat[W,I]:=0;
  1474.   end;
  1475.  
  1476.   // now get scale and shear
  1477.   for I:=0 to 2 do
  1478.   begin
  1479.     row[I,X]:=locmat[I,X];
  1480.     row[I,Y]:=locmat[I,Y];
  1481.     row[I,Z]:=locmat[I,Z];
  1482.   end;
  1483.  
  1484.   // compute X scale factor and normalize first row
  1485.   Tran[ttScaleX]:=Sqr(VectorNormalize(row[0])); // (ML) calc. optimized
  1486.  
  1487.   // compute XY shear factor and make 2nd row orthogonal to 1st
  1488.   Tran[ttShearXY]:=VectorAffineDotProduct(row[0],row[1]);
  1489.   row[1]:=VectorAffineCombine(row[1],row[0],1,-Tran[ttShearXY]);
  1490.  
  1491.   // now, compute Y scale and normalize 2nd row
  1492.   Tran[ttScaleY]:=Sqr(VectorNormalize(row[1])); // (ML) calc. optimized
  1493.   Tran[ttShearXY]:=Tran[ttShearXY]/Tran[ttScaleY];
  1494.  
  1495.   // compute XZ and YZ shears, orthogonalize 3rd row
  1496.   Tran[ttShearXZ]:=VectorAffineDotProduct(row[0],row[2]);
  1497.   row[2]:=VectorAffineCombine(row[2],row[0],1,-Tran[ttShearXZ]);
  1498.   Tran[ttShearYZ]:=VectorAffineDotProduct(row[1],row[2]);
  1499.   row[2]:=VectorAffineCombine(row[2],row[1],1,-Tran[ttShearYZ]);
  1500.  
  1501.   // next, get Z scale and normalize 3rd row
  1502.   Tran[ttScaleZ]:=Sqr(VectorNormalize(row[1])); // (ML) calc. optimized
  1503.   Tran[ttShearXZ]:=Tran[ttShearXZ]/tran[ttScaleZ];
  1504.   Tran[ttShearYZ]:=Tran[ttShearYZ]/Tran[ttScaleZ];
  1505.  
  1506.   // At this point, the matrix (in rows[]) is orthonormal.
  1507.   // Check for a coordinate system flip.  If the determinant
  1508.   // is -1, then negate the matrix and the scaling factors.
  1509.   if VectorAffineDotProduct(row[0],VectorCrossProduct(row[1],row[2])) < 0 then
  1510.     for I:=0 to 2 do
  1511.     begin
  1512.       Tran[TTransType(Ord(ttScaleX)+I)]:=-Tran[TTransType(Ord(ttScaleX)+I)];
  1513.       row[I,X]:=-row[I,X];
  1514.       row[I,Y]:=-row[I,Y];
  1515.       row[I,Z]:=-row[I,Z];
  1516.     end;
  1517.  
  1518.   // now, get the rotations out, as described in the gem
  1519.   Tran[ttRotateY]:=arcsin(-row[0,Z]);
  1520.   if cos(Tran[ttRotateY]) <> 0 then
  1521.   begin
  1522.     Tran[ttRotateX]:=arctan2(row[1,Z],row[2,Z]);
  1523.     Tran[ttRotateZ]:=arctan2(row[0,Y],row[0,X]);
  1524.   end
  1525.   else
  1526.   begin
  1527.     tran[ttRotateX]:=arctan2(row[1,X],row[1,Y]);
  1528.     tran[ttRotateZ]:=0;
  1529.   end;
  1530.   // All done!
  1531.   Result:=True;
  1532. end;
  1533.  
  1534. //------------------------------------------------------------------------------
  1535.  
  1536. function VectorDblToFlt(V: THomogeneousDblVector): THomogeneousVector; register; assembler;
  1537.  
  1538. asm
  1539.               FLD  QWORD PTR [EAX]
  1540.               FSTP DWORD PTR [EDX]
  1541.               FLD  QWORD PTR [EAX+8]
  1542.               FSTP DWORD PTR [EDX+4]
  1543.               FLD  QWORD PTR [EAX+16]
  1544.               FSTP DWORD PTR [EDX+8]
  1545.               FLD  QWORD PTR [EAX+24]
  1546.               FSTP DWORD PTR [EDX+12]
  1547. end;
  1548.  
  1549. //------------------------------------------------------------------------------
  1550.  
  1551. function VectorAffineDblToFlt(V: TAffineDblVector): TAffineVector; register; assembler;
  1552.  
  1553. asm
  1554.               FLD  QWORD PTR [EAX]
  1555.               FSTP DWORD PTR [EDX]
  1556.               FLD  QWORD PTR [EAX+8]
  1557.               FSTP DWORD PTR [EDX+4]
  1558.               FLD  QWORD PTR [EAX+16]
  1559.               FSTP DWORD PTR [EDX+8]
  1560. end;
  1561.  
  1562. //------------------------------------------------------------------------------
  1563.  
  1564. function VectorAffineFltToDbl(V: TAffineVector): TAffineDblVector; register; assembler;
  1565.  
  1566. asm
  1567.               FLD  DWORD PTR [EAX]
  1568.               FSTP QWORD PTR [EDX]
  1569.               FLD  DWORD PTR [EAX+8]
  1570.               FSTP QWORD PTR [EDX+4]
  1571.               FLD  DWORD PTR [EAX+16]
  1572.               FSTP QWORD PTR [EDX+8]
  1573. end;
  1574.  
  1575. //------------------------------------------------------------------------------
  1576.  
  1577. function VectorFltToDbl(V: TVector): THomogeneousDblVector; register; assembler;
  1578.  
  1579. asm
  1580.               FLD  DWORD PTR [EAX]
  1581.               FSTP QWORD PTR [EDX]
  1582.               FLD  DWORD PTR [EAX+8]
  1583.               FSTP QWORD PTR [EDX+4]
  1584.               FLD  DWORD PTR [EAX+16]
  1585.               FSTP QWORD PTR [EDX+8]
  1586.               FLD  DWORD PTR [EAX+24]
  1587.               FSTP QWORD PTR [EDX+12]
  1588. end;
  1589.  
  1590. //------------------------------------------------------------------------------
  1591.  
  1592. end.
  1593.