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