94f71cad |
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 | |
94f71cad |
24 | |
92efcf78 |
25 | IMPLEMENT_STANDARD_RTTIEXT(BSplSLib_Cache,Standard_Transient) |
26 | |
94f71cad |
27 | //! Converts handle of array of Standard_Real into the pointer to Standard_Real |
35c0599a |
28 | static Standard_Real* ConvertArray(const Handle(TColStd_HArray2OfReal)& theHArray) |
94f71cad |
29 | { |
30 | const TColStd_Array2OfReal& anArray = theHArray->Array2(); |
31 | return (Standard_Real*) &(anArray(anArray.LowerRow(), anArray.LowerCol())); |
32 | } |
33 | |
94f71cad |
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, |
0e14656b |
40 | const TColStd_Array2OfReal* theWeights) |
0a96e0bb |
41 | : myIsRational(theWeights != NULL), |
42 | myParamsU (theDegreeU, thePeriodicU, theFlatKnotsU), |
43 | myParamsV (theDegreeV, thePeriodicV, theFlatKnotsV) |
94f71cad |
44 | { |
0a96e0bb |
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)); |
94f71cad |
49 | } |
50 | |
94f71cad |
51 | Standard_Boolean BSplSLib_Cache::IsCacheValid(Standard_Real theParameterU, |
52 | Standard_Real theParameterV) const |
53 | { |
0a96e0bb |
54 | return myParamsU.IsCacheValid (theParameterU) && |
55 | myParamsV.IsCacheValid (theParameterV); |
94f71cad |
56 | } |
57 | |
0a96e0bb |
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, |
0e14656b |
63 | const TColStd_Array2OfReal* theWeights) |
94f71cad |
64 | { |
65 | // Normalize the parameters for periodical B-splines |
0a96e0bb |
66 | Standard_Real aNewParamU = myParamsU.PeriodicNormalization (theParameterU); |
67 | Standard_Real aNewParamV = myParamsV.PeriodicNormalization (theParameterV); |
283b833c |
68 | |
0a96e0bb |
69 | myParamsU.LocateParameter (aNewParamU, theFlatKnotsU); |
70 | myParamsV.LocateParameter (aNewParamV, theFlatKnotsV); |
283b833c |
71 | |
0a96e0bb |
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; |
94f71cad |
79 | |
80 | // Calculate new cache data |
0a96e0bb |
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()); |
94f71cad |
88 | } |
89 | |
90 | |
91 | void BSplSLib_Cache::D0(const Standard_Real& theU, |
92 | const Standard_Real& theV, |
93 | gp_Pnt& thePoint) const |
94 | { |
0a96e0bb |
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; |
94f71cad |
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(); |
0a96e0bb |
114 | Standard_Integer aMinMaxDegree[2] = {Min(myParamsU.Degree, myParamsV.Degree), |
115 | Max(myParamsU.Degree, myParamsV.Degree)}; |
94f71cad |
116 | Standard_Real aParameters[2]; |
0a96e0bb |
117 | if (myParamsU.Degree > myParamsV.Degree) |
94f71cad |
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 | { |
0a96e0bb |
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; |
94f71cad |
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(); |
0a96e0bb |
174 | Standard_Integer aMinMaxDegree[2] = {Min(myParamsU.Degree, myParamsV.Degree), |
175 | Max(myParamsU.Degree, myParamsV.Degree)}; |
94f71cad |
176 | |
177 | Standard_Real aParameters[2]; |
0a96e0bb |
178 | if (myParamsU.Degree > myParamsV.Degree) |
94f71cad |
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]); |
0a96e0bb |
212 | if (myParamsU.Degree > myParamsV.Degree) |
94f71cad |
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 | { |
0a96e0bb |
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; |
94f71cad |
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(); |
0a96e0bb |
260 | Standard_Integer aMinMaxDegree[2] = {Min(myParamsU.Degree, myParamsV.Degree), |
261 | Max(myParamsU.Degree, myParamsV.Degree)}; |
94f71cad |
262 | |
263 | Standard_Real aParameters[2]; |
0a96e0bb |
264 | if (myParamsU.Degree > myParamsV.Degree) |
94f71cad |
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]); |
0a96e0bb |
314 | if (myParamsU.Degree > myParamsV.Degree) |
94f71cad |
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 | |