0026936: Drawbacks of inlining in new type system in OCCT 7.0 -- automatic
[occt.git] / src / BSplCLib / BSplCLib_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 <BSplCLib_Cache.hxx>
15 #include <BSplCLib.hxx>
16
17 #include <NCollection_LocalArray.hxx>
18
19 #include <TColgp_HArray1OfPnt.hxx>
20 #include <TColgp_HArray1OfPnt2d.hxx>
21 #include <TColStd_HArray1OfReal.hxx>
22 #include <TColStd_HArray2OfReal.hxx>
23
24
25 IMPLEMENT_STANDARD_RTTIEXT(BSplCLib_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 BSplCLib_Cache::BSplCLib_Cache()
36 {
37   myPolesWeights.Nullify();
38   myIsRational = Standard_False;
39   mySpanStart = 0.0;
40   mySpanLength = 0.0;
41   mySpanIndex = 0;
42   myDegree = 0;
43   myFlatKnots.Nullify();
44 }
45
46 BSplCLib_Cache::BSplCLib_Cache(const Standard_Integer&        theDegree,
47                                const Standard_Boolean&        thePeriodic,
48                                const TColStd_Array1OfReal&    theFlatKnots,
49                                const TColgp_Array1OfPnt2d&    thePoles2d,
50                                const TColStd_Array1OfReal*    theWeights)
51 {
52   Standard_Real aCacheParam = theFlatKnots.Value(theFlatKnots.Lower() + theDegree);
53   BuildCache(aCacheParam, theDegree, thePeriodic, 
54              theFlatKnots, thePoles2d, theWeights);
55 }
56
57 BSplCLib_Cache::BSplCLib_Cache(const Standard_Integer&        theDegree,
58                                const Standard_Boolean&        thePeriodic,
59                                const TColStd_Array1OfReal&    theFlatKnots,
60                                const TColgp_Array1OfPnt&      thePoles,
61                                const TColStd_Array1OfReal*    theWeights)
62 {
63   Standard_Real aCacheParam = theFlatKnots.Value(theFlatKnots.Lower() + theDegree);
64   BuildCache(aCacheParam, theDegree, thePeriodic, 
65              theFlatKnots, thePoles, theWeights);
66 }
67
68
69 Standard_Boolean BSplCLib_Cache::IsCacheValid(Standard_Real theParameter) const
70 {
71   Standard_Real aNewParam = theParameter;
72   if (!myFlatKnots.IsNull())
73     PeriodicNormalization(myFlatKnots->Array1(), aNewParam);
74
75   Standard_Real aDelta = aNewParam - mySpanStart;
76   return (aDelta >= 0.0 && (aDelta < mySpanLength || mySpanIndex == mySpanIndexMax));
77 }
78
79 void BSplCLib_Cache::PeriodicNormalization(const TColStd_Array1OfReal& theFlatKnots, 
80                                            Standard_Real& theParameter) const
81 {
82   Standard_Real aPeriod = theFlatKnots.Value(theFlatKnots.Upper() - myDegree) - 
83                           theFlatKnots.Value(myDegree + 1) ;
84   if (theParameter < theFlatKnots.Value(myDegree + 1))
85   {
86     Standard_Real aScale = IntegerPart(
87         (theFlatKnots.Value(myDegree + 1) - theParameter) / aPeriod);
88     theParameter += aPeriod * (aScale + 1.0);
89   }
90   if (theParameter > theFlatKnots.Value(theFlatKnots.Upper() - myDegree))
91   {
92     Standard_Real aScale = IntegerPart(
93         (theParameter - theFlatKnots.Value(theFlatKnots.Upper() - myDegree)) / aPeriod);
94     theParameter -= aPeriod * (aScale + 1.0);
95   }
96 }
97
98
99 void BSplCLib_Cache::BuildCache(const Standard_Real&           theParameter,
100                                 const Standard_Integer&        theDegree,
101                                 const Standard_Boolean&        thePeriodic,
102                                 const TColStd_Array1OfReal&    theFlatKnots,
103                                 const TColgp_Array1OfPnt2d&    thePoles2d,
104                                 const TColStd_Array1OfReal*    theWeights)
105 {
106   // Normalize theParameter for periodical B-splines
107   Standard_Real aNewParam = theParameter;
108   if (thePeriodic)
109   {
110     PeriodicNormalization(theFlatKnots, aNewParam);
111     myFlatKnots = new TColStd_HArray1OfReal(1, theFlatKnots.Length());
112     myFlatKnots->ChangeArray1() = theFlatKnots;
113   }
114   else if (!myFlatKnots.IsNull()) // Periodical curve became non-periodical
115     myFlatKnots.Nullify();
116
117   // Change the size of cached data if needed
118   myIsRational = (theWeights != NULL);
119   Standard_Integer aPWColNumber = myIsRational ? 3 : 2;
120   if (theDegree > myDegree)
121     myPolesWeights = new TColStd_HArray2OfReal(1, theDegree + 1, 1, aPWColNumber);
122
123   myDegree = theDegree;
124   mySpanIndex = 0;
125   BSplCLib::LocateParameter(theDegree, theFlatKnots, BSplCLib::NoMults(), 
126                             aNewParam, thePeriodic, mySpanIndex, aNewParam);
127   mySpanStart  = theFlatKnots.Value(mySpanIndex);
128   mySpanLength = theFlatKnots.Value(mySpanIndex + 1) - mySpanStart;
129   mySpanIndexMax = theFlatKnots.Length() - 1 - theDegree;
130
131   // Calculate new cache data
132   BSplCLib::BuildCache(mySpanStart, mySpanLength, thePeriodic, theDegree, 
133                        theFlatKnots, thePoles2d, theWeights, 
134                        myPolesWeights->ChangeArray2());
135 }
136
137 void BSplCLib_Cache::BuildCache(const Standard_Real&           theParameter,
138                                 const Standard_Integer&        theDegree,
139                                 const Standard_Boolean&        thePeriodic,
140                                 const TColStd_Array1OfReal&    theFlatKnots,
141                                 const TColgp_Array1OfPnt&      thePoles,
142                                 const TColStd_Array1OfReal*    theWeights)
143 {
144   // Create list of knots with repetitions and normalize theParameter for periodical B-splines
145   Standard_Real aNewParam = theParameter;
146   if (thePeriodic)
147   {
148     PeriodicNormalization(theFlatKnots, aNewParam);
149     myFlatKnots = new TColStd_HArray1OfReal(1, theFlatKnots.Length());
150     myFlatKnots->ChangeArray1() = theFlatKnots;
151   }
152   else if (!myFlatKnots.IsNull()) // Periodical curve became non-periodical
153     myFlatKnots.Nullify();
154
155   // Change the size of cached data if needed
156   myIsRational = (theWeights != NULL);
157   Standard_Integer aPWColNumber = myIsRational ? 4 : 3;
158   if (theDegree > myDegree)
159     myPolesWeights = new TColStd_HArray2OfReal(1, theDegree + 1, 1, aPWColNumber);
160
161   myDegree = theDegree;
162   mySpanIndex = 0;
163   BSplCLib::LocateParameter(theDegree, theFlatKnots, BSplCLib::NoMults(), 
164                             aNewParam, thePeriodic, mySpanIndex, aNewParam);
165   mySpanStart  = theFlatKnots.Value(mySpanIndex);
166   mySpanLength = theFlatKnots.Value(mySpanIndex + 1) - mySpanStart;
167   mySpanIndexMax = theFlatKnots.Length() - 1 - theDegree;
168
169   // Calculate new cache data
170   BSplCLib::BuildCache(mySpanStart, mySpanLength, thePeriodic, theDegree, 
171                        theFlatKnots, thePoles, theWeights, 
172                        myPolesWeights->ChangeArray2());
173 }
174
175
176 void BSplCLib_Cache::CalculateDerivative(const Standard_Real&    theParameter, 
177                                          const Standard_Integer& theDerivative, 
178                                                Standard_Real&    theDerivArray) const
179 {
180   Standard_Real aNewParameter = theParameter;
181   if (!myFlatKnots.IsNull()) // B-spline is periodical
182     PeriodicNormalization(myFlatKnots->Array1(), aNewParameter);
183   aNewParameter = (aNewParameter - mySpanStart) / mySpanLength;
184
185   Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
186   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
187
188   // Temporary container. The maximal size of this container is defined by:
189   //    1) maximal derivative for cache evaluation, which is 3, plus one row for function values, 
190   //    2) and maximal dimension of the point, which is 3, plus one column for weights.
191   Standard_Real aTmpContainer[16];
192
193   // When the PLib::RationaDerivative needs to be called, use temporary container
194   Standard_Real* aPntDeriv = myIsRational ? aTmpContainer : &theDerivArray;
195
196   // When the degree of curve is lesser than the requested derivative,
197   // nullify array cells corresponding to greater derivatives
198   Standard_Integer aDerivative = theDerivative;
199   if (myDegree < theDerivative)
200   {
201     aDerivative = myDegree;
202     for (Standard_Integer ind = myDegree * aDimension; ind < (theDerivative + 1) * aDimension; ind++)
203     {
204       aPntDeriv[ind] = 0.0;
205       (&theDerivArray)[ind] = 0.0; // should be cleared separately, because aPntDeriv may look to another memory area
206     }
207   }
208
209   PLib::EvalPolynomial(aNewParameter, aDerivative, myDegree, aDimension, 
210                        aPolesArray[0], aPntDeriv[0]);
211   // Unnormalize derivatives since those are computed normalized
212   Standard_Real aFactor = 1.0;
213   for (Standard_Integer deriv = 1; deriv <= aDerivative; deriv++)
214   {
215     aFactor /= mySpanLength;
216     for (Standard_Integer ind = 0; ind < aDimension; ind++)
217       aPntDeriv[aDimension * deriv + ind] *= aFactor;
218   }
219
220   if (myIsRational) // calculate derivatives divided by weights derivatives
221     PLib::RationalDerivative(aDerivative, aDerivative, aDimension - 1, aPntDeriv[0], theDerivArray);
222 }
223
224
225 void BSplCLib_Cache::D0(const Standard_Real& theParameter, gp_Pnt2d& thePoint) const
226 {
227   Standard_Real aNewParameter = theParameter;
228   if (!myFlatKnots.IsNull()) // B-spline is periodical
229     PeriodicNormalization(myFlatKnots->Array1(), aNewParameter);
230   aNewParameter = (aNewParameter - mySpanStart) / mySpanLength;
231
232   Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
233   Standard_Real aPoint[4];
234   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
235
236   PLib::NoDerivativeEvalPolynomial(aNewParameter, myDegree,
237                                    aDimension, myDegree * aDimension,
238                                    aPolesArray[0], aPoint[0]);
239
240   thePoint.SetCoord(aPoint[0], aPoint[1]);
241   if (myIsRational)
242     thePoint.ChangeCoord().Divide(aPoint[2]);
243 }
244
245 void BSplCLib_Cache::D0(const Standard_Real& theParameter, gp_Pnt& thePoint) const
246 {
247   Standard_Real aNewParameter = theParameter;
248   if (!myFlatKnots.IsNull()) // B-spline is periodical
249     PeriodicNormalization(myFlatKnots->Array1(), aNewParameter);
250   aNewParameter = (aNewParameter - mySpanStart) / mySpanLength;
251
252   Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
253   Standard_Real aPoint[4];
254   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
255
256   PLib::NoDerivativeEvalPolynomial(aNewParameter, myDegree,
257                                    aDimension, myDegree * aDimension,
258                                    aPolesArray[0], aPoint[0]);
259
260   thePoint.SetCoord(aPoint[0], aPoint[1], aPoint[2]);
261   if (myIsRational)
262     thePoint.ChangeCoord().Divide(aPoint[3]);
263 }
264
265
266 void BSplCLib_Cache::D1(const Standard_Real& theParameter, gp_Pnt2d& thePoint, gp_Vec2d& theTangent) const
267 {
268   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
269   Standard_Real aPntDeriv[8]; // result storage (point and derivative coordinates)
270
271   this->CalculateDerivative(theParameter, 1, aPntDeriv[0]);
272   if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
273     aDimension -= 1;
274
275   thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1]);
276   theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1]);
277 }
278
279 void BSplCLib_Cache::D1(const Standard_Real& theParameter, gp_Pnt& thePoint, gp_Vec& theTangent) const
280 {
281   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
282   Standard_Real aPntDeriv[8]; // result storage (point and derivative coordinates)
283
284   this->CalculateDerivative(theParameter, 1, aPntDeriv[0]);
285   if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
286     aDimension -= 1;
287
288   thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1], aPntDeriv[2]);
289   theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1], aPntDeriv[aDimension + 2]);
290 }
291
292 void BSplCLib_Cache::D2(const Standard_Real& theParameter, gp_Pnt2d& thePoint, gp_Vec2d& theTangent, gp_Vec2d& theCurvature) const
293 {
294   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
295   Standard_Real aPntDeriv[12]; // result storage (point and derivatives coordinates)
296
297   this->CalculateDerivative(theParameter, 2, aPntDeriv[0]);
298   if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
299     aDimension -= 1;
300
301   thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1]);
302   theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1]);
303   theCurvature.SetCoord(aPntDeriv[aDimension<<1], aPntDeriv[(aDimension<<1) + 1]);
304 }
305
306 void BSplCLib_Cache::D2(const Standard_Real& theParameter, gp_Pnt& thePoint, gp_Vec& theTangent, gp_Vec& theCurvature) const
307 {
308   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
309   Standard_Real aPntDeriv[12]; // result storage (point and derivatives coordinates)
310
311   this->CalculateDerivative(theParameter, 2, aPntDeriv[0]);
312   if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
313     aDimension -= 1;
314
315   thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1], aPntDeriv[2]);
316   theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1], aPntDeriv[aDimension + 2]);
317   theCurvature.SetCoord(aPntDeriv[aDimension<<1], aPntDeriv[(aDimension<<1) + 1], aPntDeriv[(aDimension<<1) + 2]);
318 }
319
320
321 void BSplCLib_Cache::D3(const Standard_Real& theParameter, 
322                               gp_Pnt2d&      thePoint, 
323                               gp_Vec2d&      theTangent, 
324                               gp_Vec2d&      theCurvature,
325                               gp_Vec2d&      theTorsion) const
326 {
327   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
328   Standard_Real aPntDeriv[16]; // result storage (point and derivatives coordinates)
329
330   this->CalculateDerivative(theParameter, 3, aPntDeriv[0]);
331   if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
332     aDimension -= 1;
333
334   thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1]);
335   theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1]);
336   Standard_Integer aShift = aDimension << 1;
337   theCurvature.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1]);
338   aShift += aDimension;
339   theTorsion.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1]);
340 }
341
342 void BSplCLib_Cache::D3(const Standard_Real& theParameter, 
343                               gp_Pnt&        thePoint, 
344                               gp_Vec&        theTangent, 
345                               gp_Vec&        theCurvature,
346                               gp_Vec&        theTorsion) const
347 {
348   Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
349   Standard_Real aPntDeriv[16]; // result storage (point and derivatives coordinates)
350
351   this->CalculateDerivative(theParameter, 3, aPntDeriv[0]);
352   if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
353     aDimension -= 1;
354
355   thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1], aPntDeriv[2]);
356   theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1], aPntDeriv[aDimension + 2]);
357   Standard_Integer aShift = aDimension << 1;
358   theCurvature.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1], aPntDeriv[aShift + 2]);
359   aShift += aDimension;
360   theTorsion.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1], aPntDeriv[aShift + 2]);
361 }
362