1 // Copyright (c) 2014 OPEN CASCADE SAS
3 // This file is part of Open CASCADE Technology software library.
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.
11 // Alternatively, this file may be used under the terms of Open CASCADE
12 // commercial license or contractual agreement.
14 #include <BSplSLib_Cache.hxx>
15 #include <BSplSLib.hxx>
17 #include <NCollection_LocalArray.hxx>
19 #include <TColgp_HArray2OfPnt.hxx>
20 #include <TColStd_HArray1OfInteger.hxx>
21 #include <TColStd_HArray1OfReal.hxx>
22 #include <TColStd_HArray2OfReal.hxx>
24 IMPLEMENT_STANDARD_HANDLE(BSplSLib_Cache, Standard_Transient)
25 IMPLEMENT_STANDARD_RTTIEXT(BSplSLib_Cache, Standard_Transient)
27 //! Converts handle of array of Standard_Real into the pointer to Standard_Real
28 static Standard_Real* ConvertArray(const Handle_TColStd_HArray2OfReal& theHArray)
30 const TColStd_Array2OfReal& anArray = theHArray->Array2();
31 return (Standard_Real*) &(anArray(anArray.LowerRow(), anArray.LowerCol()));
35 BSplSLib_Cache::BSplSLib_Cache()
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();
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)
56 Standard_Real aU = theFlatKnotsU.Value(theFlatKnotsU.Lower() + theDegreeU);
57 Standard_Real aV = theFlatKnotsV.Value(theFlatKnotsV.Lower() + theDegreeV);
60 theDegreeU, thePeriodicU, theFlatKnotsU,
61 theDegreeV, thePeriodicV, theFlatKnotsV,
62 thePoles, theWeights);
66 Standard_Boolean BSplSLib_Cache::IsCacheValid(Standard_Real theParameterU,
67 Standard_Real theParameterV) const
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);
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]));
82 void BSplSLib_Cache::PeriodicNormalization(const Standard_Integer& theDegree,
83 const TColStd_Array1OfReal& theFlatKnots,
84 Standard_Real& theParameter) const
86 Standard_Real aPeriod = theFlatKnots.Value(theFlatKnots.Upper() - theDegree) -
87 theFlatKnots.Value(theDegree + 1) ;
88 if (theParameter < theFlatKnots.Value(theDegree + 1))
90 Standard_Real aScale = IntegerPart(
91 (theFlatKnots.Value(theDegree + 1) - theParameter) / aPeriod);
92 theParameter += aPeriod * (aScale + 1.0);
94 if (theParameter > theFlatKnots.Value(theFlatKnots.Upper() - theDegree))
96 Standard_Real aScale = IntegerPart(
97 (theParameter - theFlatKnots.Value(theFlatKnots.Upper() - theDegree)) / aPeriod);
98 theParameter -= aPeriod * (aScale + 1.0);
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)
114 // Normalize the parameters for periodical B-splines
115 Standard_Real aNewParamU = theParameterU;
118 PeriodicNormalization(theDegreeU, theFlatKnotsU, aNewParamU);
119 myFlatKnots[0] = new TColStd_HArray1OfReal(1, theFlatKnotsU.Length());
120 myFlatKnots[0]->ChangeArray1() = theFlatKnotsU;
122 else if (!myFlatKnots[0].IsNull()) // Periodical curve became non-periodical
123 myFlatKnots[0].Nullify();
125 Standard_Real aNewParamV = theParameterV;
128 PeriodicNormalization(theDegreeV, theFlatKnotsV, aNewParamV);
129 myFlatKnots[1] = new TColStd_HArray1OfReal(1, theFlatKnotsV.Length());
130 myFlatKnots[1]->ChangeArray1() = theFlatKnotsV;
132 else if (!myFlatKnots[1].IsNull()) // Periodical curve became non-periodical
133 myFlatKnots[1].Nullify();
135 Standard_Integer aMinDegree = Min(theDegreeU, theDegreeV);
136 Standard_Integer aMaxDegree = Max(theDegreeU, theDegreeV);
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));
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);
152 // Protection against Out of Range exception.
153 if (mySpanIndex[0] >= theFlatKnotsU.Length()) {
154 mySpanIndex[0] = theFlatKnotsU.Length() - 1;
157 mySpanLength[0] = (theFlatKnotsU.Value(mySpanIndex[0] + 1) - theFlatKnotsU.Value(mySpanIndex[0])) * 0.5;
158 mySpanStart[0] = theFlatKnotsU.Value(mySpanIndex[0]) + mySpanLength[0];
160 // Protection against Out of Range exception.
161 if (mySpanIndex[1] >= theFlatKnotsV.Length()) {
162 mySpanIndex[1] = theFlatKnotsV.Length() - 1;
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;
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());
181 void BSplSLib_Cache::D0(const Standard_Real& theU,
182 const Standard_Real& theV,
183 gp_Pnt& thePoint) const
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];
194 Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
195 Standard_Real aPoint[4];
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])
204 aParameters[0] = aNewV;
205 aParameters[1] = aNewU;
209 aParameters[0] = aNewU;
210 aParameters[1] = aNewV;
213 NCollection_LocalArray<Standard_Real> aTransientCoeffs(aCacheCols); // array for intermediate results
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]);
220 // Calculate total value
221 PLib::NoDerivativeEvalPolynomial(aParameters[0], aMinMaxDegree[0],
222 aDimension, aDimension * aMinMaxDegree[0],
223 aTransientCoeffs[0], aPoint[0]);
225 thePoint.SetCoord(aPoint[0], aPoint[1], aPoint[2]);
227 thePoint.ChangeCoord().Divide(aPoint[3]);
231 void BSplSLib_Cache::D1(const Standard_Real& theU,
232 const Standard_Real& theV,
235 gp_Vec& theTangentV) const
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;
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;
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])};
257 Standard_Real aParameters[2];
258 if (myDegree[0] > myDegree[1])
260 aParameters[0] = aNewV;
261 aParameters[1] = aNewU;
265 aParameters[0] = aNewU;
266 aParameters[1] = aNewV;
269 NCollection_LocalArray<Standard_Real> aTransientCoeffs(aCacheCols<<1); // array for intermediate results
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]);
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]);
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]);
282 Standard_Real* aResult = aPntDeriv;
283 Standard_Real aTempStorage[12];
284 if (myIsRational) // calculate derivatives divided by weight's derivatives
286 BSplSLib::RationalDerivative(1, 1, 1, 1, aPntDeriv[0], aTempStorage[0]);
287 aResult = aTempStorage;
291 thePoint.SetCoord(aResult[0], aResult[1], aResult[2]);
292 if (myDegree[0] > myDegree[1])
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]);
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]);
304 theTangentU.Multiply(anInvU);
305 theTangentV.Multiply(anInvV);
309 void BSplSLib_Cache::D2(const Standard_Real& theU,
310 const Standard_Real& theV,
314 gp_Vec& theCurvatureU,
315 gp_Vec& theCurvatureV,
316 gp_Vec& theCurvatureUV) const
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;
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;
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])};
338 Standard_Real aParameters[2];
339 if (myDegree[0] > myDegree[1])
341 aParameters[0] = aNewV;
342 aParameters[1] = aNewU;
346 aParameters[0] = aNewU;
347 aParameters[1] = aNewV;
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++)
357 Standard_Integer index = i * aCacheCols;
358 for (Standard_Integer j = 0; j < aCacheCols; j++)
359 aTransientCoeffs[index++] = 0.0;
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]);
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]);
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]);
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]);
379 Standard_Real* aResult = aPntDeriv;
380 Standard_Real aTempStorage[36];
381 if (myIsRational) // calculate derivatives divided by weight's derivatives
383 BSplSLib::RationalDerivative(2, 2, 2, 2, aPntDeriv[0], aTempStorage[0]);
384 aResult = aTempStorage;
388 thePoint.SetCoord(aResult[0], aResult[1], aResult[2]);
389 if (myDegree[0] > myDegree[1])
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]);
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]);
413 theTangentU.Multiply(anInvU);
414 theTangentV.Multiply(anInvV);
415 theCurvatureU.Multiply(anInvU * anInvU);
416 theCurvatureV.Multiply(anInvV * anInvV);
417 theCurvatureUV.Multiply(anInvU * anInvV);