[occt.git] / src / BSplCLib / BSplCLib_Cache.cxx

// Copyright (c) 2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

#include <BSplCLib_Cache.hxx>
#include <BSplCLib.hxx>

#include <NCollection_LocalArray.hxx>

#include <TColgp_HArray1OfPnt.hxx>
#include <TColgp_HArray1OfPnt2d.hxx>
#include <TColStd_HArray1OfReal.hxx>
#include <TColStd_HArray2OfReal.hxx>

IMPLEMENT_STANDARD_HANDLE(BSplCLib_Cache, Standard_Transient)
IMPLEMENT_STANDARD_RTTIEXT(BSplCLib_Cache, Standard_Transient)

//! Converts handle of array of Standard_Real into the pointer to Standard_Real
static Standard_Real* ConvertArray(const Handle_TColStd_HArray2OfReal& theHArray)
{
  const TColStd_Array2OfReal& anArray = theHArray->Array2();
  return (Standard_Real*) &(anArray(anArray.LowerRow(), anArray.LowerCol()));
}


BSplCLib_Cache::BSplCLib_Cache()
{
  myPolesWeights.Nullify();
  myIsRational = Standard_False;
  mySpanStart = 0.0;
  mySpanLength = 0.0;
  mySpanIndex = 0;
  myDegree = 0;
  myFlatKnots.Nullify();
}

BSplCLib_Cache::BSplCLib_Cache(const Standard_Integer&        theDegree,
                               const Standard_Boolean&        thePeriodic,
                               const TColStd_Array1OfReal&    theFlatKnots,
                               const TColgp_Array1OfPnt2d&    thePoles2d,
                               const TColStd_Array1OfReal&    theWeights)
{
  Standard_Real aCacheParam = theFlatKnots.Value(theFlatKnots.Lower() + theDegree);
  BuildCache(aCacheParam, theDegree, thePeriodic, 
             theFlatKnots, thePoles2d, theWeights);
}

BSplCLib_Cache::BSplCLib_Cache(const Standard_Integer&        theDegree,
                               const Standard_Boolean&        thePeriodic,
                               const TColStd_Array1OfReal&    theFlatKnots,
                               const TColgp_Array1OfPnt&      thePoles,
                               const TColStd_Array1OfReal&    theWeights)
{
  Standard_Real aCacheParam = theFlatKnots.Value(theFlatKnots.Lower() + theDegree);
  BuildCache(aCacheParam, theDegree, thePeriodic, 
             theFlatKnots, thePoles, theWeights);
}


Standard_Boolean BSplCLib_Cache::IsCacheValid(Standard_Real theParameter) const
{
  Standard_Real aNewParam = theParameter;
  if (!myFlatKnots.IsNull())
    PeriodicNormalization(myFlatKnots->Array1(), aNewParam);

  Standard_Real aDelta = aNewParam - mySpanStart;
  return (aDelta >= 0.0 && (aDelta < mySpanLength || mySpanIndex == mySpanIndexMax));
}

void BSplCLib_Cache::PeriodicNormalization(const TColStd_Array1OfReal& theFlatKnots, 
                                           Standard_Real& theParameter) const
{
  Standard_Real aPeriod = theFlatKnots.Value(theFlatKnots.Upper() - myDegree) - 
                          theFlatKnots.Value(myDegree + 1) ;
  if (theParameter < theFlatKnots.Value(myDegree + 1))
  {
    Standard_Real aScale = IntegerPart(
        (theFlatKnots.Value(myDegree + 1) - theParameter) / aPeriod);
    theParameter += aPeriod * (aScale + 1.0);
  }
  if (theParameter > theFlatKnots.Value(theFlatKnots.Upper() - myDegree))
  {
    Standard_Real aScale = IntegerPart(
        (theParameter - theFlatKnots.Value(theFlatKnots.Upper() - myDegree)) / aPeriod);
    theParameter -= aPeriod * (aScale + 1.0);
  }
}


void BSplCLib_Cache::BuildCache(const Standard_Real&           theParameter,
                                const Standard_Integer&        theDegree,
                                const Standard_Boolean&        thePeriodic,
                                const TColStd_Array1OfReal&    theFlatKnots,
                                const TColgp_Array1OfPnt2d&    thePoles2d,
                                const TColStd_Array1OfReal&    theWeights)
{
  // Normalize theParameter for periodical B-splines
  Standard_Real aNewParam = theParameter;
  if (thePeriodic)
  {
    PeriodicNormalization(theFlatKnots, aNewParam);
    myFlatKnots = new TColStd_HArray1OfReal(1, theFlatKnots.Length());
    myFlatKnots->ChangeArray1() = theFlatKnots;
  }
  else if (!myFlatKnots.IsNull()) // Periodical curve became non-periodical
    myFlatKnots.Nullify();

  // Change the size of cached data if needed
  myIsRational = (&theWeights != NULL);
  Standard_Integer aPWColNumber = myIsRational ? 3 : 2;
  if (theDegree > myDegree)
    myPolesWeights = new TColStd_HArray2OfReal(1, theDegree + 1, 1, aPWColNumber);

  myDegree = theDegree;
  mySpanIndex = 0;
  BSplCLib::LocateParameter(theDegree, theFlatKnots, BSplCLib::NoMults(), 
                            aNewParam, thePeriodic, mySpanIndex, aNewParam);
  mySpanStart  = theFlatKnots.Value(mySpanIndex);
  mySpanLength = theFlatKnots.Value(mySpanIndex + 1) - mySpanStart;
  mySpanIndexMax = theFlatKnots.Length() - 1 - theDegree;

  // Calculate new cache data
  BSplCLib::BuildCache(mySpanStart, mySpanLength, thePeriodic, theDegree, 
                       theFlatKnots, thePoles2d, theWeights, 
                       myPolesWeights->ChangeArray2());
}

void BSplCLib_Cache::BuildCache(const Standard_Real&           theParameter,
                                const Standard_Integer&        theDegree,
                                const Standard_Boolean&        thePeriodic,
                                const TColStd_Array1OfReal&    theFlatKnots,
                                const TColgp_Array1OfPnt&      thePoles,
                                const TColStd_Array1OfReal&    theWeights)
{
  // Create list of knots with repetitions and normalize theParameter for periodical B-splines
  Standard_Real aNewParam = theParameter;
  if (thePeriodic)
  {
    PeriodicNormalization(theFlatKnots, aNewParam);
    myFlatKnots = new TColStd_HArray1OfReal(1, theFlatKnots.Length());
    myFlatKnots->ChangeArray1() = theFlatKnots;
  }
  else if (!myFlatKnots.IsNull()) // Periodical curve became non-periodical
    myFlatKnots.Nullify();

  // Change the size of cached data if needed
  myIsRational = (&theWeights != NULL);
  Standard_Integer aPWColNumber = myIsRational ? 4 : 3;
  if (theDegree > myDegree)
    myPolesWeights = new TColStd_HArray2OfReal(1, theDegree + 1, 1, aPWColNumber);

  myDegree = theDegree;
  mySpanIndex = 0;
  BSplCLib::LocateParameter(theDegree, theFlatKnots, BSplCLib::NoMults(), 
                            aNewParam, thePeriodic, mySpanIndex, aNewParam);
  mySpanStart  = theFlatKnots.Value(mySpanIndex);
  mySpanLength = theFlatKnots.Value(mySpanIndex + 1) - mySpanStart;
  mySpanIndexMax = theFlatKnots.Length() - 1 - theDegree;

  // Calculate new cache data
  BSplCLib::BuildCache(mySpanStart, mySpanLength, thePeriodic, theDegree, 
                       theFlatKnots, thePoles, theWeights, 
                       myPolesWeights->ChangeArray2());
}


void BSplCLib_Cache::CalculateDerivative(const Standard_Real&    theParameter, 
                                         const Standard_Integer& theDerivative, 
                                               Standard_Real&    theDerivArray) const
{
  Standard_Real aNewParameter = theParameter;
  if (!myFlatKnots.IsNull()) // B-spline is periodical
    PeriodicNormalization(myFlatKnots->Array1(), aNewParameter);
  aNewParameter = (aNewParameter - mySpanStart) / mySpanLength;

  Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
  Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns

  // Temporary container. The maximal size of this container is defined by:
  //    1) maximal derivative for cache evaluation, which is 3, plus one row for function values, 
  //    2) and maximal dimension of the point, which is 3, plus one column for weights.
  Standard_Real aTmpContainer[16];

  // When the PLib::RationaDerivative needs to be called, use temporary container
  Standard_Real* aPntDeriv = myIsRational ? aTmpContainer : &theDerivArray;

  // When the degree of curve is lesser than the requested derivative,
  // nullify array cells corresponding to greater derivatives
  Standard_Integer aDerivative = theDerivative;
  if (myDegree < theDerivative)
  {
    aDerivative = myDegree;
    for (Standard_Integer ind = myDegree * aDimension; ind < (theDerivative + 1) * aDimension; ind++)
    {
      aPntDeriv[ind] = 0.0;
      (&theDerivArray)[ind] = 0.0; // should be cleared separately, because aPntDeriv may look to another memory area
    }
  }

  PLib::EvalPolynomial(aNewParameter, aDerivative, myDegree, aDimension, 
                       aPolesArray[0], aPntDeriv[0]);
  // Unnormalize derivatives since those are computed normalized
  Standard_Real aFactor = 1.0;
  for (Standard_Integer deriv = 1; deriv <= aDerivative; deriv++)
  {
    aFactor /= mySpanLength;
    for (Standard_Integer ind = 0; ind < aDimension; ind++)
      aPntDeriv[aDimension * deriv + ind] *= aFactor;
  }

  if (myIsRational) // calculate derivatives divided by weights derivatives
    PLib::RationalDerivative(aDerivative, aDerivative, aDimension - 1, aPntDeriv[0], theDerivArray);
}


void BSplCLib_Cache::D0(const Standard_Real& theParameter, gp_Pnt2d& thePoint) const
{
  Standard_Real aNewParameter = theParameter;
  if (!myFlatKnots.IsNull()) // B-spline is periodical
    PeriodicNormalization(myFlatKnots->Array1(), aNewParameter);
  aNewParameter = (aNewParameter - mySpanStart) / mySpanLength;

  Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
  Standard_Real aPoint[4];
  Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns

  PLib::NoDerivativeEvalPolynomial(aNewParameter, myDegree,
                                   aDimension, myDegree * aDimension,
                                   aPolesArray[0], aPoint[0]);

  thePoint.SetCoord(aPoint[0], aPoint[1]);
  if (myIsRational)
    thePoint.ChangeCoord().Divide(aPoint[2]);
}

void BSplCLib_Cache::D0(const Standard_Real& theParameter, gp_Pnt& thePoint) const
{
  Standard_Real aNewParameter = theParameter;
  if (!myFlatKnots.IsNull()) // B-spline is periodical
    PeriodicNormalization(myFlatKnots->Array1(), aNewParameter);
  aNewParameter = (aNewParameter - mySpanStart) / mySpanLength;

  Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
  Standard_Real aPoint[4];
  Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns

  PLib::NoDerivativeEvalPolynomial(aNewParameter, myDegree,
                                   aDimension, myDegree * aDimension,
                                   aPolesArray[0], aPoint[0]);

  thePoint.SetCoord(aPoint[0], aPoint[1], aPoint[2]);
  if (myIsRational)
    thePoint.ChangeCoord().Divide(aPoint[3]);
}


void BSplCLib_Cache::D1(const Standard_Real& theParameter, gp_Pnt2d& thePoint, gp_Vec2d& theTangent) const
{
  Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
  Standard_Real aPntDeriv[8]; // result storage (point and derivative coordinates)

  this->CalculateDerivative(theParameter, 1, aPntDeriv[0]);
  if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
    aDimension -= 1;

  thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1]);
  theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1]);
}

void BSplCLib_Cache::D1(const Standard_Real& theParameter, gp_Pnt& thePoint, gp_Vec& theTangent) const
{
  Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
  Standard_Real aPntDeriv[8]; // result storage (point and derivative coordinates)

  this->CalculateDerivative(theParameter, 1, aPntDeriv[0]);
  if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
    aDimension -= 1;

  thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1], aPntDeriv[2]);
  theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1], aPntDeriv[aDimension + 2]);
}

void BSplCLib_Cache::D2(const Standard_Real& theParameter, gp_Pnt2d& thePoint, gp_Vec2d& theTangent, gp_Vec2d& theCurvature) const
{
  Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
  Standard_Real aPntDeriv[12]; // result storage (point and derivatives coordinates)

  this->CalculateDerivative(theParameter, 2, aPntDeriv[0]);
  if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
    aDimension -= 1;

  thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1]);
  theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1]);
  theCurvature.SetCoord(aPntDeriv[aDimension<<1], aPntDeriv[(aDimension<<1) + 1]);
}

void BSplCLib_Cache::D2(const Standard_Real& theParameter, gp_Pnt& thePoint, gp_Vec& theTangent, gp_Vec& theCurvature) const
{
  Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
  Standard_Real aPntDeriv[12]; // result storage (point and derivatives coordinates)

  this->CalculateDerivative(theParameter, 2, aPntDeriv[0]);
  if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
    aDimension -= 1;

  thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1], aPntDeriv[2]);
  theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1], aPntDeriv[aDimension + 2]);
  theCurvature.SetCoord(aPntDeriv[aDimension<<1], aPntDeriv[(aDimension<<1) + 1], aPntDeriv[(aDimension<<1) + 2]);
}


void BSplCLib_Cache::D3(const Standard_Real& theParameter, 
                              gp_Pnt2d&      thePoint, 
                              gp_Vec2d&      theTangent, 
                              gp_Vec2d&      theCurvature,
                              gp_Vec2d&      theTorsion) const
{
  Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
  Standard_Real aPntDeriv[16]; // result storage (point and derivatives coordinates)

  this->CalculateDerivative(theParameter, 3, aPntDeriv[0]);
  if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
    aDimension -= 1;

  thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1]);
  theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1]);
  Standard_Integer aShift = aDimension << 1;
  theCurvature.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1]);
  aShift += aDimension;
  theTorsion.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1]);
}

void BSplCLib_Cache::D3(const Standard_Real& theParameter, 
                              gp_Pnt&        thePoint, 
                              gp_Vec&        theTangent, 
                              gp_Vec&        theCurvature,
                              gp_Vec&        theTorsion) const
{
  Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
  Standard_Real aPntDeriv[16]; // result storage (point and derivatives coordinates)

  this->CalculateDerivative(theParameter, 3, aPntDeriv[0]);
  if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
    aDimension -= 1;

  thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1], aPntDeriv[2]);
  theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1], aPntDeriv[aDimension + 2]);
  Standard_Integer aShift = aDimension << 1;
  theCurvature.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1], aPntDeriv[aShift + 2]);
  aShift += aDimension;
  theTorsion.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1], aPntDeriv[aShift + 2]);
}
Commit	Line	Data
94f71cad	1	// Copyright (c) 2014 OPEN CASCADE SAS
	2	//
	3	// This file is part of Open CASCADE Technology software library.
	4	//
	5	// This library is free software; you can redistribute it and/or modify it under
	6	// the terms of the GNU Lesser General Public License version 2.1 as published
	7	// by the Free Software Foundation, with special exception defined in the file
	8	// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
	9	// distribution for complete text of the license and disclaimer of any warranty.
	10	//
	11	// Alternatively, this file may be used under the terms of Open CASCADE
	12	// commercial license or contractual agreement.
	13
	14	#include <BSplCLib_Cache.hxx>
	15	#include <BSplCLib.hxx>
	16
	17	#include <NCollection_LocalArray.hxx>
	18
	19	#include <TColgp_HArray1OfPnt.hxx>
	20	#include <TColgp_HArray1OfPnt2d.hxx>
	21	#include <TColStd_HArray1OfReal.hxx>
	22	#include <TColStd_HArray2OfReal.hxx>
	23
	24	IMPLEMENT_STANDARD_HANDLE(BSplCLib_Cache, Standard_Transient)
	25	IMPLEMENT_STANDARD_RTTIEXT(BSplCLib_Cache, Standard_Transient)
	26
	27	//! Converts handle of array of Standard_Real into the pointer to Standard_Real
	28	static Standard_Real* ConvertArray(const Handle_TColStd_HArray2OfReal& theHArray)
	29	{
	30	const TColStd_Array2OfReal& anArray = theHArray->Array2();
	31	return (Standard_Real*) &(anArray(anArray.LowerRow(), anArray.LowerCol()));
	32	}
	33
	34
	35	BSplCLib_Cache::BSplCLib_Cache()
	36	{
	37	myPolesWeights.Nullify();
	38	myIsRational = Standard_False;
	39	mySpanStart = 0.0;
	40	mySpanLength = 0.0;
	41	mySpanIndex = 0;
	42	myDegree = 0;
	43	myFlatKnots.Nullify();
	44	}
	45
	46	BSplCLib_Cache::BSplCLib_Cache(const Standard_Integer& theDegree,
	47	const Standard_Boolean& thePeriodic,
	48	const TColStd_Array1OfReal& theFlatKnots,
	49	const TColgp_Array1OfPnt2d& thePoles2d,
	50	const TColStd_Array1OfReal& theWeights)
	51	{
	52	Standard_Real aCacheParam = theFlatKnots.Value(theFlatKnots.Lower() + theDegree);
	53	BuildCache(aCacheParam, theDegree, thePeriodic,
	54	theFlatKnots, thePoles2d, theWeights);
	55	}
	56
	57	BSplCLib_Cache::BSplCLib_Cache(const Standard_Integer& theDegree,
	58	const Standard_Boolean& thePeriodic,
	59	const TColStd_Array1OfReal& theFlatKnots,
	60	const TColgp_Array1OfPnt& thePoles,
	61	const TColStd_Array1OfReal& theWeights)
	62	{
	63	Standard_Real aCacheParam = theFlatKnots.Value(theFlatKnots.Lower() + theDegree);
	64	BuildCache(aCacheParam, theDegree, thePeriodic,
65	theFlatKnots, thePoles, theWeights);
66	}
67
68
69	Standard_Boolean BSplCLib_Cache::IsCacheValid(Standard_Real theParameter) const
70	{
71	Standard_Real aNewParam = theParameter;
72	if (!myFlatKnots.IsNull())
73	PeriodicNormalization(myFlatKnots->Array1(), aNewParam);
74
75	Standard_Real aDelta = aNewParam - mySpanStart;
76	return (aDelta >= 0.0 && (aDelta < mySpanLength \|\| mySpanIndex == mySpanIndexMax));
77	}
78
79	void BSplCLib_Cache::PeriodicNormalization(const TColStd_Array1OfReal& theFlatKnots,
80	Standard_Real& theParameter) const
81	{
82	Standard_Real aPeriod = theFlatKnots.Value(theFlatKnots.Upper() - myDegree) -
83	theFlatKnots.Value(myDegree + 1) ;
84	if (theParameter < theFlatKnots.Value(myDegree + 1))
85	{
86	Standard_Real aScale = IntegerPart(
87	(theFlatKnots.Value(myDegree + 1) - theParameter) / aPeriod);
88	theParameter += aPeriod * (aScale + 1.0);
89	}
90	if (theParameter > theFlatKnots.Value(theFlatKnots.Upper() - myDegree))
91	{
92	Standard_Real aScale = IntegerPart(
93	(theParameter - theFlatKnots.Value(theFlatKnots.Upper() - myDegree)) / aPeriod);
94	theParameter -= aPeriod * (aScale + 1.0);
95	}
96	}
97
98
99	void BSplCLib_Cache::BuildCache(const Standard_Real& theParameter,
100	const Standard_Integer& theDegree,
101	const Standard_Boolean& thePeriodic,
102	const TColStd_Array1OfReal& theFlatKnots,
103	const TColgp_Array1OfPnt2d& thePoles2d,
104	const TColStd_Array1OfReal& theWeights)
105	{
106	// Normalize theParameter for periodical B-splines
107	Standard_Real aNewParam = theParameter;
108	if (thePeriodic)
109	{
110	PeriodicNormalization(theFlatKnots, aNewParam);
111	myFlatKnots = new TColStd_HArray1OfReal(1, theFlatKnots.Length());
112	myFlatKnots->ChangeArray1() = theFlatKnots;
113	}
114	else if (!myFlatKnots.IsNull()) // Periodical curve became non-periodical
115	myFlatKnots.Nullify();
116
117	// Change the size of cached data if needed
118	myIsRational = (&theWeights != NULL);
119	Standard_Integer aPWColNumber = myIsRational ? 3 : 2;
120	if (theDegree > myDegree)
121	myPolesWeights = new TColStd_HArray2OfReal(1, theDegree + 1, 1, aPWColNumber);
122
123	myDegree = theDegree;
124	mySpanIndex = 0;
125	BSplCLib::LocateParameter(theDegree, theFlatKnots, BSplCLib::NoMults(),
126	aNewParam, thePeriodic, mySpanIndex, aNewParam);
127	mySpanStart = theFlatKnots.Value(mySpanIndex);
128	mySpanLength = theFlatKnots.Value(mySpanIndex + 1) - mySpanStart;
129	mySpanIndexMax = theFlatKnots.Length() - 1 - theDegree;
130
131	// Calculate new cache data
132	BSplCLib::BuildCache(mySpanStart, mySpanLength, thePeriodic, theDegree,
133	theFlatKnots, thePoles2d, theWeights,
134	myPolesWeights->ChangeArray2());
135	}
136
137	void BSplCLib_Cache::BuildCache(const Standard_Real& theParameter,
138	const Standard_Integer& theDegree,
139	const Standard_Boolean& thePeriodic,
140	const TColStd_Array1OfReal& theFlatKnots,
141	const TColgp_Array1OfPnt& thePoles,
142	const TColStd_Array1OfReal& theWeights)
143	{
144	// Create list of knots with repetitions and normalize theParameter for periodical B-splines
145	Standard_Real aNewParam = theParameter;
146	if (thePeriodic)
147	{
148	PeriodicNormalization(theFlatKnots, aNewParam);
149	myFlatKnots = new TColStd_HArray1OfReal(1, theFlatKnots.Length());
150	myFlatKnots->ChangeArray1() = theFlatKnots;
151	}
152	else if (!myFlatKnots.IsNull()) // Periodical curve became non-periodical
153	myFlatKnots.Nullify();
154
155	// Change the size of cached data if needed
156	myIsRational = (&theWeights != NULL);
157	Standard_Integer aPWColNumber = myIsRational ? 4 : 3;
158	if (theDegree > myDegree)
159	myPolesWeights = new TColStd_HArray2OfReal(1, theDegree + 1, 1, aPWColNumber);
160
161	myDegree = theDegree;
162	mySpanIndex = 0;
163	BSplCLib::LocateParameter(theDegree, theFlatKnots, BSplCLib::NoMults(),
164	aNewParam, thePeriodic, mySpanIndex, aNewParam);
165	mySpanStart = theFlatKnots.Value(mySpanIndex);
166	mySpanLength = theFlatKnots.Value(mySpanIndex + 1) - mySpanStart;
167	mySpanIndexMax = theFlatKnots.Length() - 1 - theDegree;
168
169	// Calculate new cache data
170	BSplCLib::BuildCache(mySpanStart, mySpanLength, thePeriodic, theDegree,
171	theFlatKnots, thePoles, theWeights,
172	myPolesWeights->ChangeArray2());
173	}
174
175
176	void BSplCLib_Cache::CalculateDerivative(const Standard_Real& theParameter,
177	const Standard_Integer& theDerivative,
178	Standard_Real& theDerivArray) const
179	{
180	Standard_Real aNewParameter = theParameter;
181	if (!myFlatKnots.IsNull()) // B-spline is periodical
182	PeriodicNormalization(myFlatKnots->Array1(), aNewParameter);
183	aNewParameter = (aNewParameter - mySpanStart) / mySpanLength;
184
185	Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
186	Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
187
188	// Temporary container. The maximal size of this container is defined by:
189	// 1) maximal derivative for cache evaluation, which is 3, plus one row for function values,
190	// 2) and maximal dimension of the point, which is 3, plus one column for weights.
191	Standard_Real aTmpContainer[16];
192
193	// When the PLib::RationaDerivative needs to be called, use temporary container
194	Standard_Real* aPntDeriv = myIsRational ? aTmpContainer : &theDerivArray;
195
196	// When the degree of curve is lesser than the requested derivative,
197	// nullify array cells corresponding to greater derivatives
198	Standard_Integer aDerivative = theDerivative;
199	if (myDegree < theDerivative)
200	{
201	aDerivative = myDegree;
202	for (Standard_Integer ind = myDegree * aDimension; ind < (theDerivative + 1) * aDimension; ind++)
203	{
204	aPntDeriv[ind] = 0.0;
205	(&theDerivArray)[ind] = 0.0; // should be cleared separately, because aPntDeriv may look to another memory area
206	}
207	}
208
209	PLib::EvalPolynomial(aNewParameter, aDerivative, myDegree, aDimension,
210	aPolesArray[0], aPntDeriv[0]);
211	// Unnormalize derivatives since those are computed normalized
212	Standard_Real aFactor = 1.0;
213	for (Standard_Integer deriv = 1; deriv <= aDerivative; deriv++)
214	{
215	aFactor /= mySpanLength;
216	for (Standard_Integer ind = 0; ind < aDimension; ind++)
217	aPntDeriv[aDimension * deriv + ind] *= aFactor;
218	}
219
220	if (myIsRational) // calculate derivatives divided by weights derivatives
221	PLib::RationalDerivative(aDerivative, aDerivative, aDimension - 1, aPntDeriv[0], theDerivArray);
222	}
223
224
225	void BSplCLib_Cache::D0(const Standard_Real& theParameter, gp_Pnt2d& thePoint) const
226	{
227	Standard_Real aNewParameter = theParameter;
228	if (!myFlatKnots.IsNull()) // B-spline is periodical
229	PeriodicNormalization(myFlatKnots->Array1(), aNewParameter);
230	aNewParameter = (aNewParameter - mySpanStart) / mySpanLength;
231
232	Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
233	Standard_Real aPoint[4];
234	Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
235
236	PLib::NoDerivativeEvalPolynomial(aNewParameter, myDegree,
237	aDimension, myDegree * aDimension,
238	aPolesArray[0], aPoint[0]);
239
240	thePoint.SetCoord(aPoint[0], aPoint[1]);
241	if (myIsRational)
242	thePoint.ChangeCoord().Divide(aPoint[2]);
243	}
244
245	void BSplCLib_Cache::D0(const Standard_Real& theParameter, gp_Pnt& thePoint) const
246	{
247	Standard_Real aNewParameter = theParameter;
248	if (!myFlatKnots.IsNull()) // B-spline is periodical
249	PeriodicNormalization(myFlatKnots->Array1(), aNewParameter);
250	aNewParameter = (aNewParameter - mySpanStart) / mySpanLength;
251
252	Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
253	Standard_Real aPoint[4];
254	Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
255
256	PLib::NoDerivativeEvalPolynomial(aNewParameter, myDegree,
257	aDimension, myDegree * aDimension,
258	aPolesArray[0], aPoint[0]);
259
260	thePoint.SetCoord(aPoint[0], aPoint[1], aPoint[2]);
261	if (myIsRational)
262	thePoint.ChangeCoord().Divide(aPoint[3]);
263	}
264
265
266	void BSplCLib_Cache::D1(const Standard_Real& theParameter, gp_Pnt2d& thePoint, gp_Vec2d& theTangent) const
267	{
268	Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
269	Standard_Real aPntDeriv[8]; // result storage (point and derivative coordinates)
270
271	this->CalculateDerivative(theParameter, 1, aPntDeriv[0]);
272	if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
273	aDimension -= 1;
274
275	thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1]);
276	theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1]);
277	}
278
279	void BSplCLib_Cache::D1(const Standard_Real& theParameter, gp_Pnt& thePoint, gp_Vec& theTangent) const
280	{
281	Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
282	Standard_Real aPntDeriv[8]; // result storage (point and derivative coordinates)
283
284	this->CalculateDerivative(theParameter, 1, aPntDeriv[0]);
285	if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
286	aDimension -= 1;
287
288	thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1], aPntDeriv[2]);
289	theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1], aPntDeriv[aDimension + 2]);
290	}
291
292	void BSplCLib_Cache::D2(const Standard_Real& theParameter, gp_Pnt2d& thePoint, gp_Vec2d& theTangent, gp_Vec2d& theCurvature) const
293	{
294	Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
295	Standard_Real aPntDeriv[12]; // result storage (point and derivatives coordinates)
296
297	this->CalculateDerivative(theParameter, 2, aPntDeriv[0]);
298	if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
299	aDimension -= 1;
300
301	thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1]);
302	theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1]);
303	theCurvature.SetCoord(aPntDeriv[aDimension<<1], aPntDeriv[(aDimension<<1) + 1]);
304	}
305
306	void BSplCLib_Cache::D2(const Standard_Real& theParameter, gp_Pnt& thePoint, gp_Vec& theTangent, gp_Vec& theCurvature) const
307	{
308	Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
309	Standard_Real aPntDeriv[12]; // result storage (point and derivatives coordinates)
310
311	this->CalculateDerivative(theParameter, 2, aPntDeriv[0]);
312	if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
313	aDimension -= 1;
314
315	thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1], aPntDeriv[2]);
316	theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1], aPntDeriv[aDimension + 2]);
317	theCurvature.SetCoord(aPntDeriv[aDimension<<1], aPntDeriv[(aDimension<<1) + 1], aPntDeriv[(aDimension<<1) + 2]);
318	}
319
320
321	void BSplCLib_Cache::D3(const Standard_Real& theParameter,
322	gp_Pnt2d& thePoint,
323	gp_Vec2d& theTangent,
324	gp_Vec2d& theCurvature,
325	gp_Vec2d& theTorsion) const
326	{
327	Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
328	Standard_Real aPntDeriv[16]; // result storage (point and derivatives coordinates)
329
330	this->CalculateDerivative(theParameter, 3, aPntDeriv[0]);
331	if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
332	aDimension -= 1;
333
334	thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1]);
335	theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1]);
336	Standard_Integer aShift = aDimension << 1;
337	theCurvature.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1]);
338	aShift += aDimension;
339	theTorsion.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1]);
340	}
341
342	void BSplCLib_Cache::D3(const Standard_Real& theParameter,
343	gp_Pnt& thePoint,
344	gp_Vec& theTangent,
345	gp_Vec& theCurvature,
346	gp_Vec& theTorsion) const
347	{
348	Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
349	Standard_Real aPntDeriv[16]; // result storage (point and derivatives coordinates)
350
351	this->CalculateDerivative(theParameter, 3, aPntDeriv[0]);
352	if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
353	aDimension -= 1;
354
355	thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1], aPntDeriv[2]);
356	theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1], aPntDeriv[aDimension + 2]);
357	Standard_Integer aShift = aDimension << 1;
358	theCurvature.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1], aPntDeriv[aShift + 2]);
359	aShift += aDimension;
360	theTorsion.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1], aPntDeriv[aShift + 2]);
361	}
362