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