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