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 | |
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, |
0e14656b |
54 | const TColStd_Array2OfReal* theWeights) |
94f71cad |
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] && (aDelta0 < mySpanLength[0] || mySpanIndex[0] == mySpanIndexMax[0]) && |
79 | aDelta1 >= -mySpanLength[1] && (aDelta1 < mySpanLength[1] || mySpanIndex[1] == mySpanIndexMax[1])); |
80 | } |
81 | |
82 | void BSplSLib_Cache::PeriodicNormalization(const Standard_Integer& theDegree, |
83 | const TColStd_Array1OfReal& theFlatKnots, |
84 | Standard_Real& theParameter) const |
85 | { |
86 | Standard_Real aPeriod = theFlatKnots.Value(theFlatKnots.Upper() - theDegree) - |
87 | theFlatKnots.Value(theDegree + 1) ; |
88 | if (theParameter < theFlatKnots.Value(theDegree + 1)) |
89 | { |
90 | Standard_Real aScale = IntegerPart( |
91 | (theFlatKnots.Value(theDegree + 1) - theParameter) / aPeriod); |
92 | theParameter += aPeriod * (aScale + 1.0); |
93 | } |
94 | if (theParameter > theFlatKnots.Value(theFlatKnots.Upper() - theDegree)) |
95 | { |
96 | Standard_Real aScale = IntegerPart( |
97 | (theParameter - theFlatKnots.Value(theFlatKnots.Upper() - theDegree)) / aPeriod); |
98 | theParameter -= aPeriod * (aScale + 1.0); |
99 | } |
100 | } |
101 | |
102 | |
103 | void BSplSLib_Cache::BuildCache(const Standard_Real& theParameterU, |
104 | const Standard_Real& theParameterV, |
105 | const Standard_Integer& theDegreeU, |
106 | const Standard_Boolean& thePeriodicU, |
107 | const TColStd_Array1OfReal& theFlatKnotsU, |
108 | const Standard_Integer& theDegreeV, |
109 | const Standard_Boolean& thePeriodicV, |
110 | const TColStd_Array1OfReal& theFlatKnotsV, |
111 | const TColgp_Array2OfPnt& thePoles, |
0e14656b |
112 | const TColStd_Array2OfReal* theWeights) |
94f71cad |
113 | { |
114 | // Normalize the parameters for periodical B-splines |
115 | Standard_Real aNewParamU = theParameterU; |
116 | if (thePeriodicU) |
117 | { |
118 | PeriodicNormalization(theDegreeU, theFlatKnotsU, aNewParamU); |
119 | myFlatKnots[0] = new TColStd_HArray1OfReal(1, theFlatKnotsU.Length()); |
120 | myFlatKnots[0]->ChangeArray1() = theFlatKnotsU; |
121 | } |
122 | else if (!myFlatKnots[0].IsNull()) // Periodical curve became non-periodical |
123 | myFlatKnots[0].Nullify(); |
124 | |
125 | Standard_Real aNewParamV = theParameterV; |
126 | if (thePeriodicV) |
127 | { |
128 | PeriodicNormalization(theDegreeV, theFlatKnotsV, aNewParamV); |
129 | myFlatKnots[1] = new TColStd_HArray1OfReal(1, theFlatKnotsV.Length()); |
130 | myFlatKnots[1]->ChangeArray1() = theFlatKnotsV; |
131 | } |
132 | else if (!myFlatKnots[1].IsNull()) // Periodical curve became non-periodical |
133 | myFlatKnots[1].Nullify(); |
134 | |
135 | Standard_Integer aMinDegree = Min(theDegreeU, theDegreeV); |
136 | Standard_Integer aMaxDegree = Max(theDegreeU, theDegreeV); |
137 | |
138 | // Change the size of cached data if needed |
0e14656b |
139 | myIsRational = (theWeights != NULL); |
94f71cad |
140 | Standard_Integer aPWColNumber = myIsRational ? 4 : 3; |
141 | if (theDegreeU > myDegree[0] || theDegreeV > myDegree[1]) |
142 | myPolesWeights = new TColStd_HArray2OfReal(1, aMaxDegree + 1, 1, aPWColNumber * (aMinDegree + 1)); |
143 | |
144 | myDegree[0] = theDegreeU; |
145 | myDegree[1] = theDegreeV; |
146 | mySpanIndex[0] = mySpanIndex[1] = 0; |
147 | BSplCLib::LocateParameter(theDegreeU, theFlatKnotsU, BSplCLib::NoMults(), aNewParamU, |
148 | thePeriodicU, mySpanIndex[0], aNewParamU); |
149 | BSplCLib::LocateParameter(theDegreeV, theFlatKnotsV, BSplCLib::NoMults(), aNewParamV, |
150 | thePeriodicV, mySpanIndex[1], aNewParamV); |
283b833c |
151 | |
152 | // Protection against Out of Range exception. |
153 | if (mySpanIndex[0] >= theFlatKnotsU.Length()) { |
154 | mySpanIndex[0] = theFlatKnotsU.Length() - 1; |
155 | } |
156 | |
94f71cad |
157 | mySpanLength[0] = (theFlatKnotsU.Value(mySpanIndex[0] + 1) - theFlatKnotsU.Value(mySpanIndex[0])) * 0.5; |
158 | mySpanStart[0] = theFlatKnotsU.Value(mySpanIndex[0]) + mySpanLength[0]; |
283b833c |
159 | |
160 | // Protection against Out of Range exception. |
161 | if (mySpanIndex[1] >= theFlatKnotsV.Length()) { |
162 | mySpanIndex[1] = theFlatKnotsV.Length() - 1; |
163 | } |
164 | |
94f71cad |
165 | mySpanLength[1] = (theFlatKnotsV.Value(mySpanIndex[1] + 1) - theFlatKnotsV.Value(mySpanIndex[1])) * 0.5; |
166 | mySpanStart[1] = theFlatKnotsV.Value(mySpanIndex[1]) + mySpanLength[1]; |
167 | mySpanIndexMax[0] = theFlatKnotsU.Length() - 1 - theDegreeU; |
168 | mySpanIndexMax[1] = theFlatKnotsV.Length() - 1 - theDegreeV; |
169 | |
170 | // Calculate new cache data |
171 | BSplSLib::BuildCache(mySpanStart[0], mySpanStart[1], |
172 | mySpanLength[0], mySpanLength[1], |
173 | thePeriodicU, thePeriodicV, |
174 | theDegreeU, theDegreeV, |
175 | mySpanIndex[0], mySpanIndex[1], |
176 | theFlatKnotsU, theFlatKnotsV, |
177 | thePoles, theWeights, myPolesWeights->ChangeArray2()); |
178 | } |
179 | |
180 | |
181 | void BSplSLib_Cache::D0(const Standard_Real& theU, |
182 | const Standard_Real& theV, |
183 | gp_Pnt& thePoint) const |
184 | { |
185 | Standard_Real aNewU = theU; |
186 | Standard_Real aNewV = theV; |
187 | if (!myFlatKnots[0].IsNull()) // B-spline is U-periodical |
188 | PeriodicNormalization(myDegree[0], myFlatKnots[0]->Array1(), aNewU); |
189 | aNewU = (aNewU - mySpanStart[0]) / mySpanLength[0]; |
190 | if (!myFlatKnots[1].IsNull()) // B-spline is V-periodical |
191 | PeriodicNormalization(myDegree[1], myFlatKnots[1]->Array1(), aNewV); |
192 | aNewV = (aNewV - mySpanStart[1]) / mySpanLength[1]; |
193 | |
194 | Standard_Real* aPolesArray = ConvertArray(myPolesWeights); |
195 | Standard_Real aPoint[4]; |
196 | |
197 | Standard_Integer aDimension = myIsRational ? 4 : 3; |
198 | Standard_Integer aCacheCols = myPolesWeights->RowLength(); |
199 | Standard_Integer aMinMaxDegree[2] = {Min(myDegree[0], myDegree[1]), |
200 | Max(myDegree[0], myDegree[1])}; |
201 | Standard_Real aParameters[2]; |
202 | if (myDegree[0] > myDegree[1]) |
203 | { |
204 | aParameters[0] = aNewV; |
205 | aParameters[1] = aNewU; |
206 | } |
207 | else |
208 | { |
209 | aParameters[0] = aNewU; |
210 | aParameters[1] = aNewV; |
211 | } |
212 | |
213 | NCollection_LocalArray<Standard_Real> aTransientCoeffs(aCacheCols); // array for intermediate results |
214 | |
215 | // Calculate intermediate value of cached polynomial along columns |
216 | PLib::NoDerivativeEvalPolynomial(aParameters[1], aMinMaxDegree[1], |
217 | aCacheCols, aMinMaxDegree[1] * aCacheCols, |
218 | aPolesArray[0], aTransientCoeffs[0]); |
219 | |
220 | // Calculate total value |
221 | PLib::NoDerivativeEvalPolynomial(aParameters[0], aMinMaxDegree[0], |
222 | aDimension, aDimension * aMinMaxDegree[0], |
223 | aTransientCoeffs[0], aPoint[0]); |
224 | |
225 | thePoint.SetCoord(aPoint[0], aPoint[1], aPoint[2]); |
226 | if (myIsRational) |
227 | thePoint.ChangeCoord().Divide(aPoint[3]); |
228 | } |
229 | |
230 | |
231 | void BSplSLib_Cache::D1(const Standard_Real& theU, |
232 | const Standard_Real& theV, |
233 | gp_Pnt& thePoint, |
234 | gp_Vec& theTangentU, |
235 | gp_Vec& theTangentV) const |
236 | { |
237 | Standard_Real aNewU = theU; |
238 | Standard_Real aNewV = theV; |
239 | Standard_Real anInvU = 1.0 / mySpanLength[0]; |
240 | Standard_Real anInvV = 1.0 / mySpanLength[1]; |
241 | if (!myFlatKnots[0].IsNull()) // B-spline is U-periodical |
242 | PeriodicNormalization(myDegree[0], myFlatKnots[0]->Array1(), aNewU); |
243 | aNewU = (aNewU - mySpanStart[0]) * anInvU; |
244 | if (!myFlatKnots[1].IsNull()) // B-spline is V-periodical |
245 | PeriodicNormalization(myDegree[1], myFlatKnots[1]->Array1(), aNewV); |
246 | aNewV = (aNewV - mySpanStart[1]) * anInvV; |
247 | |
248 | Standard_Real* aPolesArray = ConvertArray(myPolesWeights); |
249 | Standard_Real aPntDeriv[16]; // result storage (point and derivative coordinates) |
250 | for (Standard_Integer i = 0; i< 16; i++) aPntDeriv[i] = 0.0; |
251 | |
252 | Standard_Integer aDimension = myIsRational ? 4 : 3; |
253 | Standard_Integer aCacheCols = myPolesWeights->RowLength(); |
254 | Standard_Integer aMinMaxDegree[2] = {Min(myDegree[0], myDegree[1]), |
255 | Max(myDegree[0], myDegree[1])}; |
256 | |
257 | Standard_Real aParameters[2]; |
258 | if (myDegree[0] > myDegree[1]) |
259 | { |
260 | aParameters[0] = aNewV; |
261 | aParameters[1] = aNewU; |
262 | } |
263 | else |
264 | { |
265 | aParameters[0] = aNewU; |
266 | aParameters[1] = aNewV; |
267 | } |
268 | |
269 | NCollection_LocalArray<Standard_Real> aTransientCoeffs(aCacheCols<<1); // array for intermediate results |
270 | |
271 | // Calculate intermediate values and derivatives of bivariate polynomial along variable with maximal degree |
272 | PLib::EvalPolynomial(aParameters[1], 1, aMinMaxDegree[1], aCacheCols, aPolesArray[0], aTransientCoeffs[0]); |
273 | |
274 | // Calculate a point on surface and a derivative along variable with minimal degree |
275 | PLib::EvalPolynomial(aParameters[0], 1, aMinMaxDegree[0], aDimension, aTransientCoeffs[0], aPntDeriv[0]); |
276 | |
277 | // Calculate derivative along variable with maximal degree |
278 | PLib::NoDerivativeEvalPolynomial(aParameters[0], aMinMaxDegree[0], aDimension, |
279 | aMinMaxDegree[0] * aDimension, aTransientCoeffs[aCacheCols], |
280 | aPntDeriv[aDimension<<1]); |
281 | |
282 | Standard_Real* aResult = aPntDeriv; |
283 | Standard_Real aTempStorage[12]; |
284 | if (myIsRational) // calculate derivatives divided by weight's derivatives |
285 | { |
286 | BSplSLib::RationalDerivative(1, 1, 1, 1, aPntDeriv[0], aTempStorage[0]); |
287 | aResult = aTempStorage; |
288 | aDimension--; |
289 | } |
290 | |
291 | thePoint.SetCoord(aResult[0], aResult[1], aResult[2]); |
292 | if (myDegree[0] > myDegree[1]) |
293 | { |
294 | theTangentV.SetCoord(aResult[aDimension], aResult[aDimension + 1], aResult[aDimension + 2]); |
295 | Standard_Integer aShift = aDimension<<1; |
296 | theTangentU.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]); |
297 | } |
298 | else |
299 | { |
300 | theTangentU.SetCoord(aResult[aDimension], aResult[aDimension + 1], aResult[aDimension + 2]); |
301 | Standard_Integer aShift = aDimension<<1; |
302 | theTangentV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]); |
303 | } |
304 | theTangentU.Multiply(anInvU); |
305 | theTangentV.Multiply(anInvV); |
306 | } |
307 | |
308 | |
309 | void BSplSLib_Cache::D2(const Standard_Real& theU, |
310 | const Standard_Real& theV, |
311 | gp_Pnt& thePoint, |
312 | gp_Vec& theTangentU, |
313 | gp_Vec& theTangentV, |
314 | gp_Vec& theCurvatureU, |
315 | gp_Vec& theCurvatureV, |
316 | gp_Vec& theCurvatureUV) const |
317 | { |
318 | Standard_Real aNewU = theU; |
319 | Standard_Real aNewV = theV; |
320 | Standard_Real anInvU = 1.0 / mySpanLength[0]; |
321 | Standard_Real anInvV = 1.0 / mySpanLength[1]; |
322 | if (!myFlatKnots[0].IsNull()) // B-spline is U-periodical |
323 | PeriodicNormalization(myDegree[0], myFlatKnots[0]->Array1(), aNewU); |
324 | aNewU = (aNewU - mySpanStart[0]) * anInvU; |
325 | if (!myFlatKnots[1].IsNull()) // B-spline is V-periodical |
326 | PeriodicNormalization(myDegree[1], myFlatKnots[1]->Array1(), aNewV); |
327 | aNewV = (aNewV - mySpanStart[1]) * anInvV; |
328 | |
329 | Standard_Real* aPolesArray = ConvertArray(myPolesWeights); |
330 | Standard_Real aPntDeriv[36]; // result storage (point and derivative coordinates) |
331 | for (Standard_Integer i = 0; i < 36; i++) aPntDeriv[i] = 0.0; |
332 | |
333 | Standard_Integer aDimension = myIsRational ? 4 : 3; |
334 | Standard_Integer aCacheCols = myPolesWeights->RowLength(); |
335 | Standard_Integer aMinMaxDegree[2] = {Min(myDegree[0], myDegree[1]), |
336 | Max(myDegree[0], myDegree[1])}; |
337 | |
338 | Standard_Real aParameters[2]; |
339 | if (myDegree[0] > myDegree[1]) |
340 | { |
341 | aParameters[0] = aNewV; |
342 | aParameters[1] = aNewU; |
343 | } |
344 | else |
345 | { |
346 | aParameters[0] = aNewU; |
347 | aParameters[1] = aNewV; |
348 | } |
349 | |
350 | NCollection_LocalArray<Standard_Real> aTransientCoeffs(3 * aCacheCols); // array for intermediate results |
351 | // Calculating derivative to be evaluate and |
352 | // nulling transient coefficients when max or min derivative is less than 2 |
353 | Standard_Integer aMinMaxDeriv[2] = {Min(2, aMinMaxDegree[0]), |
354 | Min(2, aMinMaxDegree[1])}; |
355 | for (Standard_Integer i = aMinMaxDeriv[1] + 1; i < 3; i++) |
356 | { |
357 | Standard_Integer index = i * aCacheCols; |
358 | for (Standard_Integer j = 0; j < aCacheCols; j++) |
359 | aTransientCoeffs[index++] = 0.0; |
360 | } |
361 | |
362 | // Calculate intermediate values and derivatives of bivariate polynomial along variable with maximal degree |
363 | PLib::EvalPolynomial(aParameters[1], aMinMaxDeriv[1], aMinMaxDegree[1], |
364 | aCacheCols, aPolesArray[0], aTransientCoeffs[0]); |
365 | |
366 | // Calculate a point on surface and a derivatives along variable with minimal degree |
367 | PLib::EvalPolynomial(aParameters[0], aMinMaxDeriv[0], aMinMaxDegree[0], |
368 | aDimension, aTransientCoeffs[0], aPntDeriv[0]); |
369 | |
370 | // Calculate derivative along variable with maximal degree and mixed derivative |
371 | PLib::EvalPolynomial(aParameters[0], 1, aMinMaxDegree[0], aDimension, |
372 | aTransientCoeffs[aCacheCols], aPntDeriv[3 * aDimension]); |
373 | |
374 | // Calculate second derivative along variable with maximal degree |
375 | PLib::NoDerivativeEvalPolynomial(aParameters[0], aMinMaxDegree[0], aDimension, |
376 | aMinMaxDegree[0] * aDimension, aTransientCoeffs[aCacheCols<<1], |
377 | aPntDeriv[6 * aDimension]); |
378 | |
379 | Standard_Real* aResult = aPntDeriv; |
380 | Standard_Real aTempStorage[36]; |
381 | if (myIsRational) // calculate derivatives divided by weight's derivatives |
382 | { |
383 | BSplSLib::RationalDerivative(2, 2, 2, 2, aPntDeriv[0], aTempStorage[0]); |
384 | aResult = aTempStorage; |
385 | aDimension--; |
386 | } |
387 | |
388 | thePoint.SetCoord(aResult[0], aResult[1], aResult[2]); |
389 | if (myDegree[0] > myDegree[1]) |
390 | { |
391 | theTangentV.SetCoord(aResult[aDimension], aResult[aDimension + 1], aResult[aDimension + 2]); |
392 | Standard_Integer aShift = aDimension<<1; |
393 | theCurvatureV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]); |
394 | aShift += aDimension; |
395 | theTangentU.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]); |
396 | aShift += aDimension; |
397 | theCurvatureUV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]); |
398 | aShift += (aDimension << 1); |
399 | theCurvatureU.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]); |
400 | } |
401 | else |
402 | { |
403 | theTangentU.SetCoord(aResult[aDimension], aResult[aDimension + 1], aResult[aDimension + 2]); |
404 | Standard_Integer aShift = aDimension<<1; |
405 | theCurvatureU.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]); |
406 | aShift += aDimension; |
407 | theTangentV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]); |
408 | aShift += aDimension; |
409 | theCurvatureUV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]); |
410 | aShift += (aDimension << 1); |
411 | theCurvatureV.SetCoord(aResult[aShift], aResult[aShift + 1], aResult[aShift + 2]); |
412 | } |
413 | theTangentU.Multiply(anInvU); |
414 | theTangentV.Multiply(anInvV); |
415 | theCurvatureU.Multiply(anInvU * anInvU); |
416 | theCurvatureV.Multiply(anInvV * anInvV); |
417 | theCurvatureUV.Multiply(anInvU * anInvV); |
418 | } |
419 | |