dd9ba2eaac7c6d354409f76f639df5776960b6dc
[occt.git] / 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