geometry.pas

一个用Delphi编写的很好的屏保程序

                   FSTP ST                       // clean FPU stack
end;

//----------------------------------------------------------------------------------------------------------------------

procedure VectorRotate(var Vector: TVector4f; Axis: TVector3f; Angle: Single);

// rotates Vector about Axis with Angle radiants

var RotMatrix : TMatrix4f;

begin
  RotMatrix := CreateRotationMatrix(Axis, Angle);
  Vector := VectorTransform(Vector, RotMatrix);
end;

//----------------------------------------------------------------------------------------------------------------------

procedure VectorScale(V: array of Single; Factor: Single); assembler; register;

// returns a vector scaled by a factor
// EAX contains address of V
// EDX contains highest index in V
// Factor is located on the stack

asm
  {for I := Low(V) to High(V) do V[I] := V[I] * Factor;}

              FLD DWORD PTR [Factor]        // load factor
@@Loop:       FLD DWORD PTR [EAX + 4 * EDX] // load a component
              FMUL ST, ST(1)                // multiply it with the factor
              WAIT
              FSTP DWORD PTR [EAX + 4 * EDX] // store the result
              DEC EDX                       // do the entire array
              JNS @@Loop
              FSTP ST(0)                    // clean the FPU stack
end;

//----------------------------------------------------------------------------------------------------------------------

procedure VectorNegate(V: array of Single); assembler; register;

// returns a negated vector
// EAX contains address of V
// EDX contains highest index in V

asm
  {V[X] := -V[X];
  V[Y] := -V[Y];
  V[Z] := -V[Z];}

@@Loop:       FLD DWORD PTR [EAX + 4 * EDX]
              FCHS
              WAIT
              FSTP DWORD PTR [EAX + 4 * EDX]
              DEC EDX
              JNS @@Loop
end;

//----------------------------------------------------------------------------------------------------------------------

function VectorAdd(V1, V2: TVector): TVector; register;

// returns the sum of two vectors

begin
  Result[X] := V1[X] + V2[X];
  Result[Y] := V1[Y] + V2[Y];
  Result[Z] := V1[Z] + V2[Z];
  Result[W] := V1[W] + V2[W];
end;

//----------------------------------------------------------------------------------------------------------------------

function VectorAffineAdd(V1, V2: TAffineVector): TAffineVector; register;

// returns the sum of two vectors

begin
  Result[X] := V1[X] + V2[X];
  Result[Y] := V1[Y] + V2[Y];
  Result[Z] := V1[Z] + V2[Z];
end;

//----------------------------------------------------------------------------------------------------------------------

function VectorSubtract(V1, V2: TVector): TVector; register;

// returns the difference of two vectors

begin
  Result[X] := V1[X] - V2[X];
  Result[Y] := V1[Y] - V2[Y];
  Result[Z] := V1[Z] - V2[Z];
  Result[W] := V1[W] - V2[W];
end;

//----------------------------------------------------------------------------------------------------------------------

function VectorDotProduct(V1, V2: TVector): Single; register;

begin
  Result := V1[X] * V2[X] + V1[Y] * V2[Y] + V1[Z] * V2[Z] + V1[W] * V2[W];
end;

//----------------------------------------------------------------------------------------------------------------------

function VectorAffineDotProduct(V1, V2: TAffineVector): Single; assembler; register;

// calculates the dot product between V1 and V2
// EAX contains address of V1
// EDX contains address of V2
// result is stored in ST(0)

asm
  //Result := V1[X] * V2[X] + V1[Y] * V2[Y] + V1[Z] * V2[Z];

                   FLD DWORD PTR [EAX]
                   FMUL DWORD PTR [EDX]
                   FLD DWORD PTR [EAX + 4]
                   FMUL DWORD PTR [EDX + 4]
                   FADDP
                   FLD DWORD PTR [EAX + 8]
                   FMUL DWORD PTR [EDX + 8]
                   FADDP
end;

//----------------------------------------------------------------------------------------------------------------------

function VectorCrossProduct(V1, V2: TAffineVector): TAffineVector; 

// calculates the cross product between vector 1 and 2, Temp is necessary because
// either V1 or V2 could also be the result vector
// 
// EAX contains address of V1
// EDX contains address of V2
// ECX contains address of result

var Temp: TAffineVector;

asm
  {Temp[X] := V1[Y] * V2[Z]-V1[Z] * V2[Y];
  Temp[Y] := V1[Z] * V2[X]-V1[X] * V2[Z];
  Temp[Z] := V1[X] * V2[Y]-V1[Y] * V2[X];
  Result := Temp;}

              PUSH EBX                      // save EBX, must be restored to original value
              LEA EBX, [Temp]
              FLD DWORD PTR [EDX + 8]       // first load both vectors onto FPU register stack
              FLD DWORD PTR [EDX + 4]
              FLD DWORD PTR [EDX + 0]
              FLD DWORD PTR [EAX + 8]
              FLD DWORD PTR [EAX + 4]
              FLD DWORD PTR [EAX + 0]

              FLD ST(1)                     // ST(0) := V1[Y]
              FMUL ST, ST(6)                // ST(0) := V1[Y] * V2[Z]
              FLD ST(3)                     // ST(0) := V1[Z]
              FMUL ST, ST(6)                // ST(0) := V1[Z] * V2[Y]
              FSUBP ST(1), ST               // ST(0) := ST(1)-ST(0)
              FSTP DWORD [EBX]              // Temp[X] := ST(0)
              FLD ST(2)                     // ST(0) := V1[Z]
              FMUL ST, ST(4)                // ST(0) := V1[Z] * V2[X]
              FLD ST(1)                     // ST(0) := V1[X]
              FMUL ST, ST(7)                // ST(0) := V1[X] * V2[Z]
              FSUBP ST(1), ST               // ST(0) := ST(1)-ST(0)
              FSTP DWORD [EBX + 4]          // Temp[Y] := ST(0)
              FLD ST                        // ST(0) := V1[X]
              FMUL ST, ST(5)                // ST(0) := V1[X] * V2[Y]
              FLD ST(2)                     // ST(0) := V1[Y]
              FMUL ST, ST(5)                // ST(0) := V1[Y] * V2[X]
              FSUBP ST(1), ST               // ST(0) := ST(1)-ST(0)
              FSTP DWORD [EBX + 8]          // Temp[Z] := ST(0)
              FSTP ST(0)                    // clear FPU register stack
              FSTP ST(0)
              FSTP ST(0)
              FSTP ST(0)
              FSTP ST(0)
              FSTP ST(0)
              MOV EAX, [EBX]                // copy Temp to Result
              MOV [ECX], EAX
              MOV EAX, [EBX + 4]
              MOV [ECX + 4], EAX
              MOV EAX, [EBX + 8]
              MOV [ECX + 8], EAX
              POP EBX
end;

//----------------------------------------------------------------------------------------------------------------------

function VectorPerpendicular(V, N: TAffineVector): TAffineVector; 

// calculates a vector perpendicular to N (N is assumed to be of unit length)
// subtract out any component parallel to N

var Dot: Single;

begin
   Dot := VectorAffineDotProduct(V, N);
   Result[X] := V[X]-Dot * N[X];
   Result[Y] := V[Y]-Dot * N[Y];
   Result[Z] := V[Z]-Dot * N[Z];
end;

//----------------------------------------------------------------------------------------------------------------------

function VectorTransform(V: TVector4f; M: TMatrix): TVector4f; register;

// transforms a homogeneous vector by multiplying it with a matrix

var TV: TVector4f;

begin
  TV[X] := V[X] * M[X, X] + V[Y] * M[Y, X] + V[Z] * M[Z, X] + V[W] * M[W, X];
  TV[Y] := V[X] * M[X, Y] + V[Y] * M[Y, Y] + V[Z] * M[Z, Y] + V[W] * M[W, Y];
  TV[Z] := V[X] * M[X, Z] + V[Y] * M[Y, Z] + V[Z] * M[Z, Z] + V[W] * M[W, Z];
  TV[W] := V[X] * M[X, W] + V[Y] * M[Y, W] + V[Z] * M[Z, W] + V[W] * M[W, W];
  Result := TV
end;

//----------------------------------------------------------------------------------------------------------------------

function  VectorTransform(V: TVector3f; M: TMatrix): TVector3f;

// transforms an affine vector by multiplying it with a (homogeneous) matrix

var TV: TVector3f;

begin
  TV[X] := V[X] * M[X, X] + V[Y] * M[Y, X] + V[Z] * M[Z, X] + M[W, X];
  TV[Y] := V[X] * M[X, Y] + V[Y] * M[Y, Y] + V[Z] * M[Z, Y] + M[W, Y];
  TV[Z] := V[X] * M[X, Z] + V[Y] * M[Y, Z] + V[Z] * M[Z, Z] + M[W, Z];
  Result := TV;
end;


//----------------------------------------------------------------------------------------------------------------------

function VectorAffineTransform(V: TAffineVector; M: TAffineMatrix): TAffineVector; register;

// transforms an affine vector by multiplying it with a matrix

var TV: TAffineVector;

begin
  TV[X] := V[X] * M[X, X] + V[Y] * M[Y, X] + V[Z] * M[Z, X];
  TV[Y] := V[X] * M[X, Y] + V[Y] * M[Y, Y] + V[Z] * M[Z, Y];
  TV[Z] := V[X] * M[X, Z] + V[Y] * M[Y, Z] + V[Z] * M[Z, Z];
  Result := TV;
end;

//----------------------------------------------------------------------------------------------------------------------

function PointInPolygon(xp, yp : array of Single; x, y: Single): Boolean;

// The code below is from Wm. Randolph Franklin <[e-mail]>
// with some minor modifications for speed.  It returns 1 for strictly
// interior points, 0 for strictly exterior, and 0 or 1 for points on
// the boundary.
// This code is not yet tested!

var I, J: Integer;

begin
  Result := False;
  if High(XP) <> High(YP) then Exit;
  J := High(XP);
  for I := 0 to High(XP) do
  begin
    if ((((yp[I] <= y) and (y < yp[J])) or ((yp[J] <= y) and (y < yp[I]))) and
        (x < (xp[J] - xp[I]) * (y - yp[I]) / (yp[J] - yp[I]) + xp[I]))
    then Result := not Result;
    J := I + 1;
  end;
end;

//----------------------------------------------------------------------------------------------------------------------

function QuaternionConjugate(Q: TQuaternion): TQuaternion; assembler;

// returns the conjugate of a quaternion
// EAX contains address of Q
// EDX contains address of result

asm
              FLD DWORD PTR [EAX]
              FCHS
              WAIT
              FSTP DWORD PTR [EDX]
              FLD DWORD PTR [EAX + 4]
              FCHS
              WAIT
              FSTP DWORD PTR [EDX + 4]
              FLD DWORD PTR [EAX + 8]
              FCHS
              WAIT
              FSTP DWORD PTR [EDX + 8]
              MOV EAX, [EAX + 12]
              MOV [EDX + 12], EAX
end;

//----------------------------------------------------------------------------------------------------------------------

function QuaternionFromPoints(V1, V2: TAffineVector): TQuaternion; assembler;

// constructs a unit quaternion from two points on unit sphere
// EAX contains address of V1
// ECX contains address to result
// EDX contains address of V2

asm
  {Result.ImagPart := VectorCrossProduct(V1, V2);
   Result.RealPart :=  Sqrt((VectorAffineDotProduct(V1, V2) + 1)/2);}

              PUSH EAX
              CALL VectorCrossProduct       // determine axis to rotate about
              POP EAX
              FLD1                          // prepare next calculation
              Call VectorAffineDotProduct   // calculate cos(angle between V1 and V2)
              FADD ST, ST(1)                // transform angle to angle/2 by: cos(a/2)=sqrt((1 + cos(a))/2)
              FXCH ST(1)
              FADD ST, ST
              FDIVP ST(1), ST
              FSQRT
              FSTP DWORD PTR [ECX + 12]     // Result.RealPart := ST(0)
end;

//----------------------------------------------------------------------------------------------------------------------

function QuaternionMultiply(qL, qR: TQuaternion): TQuaternion;

// Returns quaternion product qL * qR.  Note: order is important!
// To combine rotations, use the product QuaternionMuliply(qSecond, qFirst),
// which gives the effect of rotating by qFirst then qSecond.

var Temp : TQuaternion;

begin
  Temp.RealPart := qL.RealPart * qR.RealPart - qL.ImagPart[X] * qR.ImagPart[X] -
                   qL.ImagPart[Y] * qR.ImagPart[Y] - qL.ImagPart[Z] * qR.ImagPart[Z];
  Temp.ImagPart[X] := qL.RealPart * qR.ImagPart[X] + qL.ImagPart[X] * qR.RealPart +
                      qL.ImagPart[Y] * qR.ImagPart[Z] - qL.ImagPart[Z] * qR.ImagPart[Y];
  Temp.ImagPart[Y] := qL.RealPart * qR.ImagPart[Y] + qL.ImagPart[Y] * qR.RealPart +
                      qL.ImagPart[Z] * qR.ImagPart[X] - qL.ImagPart[X] * qR.ImagPart[Z];
  Temp.ImagPart[Z] := qL.RealPart * qR.ImagPart[Z] + qL.ImagPart[Z] * qR.RealPart +
                      qL.ImagPart[X] * qR.ImagPart[Y] - qL.ImagPart[Y] * qR.ImagPart[X];
  Result := Temp;
end;

//----------------------------------------------------------------------------------------------------------------------

function QuaternionToMatrix(Q: TQuaternion): TMatrix; 

// Constructs rotation matrix from (possibly non-unit) quaternion.
// Assumes matrix is used to multiply column vector on the left:
// vnew = mat vold.  Works correctly for right-handed coordinate system
// and right-handed rotations.

// Essentially, this function is the same as CreateRotationMatrix and you can consider it as
// being for reference here.

{var Norm, S,
    XS, YS, ZS,
    WX, WY, WZ,
    XX, XY, XZ,
    YY, YZ, ZZ   : Single;

begin
  Norm := Q.Vector[X] * Q.Vector[X] + Q.Vector[Y] * Q.Vector[Y] + Q.Vector[Z] * Q.Vector[Z] + Q.RealPart * Q.RealPart;
  if Norm > 0 then S := 2 / Norm
              else S := 0;