/*======================================================================*/
/* 3DLib								*/
/*									*/
/* Spline-related functions						*/
/*									*/
/* AUTHOR:	Gabor Nagy						*/
/* DATE:	1996-Jan-10 23:14:43					*/
/*									*/
/* 3DLib(TM) Copyright (C) 1995 by Gabor Nagy. All rights reserved.	*/
/*======================================================================*/
#include <assert.h>
#include <math.h>
#include <float.h>
#include <stdio.h>

#include <EMalloc.h>
#include <EStrings.h>

#include <E3D/E3D.h>

#include <E3D/Matrix.h>

#include <E3D/Spline.h>



static E3dGeometryClass*	E3d_VectorGeometryClass=NULL;


//================================================
// Invert a Spline (reverse its direction)
//================================================
void E3d_SplineInvert(E3dSpline* LSpline)
{
 unsigned int	LC, LN=LSpline->NumOfCVs;

 switch(LSpline->SplineType)
 {
  case E3dSPL_BEZIER:
   {
    E3dBezierCV*	LNewBezierCVs=(E3dBezierCV*)EMalloc(sizeof(E3dBezierCV)*LN);

    if(LNewBezierCVs)
    {
     E3dBezierCV*	LBezierCV=(E3dBezierCV*)(LSpline->CVs);
     E3dBezierCV*	LNewBezierCV=LNewBezierCVs+LN-1;
     E3d3DPosition	LTmpPos;

     for(LC=0;LC<LN;LC++, LBezierCV++, LNewBezierCV--)
     {
      *LNewBezierCV=*LBezierCV;
      LTmpPos=LNewBezierCV->Next;
      LNewBezierCV->Next=LNewBezierCV->Previous;
      LNewBezierCV->Previous=LTmpPos;
     }
     EFree(LSpline->CVs);

     LSpline->CVs=LNewBezierCVs;
    }
   }
  break;

  default:
   {
    E3dSplineCV*	LNewSplineCVs=(E3dSplineCV*)EMalloc(sizeof(E3dSplineCV)*LN);

    if(LNewSplineCVs)
    {
     E3dSplineCV*	LSplineCV=LSpline->CVs;
     E3dSplineCV*	LNewSplineCV=LNewSplineCVs+LN-1;

     for(LC=0;LC<LN;LC++, LSplineCV++, LNewSplineCV--) *LNewSplineCV=*LSplineCV;
    }
    EFree(LSpline->CVs);
    LSpline->CVs=LNewSplineCVs;
   }
  break;
 }
}


//================================================
// Count tagged CVs on a Spline
//================================================
int E3d_SplineNumOfTaggedCVs(E3dSpline* LSpline)
{
 unsigned int	LC, LN=LSpline->NumOfCVs, LNumOfTaggedCVs=0;


 switch(LSpline->SplineType)
 {
  case E3dSPL_BEZIER:
   {
    E3dBezierCV*	LBezierCV=(E3dBezierCV*)(LSpline->CVs);

    for(LC=0;LC<LN;LC++, LBezierCV++) if(LBezierCV->Flags&E3dKEYF_TAGGED) LNumOfTaggedCVs++;
   }
  break;

  default:
   {
    E3dSplineCV*	LSplineCV=(E3dSplineCV*)(LSpline->CVs);

    for(LC=0;LC<LN;LC++, LSplineCV++) if(LSplineCV->Flags&E3dKEYF_TAGGED) LNumOfTaggedCVs++;
   }
  break;
 }

 return(LNumOfTaggedCVs);
}


//================================================
// Get the local bounding box of a Spline
//================================================
EBool E3d_SplineGetBoundingBox(E3dSpline* LSpline, E3d3DPosition* LBBMin, E3d3DPosition* LBBMax, unsigned int LWhatToInclude)
{
#ifdef USEOpenGL
 E3dGLCoordinate*	LPoint=LSpline->GLLinSegments;
 unsigned int		LC, LN;


 if(LPoint)
 {

// Must initialize BBox on first linear segment start
//
  if((LN=LSpline->NumOfLinSegments)>0)
  {
   LBBMin->X = LBBMax->X = LPoint[E3dX];
   LBBMin->Y = LBBMax->Y = LPoint[E3dY];
   LBBMin->Z = LBBMax->Z = LPoint[E3dZ];
   LPoint+=3;
  }

  for(LC=1;LC<LN;LC++)
  {
   if(LPoint[E3dX] < LBBMin->X) LBBMin->X = LPoint[E3dX];
   if(LPoint[E3dX] > LBBMax->X) LBBMax->X = LPoint[E3dX];

   if(LPoint[E3dY] < LBBMin->Y) LBBMin->Y = LPoint[E3dY];
   if(LPoint[E3dY] > LBBMax->Y) LBBMax->Y = LPoint[E3dY];

   if(LPoint[E3dZ] < LBBMin->Z) LBBMin->Z = LPoint[E3dZ];
   if(LPoint[E3dZ] > LBBMax->Z) LBBMax->Z = LPoint[E3dZ];
   LPoint+=3;
  }
  return(TRUE);
 }
#endif // USEOpenGL

 return(FALSE);
}


//========================================================================================
// Get the bounding box of a Spline aligned to an arbitrary coordinate system
//
// Arguments
//  E3dSpline*     LSpline     Pointer to the Spline structure
//  E3dMatrix      LMatrix     The Matrix that defines the coordinate system
//  E3d3DPosition* LBBMin      Returned bounding box point 0
//  E3d3DPosition* LBBMax      Returned bounding box point 1
//
// Description
//  This function transforms each vertex of the given Spline with the given Matrix and
//  returns the minimum and maximum XYZ values of the resulting bounding box.
//
// Example
//  E3d_SplineGetTransformedBoundingBox(LSpline, LModel->LocalToWorldMatrix, LBBMin, LBBMax);
//  This will return the "world-aligned" bounding box of the LSpline (assuming that
//  LSpline is a Geometry of LModel).
//
// Return value
//  TRUE if successful, FALSE in case of an error.
//
// See also
//  E3d_SplineGetBoundingBox
//========================================================================================
EBool E3d_SplineGetTransformedBoundingBox(E3dSpline* LSpline,
                               E3dMatrix LMatrix,
                               E3d3DPosition* LBBMin, E3d3DPosition* LBBMax,
                               unsigned int LWhatToInclude)
{
#ifdef USEOpenGL
 E3dGLCoordinate*	LPoint=LSpline->GLLinSegments;
 unsigned int		LC, LN;


 if(LPoint)
 {
  E3dCoordinate	mX, mY, mZ, LFX, LFY, LFZ;

// Must initialize BBox on first linear segment start
//
  if((LN=LSpline->NumOfLinSegments)>0)
  {
   mX=LPoint[E3dX];mY=LPoint[E3dY];mZ=LPoint[E3dZ];
   E3dM_MatrixTransform3x4(LMatrix, LBBMax->X, LBBMax->Y, LBBMax->Z);
   LBBMin->X = LBBMax->X;
   LBBMin->Y = LBBMax->Y;
   LBBMin->Z = LBBMax->Z;
   LPoint+=3;
  }

  for(LC=1;LC<LN;LC++)
  {
   mX=LPoint[E3dX];mY=LPoint[E3dY];mZ=LPoint[E3dZ];
   E3dM_MatrixTransform3x4(LMatrix, LFX, LFY, LFZ);
   if(LFX < LBBMin->X) LBBMin->X = LFX;
   if(LFX > LBBMax->X) LBBMax->X = LFX;

   if(LFY < LBBMin->Y) LBBMin->Y = LFY;
   if(LFY > LBBMax->Y) LBBMax->Y = LFY;

   if(LFZ < LBBMin->Z) LBBMin->Z = LFZ;
   if(LFZ > LBBMax->Z) LBBMax->Z = LFZ;
   LPoint+=3;
  }
  return(TRUE);
 }
#endif // USEOpenGL

 return(FALSE);
}


//================================================================
// Compute the bounding box of a Spline
//
// LMaxProjectedXYAbsValPoint will have the non-transformed XYZ
// of the Vertex with largest projected absolute value of
// X or Y.
// The W will have that maximum X or Y absolute value.
//==============================================================*/
EBool E3d_SplineGetPerspectiveProjectedBoundingBox(E3dSpline* LSpline, E3dMatrix LMatrix, E3d2DPosition* LBBMin, E3d2DPosition* LBBMax, E3dHPosition* LMaxProjectedXYAbsValPoint, unsigned int LWhatToInclude)
{
 E3dGLCoordinate*	LGLLinSegment=LSpline->GLLinSegments;


 if(LGLLinSegment)
 {
  E3dCoordinate	mX, mY, mZ, LXF, LYF, LZF, LWF,
		LXOfMaxXYAbsVal, LYOfMaxXYAbsVal, LZOfMaxXYAbsVal, LMaxXYAbsVal;
  unsigned int	LC, LN;


// Must initialize BBox on first Point
//
  if((LN=LSpline->NumOfLinSegments)>0)
  {
   mX=LGLLinSegment[E3dX];mY=LGLLinSegment[E3dY];mZ=LGLLinSegment[E3dZ];

   E3dM_MatrixTransformHPosition(LMatrix, LXF, LYF, LZF, LWF);LWF=1.0/LWF;
   LXF*=LWF;LYF*=LWF;

   LBBMin->X = LBBMax->X = LXF;
   LBBMin->Y = LBBMax->Y = LYF;

   LXOfMaxXYAbsVal = LXF;
   LYOfMaxXYAbsVal = LYF;
   LZOfMaxXYAbsVal = LZF;
   if(E3dM_ABS(LXF)>E3dM_ABS(LYF)) LMaxXYAbsVal=E3dM_ABS(LXF);
   else LMaxXYAbsVal=E3dM_ABS(LYF);

   LGLLinSegment+=3;

   for(LC=1;LC<LN;LC++)
   {
    mX=LGLLinSegment[E3dX];mY=LGLLinSegment[E3dY];mZ=LGLLinSegment[E3dZ];

    E3dM_MatrixTransformHPosition(LMatrix, LXF, LYF, LZF, LWF);LWF=1.0/LWF;
    LXF*=LWF;LYF*=LWF;

    if(LXF < LBBMin->X) LBBMin->X = LXF;
    if(LXF > LBBMax->X) LBBMax->X = LXF;

    if(LYF < LBBMin->Y) LBBMin->Y = LYF;
    if(LYF > LBBMax->Y) LBBMax->Y = LYF;

    if(E3dM_ABS(LXF)>E3dM_ABS(LYF))
    {
     if(E3dM_ABS(LXF)>LMaxXYAbsVal)
     {
      LXOfMaxXYAbsVal = LXF;
      LYOfMaxXYAbsVal = LYF;
      LZOfMaxXYAbsVal = LZF;
      LMaxXYAbsVal=E3dM_ABS(LXF);
     }
    }
    else
    {
     if(E3dM_ABS(LYF)>LMaxXYAbsVal)
     {
      LXOfMaxXYAbsVal = LXF;
      LYOfMaxXYAbsVal = LYF;
      LZOfMaxXYAbsVal = LZF;
      LMaxXYAbsVal=E3dM_ABS(LYF);
     }
    }

    LGLLinSegment+=3;
   }

   LMaxProjectedXYAbsValPoint->X = LXOfMaxXYAbsVal;
   LMaxProjectedXYAbsValPoint->Y = LYOfMaxXYAbsVal;
   LMaxProjectedXYAbsValPoint->Z = LZOfMaxXYAbsVal;
   LMaxProjectedXYAbsValPoint->W = LMaxXYAbsVal;

   return(TRUE);
  }
 }
 return(FALSE);
}


//========================================================================================
// Allocate array of independent Edges (line segments)
//
// Argument
//  unsigned int LEdgeNum          Number of Edges to allocate
//
// Description
//  This function allocates memory for an array of edges and initializes all the Edges.
//
// Return value
//  Pointer to the array of Edges or NULL in case of an error.
//========================================================================================
E3dEdge* E3d_EdgesAllocate(unsigned int LEdgeNum)
{ 
 if(LEdgeNum)
 {
  E3dEdge*	LEdges=(E3dEdge*)EMalloc(sizeof(E3dEdge)*LEdgeNum);

  if(LEdges)
  {
   E3dEdge*	LEdge=LEdges;
   unsigned int	LECnt;

   for(LECnt=0;LECnt<LEdgeNum;LECnt++, LEdge++)
   {
    LEdge->Start=0;
    LEdge->End=0;
    LEdge->Flags=0;
   }
  }
  return(LEdges);
 }
 return(NULL);
}


//========================================================================================
// Allocate array of PolyLines (connected lines)
//
// Argument
//  unsigned int LPolyLineNum      Number of PolyLines to allocate
//
// Description
//  This function allocates memory for an array of points for PolyLines and initializes
//  all the elements of the array.
//
// Return value
//  Pointer to the array or NULL in case of an error.
//========================================================================================
E3dPolyLine* E3d_PolyLinesAllocate(unsigned int LPolyLineNum)
{ 
 E3dPolyLine*	LPolyLines=NULL;
 E3dPolyLine*	LPolyLine;
 unsigned int	LLECnt;


 if(LPolyLineNum>0) LPolyLines=(E3dPolyLine*)EMalloc(sizeof(E3dPolyLine)*LPolyLineNum);
 {
  if(LPolyLines)
  {
   LPolyLine=LPolyLines;
   for(LLECnt=0;LLECnt<LPolyLineNum;LLECnt++, LPolyLine++)
   {
    LPolyLine->NumOfVertices=0;
    LPolyLine->VertexIDs=NULL;

#ifdef USEOpenGL
    LPolyLine->GLVertices=NULL;
#endif // USEOpenGL

    LPolyLine->ColorIndex=1;
   }
  }
  return(LPolyLines);
 }
}


//================================================
// Allocate VertexID array for a linked edge
//================================================
void E3d_PolyLineAllocateVertices(E3dPolyLine* LIPolyLine, unsigned int LVertexNum, int LFlags)
{
 if(LVertexNum>0)
 {
  LIPolyLine->VertexIDs=(unsigned int*)EMalloc(sizeof(unsigned int)*LVertexNum);
#ifdef USEOpenGL
  if(LFlags&E3dALLOCFORGL) LIPolyLine->GLVertices=(E3dGLType*)EMalloc(sizeof(E3dGLType)*LVertexNum*3);
#endif // USEOpenGL
 }
}


//========================================
// Free the vertices of a PolyLine
//========================================
void E3d_PolyLineFreeVertices(E3dPolyLine* LPolyLine)
{
 if(LPolyLine->VertexIDs) EFree(LPolyLine->VertexIDs);
#ifdef USEOpenGL
 if(LPolyLine->GLVertices) EFree(LPolyLine->GLVertices);
#endif // USEOpenGL
}


//================================================================================================
// Allocate a Spline of the given type
//
// Argument
//  int    LType        Type of Spline to allocate
//
// Description
//  This function allocates and initializes an E3dSpline structure and sets the type to LType.
//  The currently supported types are:
//  - E3dSPL_LINEAR    "a.k.a. PolyLine"; all segments are linear
//  - E3dSPL_BEZIER    segment and key types can vary (you can mix linear and curved segments)
//  - E3dSPL_BSPLINE
//  - E3dSPL_CARDINAL
//  Note that the Spline type can be changed later, but the old spline keys have to be removed
//  and reallocated as they might be incompatible with the new Spline type (e.g. Bezier Spline
//  keys have a different format than B-Spline keys).
//
// Return value
//  Pointer to the allocated E3dSpline structure or NULL in case of an error.
//
// See also
//  E3d_SplineFree
//================================================================================================
E3dSpline* E3d_SplineAllocate(int LType)
{
 E3dSpline*	LSpline=(E3dSpline*)EMalloc(sizeof(E3dSpline));

 if(LSpline)
 {
  E3dM_GeometryDefault((E3dGeometry*)LSpline);

  LSpline->GeoType=E3dGEO_SPLINE;

  LSpline->SplineType=LType;
  LSpline->NumOfCVs=0;
  LSpline->Length=-1.0;		// "Invalid" (not computed yet)

  LSpline->CVs=NULL;
  LSpline->NumOfDrawSteps=8;
  LSpline->Phantom=0;
  LSpline->Tension=0.0;
  LSpline->Constraint=E3dBEZ_C0;
  LSpline->Segment=E3dBEZ_CURVE;

// This is per segment...
// The arc-length parametrization function does a binary search, so powers of 2 seem nicer here
//
  LSpline->NumOfStepsForLength=2048;

// Drawing steps (linear segments) for the whole Spline
//
  LSpline->NumOfLinSegments=0;
  LSpline->GLLinSegments=NULL;

  LSpline->UniformSteps=TRUE;
  LSpline->Closed=FALSE;
 }
 return(LSpline);
}


//========================================================================================
// Free a Spline's rendering-related data, such as OpenGL vertex arrays etc.
// we usually do this when a Spline is removed from the Scene, but stored for undo.
//========================================================================================
void E3d_SplineFreeRenderData(E3dSpline* LSpline)
{
 unsigned int	LC, LN=LSpline->NumOfCVs;

 if(LSpline->GLLinSegments) { EFree(LSpline->GLLinSegments);LSpline->GLLinSegments=NULL;LSpline->NumOfLinSegments=0; }

 switch(LSpline->SplineType)
 {
  case E3dSPL_BEZIER:
   {
    E3dBezierCV*	LBezierCV=(E3dBezierCV*)(LSpline->CVs);

    for(LC=0;LC<LN;LC++, LBezierCV++)
    {
     if(LBezierCV->StepsForLength) { EFree(LBezierCV->StepsForLength);LBezierCV->StepsForLength=NULL;LBezierCV->NumOfStepsForLength=0; }
    }
   }
  break;

  default:
   {
    E3dSplineCV*	LSplineCV=(E3dSplineCV*)(LSpline->CVs);

    for(LC=0;LC<LN;LC++, LSplineCV++)
    {
     if(LSplineCV->StepsForLength) { EFree(LSplineCV->StepsForLength);LSplineCV->StepsForLength=NULL;LSplineCV->NumOfStepsForLength=0; }
    }
   }
  break;
 }
}


/*======================================================================================*/
// Free Control Vertices of a Spline
//
// Argument
//  E3dSpline* LSpline         Pointer to the Spline structure
//
// Description
//  This function frees all the CVs of the given Spline.
//
// See also
//  E3d_SplineAllocate
/*======================================================================================*/
void E3d_SplineFreeCVs(E3dSpline* LSpline)
{
 if(LSpline->CVs)
 {
  unsigned int	LC, LN=LSpline->NumOfCVs;


  switch(LSpline->SplineType)
  {
   case E3dSPL_BEZIER:
    {
     E3dBezierCV*	LBezierCV=(E3dBezierCV*)(LSpline->CVs);

     for(LC=0;LC<LN;LC++, LBezierCV++)
     {
      if(LBezierCV->StepsForLength) { EFree(LBezierCV->StepsForLength);LBezierCV->StepsForLength=NULL;LBezierCV->NumOfStepsForLength=0; }
     }
    }
   break;

   default:
    {
     E3dSplineCV*	LSplineCV=(E3dSplineCV*)(LSpline->CVs);

     for(LC=0;LC<LN;LC++, LSplineCV++)
     {
      if(LSplineCV->StepsForLength) { EFree(LSplineCV->StepsForLength);LSplineCV->StepsForLength=NULL;LSplineCV->NumOfStepsForLength=0; }
     }
    }
   break;
  }

  EFree(LSpline->CVs);LSpline->CVs=NULL;
 }

 if(LSpline->GLLinSegments) { EFree(LSpline->GLLinSegments);LSpline->GLLinSegments=NULL; }
}


//========================================================================================
// Free a Spline structure
//
// Argument
//  E3dSpline* LSpline         Pointer to the Spline structure to be freed
//
// Description
//  This function first decrements the reference count (RefCnt) of the given Spline
//  structure (if it's not already zero).
//  After that, if RefCnt is still greater than zero, it means that this Spline is
//  bein referred to somewhere else, so E3d_SplineFree() will simply return.
//  If RefCnt is zero, E3d_SplineFree() will free all the keys of the Spline and then
//  the Spline structure itself.
//
// See also
//  E3d_SplineAllocate
//========================================================================================
void E3d_SplineFree(E3dSpline* LSpline)
{
//printf(" E3d_SplineFree %p LSpline->RefCnt %d\n", LSpline, LSpline->RefCnt);fflush(stdout);

 if(LSpline->RefCnt>0) LSpline->RefCnt-=1;
 if(LSpline->RefCnt>0) return;


 E3d_SplineFreeCVs(LSpline);

 EFree(LSpline);
}


//================================================================================
// Duplicate a Spline
//
// Argument
//  E3dSpline* LSpline         Pointer to the Spline structure to be cloned
//
// Description
//  This function creates and exact copy of the given Spline.
//
// Return value
//  Pointer to the new Spline structure, or NULL in case of an error
//
// See also
//  E3d_SplineAllocate, E3d_SplineFree
//================================================================================
E3dSpline* E3d_SplineClone(E3dSpline* LSpline)
{
 E3dSpline*	LDSpline=E3d_SplineAllocate(LSpline->SplineType);

 if(LDSpline)
 {
  int	LC, LN=LSpline->NumOfCVs;

  if((LDSpline->CVs=E3d_SplineCVsAllocate(LSpline->SplineType, LSpline->NumOfCVs))!=NULL)
  {
   switch(LSpline->SplineType)
   {
    case E3dSPL_BEZIER:
     {
      E3dBezierCV*	LBezierCV=(E3dBezierCV*)(LDSpline->CVs);

      memcpy(LDSpline->CVs, LSpline->CVs, sizeof(E3dBezierCV)*LSpline->NumOfCVs);
      for(LC=0;LC<LN;LC++, LBezierCV++)
      {
       LBezierCV->NumOfStepsForLength=0;
       LBezierCV->StepsForLength=NULL;
      }
     }
    break;

    default:
     {
      E3dSplineCV*	LSplineCV=(E3dSplineCV*)(LDSpline->CVs);

      memcpy(LDSpline->CVs, LSpline->CVs, sizeof(E3dSplineCV)*LSpline->NumOfCVs);
      for(LC=0;LC<LN;LC++, LSplineCV++)
      {
       LSplineCV->NumOfStepsForLength=0;
       LSplineCV->StepsForLength=NULL;
      }
     }
    break;
   }

   LDSpline->NumOfCVs=LSpline->NumOfCVs;
   LDSpline->Length=LSpline->Length;
   LDSpline->NumOfDrawSteps=LSpline->NumOfDrawSteps;
   LDSpline->Phantom=LSpline->Phantom;
   LDSpline->Tension=LSpline->Tension;
   LDSpline->Constraint=LSpline->Constraint;
   LDSpline->Segment=LSpline->Segment;
   LDSpline->NumOfStepsForLength=LSpline->NumOfStepsForLength;

   LDSpline->UniformSteps=LSpline->UniformSteps;
   LDSpline->Closed=LSpline->Closed;

   return(LDSpline);
  }
 }
 return(NULL);
}


//================================================================================
// Set a Spline Control Vertex to a default
//
// Arguments
//  int          LSplineType   Type of the Spline the CVs will be used with
//  void*        LCV           Pointer to the Spline CV
//
// Description
//  This function initializes the given Spline key (CV), depending on the type
//  of Spline it is assigned to.
//  Bezier Spline CVs are treated differently from other types as they have
//  more information.
//
// See also
//  E3d_SplineCVsAllocate
//================================================================================
void E3d_SplineCVDefault(int LSplineType, void* LCV)
{
 switch(LSplineType)
 {
  case E3dSPL_BEZIER:
   {
    E3dBezierCV*	LBezierCV=(E3dBezierCV*)LCV;

    LBezierCV->Previous.X=0.0;LBezierCV->Previous.Y=0.0;LBezierCV->Previous.Z=0.0;
    LBezierCV->Position.X=0.0;LBezierCV->Position.Y=0.0;LBezierCV->Position.Z=0.0;
    LBezierCV->Next.X=0.0;LBezierCV->Next.Y=0.0;LBezierCV->Next.Z=0.0;
    E3dM_InitBezierCV(LBezierCV);
   }
  break;

  default:
   {
    E3dSplineCV*	LSplineCV=(E3dSplineCV*)LCV;

    LSplineCV->Position.X=0.0;LSplineCV->Position.Y=0.0;LSplineCV->Position.Z=0.0;

    E3dM_InitSplineCV(LSplineCV);
   }
  break;
 }
}


//================================================================================
// Allocate array of Spline Control Vertices (CVs)
//
// Arguments
//  int          LSplineType   Type of the spline the keys will be used with
//  unsigned int LNumOfCVs     Number of Spline CVs to allocate
//
// Description
//  This function allocates a number of Spline key structures in an array
//  and initializes them.
//  The type of the structures allocated depends on LSplineType.
//
// Return value
//  Pointer to the array of allocated spline key or NULL in case of an error.
//
// See also
//  E3d_SplineAllocate
//================================================================================
void* E3d_SplineCVsAllocate(int LSplineType, unsigned int LNumOfCVs)
{
 E3dSplineCV*	LSplineCVs;
 E3dSplineCV*	LSplineCV;
 E3dBezierCV*	LBezierCVs;
 E3dBezierCV*	LBezierCV;
 unsigned int	LCnt;


 switch(LSplineType)
 {
  case E3dSPL_BEZIER:
   if((LBezierCVs=(E3dBezierCV*)EMalloc(sizeof(E3dBezierCV)*LNumOfCVs))!=NULL)
   {
    LBezierCV=LBezierCVs;
    for(LCnt=0;LCnt<LNumOfCVs;LCnt++, LBezierCV++)
    {
     LBezierCV->Previous.X=0.0;LBezierCV->Previous.Y=0.0;LBezierCV->Previous.Z=0.0;
     LBezierCV->Position.X=0.0;LBezierCV->Position.Y=0.0;LBezierCV->Position.Z=0.0;
     LBezierCV->Next.X=0.0;LBezierCV->Next.Y=0.0;LBezierCV->Next.Z=0.0;
     E3dM_InitBezierCV(LBezierCV);
    }
    return((void*)LBezierCVs);
   }
  break;

  default:
   if((LSplineCVs=(E3dSplineCV*)EMalloc(sizeof(E3dSplineCV)*LNumOfCVs))!=NULL)
   {
    LSplineCV=LSplineCVs;
    for(LCnt=0;LCnt<LNumOfCVs;LCnt++, LSplineCV++)
    {
     LSplineCV->Position.X=0.0;LSplineCV->Position.Y=0.0;LSplineCV->Position.Z=0.0;
     E3dM_InitSplineCV(LSplineCV);
    }
   }
   return((void*)LSplineCVs);
 }
 return(NULL);
}


//========================================
// Add a Spline CV to CV array
//========================================
void* E3d_SplineAddCV(E3dSpline* LSpline)
{
 void*		LSplineCVsT;

 if(LSpline->CVs==NULL)
 {
  if((LSpline->CVs=E3d_SplineCVsAllocate(LSpline->SplineType, 1))!=NULL)
  {
   LSpline->NumOfCVs=1;
   return(LSpline->CVs);
  }
  else return(NULL);
 }
 else
 {
  unsigned int	LCVSize;

  switch(LSpline->SplineType)
  {
   case E3dSPL_BEZIER:
    LCVSize=sizeof(E3dBezierCV);
   break;

   default:
    LCVSize=sizeof(E3dSplineCV);
   break;
  }

  LSplineCVsT=ERealloc(LSpline->CVs, LCVSize*(LSpline->NumOfCVs+1));
  if(LSplineCVsT)
  {
   LSpline->CVs=LSplineCVsT;
   LSpline->NumOfCVs+=1;
   return((void*)(((unsigned char*)(LSpline->CVs))+LCVSize*(LSpline->NumOfCVs-1)));
  }
  else return(NULL);
 }
}


//================================================
// Transform a Spline by a Matrix
//================================================
void E3d_SplineTransform(E3dSpline* LSpline, E3dMatrix LMatrix)
{
 E3dCoordinate	mX, mY, mZ;
 unsigned int	LVC, LVN=LSpline->NumOfCVs;

 if(LVN==0) return;

 switch(LSpline->SplineType)
 {
  case E3dSPL_BEZIER:
   {
    E3dBezierCV*	LBezierCV=(E3dBezierCV*)(LSpline->CVs);

// Transform the CVs with LMatrix
//
    for(LVC=0;LVC<LVN;LVC++, LBezierCV++)
    {
     mX=LBezierCV->Position.X;mY=LBezierCV->Position.Y;mZ=LBezierCV->Position.Z;
     E3dM_MatrixTransform3x4(LMatrix, LBezierCV->Position.X, LBezierCV->Position.Y, LBezierCV->Position.Z);
     mX=LBezierCV->Previous.X;mY=LBezierCV->Previous.Y;mZ=LBezierCV->Previous.Z;
     E3dM_MatrixTransform3x4(LMatrix, LBezierCV->Previous.X, LBezierCV->Previous.Y, LBezierCV->Previous.Z);
     mX=LBezierCV->Next.X;mY=LBezierCV->Next.Y;mZ=LBezierCV->Next.Z;
     E3dM_MatrixTransform3x4(LMatrix, LBezierCV->Next.X, LBezierCV->Next.Y, LBezierCV->Next.Z);
    }
   }
  break;

  default:
   {
    E3dSplineCV*	LSplineCV=(E3dSplineCV*)(LSpline->CVs);

// Transform the CVs with LMatrix
//
    for(LVC=0;LVC<LVN;LVC++, LSplineCV++)
    {
     mX=LSplineCV->Position.X;mY=LSplineCV->Position.Y;mZ=LSplineCV->Position.Z;
     E3dM_MatrixTransform3x4(LMatrix, LSplineCV->Position.X, LSplineCV->Position.Y, LSplineCV->Position.Z);
    }
   }
  break;
 }
}



E3dMatrix E3d_BezierMatrix=
{
 -1.0,  3.0, -3.0,  1.0,
  3.0, -6.0,  3.0,  0.0,
 -3.0,  3.0,  0.0,  0.0,
  1.0,  0.0,  0.0,  0.0
};

E3dMatrix E3d_CardinalMatrix=
{
 -0.5,  1.5, -1.5,  0.5,
  1.0, -2.5,  2.0, -0.5,
 -0.5,  0.0,  0.5,  0.0,
  0.0,  1.0,  0.0,  0.0
};

E3dMatrix E3d_BSplineMatrix=
{
 -1.0/6.0,  3.0/6.0, -3.0/6.0, 1.0/6.0,
  3.0/6.0, -6.0/6.0,  3.0/6.0, 0.0,
 -3.0/6.0,  0.0,      3.0/6.0, 0.0,
  1.0/6.0,  4.0/6.0,  1.0/6.0, 0.0
};


//========================================================
// Tesselate a BSpline or a Cardinal Spline for drawing
//========================================================
static EBool _TesselateSplineForDrawing(E3dMatrix LBaseMatrix, E3dSpline* LSpline)
{
 E3dMatrix		LPrecMatrix, LGMatrix;				// Precision Matrix, Geometry Matrix (CVs)
 E3dSplineCV*		LSplineCVs=LSpline->CVs;
 E3dSplineCV*		LSplineCV;
 E3dGLCoordinate*	LGLSegments;
 E3dCoordinate		LStepsF, LStepsF2, LStepsF3;
 unsigned int		LCnt, LGLCnt, LSCnt, LNumOfCVs=LSpline->NumOfCVs, LNCVs,
			LLinSteps=0, LNumOfLinearSegments;
 EBool			LPerSegmentSteps;


 if(LNumOfCVs<2) return(FALSE);

 LNCVs=LNumOfCVs-2;

 if(LSpline->UniformSteps)
 {
  LLinSteps=LSpline->NumOfDrawSteps;
  if(LSpline->Closed) LNumOfLinearSegments=LNumOfCVs*LLinSteps;
  else
  {
   if(LNumOfCVs<3) return(FALSE);
   LNumOfLinearSegments=(LNumOfCVs-3)*LLinSteps+1;
  }
  LPerSegmentSteps=FALSE;
  E3dM_CurvePrecision(LLinSteps);
 }
 else
 {
  LNumOfLinearSegments=0;
  LPerSegmentSteps=TRUE;
  LSplineCV=LSplineCVs+1;
  for(LCnt=1, LGLCnt=0;LCnt<LNCVs;LCnt++, LSplineCV++) LNumOfLinearSegments+=LSplineCV->NumOfDrawSteps;
  if(LSpline->Closed)
  {
   LNumOfLinearSegments+=LSplineCV->NumOfDrawSteps;LSplineCV++;
   LNumOfLinearSegments+=LSplineCV->NumOfDrawSteps;
   LNumOfLinearSegments+=LSplineCVs->NumOfDrawSteps;
  }
  else LNumOfLinearSegments++;
 }

 if(LNumOfLinearSegments<1) return(FALSE);

 if(LSpline->NumOfLinSegments==LNumOfLinearSegments) LGLSegments=LSpline->GLLinSegments;
 else { if(LSpline->GLLinSegments) EFree(LSpline->GLLinSegments);LGLSegments=LSpline->GLLinSegments=(E3dGLCoordinate*)EMalloc(sizeof(E3dGLCoordinate)*3*LNumOfLinearSegments); }

 if(LGLSegments)
 {
  LSplineCV=LSplineCVs+1;
  for(LCnt=1, LGLCnt=0;LCnt<LNCVs;LCnt++, LSplineCV++)
  {
   if(LPerSegmentSteps)
   {
    LLinSteps=LSplineCV->NumOfDrawSteps;
    E3dM_CurvePrecision(LLinSteps);
   }

   E3dM_CurveBasis(LBaseMatrix, LSplineCVs[LCnt-1].Position, LSplineCV->Position, LSplineCVs[LCnt+1].Position, LSplineCVs[LCnt+2].Position);

   LGLSegments[LGLCnt]=LGMatrix[M30];
   LGLSegments[LGLCnt+1]=LGMatrix[M31];
   LGLSegments[LGLCnt+2]=LGMatrix[M32];
   LGLCnt+=3;
   for(LSCnt=1;LSCnt<LLinSteps;LSCnt++)
   {
    LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
    LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
    LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
    LGLSegments[LGLCnt]=LGMatrix[M30];
    LGLSegments[LGLCnt+1]=LGMatrix[M31];
    LGLSegments[LGLCnt+2]=LGMatrix[M32];
    LGLCnt+=3;
   }
  }

  if(LSpline->Closed)
  {
   if(LPerSegmentSteps)
   {
    LLinSteps=LSplineCV->NumOfDrawSteps;
    E3dM_CurvePrecision(LLinSteps);
   }
   E3dM_CurveBasis(LBaseMatrix, LSplineCVs[LNumOfCVs-3].Position, LSplineCV->Position, LSplineCVs[LNumOfCVs-1].Position, LSplineCVs[0].Position);
   LGLSegments[LGLCnt]=LGMatrix[M30];
   LGLSegments[LGLCnt+1]=LGMatrix[M31];
   LGLSegments[LGLCnt+2]=LGMatrix[M32];
   LGLCnt+=3;
   for(LSCnt=1;LSCnt<LLinSteps;LSCnt++)
   {
    LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
    LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
    LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
    LGLSegments[LGLCnt]=LGMatrix[M30];
    LGLSegments[LGLCnt+1]=LGMatrix[M31];
    LGLSegments[LGLCnt+2]=LGMatrix[M32];
    LGLCnt+=3;
   }

   LSplineCV++;
   if(LPerSegmentSteps)
   {
    LLinSteps=LSplineCV->NumOfDrawSteps;
    E3dM_CurvePrecision(LLinSteps);
   }
   E3dM_CurveBasis(LBaseMatrix, LSplineCVs[LNumOfCVs-2].Position, LSplineCV->Position, LSplineCVs[0].Position, LSplineCVs[1].Position);
   LGLSegments[LGLCnt]=LGMatrix[M30];
   LGLSegments[LGLCnt+1]=LGMatrix[M31];
   LGLSegments[LGLCnt+2]=LGMatrix[M32];
   LGLCnt+=3;
   for(LSCnt=1;LSCnt<LLinSteps;LSCnt++)
   {
    LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
    LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
    LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
    LGLSegments[LGLCnt]=LGMatrix[M30];
    LGLSegments[LGLCnt+1]=LGMatrix[M31];
    LGLSegments[LGLCnt+2]=LGMatrix[M32];
    LGLCnt+=3;
   }

   LSplineCV=LSplineCVs;
   if(LPerSegmentSteps)
   {
    LLinSteps=LSplineCV->NumOfDrawSteps;
    E3dM_CurvePrecision(LLinSteps);
   }
   E3dM_CurveBasis(LBaseMatrix, LSplineCVs[LNumOfCVs-1].Position, LSplineCV->Position, LSplineCVs[1].Position, LSplineCVs[2].Position);
   LGLSegments[LGLCnt]=LGMatrix[M30];
   LGLSegments[LGLCnt+1]=LGMatrix[M31];
   LGLSegments[LGLCnt+2]=LGMatrix[M32];
   LGLCnt+=3;
   for(LSCnt=1;LSCnt<LLinSteps;LSCnt++)
   {
    LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
    LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
    LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
    LGLSegments[LGLCnt]=LGMatrix[M30];
    LGLSegments[LGLCnt+1]=LGMatrix[M31];
    LGLSegments[LGLCnt+2]=LGMatrix[M32];
    LGLCnt+=3;
   }
//printf("Closed Cardinal LinSeg: %d <-> %d\n", LNumOfLinearSegments, LGLCnt/3);fflush(stdout);
  }
  else
  {
   switch(LSpline->SplineType)
   {
    case E3dSPL_BSPLINE:
     E3dM_CurveBasis(LBaseMatrix, LSplineCVs[LNCVs-1].Position, LSplineCVs[LNCVs].Position, LSplineCVs[LNCVs+1].Position, LSplineCVs[0].Position);
     LGLSegments[LGLCnt]=LGMatrix[M30];
     LGLSegments[LGLCnt+1]=LGMatrix[M31];
     LGLSegments[LGLCnt+2]=LGMatrix[M32];
     LGLCnt+=3;
    break;

    case E3dSPL_CARDINAL:
     LGLSegments[LGLCnt]=LSplineCVs[LCnt].Position.X;
     LGLSegments[LGLCnt+1]=LSplineCVs[LCnt].Position.Y;
     LGLSegments[LGLCnt+2]=LSplineCVs[LCnt].Position.Z;
    break;
   }
//printf("Open Cardinal LinSeg: %d <-> %d\n", LNumOfLinearSegments, LGLCnt/3+1);fflush(stdout);
  }
  LSpline->NumOfLinSegments=LNumOfLinearSegments;
  return(TRUE);
 }
 else LSpline->NumOfLinSegments=0;

 return(FALSE);
}


//================================================================
// Create linear approximation of spline for drawing
//================================================================
EBool E3d_SplineCreateLinearSegments(E3dSpline* LSpline)
{
 E3dGLCoordinate*	LGLSegments;
 E3dMatrix		LPrecMatrix, LGMatrix;	// Precision matrix, Geometry matrix (control points)
 E3dGLCoordinate	LStepsF, LStepsF2, LStepsF3;
 unsigned int		LCnt, LSCnt, LGLCnt, LNumOfCVs=LSpline->NumOfCVs,
			LNCVs, LNumOfLinearSegments;


 if(LNumOfCVs)
 {
  switch(LSpline->SplineType)
  {
   case E3dSPL_LINEAR:
    {
     E3dSplineCV*	LSplineCVs;

     LSplineCVs=LSpline->CVs;
     LNumOfLinearSegments=LNumOfCVs;
     if(LSpline->NumOfLinSegments==LNumOfLinearSegments) LGLSegments=LSpline->GLLinSegments;
     else
     {
      if(LSpline->GLLinSegments) EFree(LSpline->GLLinSegments);
      LGLSegments=LSpline->GLLinSegments=(E3dGLCoordinate*)EMalloc(sizeof(E3dGLCoordinate)*3*LNumOfLinearSegments);
     }

     if(LGLSegments)
     {
      for(LCnt=0, LGLCnt=0;LCnt<LNumOfCVs;LCnt++, LGLCnt+=3)
      {
       LGLSegments[LGLCnt]=LSplineCVs[LCnt].Position.X;
       LGLSegments[LGLCnt+1]=LSplineCVs[LCnt].Position.Y;
       LGLSegments[LGLCnt+2]=LSplineCVs[LCnt].Position.Z;
      }
      LSpline->NumOfLinSegments=LNumOfCVs;
      return(TRUE);
     }
    }
   break;

   case E3dSPL_BEZIER:
    {
     E3dBezierCV*	LBezierCVs=(E3dBezierCV*)(LSpline->CVs);
     E3dBezierCV*	LBezierCV=LBezierCVs;
     E3dBezierCV*	LNextBezierCV=LNextBezierCV=LBezierCV+1;
     unsigned int	LLinSteps=0;
     EBool		LPerSegmentSteps;


     LNCVs=LNumOfCVs-1;

     if(LSpline->UniformSteps)
     {
      LLinSteps=LSpline->NumOfDrawSteps;E3dM_CurvePrecision(LLinSteps); 
      LPerSegmentSteps=FALSE;
     }
     else LPerSegmentSteps=TRUE;

// Calculate number of segments needed
//
     LNumOfLinearSegments=0;
     for(LCnt=0;LCnt<LNCVs;LCnt++, LBezierCV++)
     {
      switch(LBezierCV->Segment)
      {
       case E3dBEZ_CURVE:
        if(LPerSegmentSteps) LNumOfLinearSegments+=LBezierCV->NumOfDrawSteps;
        else LNumOfLinearSegments+=LLinSteps;
       break;

       case E3dBEZ_LINEAR:	LNumOfLinearSegments++;break;
      }
     }
     if(LSpline->Closed)
     {
      switch(LBezierCV->Segment)
      {
       case E3dBEZ_CURVE:
        if(LPerSegmentSteps) LNumOfLinearSegments+=LBezierCV->NumOfDrawSteps;
        else LNumOfLinearSegments+=LLinSteps;
       break;

       case E3dBEZ_LINEAR:	LNumOfLinearSegments++;break;
      }
     }

     if(!(LSpline->Closed)) LNumOfLinearSegments++;		// Starting point of Spline

     LBezierCV=LBezierCVs;
     if(LSpline->NumOfLinSegments==LNumOfLinearSegments) LGLSegments=LSpline->GLLinSegments;
     else
     {
      if(LSpline->GLLinSegments) EFree(LSpline->GLLinSegments);
      LGLSegments=LSpline->GLLinSegments=(E3dGLCoordinate*)EMalloc(sizeof(E3dGLCoordinate)*3*LNumOfLinearSegments);
     }

     if(LGLSegments)
     {
      for(LCnt=0, LGLCnt=0;LCnt<LNCVs;LCnt++, LBezierCV++, LNextBezierCV++)
      {
       if(LPerSegmentSteps) { LLinSteps=LBezierCV->NumOfDrawSteps;E3dM_CurvePrecision(LLinSteps); }

       switch(LBezierCV->Segment)
       {
	case E3dBEZ_CURVE:
	 E3dM_CurveBasis(E3d_BezierMatrix, LBezierCV->Position, LBezierCV->Next, LNextBezierCV->Previous, LNextBezierCV->Position);

	 LGLSegments[LGLCnt]=LGMatrix[M30];
	 LGLSegments[LGLCnt+1]=LGMatrix[M31];
	 LGLSegments[LGLCnt+2]=LGMatrix[M32];

	 for(LSCnt=0;LSCnt<LLinSteps;LSCnt++)
	 {
	  LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
	  LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
	  LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];

	  LGLCnt+=3;
	  LGLSegments[LGLCnt]=LGMatrix[M30];
	  LGLSegments[LGLCnt+1]=LGMatrix[M31];
	  LGLSegments[LGLCnt+2]=LGMatrix[M32];
	 }
	break;

	case E3dBEZ_LINEAR:
	 LGLSegments[LGLCnt]=LBezierCV->Position.X;
	 LGLSegments[LGLCnt+1]=LBezierCV->Position.Y;
	 LGLSegments[LGLCnt+2]=LBezierCV->Position.Z;
	 LGLCnt+=3;
	break;
       }
      }

      if(LSpline->Closed)
      {
       switch(LBezierCV->Segment)
       {
	case E3dBEZ_CURVE:
         if(LPerSegmentSteps) { LLinSteps=LBezierCV->NumOfDrawSteps;E3dM_CurvePrecision(LLinSteps); }

	 E3dM_CurveBasis(E3d_BezierMatrix, LBezierCV->Position, LBezierCV->Next, LBezierCVs[0].Previous, LBezierCVs[0].Position);

	 LGLSegments[LGLCnt]=LGMatrix[M30];
	 LGLSegments[LGLCnt+1]=LGMatrix[M31];
	 LGLSegments[LGLCnt+2]=LGMatrix[M32];

	 for(LSCnt=1;LSCnt<LLinSteps;LSCnt++)
	 {
	  LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
	  LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
	  LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];

	  LGLCnt+=3;

// Back to the beginning...
//
	  LGLSegments[LGLCnt]=LGMatrix[M30];
	  LGLSegments[LGLCnt+1]=LGMatrix[M31];
	  LGLSegments[LGLCnt+2]=LGMatrix[M32];
	 }
	break;

	case E3dBEZ_LINEAR:
	 LGLSegments[LGLCnt]=LBezierCV->Position.X;
	 LGLSegments[LGLCnt+1]=LBezierCV->Position.Y;
	 LGLSegments[LGLCnt+2]=LBezierCV->Position.Z;
	break;
       }
      }
      else
      {
       LGLSegments[LGLCnt]=LBezierCV->Position.X;
       LGLSegments[LGLCnt+1]=LBezierCV->Position.Y;
       LGLSegments[LGLCnt+2]=LBezierCV->Position.Z;
      }
      LSpline->NumOfLinSegments=LNumOfLinearSegments;


//Printf("Bez-LinSeg: %d <-> %d", LNumOfLinearSegments, LGLCnt/3);fflush(stdout);
      return(TRUE);
     }
     else LSpline->NumOfLinSegments=0;
    }
   break;

   case E3dSPL_CARDINAL:
   return(_TesselateSplineForDrawing(E3d_CardinalMatrix, LSpline));

   case E3dSPL_BSPLINE:
   return(_TesselateSplineForDrawing(E3d_BSplineMatrix, LSpline));
  }
 }
 return(FALSE);
}


#define E3dM_LinearSegmentLength()\
{\
 LDX=LNextCV->Position.X-LCurrentCV->Position.X;\
 LDY=LNextCV->Position.Y-LCurrentCV->Position.Y;\
 LDZ=LNextCV->Position.Z-LCurrentCV->Position.Z;\
\
 LCurrentCV->SegmentLength=sqrt(LDX*LDX+LDY*LDY+LDZ*LDZ);\
}


//========================================================================
// Tesselate a BSpline or a Cardinal Spline for arc-length computation
//========================================================================
static E3dCoordinate _TesselateSplineForLength(E3dMatrix LBaseMatrix, E3dSpline* LSpline, unsigned int LLinSteps)
{
 E3dMatrix		LPrecMatrix, LGMatrix;				// Precision matrix, Geometry Matrix (CVs)
 E3dSplineCV*		LSplineCVs=LSpline->CVs;
 E3dSplineCV*		LSplineCV;
 E3dVector*		LLinSegmentForLength;
 E3dCoordinate		LStepsF, LStepsF2, LStepsF3;
 E3dCoordinate		LDX, LDY, LDZ, LTotalLength=0.0;
 E3dCoordinate		LCurrentX, LCurrentY, LCurrentZ, LPrevX, LPrevY, LPrevZ,
			LSegmentLength;
 unsigned int		LCnt, LSCnt, LNumOfCVs=LSpline->NumOfCVs, LNCVs;


 LNCVs=LNumOfCVs-2;

 if(LLinSteps==0) LLinSteps=LSpline->NumOfStepsForLength;
 
 E3dM_CurvePrecision(LLinSteps);

 LSplineCV=LSplineCVs+1;
 for(LCnt=1;LCnt<LNCVs;LCnt++, LSplineCV++)
 {
  if(LSplineCV->NumOfStepsForLength!=(LLinSteps+1))
  {
   if(LSplineCV->StepsForLength) { EFree(LSplineCV->StepsForLength);LSplineCV->StepsForLength=NULL; }
  }
  if(LSplineCV->StepsForLength==NULL) LSplineCV->StepsForLength=(E3dVector*)EMalloc(sizeof(E3dVector)*(LLinSteps+1));
  LSplineCV->NumOfStepsForLength=LLinSteps;
  LLinSegmentForLength=LSplineCV->StepsForLength;


  LSegmentLength=0.0;
  E3dM_CurveBasis(LBaseMatrix, LSplineCVs[LCnt-1].Position, LSplineCV->Position, LSplineCVs[LCnt+1].Position, LSplineCVs[LCnt+2].Position);

  LPrevX=LGMatrix[M30];LPrevY=LGMatrix[M31];LPrevZ=LGMatrix[M32];
  LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
  LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
  LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
  LCurrentX=LGMatrix[M30];LCurrentY=LGMatrix[M31];LCurrentZ=LGMatrix[M32];

  for(LSCnt=1;LSCnt<LLinSteps;LSCnt++)
  {
   LDX=LCurrentX-LPrevX;LDY=LCurrentY-LPrevY;LDZ=LCurrentZ-LPrevZ;
   LSegmentLength+=sqrt(LDX*LDX+LDY*LDY+LDZ*LDZ);

   LPrevX=LCurrentX;LPrevY=LCurrentY;LPrevZ=LCurrentZ;

   LLinSegmentForLength->X=LCurrentX;LLinSegmentForLength->Y=LCurrentY;LLinSegmentForLength->Z=LCurrentZ;
   LLinSegmentForLength->Length=LSegmentLength;
   LLinSegmentForLength++;

   LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
   LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
   LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
   LCurrentX=LGMatrix[M30];LCurrentY=LGMatrix[M31];LCurrentZ=LGMatrix[M32];
  }

  LSplineCV->SegmentLength=LSegmentLength;
  LTotalLength+=LSegmentLength;
 }

 if(LSpline->Closed)
 {
  if(LSplineCV->NumOfStepsForLength!=LLinSteps)
  {
   if(LSplineCV->StepsForLength) { EFree(LSplineCV->StepsForLength);LSplineCV->StepsForLength=NULL; }
  }
  if(LSplineCV->StepsForLength==NULL) LSplineCV->StepsForLength=(E3dVector*)EMalloc(sizeof(E3dVector)*LLinSteps);
  LSplineCV->NumOfStepsForLength=LLinSteps;
  LLinSegmentForLength=LSplineCV->StepsForLength;

  LSegmentLength=0.0;
  E3dM_CurveBasis(LBaseMatrix, LSplineCVs[LNumOfCVs-3].Position, LSplineCV->Position, LSplineCVs[LNumOfCVs-1].Position, LSplineCVs[0].Position);

  LPrevX=LGMatrix[M30];LPrevY=LGMatrix[M31];LPrevZ=LGMatrix[M32];
  LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
  LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
  LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
  LCurrentX=LGMatrix[M30];LCurrentY=LGMatrix[M31];LCurrentZ=LGMatrix[M32];

  for(LSCnt=1;LSCnt<LLinSteps;LSCnt++)
  {
   LDX=LCurrentX-LPrevX;LDY=LCurrentY-LPrevY;LDZ=LCurrentZ-LPrevZ;
   LSegmentLength+=sqrt(LDX*LDX+LDY*LDY+LDZ*LDZ);

   LPrevX=LCurrentX;LPrevY=LCurrentY;LPrevZ=LCurrentZ;

   LLinSegmentForLength->X=LCurrentX;LLinSegmentForLength->Y=LCurrentY;LLinSegmentForLength->Z=LCurrentZ;
   LLinSegmentForLength->Length=LSegmentLength;
   LLinSegmentForLength++;

   LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
   LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
   LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
   LCurrentX=LGMatrix[M30];LCurrentY=LGMatrix[M31];LCurrentZ=LGMatrix[M32];
  }
  LSplineCV->SegmentLength=LSegmentLength;
  LTotalLength+=LSegmentLength;


  LSplineCV++;
  if(LSplineCV->NumOfStepsForLength!=LLinSteps)
  {
   if(LSplineCV->StepsForLength) { EFree(LSplineCV->StepsForLength);LSplineCV->StepsForLength=NULL; }
  }
  if(LSplineCV->StepsForLength==NULL) LSplineCV->StepsForLength=(E3dVector*)EMalloc(sizeof(E3dVector)*LLinSteps);
  LSplineCV->NumOfStepsForLength=LLinSteps;
  LLinSegmentForLength=LSplineCV->StepsForLength;

  LSegmentLength=0.0;
  E3dM_CurveBasis(LBaseMatrix, LSplineCVs[LNumOfCVs-2].Position, LSplineCV->Position, LSplineCVs[0].Position, LSplineCVs[1].Position);

  LPrevX=LGMatrix[M30];LPrevY=LGMatrix[M31];LPrevZ=LGMatrix[M32];
  LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
  LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
  LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
  LCurrentX=LGMatrix[M30];LCurrentY=LGMatrix[M31];LCurrentZ=LGMatrix[M32];

  for(LSCnt=1;LSCnt<LLinSteps;LSCnt++)
  {
   LDX=LCurrentX-LPrevX;LDY=LCurrentY-LPrevY;LDZ=LCurrentZ-LPrevZ;
   LSegmentLength+=sqrt(LDX*LDX+LDY*LDY+LDZ*LDZ);

   LPrevX=LCurrentX;LPrevY=LCurrentY;LPrevZ=LCurrentZ;

   LLinSegmentForLength->X=LCurrentX;LLinSegmentForLength->Y=LCurrentY;LLinSegmentForLength->Z=LCurrentZ;
   LLinSegmentForLength->Length=LSegmentLength;
   LLinSegmentForLength++;

   LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
   LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
   LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
   LCurrentX=LGMatrix[M30];LCurrentY=LGMatrix[M31];LCurrentZ=LGMatrix[M32];
  }
  LSplineCV->SegmentLength=LSegmentLength;
  LTotalLength+=LSegmentLength;


  LSplineCV=LSplineCVs;
  if(LSplineCV->NumOfStepsForLength!=LLinSteps)
  {
   if(LSplineCV->StepsForLength) { EFree(LSplineCV->StepsForLength);LSplineCV->StepsForLength=NULL; }
  }
  if(LSplineCV->StepsForLength==NULL) LSplineCV->StepsForLength=(E3dVector*)EMalloc(sizeof(E3dVector)*LLinSteps);
  LSplineCV->NumOfStepsForLength=LLinSteps;
  LLinSegmentForLength=LSplineCV->StepsForLength;

  LSegmentLength=0.0;
  E3dM_CurveBasis(LBaseMatrix, LSplineCVs[LNumOfCVs-1].Position, LSplineCV->Position, LSplineCVs[1].Position, LSplineCVs[2].Position);

  LPrevX=LGMatrix[M30];LPrevY=LGMatrix[M31];LPrevZ=LGMatrix[M32];
  LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
  LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
  LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
  LCurrentX=LGMatrix[M30];LCurrentY=LGMatrix[M31];LCurrentZ=LGMatrix[M32];

  for(LSCnt=1;LSCnt<LLinSteps;LSCnt++)
  {
   LDX=LCurrentX-LPrevX;LDY=LCurrentY-LPrevY;LDZ=LCurrentZ-LPrevZ;
   LSegmentLength+=sqrt(LDX*LDX+LDY*LDY+LDZ*LDZ);

   LPrevX=LCurrentX;LPrevY=LCurrentY;LPrevZ=LCurrentZ;

   LLinSegmentForLength->X=LCurrentX;LLinSegmentForLength->Y=LCurrentY;LLinSegmentForLength->Z=LCurrentZ;
   LLinSegmentForLength->Length=LSegmentLength;
   LLinSegmentForLength++;

   LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
   LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
   LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
   LCurrentX=LGMatrix[M30];LCurrentY=LGMatrix[M31];LCurrentZ=LGMatrix[M32];
  }
  LSplineCV->SegmentLength=LSegmentLength;
  LTotalLength+=LSegmentLength;
 }
 else
 {
// If the Spline is not closed, CV 0, n-2 and n-1 don't have segments associated with them
//
  LSplineCV=LSplineCVs;
  if(LSplineCV->StepsForLength) { EFree(LSplineCV->StepsForLength);LSplineCV->StepsForLength=NULL; }
  LSplineCV->SegmentLength=0.0;

  LSplineCV=LSplineCVs+LNumOfCVs-2;
  if(LSplineCV->StepsForLength) { EFree(LSplineCV->StepsForLength);LSplineCV->StepsForLength=NULL; }
  LSplineCV->SegmentLength=0.0;

  LSplineCV=LSplineCVs+LNumOfCVs-1;
  if(LSplineCV->StepsForLength) { EFree(LSplineCV->StepsForLength);LSplineCV->StepsForLength=NULL; }
  LSplineCV->SegmentLength=0.0;
 }
 return(LTotalLength);
}


//========================================================================================
// (Re)compute segment lengths and the total length of a Spline
//
// Argument
//  E3dSpline*    LSpline          Pointer to the Spline structure
//  unsigned int  LLinSteps        Number of steps for curved segments
//
// Description
//  This function recomputes the length of all the segments in the given Spline and
//  stores the length value on each Spline CV (CV).
//  As any other segment-related information, the segment length is stored on the
//  segment's <B>starting key</B>.
//  Lengths of curved segments are computed with implicit integration. The LLinSteps
//  argument determines the number of steps (samples) within each curved segment.
//  The higher this number, the higher the precision of the computation will be.
//
// Return value	
//  The total length of the Spline.
//========================================================================================
E3dCoordinate E3d_SplineUpdateSegmentLengths(E3dSpline* LSpline, unsigned int LLinSteps)
{
 E3dCoordinate	LDX, LDY, LDZ, LTotalLength;
 unsigned int	LC, LNumOfCVs;


 LTotalLength=0.0;
 if(LSpline->NumOfCVs<2) return(0.0);
 switch(LSpline->SplineType)
 {
  case E3dSPL_LINEAR:
   {
    E3dSplineCV*	LSplineCVs;
    E3dSplineCV*	LCurrentCV;
    E3dSplineCV*	LNextCV;

    LSplineCVs=LSpline->CVs;
    LNumOfCVs=LSpline->NumOfCVs-1;

    LCurrentCV=LSplineCVs;LNextCV=LCurrentCV+1;
    for(LC=0;LC<LNumOfCVs;LC++, LCurrentCV++, LNextCV++)
    {
     E3dM_LinearSegmentLength();
     LTotalLength+=LCurrentCV->SegmentLength;
    }

    if(LSpline->Closed)
    {
     LCurrentCV=LSplineCVs+LNumOfCVs;
     LNextCV=LSplineCVs;

     E3dM_LinearSegmentLength();

     LTotalLength+=LCurrentCV->SegmentLength;
    }
   }
  break;

  case E3dSPL_BEZIER:
   {
    register unsigned int	LCnt, LSCnt;
    register E3dBezierCV*	LBezierCVs=(E3dBezierCV*)(LSpline->CVs);
    register E3dBezierCV*	LBezierCV=LBezierCVs;
    register E3dBezierCV*	LNextBezierCV=LBezierCV+1;
    E3dCoordinate		LCurrentX, LCurrentY, LCurrentZ, LPrevX, LPrevY, LPrevZ,
				LSegmentLength;
    E3dMatrix			LPrecMatrix, LGMatrix;		// Precision matrix, Geometry matrix (control points)
    E3dVector*			LLinSegmentForLength;
    float			LStepsF, LStepsF2, LStepsF3;
    unsigned int		LNCVs;


    LNumOfCVs=LSpline->NumOfCVs;
    LNCVs=LNumOfCVs-1;

    if(LLinSteps==0) LLinSteps=LSpline->NumOfStepsForLength;

    E3dM_CurvePrecision(LLinSteps);

    for(LCnt=0;LCnt<LNCVs;LCnt++, LBezierCV++, LNextBezierCV++)
    {
     if(LBezierCV->NumOfStepsForLength!=LLinSteps)
     {
      if(LBezierCV->StepsForLength) { EFree(LBezierCV->StepsForLength);LBezierCV->StepsForLength=NULL; }
     }
     if(LBezierCV->StepsForLength==NULL) LBezierCV->StepsForLength=(E3dVector*)EMalloc(sizeof(E3dVector)*LLinSteps);
     LBezierCV->NumOfStepsForLength=LLinSteps;
     LLinSegmentForLength=LBezierCV->StepsForLength;

     switch(LBezierCV->Segment)
     {
      case E3dBEZ_CURVE:
       LSegmentLength=0.0;
       E3dM_CurveBasis(E3d_BezierMatrix, LBezierCV->Position, LBezierCV->Next, LNextBezierCV->Previous, LNextBezierCV->Position);

       LPrevX=LGMatrix[M30];LPrevY=LGMatrix[M31];LPrevZ=LGMatrix[M32];
       LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
       LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
       LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
       LCurrentX=LGMatrix[M30];LCurrentY=LGMatrix[M31];LCurrentZ=LGMatrix[M32];

       for(LSCnt=1;LSCnt<LLinSteps;LSCnt++)
       {
	LDX=LCurrentX-LPrevX;LDY=LCurrentY-LPrevY;LDZ=LCurrentZ-LPrevZ;
	LSegmentLength+=sqrt(LDX*LDX+LDY*LDY+LDZ*LDZ);

	LPrevX=LCurrentX;LPrevY=LCurrentY;LPrevZ=LCurrentZ;

        LLinSegmentForLength->X=LCurrentX;LLinSegmentForLength->Y=LCurrentY;LLinSegmentForLength->Z=LCurrentZ;
        LLinSegmentForLength->Length=LSegmentLength;
        LLinSegmentForLength++;

	LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
	LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
	LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
	LCurrentX=LGMatrix[M30];LCurrentY=LGMatrix[M31];LCurrentZ=LGMatrix[M32];
       }
       LBezierCV->SegmentLength=LSegmentLength;
       LTotalLength+=LSegmentLength;
      break;

      case E3dBEZ_LINEAR:
       LDX=LNextBezierCV->Position.X-LBezierCV->Position.X;
       LDY=LNextBezierCV->Position.Y-LBezierCV->Position.Y;
       LDZ=LNextBezierCV->Position.Z-LBezierCV->Position.Z;

       LBezierCV->SegmentLength=sqrt(LDX*LDX+LDY*LDY+LDZ*LDZ);

       LTotalLength+=LBezierCV->SegmentLength;
      break;
     }
    }

    if(LSpline->Closed)
    {
     LBezierCV=LBezierCVs+LNCVs;

     LNextBezierCV=LBezierCVs;
     switch(LNextBezierCV->Segment)
     {
      case E3dBEZ_CURVE:
       if(LBezierCV->NumOfStepsForLength!=LLinSteps)
       {
	if(LBezierCV->StepsForLength) { EFree(LBezierCV->StepsForLength);LBezierCV->StepsForLength=NULL; }
       }
       if(LBezierCV->StepsForLength==NULL) LBezierCV->StepsForLength=(E3dVector*)EMalloc(sizeof(E3dVector)*LLinSteps);
       LBezierCV->NumOfStepsForLength=LLinSteps;
       LLinSegmentForLength=LBezierCV->StepsForLength;

       LSegmentLength=0.0;
       E3dM_CurveBasis(E3d_BezierMatrix, LBezierCV->Position, LBezierCV->Next, LNextBezierCV->Previous, LNextBezierCV->Position);

       LPrevX=LGMatrix[M30];LPrevY=LGMatrix[M31];LPrevZ=LGMatrix[M32];
       LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
       LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
       LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
       LCurrentX=LGMatrix[M30];LCurrentY=LGMatrix[M31];LCurrentZ=LGMatrix[M32];

       for(LSCnt=1;LSCnt<LLinSteps;LSCnt++)
       {
	LDX=LCurrentX-LPrevX;LDY=LCurrentY-LPrevY;LDZ=LCurrentZ-LPrevZ;
	LSegmentLength+=sqrt(LDX*LDX+LDY*LDY+LDZ*LDZ);

	LPrevX=LCurrentX;LPrevY=LCurrentY;LPrevZ=LCurrentZ;

        LLinSegmentForLength->X=LCurrentX;LLinSegmentForLength->Y=LCurrentY;LLinSegmentForLength->Z=LCurrentZ;
        LLinSegmentForLength->Length=LSegmentLength;
        LLinSegmentForLength++;

	LGMatrix[M30]+=LGMatrix[M20];LGMatrix[M31]+=LGMatrix[M21];LGMatrix[M32]+=LGMatrix[M22];LGMatrix[M33]+=LGMatrix[M23];
	LGMatrix[M20]+=LGMatrix[M10];LGMatrix[M21]+=LGMatrix[M11];LGMatrix[M22]+=LGMatrix[M12];LGMatrix[M23]+=LGMatrix[M13];
	LGMatrix[M10]+=LGMatrix[M00];LGMatrix[M11]+=LGMatrix[M01];LGMatrix[M12]+=LGMatrix[M02];LGMatrix[M13]+=LGMatrix[M03];
	LCurrentX=LGMatrix[M30];LCurrentY=LGMatrix[M31];LCurrentZ=LGMatrix[M32];
       }
       LBezierCV->SegmentLength=LSegmentLength;
       LTotalLength+=LSegmentLength;
      break;

      case E3dBEZ_LINEAR:
       LDX=LNextBezierCV->Position.X-LBezierCV->Position.X;
       LDY=LNextBezierCV->Position.Y-LBezierCV->Position.Y;
       LDZ=LNextBezierCV->Position.Z-LBezierCV->Position.Z;

       LBezierCV->SegmentLength=sqrt(LDX*LDX+LDY*LDY+LDZ*LDZ);

       LTotalLength+=LBezierCV->SegmentLength;
      break;
     }
    }
   }
  break;

  case E3dSPL_BSPLINE:
   LTotalLength=_TesselateSplineForLength(E3d_BSplineMatrix, LSpline, LLinSteps);
  break;

  case E3dSPL_CARDINAL:
   LTotalLength=_TesselateSplineForLength(E3d_CardinalMatrix, LSpline, LLinSteps);
  break;
 }
 LSpline->Length=LTotalLength;
 return(LTotalLength);
}



#define E3dM_ArcLenCmp(mIdx)	LStepsForLength[mIdx].Length<LPosOnSegment


//========================================================================================
// Compute 3D position along a Spline
//
// Arguments
//  E3dSpline*     LSpline         Pointer to the Spline structure
//  double         LTime	   Parameter (0.0...1.0)
//  E3d3DPosition* L3DPosition	   Pointer to structure for the return value
//
// Description
//  This function computes the 3D position along the given Spline, defined by the
//  "Time" parameter.
//  The parameter value 0.0 represents the starting point and 1.0 represents the end
//  point of the Spline, while values in-between refer to a point along the Spline.
//  The segment lengths and the total length of the Spline must be up-to-date when
//  calling this function.
//  The postion returned is in the Spline's local coordinate system.
//
// Return value	
//  0 if all went well.
//
// See also
//  E3d_SplineUpdateSegmentLengths
//========================================================================================
int E3d_SplineGet3DPosition(E3dSpline* LSpline, double LTime, E3d3DPosition* L3DPosition)
{
 double		LFP, LNextPos, LPosOnSegment;
 unsigned int	LNumOfCVs, LNCVs,
		LC;


 if(LSpline->Closed)
 {
  if(LTime>1.0)
  {
   int	LTI=LTime;

   LTime=LTime-LTI;	// Get fraction part...
  }
  else if(LTime<0.0)
       {
	int	LTI=LTime;

	LTime=1.0+(LTime-LTI);
       }
 }
 else
 {
  if(LTime>1.0) LTime=1.0;
  else if(LTime<0.0) LTime=0.0;
 }

 if(LSpline->Length<0.0) E3d_SplineUpdateSegmentLengths(LSpline, LSpline->NumOfStepsForLength);


 LFP=LTime*LSpline->Length;
 LNumOfCVs=LSpline->NumOfCVs;
 LNCVs=LNumOfCVs-1;
 switch(LSpline->SplineType)
 {
  case E3dSPL_LINEAR:
   {
    E3dSplineCV*	LSplineCVs;
    E3dSplineCV*	LCurrentCV;
    E3dSplineCV*	LNextCV;

    LCurrentCV=LSplineCVs=LSpline->CVs;
    LNextCV=LCurrentCV+1;

// Find segment that contains this position
//
    LNextPos=LCurrentCV->SegmentLength;
    for(LC=0;LC<LNCVs;)
    {
     if(LNextPos>LFP) break;
     else
     {
      printf("Segment: %f\n", LCurrentCV->SegmentLength);fflush(stdout);

      LNextPos+=LNextCV->SegmentLength;
      LC++;LCurrentCV++;LNextCV++;
     }
    }

    if(LSpline->Closed&&(LC==LNCVs))
    {
     LNextCV=LSplineCVs;
    }

    LPosOnSegment=(LFP-(LNextPos-LCurrentCV->SegmentLength))/LCurrentCV->SegmentLength;
    L3DPosition->X=(LNextCV->Position.X-LCurrentCV->Position.X)*LPosOnSegment+LCurrentCV->Position.X;
    L3DPosition->Y=(LNextCV->Position.Y-LCurrentCV->Position.Y)*LPosOnSegment+LCurrentCV->Position.Y;
    L3DPosition->Z=(LNextCV->Position.Z-LCurrentCV->Position.Z)*LPosOnSegment+LCurrentCV->Position.Z;
   }
  break;

  case E3dSPL_BEZIER:
   {
    register E3dBezierCV*	LBezierCVs=(E3dBezierCV*)(LSpline->CVs);
    register E3dBezierCV*	LBezierCV=LBezierCVs;
    register E3dBezierCV*	LNextBezierCV=LBezierCV+1;
//    double			LNPosOnSegment, LUSq, LUCb, LNUSq, LNUCb;


// Find the segment that contains this position
//
    LNextPos=LBezierCV->SegmentLength;
    for(LC=0;LC<LNCVs;)
    {
     if(LNextPos>LFP) break;
     else
     {
      LNextPos+=LNextBezierCV->SegmentLength;
      LC++;LBezierCV++;LNextBezierCV++;
     }
    }

    if(LSpline->Closed)
    {
     if(LC==LNCVs) LNextBezierCV=LBezierCVs;
    }
    else
    {
     if(LTime>=1.0)
     {
      LBezierCV=LBezierCVs+LSpline->NumOfCVs-1;
      L3DPosition->X=LBezierCV->Position.X;
      L3DPosition->Y=LBezierCV->Position.Y;
      L3DPosition->Z=LBezierCV->Position.Z;
      return(0);
     }
    }


    switch(LBezierCV->Segment)
    {
     case E3dBEZ_CURVE:
      {
       E3dVector*	LStepsForLength=LBezierCV->StepsForLength;
       E3dVector*	LLinSegment;
       E3dVector*	LNextLinSegment;
       E3dCoordinate	LPos;
       int		LIndex;

// Do binary search to find linear segment closest to arc-length parameter
//
// Get position within segment
//
       LPosOnSegment=(LFP-(LNextPos-LBezierCV->SegmentLength));

//printf("LPosOnSegment %f  %d\n", LPosOnSegment, LBezierCV-LBezierCVs);fflush(stdout);

       EM_BinSearch(LBezierCV->NumOfStepsForLength, E3dM_ArcLenCmp);

       LLinSegment=LStepsForLength+LIndex;
       LNextLinSegment=LStepsForLength+LIndex+1;

       LPos=(LPosOnSegment-LLinSegment->Length)/(LNextLinSegment->Length-LLinSegment->Length);


//printf("Pos %f\n", LPos);fflush(stdout);

       L3DPosition->X=(LNextLinSegment->X-LLinSegment->X)*LPos+LLinSegment->X;
       L3DPosition->Y=(LNextLinSegment->Y-LLinSegment->Y)*LPos+LLinSegment->Y;
       L3DPosition->Z=(LNextLinSegment->Z-LLinSegment->Z)*LPos+LLinSegment->Z;
      }

/*
      LNPosOnSegment=1.0-LPosOnSegment;
      LUSq=LPosOnSegment*LPosOnSegment;
      LUCb=LUSq*LPosOnSegment;
      LNUSq=LNPosOnSegment*LNPosOnSegment;
      LNUCb=LNUSq*LNPosOnSegment;

      L3DPosition->X=LNUCb*LBezierCV->Position.X + 3.0*LPosOnSegment*LNUSq*LBezierCV->Next.X + 3.0*LUSq*LNPosOnSegment*LNextBezierCV->Previous.X + LUCb*LNextBezierCV->Position.X;
      L3DPosition->Y=LNUCb*LBezierCV->Position.Y + 3.0*LPosOnSegment*LNUSq*LBezierCV->Next.Y + 3.0*LUSq*LNPosOnSegment*LNextBezierCV->Previous.Y + LUCb*LNextBezierCV->Position.Y;
      L3DPosition->Z=LNUCb*LBezierCV->Position.Z + 3.0*LPosOnSegment*LNUSq*LBezierCV->Next.Z + 3.0*LUSq*LNPosOnSegment*LNextBezierCV->Previous.Z + LUCb*LNextBezierCV->Position.Z;
*/
     break;

     case E3dBEZ_LINEAR:
      LPosOnSegment=(LFP-(LNextPos-LBezierCV->SegmentLength))/LBezierCV->SegmentLength;
      L3DPosition->X=(LNextBezierCV->Position.X-LBezierCV->Position.X)*LPosOnSegment+LBezierCV->Position.X;
      L3DPosition->Y=(LNextBezierCV->Position.Y-LBezierCV->Position.Y)*LPosOnSegment+LBezierCV->Position.Y;
      L3DPosition->Z=(LNextBezierCV->Position.Z-LBezierCV->Position.Z)*LPosOnSegment+LBezierCV->Position.Z;
     break;
    }
   }
  break;

  case E3dSPL_BSPLINE:
  case E3dSPL_CARDINAL:
   {
    register E3dSplineCV*	LSplineCVs=(E3dSplineCV*)(LSpline->CVs);
    register E3dSplineCV*	LSplineCV=LSplineCVs;
    register E3dSplineCV*	LNextSplineCV=LSplineCV+1;


// Find the segment that contains this position
//
    LNextPos=LSplineCV->SegmentLength;
    for(LC=0;LC<LNCVs;)
    {
     if(LNextPos>LFP) break;
     else
     {
      LNextPos+=LNextSplineCV->SegmentLength;
      LC++;LSplineCV++;LNextSplineCV++;
     }
    }

    if(LSpline->Closed)
    {
     if(LC==LNCVs) LNextSplineCV=LSplineCVs;
    }
    else
    {
     if(LTime>=1.0)
     {
      if(LSpline->SplineType==E3dSPL_CARDINAL)
      {
       LSplineCV=LSplineCVs+LSpline->NumOfCVs-2;
       L3DPosition->X=LSplineCV->Position.X;
       L3DPosition->Y=LSplineCV->Position.Y;
       L3DPosition->Z=LSplineCV->Position.Z;
      }
      else
      {
       LSplineCV=LSplineCVs+LSpline->NumOfCVs-3;
       if(LSplineCV->StepsForLength)
       {
        E3dVector*	LLastPos=LSplineCV->StepsForLength+LSplineCV->NumOfStepsForLength-2;

        L3DPosition->X=LLastPos->X;
        L3DPosition->Y=LLastPos->Y;
        L3DPosition->Z=LLastPos->Z;
       }
      }
      return(0);
     }
    }


    {
     E3dVector*		LStepsForLength=LSplineCV->StepsForLength;
     E3dVector*		LLinSegment;
     E3dVector*		LNextLinSegment;
     E3dCoordinate	LPos;
     int		LIndex;

// Do binary search to find linear segment closest to arc-length parameter
//
// Get position within segment
//
     LPosOnSegment=(LFP-(LNextPos-LSplineCV->SegmentLength));

     EM_BinSearch(LSplineCV->NumOfStepsForLength, E3dM_ArcLenCmp);

     LLinSegment=LStepsForLength+LIndex;
     LNextLinSegment=LStepsForLength+LIndex+1;

     LPos=(LPosOnSegment-LLinSegment->Length)/(LNextLinSegment->Length-LLinSegment->Length);


//printf("Pos %f\n", LPos);fflush(stdout);

     L3DPosition->X=(LNextLinSegment->X-LLinSegment->X)*LPos+LLinSegment->X;
     L3DPosition->Y=(LNextLinSegment->Y-LLinSegment->Y)*LPos+LLinSegment->Y;
     L3DPosition->Z=(LNextLinSegment->Z-LLinSegment->Z)*LPos+LLinSegment->Z;
    }
   }
  break;
 }

 return(0);
}



//========================================
// Unselect components of a Spline
//========================================
void E3d_SplineUnselect(E3dSpline* LSpline, int LFlags)
{
 LSpline->Selection=E3dGSEL_NONE;
 if(LFlags&E3dSF_SPLINE_SEGMENTS)
 {
  unsigned int	LC, LN=LSpline->NumOfCVs;

  switch(LSpline->SplineType)
  {
   case E3dSPL_BEZIER:
    {
     E3dBezierCV*	LBezierCV=(E3dBezierCV*)(LSpline->CVs);

     for(LC=0;LC<LN;LC++, LBezierCV++)
     {
      LBezierCV->Flags&=(E3dKEYF_ALL-(E3dKEYF_SEGMENT_SELECTED|E3dKEYF_ACTIVE));
     }
    }
   break;

   default:
    {
     E3dSplineCV*	LSplineCV=(E3dSplineCV*)(LSpline->CVs);

     for(LC=0;LC<LN;LC++, LSplineCV++)
     {
      LSplineCV->Flags&=(E3dKEYF_ALL-(E3dKEYF_SEGMENT_SELECTED|E3dKEYF_ACTIVE));
     }
     LSpline->Selection=E3dGSEL_NONE;
    }
   break;
  }
 }
}



#define M3M_RotV(lxv, lyv, lzv)	\
 if((LDZ==0.0)||(LDXP==0.0))	\
 {				\
  LX=(lxv)*LDX-(lyv)*LDY,	\
  LY=(lxv)*LDY+(lyv)*LDX,	\
  LXR=LX+LX1;LYR=LY+LY1;LZR=(lzv)*LLen+LZ1;\
 }				\
 else				\
 {				\
  LX=(lxv)*LDXP-(lyv)*LDY,	\
  LY=(lxv)*LDY+(lyv)*LDXP,	\
  LZ=(lzv)*LLen,		\
  LXR=LX*LZRCos-LZ*LZRSin+LX1,	\
  LYR=LY+LY1,			\
  LZR=LX*LZRSin+LZ*LZRCos+LZ1;	\
 }


//========================================================================
// Vector primitive
//========================================================================

//========================================================
// Create the arrow on a vector primitive Spline
//========================================================
void E3d_CreateVectorArrow(E3dSplineCV* LCV, double LX2, double LY2, double LZ2)
{
 double		LX, LY, LZ, LX1, LY1, LZ1, LDX, LDY, LDZ, LDXP, LXR, LYR, LZR,
		LZRSin=0.0, LZRCos=0.0,
		LLen;


 LDX=LX2-(LX1=LCV[0].Position.X);
 LDY=LY2-(LY1=LCV[0].Position.Y);
 LDZ=LZ2-(LZ1=LCV[0].Position.Z);
 LLen=sqrt(LDX*LDX+LDY*LDY+LDZ*LDZ);
 LCV[1].Position.X=LCV[5].Position.X=LCV[7].Position.X=LCV[9].Position.X=0.92*LDX+LX1;
 LCV[1].Position.Y=LCV[5].Position.Y=LCV[7].Position.Y=LCV[9].Position.Y=0.92*LDY+LY1;
 LCV[1].Position.Z=LCV[5].Position.Z=LCV[7].Position.Z=LCV[9].Position.Z=0.92*LDZ+LZ1;

 if((LDXP=LDX*LDX+LDZ*LDZ)!=0.0)
 {
  LDXP=sqrt(LDXP);if(LDX<0.0) LDXP*=-1.0;
  LZRSin=LDZ/LDXP;LZRCos=LDX/LDXP;
 }

 M3M_RotV(0.9, -0.04, 0.0);LCV[2].Position.X=LXR;LCV[2].Position.Y=LYR;LCV[2].Position.Z=LZR;

 LCV[3].Position.X=LCV[7].Position.X=LX2;
 LCV[3].Position.Y=LCV[7].Position.Y=LY2;
 LCV[3].Position.Z=LCV[7].Position.Z=LZ2;

 M3M_RotV(0.9, 0.04, 0.0);LCV[4].Position.X=LXR;LCV[4].Position.Y=LYR;LCV[4].Position.Z=LZR;

 M3M_RotV(0.9, 0.0, 0.04);LCV[6].Position.X=LXR;LCV[6].Position.Y=LYR;LCV[6].Position.Z=LZR;

 M3M_RotV(0.9, 0.0, -0.04);LCV[8].Position.X=LXR;LCV[8].Position.Y=LYR;LCV[8].Position.Z=LZR;
}


//========================================================
// Create the Model for the vector (linear Spline)
//========================================================
E3dModel* E3d_VectorTest(double LX1, double LY1, double LZ1, double LX2, double LY2, double LZ2, EBool LCreateArrow)
{
 E3dModel*	LModel;


 if(E3d_VectorGeometryClass==NULL)
 {
  E3dGeometryClass	LGeoClass=
   {
    "VectorTest",			// Name
    sizeof(E3dVectorInfo),		// StructSize
    NULL,				// EditProc
    NULL,				// DestroyProc
    NULL				// Resources
   };

  if((E3d_VectorGeometryClass=E3d_GeometryClassRegister(&LGeoClass))==NULL) return(NULL);
 }

 if((LModel=E3d_ModelAllocate("Vector"))!=NULL)
 {
  E3dVectorInfo*	LVectorInfo;
  E3dSpline*		LSpline=E3d_SplineAllocate(E3dSPL_LINEAR);
  E3dSplineCV*		LCV;


  if(LSpline)
  {
   if(E3d_ModelAppendGeometry(LModel, (E3dGeometry*)LSpline)==NULL) { E3d_ModelHrcFree(LModel, TRUE);return(NULL); }
  }
  else { E3d_ModelHrcFree(LModel, TRUE);return(NULL); }

  LSpline->LockCount=1;

  LVectorInfo=(E3dVectorInfo*)E3d_GeometryInfoAdd((E3dGeometry*)LSpline, E3d_VectorGeometryClass);
  if(LVectorInfo)
  {
   E3d_VectorInit(&(LVectorInfo->Vector), LX2-LX1, LY2-LY1, LZ2-LZ1);
   E3d_3DPositionInit(&(LVectorInfo->Position), LX1, LY1, LZ1);
   LVectorInfo->ArrowOn=TRUE;
   LVectorInfo->CoordinateMode=E3dVTM_P1_P2;
  }

  if((LSpline->CVs=E3d_SplineCVsAllocate(E3dSPL_LINEAR, 10))!=NULL)
  {
   LCV=LSpline->CVs;
   LCV[0].Position.X=LX1;LCV[0].Position.Y=LY1;LCV[0].Position.Z=LZ1;
   if(LCreateArrow) { LSpline->NumOfCVs=10;E3d_CreateVectorArrow(LCV, LX2, LY2, LZ2); }
   else
   {
    LSpline->NumOfCVs=2;
    LCV[1].Position.X=LX2;LCV[1].Position.Y=LY2;LCV[1].Position.Z=LZ2;
   }
   E3d_SplineCreateLinearSegments(LSpline);
  }
 }
 return(LModel);
}
