src/BSplSLib/BSplSLib_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 <BSplSLib_Cache.hxx>
  15 #include <BSplSLib.hxx>
  16
  17 #include <NCollection_LocalArray.hxx>
  18
  19 #include <TColgp_HArray2OfPnt.hxx>
  20 #include <TColStd_HArray1OfInteger.hxx>
  21 #include <TColStd_HArray1OfReal.hxx>
  22 #include <TColStd_HArray2OfReal.hxx>
  23
  24 IMPLEMENT_STANDARD_HANDLE(BSplSLib_Cache, Standard_Transient)
  25 IMPLEMENT_STANDARD_RTTIEXT(BSplSLib_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 BSplSLib_Cache::BSplSLib_Cache()
  36 {
  37   myPolesWeights.Nullify();
  38   myIsRational = Standard_False;
  39   mySpanStart[0]  = mySpanStart[1]  = 0.0;
  40   mySpanLength[0] = mySpanLength[1] = 0.0;
  41   mySpanIndex[0]  = mySpanIndex[1]  = 0;
  42   myDegree[0]     = myDegree[1]     = 0;
  43   myFlatKnots[0].Nullify();
  44   myFlatKnots[1].Nullify();
  45 }
  46
  47 BSplSLib_Cache::BSplSLib_Cache(const Standard_Integer&        theDegreeU,
  48                                const Standard_Boolean&        thePeriodicU,
  49                                const TColStd_Array1OfReal&    theFlatKnotsU,
  50                                const Standard_Integer&        theDegreeV,
  51                                const Standard_Boolean&        thePeriodicV,
  52                                const TColStd_Array1OfReal&    theFlatKnotsV,
  53                                const TColgp_Array2OfPnt&      thePoles,
  54                                const TColStd_Array2OfReal&    theWeights)
  55 {
  56   Standard_Real aU = theFlatKnotsU.Value(theFlatKnotsU.Lower() + theDegreeU);
  57   Standard_Real aV = theFlatKnotsV.Value(theFlatKnotsV.Lower() + theDegreeV);
  58
  59   BuildCache(aU, aV,
  60              theDegreeU, thePeriodicU, theFlatKnotsU,
  61              theDegreeV, thePeriodicV, theFlatKnotsV,
  62              thePoles, theWeights);
  63 }
  64
  65
  66 Standard_Boolean BSplSLib_Cache::IsCacheValid(Standard_Real theParameterU,
  67                                               Standard_Real theParameterV) const
  68 {
  69   Standard_Real aNewU = theParameterU;
  70   Standard_Real aNewV = theParameterV;
  71   if (!myFlatKnots[0].IsNull())
  72     PeriodicNormalization(myDegree[0], myFlatKnots[0]->Array1(), aNewU);
  73   if (!myFlatKnots[1].IsNull())
  74     PeriodicNormalization(myDegree[1], myFlatKnots[1]->Array1(), aNewV);
  75
  76   Standard_Real aDelta0 = aNewU - mySpanStart[0];
  77   Standard_Real aDelta1 = aNewV - mySpanStart[1];
  78   return (aDelta0 >= -mySpanLength[0] && (aDelta0 < mySpanLength[0] || mySpanIndex[0] == mySpanIndexMax[0]) &&
  79           aDelta1 >= -mySpanLength[1] && (aDelta1 < mySpanLength[1] || mySpanIndex[1] == mySpanIndexMax[1]));
  80 }
  81
  82 void BSplSLib_Cache::PeriodicNormalization(const Standard_Integer& theDegree,
  83                                            const TColStd_Array1OfReal& theFlatKnots,
  84                                            Standard_Real& theParameter) const
  85 {
  86   Standard_Real aPeriod = theFlatKnots.Value(theFlatKnots.Upper() - theDegree) -
  87                           theFlatKnots.Value(theDegree + 1) ;
  88   if (theParameter < theFlatKnots.Value(theDegree + 1))
  89   {
  90     Standard_Real aScale = IntegerPart(
  91         (theFlatKnots.Value(theDegree + 1) - theParameter) / aPeriod);
  92     theParameter += aPeriod * (aScale + 1.0);
  93   }
  94   if (theParameter > theFlatKnots.Value(theFlatKnots.Upper() - theDegree))
  95   {
  96     Standard_Real aScale = IntegerPart(
  97         (theParameter - theFlatKnots.Value(theFlatKnots.Upper() - theDegree)) / aPeriod);
  98     theParameter -= aPeriod * (aScale + 1.0);
  99   }
 100 }
 101
 102
 103 void BSplSLib_Cache::BuildCache(const Standard_Real&           theParameterU,
 104                                 const Standard_Real&           theParameterV,
 105                                 const Standard_Integer&        theDegreeU,
 106                                 const Standard_Boolean&        thePeriodicU,
 107                                 const TColStd_Array1OfReal&    theFlatKnotsU,
 108                                 const Standard_Integer&        theDegreeV,
 109                                 const Standard_Boolean&        thePeriodicV,
 110                                 const TColStd_Array1OfReal&    theFlatKnotsV,
 111                                 const TColgp_Array2OfPnt&      thePoles,
 112                                 const TColStd_Array2OfReal&    theWeights)
 113 {
 114   // Normalize the parameters for periodical B-splines
 115   Standard_Real aNewParamU = theParameterU;
 116   if (thePeriodicU)
 117   {
 118     PeriodicNormalization(theDegreeU, theFlatKnotsU, aNewParamU);
 119     myFlatKnots[0] = new TColStd_HArray1OfReal(1, theFlatKnotsU.Length());
 120     myFlatKnots[0]->ChangeArray1() = theFlatKnotsU;
 121   }
 122   else if (!myFlatKnots[0].IsNull()) // Periodical curve became non-periodical
 123     myFlatKnots[0].Nullify();
 124
 125   Standard_Real aNewParamV = theParameterV;
 126   if (thePeriodicV)
 127   {
 128     PeriodicNormalization(theDegreeV, theFlatKnotsV, aNewParamV);
 129     myFlatKnots[1] = new TColStd_HArray1OfReal(1, theFlatKnotsV.Length());
 130     myFlatKnots[1]->ChangeArray1() = theFlatKnotsV;
 131   }
 132   else if (!myFlatKnots[1].IsNull()) // Periodical curve became non-periodical
 133     myFlatKnots[1].Nullify();
 134
 135   Standard_Integer aMinDegree = Min(theDegreeU, theDegreeV);
 136   Standard_Integer aMaxDegree = Max(theDegreeU, theDegreeV);
 137
 138   // Change the size of cached data if needed
 139   myIsRational = (&theWeights != NULL);
 140   Standard_Integer aPWColNumber = myIsRational ? 4 : 3;
 141   if (theDegreeU > myDegree[0] || theDegreeV > myDegree[1])
 142     myPolesWeights = new TColStd_HArray2OfReal(1, aMaxDegree + 1, 1, aPWColNumber * (aMinDegree + 1));
 143
 144   myDegree[0] = theDegreeU;
 145   myDegree[1] = theDegreeV;
 146   mySpanIndex[0] = mySpanIndex[1] = 0;
 147   BSplCLib::LocateParameter(theDegreeU, theFlatKnotsU, BSplCLib::NoMults(), aNewParamU,
 148                             thePeriodicU, mySpanIndex[0], aNewParamU);
 149   BSplCLib::LocateParameter(theDegreeV, theFlatKnotsV, BSplCLib::NoMults(), aNewParamV,
 150                             thePeriodicV, mySpanIndex[1], aNewParamV);
 151
 152   // Protection against Out of Range exception.
 153   if (mySpanIndex[0] >= theFlatKnotsU.Length()) {
 154     mySpanIndex[0] = theFlatKnotsU.Length() - 1;
 155   }
 156
 157   mySpanLength[0] = (theFlatKnotsU.Value(mySpanIndex[0] + 1) - theFlatKnotsU.Value(mySpanIndex[0])) * 0.5;
 158   mySpanStart[0]  = theFlatKnotsU.Value(mySpanIndex[0]) + mySpanLength[0];
 159
 160   // Protection against Out of Range exception.
 161   if (mySpanIndex[1] >= theFlatKnotsV.Length()) {
 162     mySpanIndex[1] = theFlatKnotsV.Length() - 1;
 163   }
 164
 165   mySpanLength[1] = (theFlatKnotsV.Value(mySpanIndex[1] + 1) - theFlatKnotsV.Value(mySpanIndex[1])) * 0.5;
 166   mySpanStart[1]  = theFlatKnotsV.Value(mySpanIndex[1]) + mySpanLength[1];
 167   mySpanIndexMax[0] = theFlatKnotsU.Length() - 1 - theDegreeU;
 168   mySpanIndexMax[1] = theFlatKnotsV.Length() - 1 - theDegreeV;
 169
 170   // Calculate new cache data
 171   BSplSLib::BuildCache(mySpanStart[0],  mySpanStart[1],
 172                        mySpanLength[0], mySpanLength[1],
 173                        thePeriodicU,    thePeriodicV,
 174                        theDegreeU,      theDegreeV,
 175                        mySpanIndex[0],  mySpanIndex[1],
 176                        theFlatKnotsU,   theFlatKnotsV,
 177                        thePoles, theWeights, myPolesWeights->ChangeArray2());
 178 }
 179
 180
 181 void BSplSLib_Cache::D0(const Standard_Real& theU,
 182                         const Standard_Real& theV,
 183                               gp_Pnt&        thePoint) const
 184 {
 185   Standard_Real aNewU = theU;
 186   Standard_Real aNewV = theV;
 187   if (!myFlatKnots[0].IsNull()) // B-spline is U-periodical
 188     PeriodicNormalization(myDegree[0], myFlatKnots[0]->Array1(), aNewU);
 189   aNewU = (aNewU - mySpanStart[0]) / mySpanLength[0];
 190   if (!myFlatKnots[1].IsNull()) // B-spline is V-periodical
 191     PeriodicNormalization(myDegree[1], myFlatKnots[1]->Array1(), aNewV);
 192   aNewV = (aNewV - mySpanStart[1]) / mySpanLength[1];
 193
 194   Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
 195   Standard_Real aPoint[4];
 196
 197   Standard_Integer aDimension = myIsRational ? 4 : 3;
 198   Standard_Integer aCacheCols = myPolesWeights->RowLength();
 199   Standard_Integer aMinMaxDegree[2] = {Min(myDegree[0], myDegree[1]),
 200                                        Max(myDegree[0], myDegree[1])};
 201   Standard_Real aParameters[2];
 202   if (myDegree[0] > myDegree[1])
 203   {
 204     aParameters[0] = aNewV;
 205     aParameters[1] = aNewU;
 206   }
 207   else
 208   {
 209     aParameters[0] = aNewU;
 210     aParameters[1] = aNewV;
 211   }
 212
 213   NCollection_LocalArray<Standard_Real> aTransientCoeffs(aCacheCols); // array for intermediate results
 214
 215   // Calculate intermediate value of cached polynomial along columns
 216   PLib::NoDerivativeEvalPolynomial(aParameters[1], aMinMaxDegree[1],
 217                                    aCacheCols, aMinMaxDegree[1] * aCacheCols,
 218                                    aPolesArray[0], aTransientCoeffs[0]);
 219
 220   // Calculate total value
 221   PLib::NoDerivativeEvalPolynomial(aParameters[0], aMinMaxDegree[0],
 222                                    aDimension, aDimension * aMinMaxDegree[0],
 223                                    aTransientCoeffs[0], aPoint[0]);
 224
 225   thePoint.SetCoord(aPoint[0], aPoint[1], aPoint[2]);
 226   if (myIsRational)
 227     thePoint.ChangeCoord().Divide(aPoint[3]);
 228 }
 229
 230
 231 void BSplSLib_Cache::D1(const Standard_Real& theU,
 232                         const Standard_Real& theV,
 233                               gp_Pnt&        thePoint,
 234                               gp_Vec&        theTangentU,
 235                               gp_Vec&        theTangentV) const
 236 {
 237   Standard_Real aNewU = theU;
 238   Standard_Real aNewV = theV;
 239   Standard_Real anInvU = 1.0 / mySpanLength[0];
 240   Standard_Real anInvV = 1.0 / mySpanLength[1];
 241   if (!myFlatKnots[0].IsNull()) // B-spline is U-periodical
 242     PeriodicNormalization(myDegree[0], myFlatKnots[0]->Array1(), aNewU);
 243   aNewU = (aNewU - mySpanStart[0]) * anInvU;
 244   if (!myFlatKnots[1].IsNull()) // B-spline is V-periodical
 245     PeriodicNormalization(myDegree[1], myFlatKnots[1]->Array1(), aNewV);
 246   aNewV = (aNewV - mySpanStart[1]) * anInvV;
 247
 248   Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
 249   Standard_Real aPntDeriv[16]; // result storage (point and derivative coordinates)
 250   for (Standard_Integer i = 0; i< 16; i++) aPntDeriv[i] = 0.0;
 251
 252   Standard_Integer aDimension = myIsRational ? 4 : 3;
 253   Standard_Integer aCacheCols = myPolesWeights->RowLength();
 254   Standard_Integer aMinMaxDegree[2] = {Min(myDegree[0], myDegree[1]),
 255                                        Max(myDegree[0], myDegree[1])};
 256
 257   Standard_Real aParameters[2];
 258   if (myDegree[0] > myDegree[1])
 259   {
 260     aParameters[0] = aNewV;
 261     aParameters[1] = aNewU;
 262   }
 263   else
 264   {
 265     aParameters[0] = aNewU;
 266     aParameters[1] = aNewV;
 267   }
 268
 269   NCollection_LocalArray<Standard_Real> aTransientCoeffs(aCacheCols<<1); // array for intermediate results
 270
 271   // Calculate intermediate values and derivatives of bivariate polynomial along variable with maximal degree
 272   PLib::EvalPolynomial(aParameters[1], 1, aMinMaxDegree[1], aCacheCols, aPolesArray[0], aTransientCoeffs[0]);
 273
 274   // Calculate a point on surface and a derivative along variable with minimal degree
 275   PLib::EvalPolynomial(aParameters[0], 1, aMinMaxDegree[0], aDimension, aTransientCoeffs[0], aPntDeriv[0]);
 276
 277   // Calculate derivative along variable with maximal degree
 278   PLib::NoDerivativeEvalPolynomial(aParameters[0], aMinMaxDegree[0], aDimension,
 279                                    aMinMaxDegree[0] * aDimension, aTransientCoeffs[aCacheCols],
 280                                    aPntDeriv[aDimension<<1]);
 281
 282   Standard_Real* aResult = aPntDeriv;
 283   Standard_Real aTempStorage[12];
 284   if (myIsRational) // calculate derivatives divided by weight's derivatives
 285   {
 286     BSplSLib::RationalDerivative(1, 1, 1, 1, aPntDeriv[0], aTempStorage[0]);
 287     aResult = aTempStorage;
 288     aDimension--;
 289   }
 290
 291   thePoint.SetCoord(aResult[0], aResult[1], aResult[2]);
 292   if (myDegree[0] > myDegree[1])
 293   {
 294     theTangentV.SetCoord(aResult[aDimension], aResult[aDimension + 1], aResult[aDimension + 2]);
 295     Standard_Integer aShift = aDimension<<1;
 296     theTangentU.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 297   }
 298   else
 299   {
 300     theTangentU.SetCoord(aResult[aDimension], aResult[aDimension + 1], aResult[aDimension + 2]);
 301     Standard_Integer aShift = aDimension<<1;
 302     theTangentV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 303   }
 304   theTangentU.Multiply(anInvU);
 305   theTangentV.Multiply(anInvV);
 306 }
 307
 308
 309 void BSplSLib_Cache::D2(const Standard_Real& theU,
 310                         const Standard_Real& theV,
 311                               gp_Pnt&        thePoint,
 312                               gp_Vec&        theTangentU,
 313                               gp_Vec&        theTangentV,
 314                               gp_Vec&        theCurvatureU,
 315                               gp_Vec&        theCurvatureV,
 316                               gp_Vec&        theCurvatureUV) const
 317 {
 318   Standard_Real aNewU = theU;
 319   Standard_Real aNewV = theV;
 320   Standard_Real anInvU = 1.0 / mySpanLength[0];
 321   Standard_Real anInvV = 1.0 / mySpanLength[1];
 322   if (!myFlatKnots[0].IsNull()) // B-spline is U-periodical
 323     PeriodicNormalization(myDegree[0], myFlatKnots[0]->Array1(), aNewU);
 324   aNewU = (aNewU - mySpanStart[0]) * anInvU;
 325   if (!myFlatKnots[1].IsNull()) // B-spline is V-periodical
 326     PeriodicNormalization(myDegree[1], myFlatKnots[1]->Array1(), aNewV);
 327   aNewV = (aNewV - mySpanStart[1]) * anInvV;
 328
 329   Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
 330   Standard_Real aPntDeriv[36]; // result storage (point and derivative coordinates)
 331   for (Standard_Integer i = 0; i < 36; i++) aPntDeriv[i] = 0.0;
 332
 333   Standard_Integer aDimension = myIsRational ? 4 : 3;
 334   Standard_Integer aCacheCols = myPolesWeights->RowLength();
 335   Standard_Integer aMinMaxDegree[2] = {Min(myDegree[0], myDegree[1]),
 336                                        Max(myDegree[0], myDegree[1])};
 337
 338   Standard_Real aParameters[2];
 339   if (myDegree[0] > myDegree[1])
 340   {
 341     aParameters[0] = aNewV;
 342     aParameters[1] = aNewU;
 343   }
 344   else
 345   {
 346     aParameters[0] = aNewU;
 347     aParameters[1] = aNewV;
 348   }
 349
 350   NCollection_LocalArray<Standard_Real> aTransientCoeffs(3 * aCacheCols); // array for intermediate results
 351   // Calculating derivative to be evaluate and
 352   // nulling transient coefficients when max or min derivative is less than 2
 353   Standard_Integer aMinMaxDeriv[2] = {Min(2, aMinMaxDegree[0]),
 354                                       Min(2, aMinMaxDegree[1])};
 355   for (Standard_Integer i = aMinMaxDeriv[1] + 1; i < 3; i++)
 356   {
 357     Standard_Integer index = i * aCacheCols;
 358     for (Standard_Integer j = 0; j < aCacheCols; j++)
 359       aTransientCoeffs[index++] = 0.0;
 360   }
 361
 362   // Calculate intermediate values and derivatives of bivariate polynomial along variable with maximal degree
 363   PLib::EvalPolynomial(aParameters[1], aMinMaxDeriv[1], aMinMaxDegree[1],
 364                        aCacheCols, aPolesArray[0], aTransientCoeffs[0]);
 365
 366   // Calculate a point on surface and a derivatives along variable with minimal degree
 367   PLib::EvalPolynomial(aParameters[0], aMinMaxDeriv[0], aMinMaxDegree[0],
 368                        aDimension, aTransientCoeffs[0], aPntDeriv[0]);
 369
 370   // Calculate derivative along variable with maximal degree and mixed derivative
 371   PLib::EvalPolynomial(aParameters[0], 1, aMinMaxDegree[0], aDimension,
 372                        aTransientCoeffs[aCacheCols], aPntDeriv[3 * aDimension]);
 373
 374   // Calculate second derivative along variable with maximal degree
 375   PLib::NoDerivativeEvalPolynomial(aParameters[0], aMinMaxDegree[0], aDimension,
 376                                    aMinMaxDegree[0] * aDimension, aTransientCoeffs[aCacheCols<<1],
 377                                    aPntDeriv[6 * aDimension]);
 378
 379   Standard_Real* aResult = aPntDeriv;
 380   Standard_Real aTempStorage[36];
 381   if (myIsRational) // calculate derivatives divided by weight's derivatives
 382   {
 383     BSplSLib::RationalDerivative(2, 2, 2, 2, aPntDeriv[0], aTempStorage[0]);
 384     aResult = aTempStorage;
 385     aDimension--;
 386   }
 387
 388   thePoint.SetCoord(aResult[0], aResult[1], aResult[2]);
 389   if (myDegree[0] > myDegree[1])
 390   {
 391     theTangentV.SetCoord(aResult[aDimension], aResult[aDimension + 1], aResult[aDimension + 2]);
 392     Standard_Integer aShift = aDimension<<1;
 393     theCurvatureV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 394     aShift += aDimension;
 395     theTangentU.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 396     aShift += aDimension;
 397     theCurvatureUV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 398     aShift += (aDimension << 1);
 399     theCurvatureU.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 400   }
 401   else
 402   {
 403     theTangentU.SetCoord(aResult[aDimension], aResult[aDimension + 1], aResult[aDimension + 2]);
 404     Standard_Integer aShift = aDimension<<1;
 405     theCurvatureU.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 406     aShift += aDimension;
 407     theTangentV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 408     aShift += aDimension;
 409     theCurvatureUV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 410     aShift += (aDimension << 1);
 411     theCurvatureV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 412   }
 413   theTangentU.Multiply(anInvU);
 414   theTangentV.Multiply(anInvV);
 415   theCurvatureU.Multiply(anInvU * anInvU);
 416   theCurvatureV.Multiply(anInvV * anInvV);
 417   theCurvatureUV.Multiply(anInvU * anInvV);
 418 }
 419