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 BSplSLib_Cache::BSplSLib_Cache(const Standard_Integer&        theDegreeU,
  35                                const Standard_Boolean&        thePeriodicU,
  36                                const TColStd_Array1OfReal&    theFlatKnotsU,
  37                                const Standard_Integer&        theDegreeV,
  38                                const Standard_Boolean&        thePeriodicV,
  39                                const TColStd_Array1OfReal&    theFlatKnotsV,
  40                                const TColStd_Array2OfReal*    theWeights)
  41 : myIsRational(theWeights != NULL),
  42   myParamsU (theDegreeU, thePeriodicU, theFlatKnotsU),
  43   myParamsV (theDegreeV, thePeriodicV, theFlatKnotsV)
  44 {
  45   Standard_Integer aMinDegree = Min (theDegreeU, theDegreeV);
  46   Standard_Integer aMaxDegree = Max (theDegreeU, theDegreeV);
  47   Standard_Integer aPWColNumber = (myIsRational ? 4 : 3);
  48   myPolesWeights = new TColStd_HArray2OfReal(1, aMaxDegree + 1, 1, aPWColNumber * (aMinDegree + 1));
  49 }
  50
  51 Standard_Boolean BSplSLib_Cache::IsCacheValid(Standard_Real theParameterU,
  52                                               Standard_Real theParameterV) const
  53 {
  54   return myParamsU.IsCacheValid (theParameterU) &&
  55          myParamsV.IsCacheValid (theParameterV);
  56 }
  57
  58 void BSplSLib_Cache::BuildCache(const Standard_Real&           theParameterU,
  59                                 const Standard_Real&           theParameterV,
  60                                 const TColStd_Array1OfReal&    theFlatKnotsU,
  61                                 const TColStd_Array1OfReal&    theFlatKnotsV,
  62                                 const TColgp_Array2OfPnt&      thePoles,
  63                                 const TColStd_Array2OfReal*    theWeights)
  64 {
  65   // Normalize the parameters for periodical B-splines
  66   Standard_Real aNewParamU = myParamsU.PeriodicNormalization (theParameterU);
  67   Standard_Real aNewParamV = myParamsV.PeriodicNormalization (theParameterV);
  68
  69   myParamsU.LocateParameter (aNewParamU, theFlatKnotsU);
  70   myParamsV.LocateParameter (aNewParamV, theFlatKnotsV);
  71
  72   // BSplSLib uses different convention for span parameters than BSplCLib
  73   // (Start is in the middle of the span and length is half-span),
  74   // thus we need to amend them here
  75   Standard_Real aSpanLengthU = 0.5 * myParamsU.SpanLength;
  76   Standard_Real aSpanStartU = myParamsU.SpanStart + aSpanLengthU;
  77   Standard_Real aSpanLengthV = 0.5 * myParamsV.SpanLength;
  78   Standard_Real aSpanStartV = myParamsV.SpanStart + aSpanLengthV;
  79
  80   // Calculate new cache data
  81   BSplSLib::BuildCache (aSpanStartU,  aSpanStartV,
  82                         aSpanLengthU, aSpanLengthV,
  83                         myParamsU.IsPeriodic, myParamsV.IsPeriodic,
  84                         myParamsU.Degree,     myParamsV.Degree,
  85                         myParamsU.SpanIndex,  myParamsV.SpanIndex,
  86                         theFlatKnotsU,        theFlatKnotsV,
  87                         thePoles, theWeights, myPolesWeights->ChangeArray2());
  88 }
  89
  90
  91 void BSplSLib_Cache::D0(const Standard_Real& theU,
  92                         const Standard_Real& theV,
  93                               gp_Pnt&        thePoint) const
  94 {
  95   Standard_Real aNewU = myParamsU.PeriodicNormalization (theU);
  96   Standard_Real aNewV = myParamsV.PeriodicNormalization (theV);
  97
  98   // BSplSLib uses different convention for span parameters than BSplCLib
  99   // (Start is in the middle of the span and length is half-span),
 100   // thus we need to amend them here
 101   Standard_Real aSpanLengthU = 0.5 * myParamsU.SpanLength;
 102   Standard_Real aSpanStartU = myParamsU.SpanStart + aSpanLengthU;
 103   Standard_Real aSpanLengthV = 0.5 * myParamsV.SpanLength;
 104   Standard_Real aSpanStartV = myParamsV.SpanStart + aSpanLengthV;
 105
 106   aNewU = (aNewU - aSpanStartU) / aSpanLengthU;
 107   aNewV = (aNewV - aSpanStartV) / aSpanLengthV;
 108
 109   Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
 110   Standard_Real aPoint[4];
 111
 112   Standard_Integer aDimension = myIsRational ? 4 : 3;
 113   Standard_Integer aCacheCols = myPolesWeights->RowLength();
 114   Standard_Integer aMinMaxDegree[2] = {Min(myParamsU.Degree, myParamsV.Degree),
 115                                        Max(myParamsU.Degree, myParamsV.Degree)};
 116   Standard_Real aParameters[2];
 117   if (myParamsU.Degree > myParamsV.Degree)
 118   {
 119     aParameters[0] = aNewV;
 120     aParameters[1] = aNewU;
 121   }
 122   else
 123   {
 124     aParameters[0] = aNewU;
 125     aParameters[1] = aNewV;
 126   }
 127
 128   NCollection_LocalArray<Standard_Real> aTransientCoeffs(aCacheCols); // array for intermediate results
 129
 130   // Calculate intermediate value of cached polynomial along columns
 131   PLib::NoDerivativeEvalPolynomial(aParameters[1], aMinMaxDegree[1],
 132                                    aCacheCols, aMinMaxDegree[1] * aCacheCols,
 133                                    aPolesArray[0], aTransientCoeffs[0]);
 134
 135   // Calculate total value
 136   PLib::NoDerivativeEvalPolynomial(aParameters[0], aMinMaxDegree[0],
 137                                    aDimension, aDimension * aMinMaxDegree[0],
 138                                    aTransientCoeffs[0], aPoint[0]);
 139
 140   thePoint.SetCoord(aPoint[0], aPoint[1], aPoint[2]);
 141   if (myIsRational)
 142     thePoint.ChangeCoord().Divide(aPoint[3]);
 143 }
 144
 145
 146 void BSplSLib_Cache::D1(const Standard_Real& theU,
 147                         const Standard_Real& theV,
 148                               gp_Pnt&        thePoint,
 149                               gp_Vec&        theTangentU,
 150                               gp_Vec&        theTangentV) const
 151 {
 152   Standard_Real aNewU = myParamsU.PeriodicNormalization (theU);
 153   Standard_Real aNewV = myParamsV.PeriodicNormalization (theV);
 154
 155   // BSplSLib uses different convention for span parameters than BSplCLib
 156   // (Start is in the middle of the span and length is half-span),
 157   // thus we need to amend them here
 158   Standard_Real aSpanLengthU = 0.5 * myParamsU.SpanLength;
 159   Standard_Real aSpanStartU = myParamsU.SpanStart + aSpanLengthU;
 160   Standard_Real aSpanLengthV = 0.5 * myParamsV.SpanLength;
 161   Standard_Real aSpanStartV = myParamsV.SpanStart + aSpanLengthV;
 162
 163   Standard_Real anInvU = 1.0 / aSpanLengthU;
 164   Standard_Real anInvV = 1.0 / aSpanLengthV;
 165   aNewU = (aNewU - aSpanStartU) * anInvU;
 166   aNewV = (aNewV - aSpanStartV) * anInvV;
 167
 168   Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
 169   Standard_Real aPntDeriv[16]; // result storage (point and derivative coordinates)
 170   for (Standard_Integer i = 0; i< 16; i++) aPntDeriv[i] = 0.0;
 171
 172   Standard_Integer aDimension = myIsRational ? 4 : 3;
 173   Standard_Integer aCacheCols = myPolesWeights->RowLength();
 174   Standard_Integer aMinMaxDegree[2] = {Min(myParamsU.Degree, myParamsV.Degree),
 175                                        Max(myParamsU.Degree, myParamsV.Degree)};
 176
 177   Standard_Real aParameters[2];
 178   if (myParamsU.Degree > myParamsV.Degree)
 179   {
 180     aParameters[0] = aNewV;
 181     aParameters[1] = aNewU;
 182   }
 183   else
 184   {
 185     aParameters[0] = aNewU;
 186     aParameters[1] = aNewV;
 187   }
 188
 189   NCollection_LocalArray<Standard_Real> aTransientCoeffs(aCacheCols<<1); // array for intermediate results
 190
 191   // Calculate intermediate values and derivatives of bivariate polynomial along variable with maximal degree
 192   PLib::EvalPolynomial(aParameters[1], 1, aMinMaxDegree[1], aCacheCols, aPolesArray[0], aTransientCoeffs[0]);
 193
 194   // Calculate a point on surface and a derivative along variable with minimal degree
 195   PLib::EvalPolynomial(aParameters[0], 1, aMinMaxDegree[0], aDimension, aTransientCoeffs[0], aPntDeriv[0]);
 196
 197   // Calculate derivative along variable with maximal degree
 198   PLib::NoDerivativeEvalPolynomial(aParameters[0], aMinMaxDegree[0], aDimension,
 199                                    aMinMaxDegree[0] * aDimension, aTransientCoeffs[aCacheCols],
 200                                    aPntDeriv[aDimension<<1]);
 201
 202   Standard_Real* aResult = aPntDeriv;
 203   Standard_Real aTempStorage[12];
 204   if (myIsRational) // calculate derivatives divided by weight's derivatives
 205   {
 206     BSplSLib::RationalDerivative(1, 1, 1, 1, aPntDeriv[0], aTempStorage[0]);
 207     aResult = aTempStorage;
 208     aDimension--;
 209   }
 210
 211   thePoint.SetCoord(aResult[0], aResult[1], aResult[2]);
 212   if (myParamsU.Degree > myParamsV.Degree)
 213   {
 214     theTangentV.SetCoord(aResult[aDimension], aResult[aDimension + 1], aResult[aDimension + 2]);
 215     Standard_Integer aShift = aDimension<<1;
 216     theTangentU.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 217   }
 218   else
 219   {
 220     theTangentU.SetCoord(aResult[aDimension], aResult[aDimension + 1], aResult[aDimension + 2]);
 221     Standard_Integer aShift = aDimension<<1;
 222     theTangentV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 223   }
 224   theTangentU.Multiply(anInvU);
 225   theTangentV.Multiply(anInvV);
 226 }
 227
 228
 229 void BSplSLib_Cache::D2(const Standard_Real& theU,
 230                         const Standard_Real& theV,
 231                               gp_Pnt&        thePoint,
 232                               gp_Vec&        theTangentU,
 233                               gp_Vec&        theTangentV,
 234                               gp_Vec&        theCurvatureU,
 235                               gp_Vec&        theCurvatureV,
 236                               gp_Vec&        theCurvatureUV) const
 237 {
 238   Standard_Real aNewU = myParamsU.PeriodicNormalization (theU);
 239   Standard_Real aNewV = myParamsV.PeriodicNormalization (theV);
 240
 241   // BSplSLib uses different convention for span parameters than BSplCLib
 242   // (Start is in the middle of the span and length is half-span),
 243   // thus we need to amend them here
 244   Standard_Real aSpanLengthU = 0.5 * myParamsU.SpanLength;
 245   Standard_Real aSpanStartU = myParamsU.SpanStart + aSpanLengthU;
 246   Standard_Real aSpanLengthV = 0.5 * myParamsV.SpanLength;
 247   Standard_Real aSpanStartV = myParamsV.SpanStart + aSpanLengthV;
 248
 249   Standard_Real anInvU = 1.0 / aSpanLengthU;
 250   Standard_Real anInvV = 1.0 / aSpanLengthV;
 251   aNewU = (aNewU - aSpanStartU) * anInvU;
 252   aNewV = (aNewV - aSpanStartV) * anInvV;
 253
 254   Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
 255   Standard_Real aPntDeriv[36]; // result storage (point and derivative coordinates)
 256   for (Standard_Integer i = 0; i < 36; i++) aPntDeriv[i] = 0.0;
 257
 258   Standard_Integer aDimension = myIsRational ? 4 : 3;
 259   Standard_Integer aCacheCols = myPolesWeights->RowLength();
 260   Standard_Integer aMinMaxDegree[2] = {Min(myParamsU.Degree, myParamsV.Degree),
 261                                        Max(myParamsU.Degree, myParamsV.Degree)};
 262
 263   Standard_Real aParameters[2];
 264   if (myParamsU.Degree > myParamsV.Degree)
 265   {
 266     aParameters[0] = aNewV;
 267     aParameters[1] = aNewU;
 268   }
 269   else
 270   {
 271     aParameters[0] = aNewU;
 272     aParameters[1] = aNewV;
 273   }
 274
 275   NCollection_LocalArray<Standard_Real> aTransientCoeffs(3 * aCacheCols); // array for intermediate results
 276   // Calculating derivative to be evaluate and
 277   // nulling transient coefficients when max or min derivative is less than 2
 278   Standard_Integer aMinMaxDeriv[2] = {Min(2, aMinMaxDegree[0]),
 279                                       Min(2, aMinMaxDegree[1])};
 280   for (Standard_Integer i = aMinMaxDeriv[1] + 1; i < 3; i++)
 281   {
 282     Standard_Integer index = i * aCacheCols;
 283     for (Standard_Integer j = 0; j < aCacheCols; j++)
 284       aTransientCoeffs[index++] = 0.0;
 285   }
 286
 287   // Calculate intermediate values and derivatives of bivariate polynomial along variable with maximal degree
 288   PLib::EvalPolynomial(aParameters[1], aMinMaxDeriv[1], aMinMaxDegree[1],
 289                        aCacheCols, aPolesArray[0], aTransientCoeffs[0]);
 290
 291   // Calculate a point on surface and a derivatives along variable with minimal degree
 292   PLib::EvalPolynomial(aParameters[0], aMinMaxDeriv[0], aMinMaxDegree[0],
 293                        aDimension, aTransientCoeffs[0], aPntDeriv[0]);
 294
 295   // Calculate derivative along variable with maximal degree and mixed derivative
 296   PLib::EvalPolynomial(aParameters[0], 1, aMinMaxDegree[0], aDimension,
 297                        aTransientCoeffs[aCacheCols], aPntDeriv[3 * aDimension]);
 298
 299   // Calculate second derivative along variable with maximal degree
 300   PLib::NoDerivativeEvalPolynomial(aParameters[0], aMinMaxDegree[0], aDimension,
 301                                    aMinMaxDegree[0] * aDimension, aTransientCoeffs[aCacheCols<<1],
 302                                    aPntDeriv[6 * aDimension]);
 303
 304   Standard_Real* aResult = aPntDeriv;
 305   Standard_Real aTempStorage[36];
 306   if (myIsRational) // calculate derivatives divided by weight's derivatives
 307   {
 308     BSplSLib::RationalDerivative(2, 2, 2, 2, aPntDeriv[0], aTempStorage[0]);
 309     aResult = aTempStorage;
 310     aDimension--;
 311   }
 312
 313   thePoint.SetCoord(aResult[0], aResult[1], aResult[2]);
 314   if (myParamsU.Degree > myParamsV.Degree)
 315   {
 316     theTangentV.SetCoord(aResult[aDimension], aResult[aDimension + 1], aResult[aDimension + 2]);
 317     Standard_Integer aShift = aDimension<<1;
 318     theCurvatureV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 319     aShift += aDimension;
 320     theTangentU.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 321     aShift += aDimension;
 322     theCurvatureUV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 323     aShift += (aDimension << 1);
 324     theCurvatureU.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 325   }
 326   else
 327   {
 328     theTangentU.SetCoord(aResult[aDimension], aResult[aDimension + 1], aResult[aDimension + 2]);
 329     Standard_Integer aShift = aDimension<<1;
 330     theCurvatureU.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 331     aShift += aDimension;
 332     theTangentV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 333     aShift += aDimension;
 334     theCurvatureUV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 335     aShift += (aDimension << 1);
 336     theCurvatureV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]);
 337   }
 338   theTangentU.Multiply(anInvU);
 339   theTangentV.Multiply(anInvV);
 340   theCurvatureU.Multiply(anInvU * anInvU);
 341   theCurvatureV.Multiply(anInvV * anInvV);
 342   theCurvatureUV.Multiply(anInvU * anInvV);
 343 }
 344