src/BSplCLib/BSplCLib_Cache.cxx

   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_HArray2OfReal.hxx>
  22
  23
  24 IMPLEMENT_STANDARD_RTTIEXT(BSplCLib_Cache,Standard_Transient)
  25
  26 //! Converts handle of array of Standard_Real into the pointer to Standard_Real
  27 static Standard_Real* ConvertArray(const Handle(TColStd_HArray2OfReal)& theHArray)
  28 {
  29   const TColStd_Array2OfReal& anArray = theHArray->Array2();
  30   return (Standard_Real*) &(anArray(anArray.LowerRow(), anArray.LowerCol()));
  31 }
  32
  33 BSplCLib_Cache::BSplCLib_Cache(const Standard_Integer&        theDegree,
  34                                const Standard_Boolean&        thePeriodic,
  35                                const TColStd_Array1OfReal&    theFlatKnots,
  36                                const TColgp_Array1OfPnt2d&    /* only used to distinguish from 3d variant */,
  37                                const TColStd_Array1OfReal*    theWeights)
  38 : myIsRational(theWeights != NULL),
  39   myParams (theDegree, thePeriodic, theFlatKnots)
  40 {
  41   Standard_Integer aPWColNumber = (myIsRational ? 3 : 2);
  42   myPolesWeights = new TColStd_HArray2OfReal (1, theDegree + 1, 1, aPWColNumber);
  43 }
  44
  45 BSplCLib_Cache::BSplCLib_Cache(const Standard_Integer&        theDegree,
  46                                const Standard_Boolean&        thePeriodic,
  47                                const TColStd_Array1OfReal&    theFlatKnots,
  48                                const TColgp_Array1OfPnt&      /* only used to distinguish from 2d variant */,
  49                                const TColStd_Array1OfReal*    theWeights)
  50 : myIsRational(theWeights != NULL),
  51   myParams (theDegree, thePeriodic, theFlatKnots)
  52 {
  53   Standard_Integer aPWColNumber = (myIsRational ? 4 : 3);
  54   myPolesWeights = new TColStd_HArray2OfReal (1, theDegree + 1, 1, aPWColNumber);
  55 }
  56
  57 Standard_Boolean BSplCLib_Cache::IsCacheValid(Standard_Real theParameter) const
  58 {
  59   return myParams.IsCacheValid (theParameter);
  60 }
  61
  62 void BSplCLib_Cache::BuildCache(const Standard_Real&           theParameter,
  63                                 const TColStd_Array1OfReal&    theFlatKnots,
  64                                 const TColgp_Array1OfPnt2d&    thePoles2d,
  65                                 const TColStd_Array1OfReal*    theWeights)
  66 {
  67   // Normalize theParameter for periodical B-splines
  68   Standard_Real aNewParam = myParams.PeriodicNormalization (theParameter);
  69   myParams.LocateParameter (aNewParam, theFlatKnots);
  70
  71   // Calculate new cache data
  72   BSplCLib::BuildCache (myParams.SpanStart, myParams.SpanLength, myParams.IsPeriodic,
  73                         myParams.Degree, myParams.SpanIndex, theFlatKnots, thePoles2d,
  74                         theWeights, myPolesWeights->ChangeArray2());
  75 }
  76
  77 void BSplCLib_Cache::BuildCache(const Standard_Real&           theParameter,
  78                                 const TColStd_Array1OfReal&    theFlatKnots,
  79                                 const TColgp_Array1OfPnt&      thePoles,
  80                                 const TColStd_Array1OfReal*    theWeights)
  81 {
  82   // Create list of knots with repetitions and normalize theParameter for periodical B-splines
  83   Standard_Real aNewParam = myParams.PeriodicNormalization (theParameter);
  84   myParams.LocateParameter (aNewParam, theFlatKnots);
  85
  86   // Calculate new cache data
  87   BSplCLib::BuildCache (myParams.SpanStart, myParams.SpanLength, myParams.IsPeriodic,
  88                         myParams.Degree, myParams.SpanIndex, theFlatKnots, thePoles,
  89                         theWeights, myPolesWeights->ChangeArray2());
  90 }
  91
  92 void BSplCLib_Cache::CalculateDerivative(const Standard_Real&    theParameter,
  93                                          const Standard_Integer& theDerivative,
  94                                                Standard_Real&    theDerivArray) const
  95 {
  96   Standard_Real aNewParameter = myParams.PeriodicNormalization (theParameter);
  97   aNewParameter = (aNewParameter - myParams.SpanStart) / myParams.SpanLength;
  98
  99   Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
 100   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
 101
 102   // Temporary container. The maximal size of this container is defined by:
 103   //    1) maximal derivative for cache evaluation, which is 3, plus one row for function values,
 104   //    2) and maximal dimension of the point, which is 3, plus one column for weights.
 105   Standard_Real aTmpContainer[16];
 106
 107   // When the PLib::RationaDerivative needs to be called, use temporary container
 108   Standard_Real* aPntDeriv = myIsRational ? aTmpContainer : &theDerivArray;
 109
 110   // When the degree of curve is lesser than the requested derivative,
 111   // nullify array cells corresponding to greater derivatives
 112   Standard_Integer aDerivative = theDerivative;
 113   if (!myIsRational && myParams.Degree < theDerivative)
 114   {
 115     aDerivative = myParams.Degree;
 116     for (Standard_Integer ind = myParams.Degree * aDimension; ind < (theDerivative + 1) * aDimension; ind++)
 117     {
 118       aPntDeriv[ind] = 0.0;
 119       (&theDerivArray)[ind] = 0.0; // should be cleared separately, because aPntDeriv may look to another memory area
 120     }
 121   }
 122
 123   PLib::EvalPolynomial(aNewParameter, aDerivative, myParams.Degree, aDimension,
 124                        aPolesArray[0], aPntDeriv[0]);
 125   // Unnormalize derivatives since those are computed normalized
 126   Standard_Real aFactor = 1.0;
 127   for (Standard_Integer deriv = 1; deriv <= aDerivative; deriv++)
 128   {
 129     aFactor /= myParams.SpanLength;
 130     for (Standard_Integer ind = 0; ind < aDimension; ind++)
 131       aPntDeriv[aDimension * deriv + ind] *= aFactor;
 132   }
 133
 134   if (myIsRational) // calculate derivatives divided by weights derivatives
 135     PLib::RationalDerivative(aDerivative, aDerivative, aDimension - 1, aPntDeriv[0], theDerivArray);
 136 }
 137
 138
 139 void BSplCLib_Cache::D0(const Standard_Real& theParameter, gp_Pnt2d& thePoint) const
 140 {
 141   Standard_Real aNewParameter = myParams.PeriodicNormalization (theParameter);
 142   aNewParameter = (aNewParameter - myParams.SpanStart) / myParams.SpanLength;
 143
 144   Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
 145   Standard_Real aPoint[4];
 146   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
 147
 148   PLib::NoDerivativeEvalPolynomial(aNewParameter, myParams.Degree,
 149                                    aDimension, myParams.Degree * aDimension,
 150                                    aPolesArray[0], aPoint[0]);
 151
 152   thePoint.SetCoord(aPoint[0], aPoint[1]);
 153   if (myIsRational)
 154     thePoint.ChangeCoord().Divide(aPoint[2]);
 155 }
 156
 157 void BSplCLib_Cache::D0(const Standard_Real& theParameter, gp_Pnt& thePoint) const
 158 {
 159   Standard_Real aNewParameter = myParams.PeriodicNormalization (theParameter);
 160   aNewParameter = (aNewParameter - myParams.SpanStart) / myParams.SpanLength;
 161
 162   Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
 163   Standard_Real aPoint[4];
 164   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
 165
 166   PLib::NoDerivativeEvalPolynomial(aNewParameter, myParams.Degree,
 167                                    aDimension, myParams.Degree * aDimension,
 168                                    aPolesArray[0], aPoint[0]);
 169
 170   thePoint.SetCoord(aPoint[0], aPoint[1], aPoint[2]);
 171   if (myIsRational)
 172     thePoint.ChangeCoord().Divide(aPoint[3]);
 173 }
 174
 175
 176 void BSplCLib_Cache::D1(const Standard_Real& theParameter, gp_Pnt2d& thePoint, gp_Vec2d& theTangent) const
 177 {
 178   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
 179   Standard_Real aPntDeriv[8]; // result storage (point and derivative coordinates)
 180
 181   this->CalculateDerivative(theParameter, 1, aPntDeriv[0]);
 182   if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
 183     aDimension -= 1;
 184
 185   thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1]);
 186   theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1]);
 187 }
 188
 189 void BSplCLib_Cache::D1(const Standard_Real& theParameter, gp_Pnt& thePoint, gp_Vec& theTangent) const
 190 {
 191   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
 192   Standard_Real aPntDeriv[8]; // result storage (point and derivative coordinates)
 193
 194   this->CalculateDerivative(theParameter, 1, aPntDeriv[0]);
 195   if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
 196     aDimension -= 1;
 197
 198   thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1], aPntDeriv[2]);
 199   theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1], aPntDeriv[aDimension + 2]);
 200 }
 201
 202 void BSplCLib_Cache::D2(const Standard_Real& theParameter, gp_Pnt2d& thePoint, gp_Vec2d& theTangent, gp_Vec2d& theCurvature) const
 203 {
 204   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
 205   Standard_Real aPntDeriv[12]; // result storage (point and derivatives coordinates)
 206
 207   this->CalculateDerivative(theParameter, 2, aPntDeriv[0]);
 208   if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
 209     aDimension -= 1;
 210
 211   thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1]);
 212   theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1]);
 213   theCurvature.SetCoord(aPntDeriv[aDimension<<1], aPntDeriv[(aDimension<<1) + 1]);
 214 }
 215
 216 void BSplCLib_Cache::D2(const Standard_Real& theParameter, gp_Pnt& thePoint, gp_Vec& theTangent, gp_Vec& theCurvature) const
 217 {
 218   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
 219   Standard_Real aPntDeriv[12]; // result storage (point and derivatives coordinates)
 220
 221   this->CalculateDerivative(theParameter, 2, aPntDeriv[0]);
 222   if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
 223     aDimension -= 1;
 224
 225   thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1], aPntDeriv[2]);
 226   theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1], aPntDeriv[aDimension + 2]);
 227   theCurvature.SetCoord(aPntDeriv[aDimension<<1], aPntDeriv[(aDimension<<1) + 1], aPntDeriv[(aDimension<<1) + 2]);
 228 }
 229
 230
 231 void BSplCLib_Cache::D3(const Standard_Real& theParameter,
 232                               gp_Pnt2d&      thePoint,
 233                               gp_Vec2d&      theTangent,
 234                               gp_Vec2d&      theCurvature,
 235                               gp_Vec2d&      theTorsion) const
 236 {
 237   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
 238   Standard_Real aPntDeriv[16]; // result storage (point and derivatives coordinates)
 239
 240   this->CalculateDerivative(theParameter, 3, aPntDeriv[0]);
 241   if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
 242     aDimension -= 1;
 243
 244   thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1]);
 245   theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1]);
 246   Standard_Integer aShift = aDimension << 1;
 247   theCurvature.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1]);
 248   aShift += aDimension;
 249   theTorsion.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1]);
 250 }
 251
 252 void BSplCLib_Cache::D3(const Standard_Real& theParameter,
 253                               gp_Pnt&        thePoint,
 254                               gp_Vec&        theTangent,
 255                               gp_Vec&        theCurvature,
 256                               gp_Vec&        theTorsion) const
 257 {
 258   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
 259   Standard_Real aPntDeriv[16]; // result storage (point and derivatives coordinates)
 260
 261   this->CalculateDerivative(theParameter, 3, aPntDeriv[0]);
 262   if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
 263     aDimension -= 1;
 264
 265   thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1], aPntDeriv[2]);
 266   theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1], aPntDeriv[aDimension + 2]);
 267   Standard_Integer aShift = aDimension << 1;
 268   theCurvature.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1], aPntDeriv[aShift + 2]);
 269   aShift += aDimension;
 270   theTorsion.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1], aPntDeriv[aShift + 2]);
 271 }
 272