0029769: Uninitialized data with BSplCLib_Cache, BSplSLib_Cache
[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 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