0031381: Foundation Classes -wrong evaluations for rational BSpline curves using...
[occt.git] / src / BSplCLib / BSplCLib_Cache.cxx
CommitLineData
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 <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
94f71cad 24
92efcf78 25IMPLEMENT_STANDARD_RTTIEXT(BSplCLib_Cache,Standard_Transient)
26
94f71cad 27//! Converts handle of array of Standard_Real into the pointer to Standard_Real
35c0599a 28static 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 34BSplCLib_Cache::BSplCLib_Cache(const Standard_Integer& theDegree,
35 const Standard_Boolean& thePeriodic,
36 const TColStd_Array1OfReal& theFlatKnots,
0a96e0bb 37 const TColgp_Array1OfPnt2d& /* only used to distinguish from 3d variant */,
0e14656b 38 const TColStd_Array1OfReal* theWeights)
0a96e0bb 39: myIsRational(theWeights != NULL),
40 myParams (theDegree, thePeriodic, theFlatKnots)
94f71cad 41{
0a96e0bb 42 Standard_Integer aPWColNumber = (myIsRational ? 3 : 2);
43 myPolesWeights = new TColStd_HArray2OfReal (1, theDegree + 1, 1, aPWColNumber);
94f71cad 44}
45
46BSplCLib_Cache::BSplCLib_Cache(const Standard_Integer& theDegree,
47 const Standard_Boolean& thePeriodic,
48 const TColStd_Array1OfReal& theFlatKnots,
0a96e0bb 49 const TColgp_Array1OfPnt& /* only used to distinguish from 2d variant */,
0e14656b 50 const TColStd_Array1OfReal* theWeights)
0a96e0bb 51: myIsRational(theWeights != NULL),
52 myParams (theDegree, thePeriodic, theFlatKnots)
94f71cad 53{
0a96e0bb 54 Standard_Integer aPWColNumber = (myIsRational ? 4 : 3);
55 myPolesWeights = new TColStd_HArray2OfReal (1, theDegree + 1, 1, aPWColNumber);
94f71cad 56}
57
94f71cad 58Standard_Boolean BSplCLib_Cache::IsCacheValid(Standard_Real theParameter) const
59{
0a96e0bb 60 return myParams.IsCacheValid (theParameter);
94f71cad 61}
62
94f71cad 63void BSplCLib_Cache::BuildCache(const Standard_Real& theParameter,
94f71cad 64 const TColStd_Array1OfReal& theFlatKnots,
65 const TColgp_Array1OfPnt2d& thePoles2d,
0e14656b 66 const TColStd_Array1OfReal* theWeights)
94f71cad 67{
68 // Normalize theParameter for periodical B-splines
0a96e0bb 69 Standard_Real aNewParam = myParams.PeriodicNormalization (theParameter);
70 myParams.LocateParameter (aNewParam, theFlatKnots);
94f71cad 71
72 // Calculate new cache data
0a96e0bb 73 BSplCLib::BuildCache (myParams.SpanStart, myParams.SpanLength, myParams.IsPeriodic,
74 myParams.Degree, myParams.SpanIndex, theFlatKnots, thePoles2d,
75 theWeights, myPolesWeights->ChangeArray2());
94f71cad 76}
77
78void BSplCLib_Cache::BuildCache(const Standard_Real& theParameter,
94f71cad 79 const TColStd_Array1OfReal& theFlatKnots,
80 const TColgp_Array1OfPnt& thePoles,
0e14656b 81 const TColStd_Array1OfReal* theWeights)
94f71cad 82{
83 // Create list of knots with repetitions and normalize theParameter for periodical B-splines
0a96e0bb 84 Standard_Real aNewParam = myParams.PeriodicNormalization (theParameter);
85 myParams.LocateParameter (aNewParam, theFlatKnots);
94f71cad 86
87 // Calculate new cache data
0a96e0bb 88 BSplCLib::BuildCache (myParams.SpanStart, myParams.SpanLength, myParams.IsPeriodic,
89 myParams.Degree, myParams.SpanIndex, theFlatKnots, thePoles,
90 theWeights, myPolesWeights->ChangeArray2());
94f71cad 91}
92
94f71cad 93void BSplCLib_Cache::CalculateDerivative(const Standard_Real& theParameter,
94 const Standard_Integer& theDerivative,
95 Standard_Real& theDerivArray) const
96{
0a96e0bb 97 Standard_Real aNewParameter = myParams.PeriodicNormalization (theParameter);
98 aNewParameter = (aNewParameter - myParams.SpanStart) / myParams.SpanLength;
94f71cad 99
100 Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
101 Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
102
103 // Temporary container. The maximal size of this container is defined by:
104 // 1) maximal derivative for cache evaluation, which is 3, plus one row for function values,
105 // 2) and maximal dimension of the point, which is 3, plus one column for weights.
106 Standard_Real aTmpContainer[16];
107
108 // When the PLib::RationaDerivative needs to be called, use temporary container
109 Standard_Real* aPntDeriv = myIsRational ? aTmpContainer : &theDerivArray;
110
111 // When the degree of curve is lesser than the requested derivative,
112 // nullify array cells corresponding to greater derivatives
113 Standard_Integer aDerivative = theDerivative;
f732ea1a 114 if (!myIsRational && myParams.Degree < theDerivative)
94f71cad 115 {
0a96e0bb 116 aDerivative = myParams.Degree;
117 for (Standard_Integer ind = myParams.Degree * aDimension; ind < (theDerivative + 1) * aDimension; ind++)
94f71cad 118 {
119 aPntDeriv[ind] = 0.0;
120 (&theDerivArray)[ind] = 0.0; // should be cleared separately, because aPntDeriv may look to another memory area
121 }
122 }
123
0a96e0bb 124 PLib::EvalPolynomial(aNewParameter, aDerivative, myParams.Degree, aDimension,
94f71cad 125 aPolesArray[0], aPntDeriv[0]);
126 // Unnormalize derivatives since those are computed normalized
127 Standard_Real aFactor = 1.0;
128 for (Standard_Integer deriv = 1; deriv <= aDerivative; deriv++)
129 {
0a96e0bb 130 aFactor /= myParams.SpanLength;
94f71cad 131 for (Standard_Integer ind = 0; ind < aDimension; ind++)
132 aPntDeriv[aDimension * deriv + ind] *= aFactor;
133 }
134
135 if (myIsRational) // calculate derivatives divided by weights derivatives
136 PLib::RationalDerivative(aDerivative, aDerivative, aDimension - 1, aPntDeriv[0], theDerivArray);
137}
138
139
140void BSplCLib_Cache::D0(const Standard_Real& theParameter, gp_Pnt2d& thePoint) const
141{
0a96e0bb 142 Standard_Real aNewParameter = myParams.PeriodicNormalization (theParameter);
143 aNewParameter = (aNewParameter - myParams.SpanStart) / myParams.SpanLength;
94f71cad 144
145 Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
146 Standard_Real aPoint[4];
147 Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
148
0a96e0bb 149 PLib::NoDerivativeEvalPolynomial(aNewParameter, myParams.Degree,
150 aDimension, myParams.Degree * aDimension,
94f71cad 151 aPolesArray[0], aPoint[0]);
152
153 thePoint.SetCoord(aPoint[0], aPoint[1]);
154 if (myIsRational)
155 thePoint.ChangeCoord().Divide(aPoint[2]);
156}
157
158void BSplCLib_Cache::D0(const Standard_Real& theParameter, gp_Pnt& thePoint) const
159{
0a96e0bb 160 Standard_Real aNewParameter = myParams.PeriodicNormalization (theParameter);
161 aNewParameter = (aNewParameter - myParams.SpanStart) / myParams.SpanLength;
94f71cad 162
163 Standard_Real* aPolesArray = ConvertArray(myPolesWeights);
164 Standard_Real aPoint[4];
165 Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
166
0a96e0bb 167 PLib::NoDerivativeEvalPolynomial(aNewParameter, myParams.Degree,
168 aDimension, myParams.Degree * aDimension,
94f71cad 169 aPolesArray[0], aPoint[0]);
170
171 thePoint.SetCoord(aPoint[0], aPoint[1], aPoint[2]);
172 if (myIsRational)
173 thePoint.ChangeCoord().Divide(aPoint[3]);
174}
175
176
177void BSplCLib_Cache::D1(const Standard_Real& theParameter, gp_Pnt2d& thePoint, gp_Vec2d& theTangent) const
178{
179 Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
180 Standard_Real aPntDeriv[8]; // result storage (point and derivative coordinates)
181
182 this->CalculateDerivative(theParameter, 1, aPntDeriv[0]);
183 if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
184 aDimension -= 1;
185
186 thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1]);
187 theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1]);
188}
189
190void BSplCLib_Cache::D1(const Standard_Real& theParameter, gp_Pnt& thePoint, gp_Vec& theTangent) const
191{
192 Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
193 Standard_Real aPntDeriv[8]; // result storage (point and derivative coordinates)
194
195 this->CalculateDerivative(theParameter, 1, aPntDeriv[0]);
196 if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
197 aDimension -= 1;
198
199 thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1], aPntDeriv[2]);
200 theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1], aPntDeriv[aDimension + 2]);
201}
202
203void BSplCLib_Cache::D2(const Standard_Real& theParameter, gp_Pnt2d& thePoint, gp_Vec2d& theTangent, gp_Vec2d& theCurvature) const
204{
205 Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
206 Standard_Real aPntDeriv[12]; // result storage (point and derivatives coordinates)
207
208 this->CalculateDerivative(theParameter, 2, aPntDeriv[0]);
209 if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
210 aDimension -= 1;
211
212 thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1]);
213 theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1]);
214 theCurvature.SetCoord(aPntDeriv[aDimension<<1], aPntDeriv[(aDimension<<1) + 1]);
215}
216
217void BSplCLib_Cache::D2(const Standard_Real& theParameter, gp_Pnt& thePoint, gp_Vec& theTangent, gp_Vec& theCurvature) const
218{
219 Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
220 Standard_Real aPntDeriv[12]; // result storage (point and derivatives coordinates)
221
222 this->CalculateDerivative(theParameter, 2, aPntDeriv[0]);
223 if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
224 aDimension -= 1;
225
226 thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1], aPntDeriv[2]);
227 theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1], aPntDeriv[aDimension + 2]);
228 theCurvature.SetCoord(aPntDeriv[aDimension<<1], aPntDeriv[(aDimension<<1) + 1], aPntDeriv[(aDimension<<1) + 2]);
229}
230
231
232void BSplCLib_Cache::D3(const Standard_Real& theParameter,
233 gp_Pnt2d& thePoint,
234 gp_Vec2d& theTangent,
235 gp_Vec2d& theCurvature,
236 gp_Vec2d& theTorsion) const
237{
238 Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
239 Standard_Real aPntDeriv[16]; // result storage (point and derivatives coordinates)
240
241 this->CalculateDerivative(theParameter, 3, aPntDeriv[0]);
242 if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
243 aDimension -= 1;
244
245 thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1]);
246 theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1]);
247 Standard_Integer aShift = aDimension << 1;
248 theCurvature.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1]);
249 aShift += aDimension;
250 theTorsion.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1]);
251}
252
253void BSplCLib_Cache::D3(const Standard_Real& theParameter,
254 gp_Pnt& thePoint,
255 gp_Vec& theTangent,
256 gp_Vec& theCurvature,
257 gp_Vec& theTorsion) const
258{
259 Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns
260 Standard_Real aPntDeriv[16]; // result storage (point and derivatives coordinates)
261
262 this->CalculateDerivative(theParameter, 3, aPntDeriv[0]);
263 if (myIsRational) // the size of aPntDeriv was changed by PLib::RationalDerivative
264 aDimension -= 1;
265
266 thePoint.SetCoord(aPntDeriv[0], aPntDeriv[1], aPntDeriv[2]);
267 theTangent.SetCoord(aPntDeriv[aDimension], aPntDeriv[aDimension + 1], aPntDeriv[aDimension + 2]);
268 Standard_Integer aShift = aDimension << 1;
269 theCurvature.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1], aPntDeriv[aShift + 2]);
270 aShift += aDimension;
271 theTorsion.SetCoord(aPntDeriv[aShift], aPntDeriv[aShift + 1], aPntDeriv[aShift + 2]);
272}
273