0030582: Coding - avoid defining references to properties of NULL objects
[occt.git] / src / BSplCLib / BSplCLib_CurveComputation.gxx
CommitLineData
b311480e 1// Copyright (c) 1995-1999 Matra Datavision
973c2be1 2// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 3//
973c2be1 4// This file is part of Open CASCADE Technology software library.
b311480e 5//
d5f74e42 6// This library is free software; you can redistribute it and/or modify it under
7// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1 8// by the Free Software Foundation, with special exception defined in the file
9// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10// distribution for complete text of the license and disclaimer of any warranty.
b311480e 11//
973c2be1 12// Alternatively, this file may be used under the terms of Open CASCADE
13// commercial license or contractual agreement.
b311480e 14
7fd59977 15// xab : modified 15-Mar-94 added EvalBSplineMatrix, FactorBandedMatrix, SolveBandedSystem
16// Eval, BuildCache, cacheD0, cacheD1, cacheD2, cacheD3, LocalMatrix,
17// EvalPolynomial, RationalDerivatives.
18// xab : 22-Mar-94 : fixed rational problem in BuildCache
19// xab : 12-Mar-96 : added MovePointAndTangent
7fd59977 20#include <TColStd_Array1OfInteger.hxx>
21#include <TColStd_Array1OfReal.hxx>
94f71cad 22#include <TColStd_Array2OfReal.hxx>
7fd59977 23#include <gp_Vec2d.hxx>
24#include <Standard_ConstructionError.hxx>
25#include <PLib.hxx>
26#include <math_Matrix.hxx>
27
7fd59977 28//=======================================================================
29//struct : BSplCLib_DataContainer
30//purpose: Auxiliary structure providing buffers for poles and knots used in
31// evaluation of bspline (allocated in the stack)
32//=======================================================================
33
34struct BSplCLib_DataContainer
35{
36 BSplCLib_DataContainer(Standard_Integer Degree)
37 {
105aae76 38 (void)Degree; // avoid compiler warning
41194117
K
39 Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree() ||
40 BSplCLib::MaxDegree() > 25,
41 "BSplCLib: bspline degree is greater than maximum supported");
7fd59977 42 }
43
44 Standard_Real poles[(25+1)*(Dimension_gen+1)];
45 Standard_Real knots[2*25];
46 Standard_Real ders[Dimension_gen*4];
47};
48
49//=======================================================================
50//function : Reverse
51//purpose :
52//=======================================================================
53
54void BSplCLib::Reverse(Array1OfPoints& Poles,
55 const Standard_Integer L)
56{
57 Standard_Integer i, l = L;
58 l = Poles.Lower() + (l-Poles.Lower()) % (Poles.Upper()-Poles.Lower()+1);
59
60 Array1OfPoints temp(0,Poles.Length()-1);
61
62 for (i = Poles.Lower(); i <= l; i++)
63 temp(l-i) = Poles(i);
64
65 for (i = l+1; i <= Poles.Upper(); i++)
66 temp(l-Poles.Lower()+Poles.Upper()-i+1) = Poles(i);
67
68 for (i = Poles.Lower(); i <= Poles.Upper(); i++)
69 Poles(i) = temp(i-Poles.Lower());
70}
71
72//
73// CURVES MODIFICATIONS
74//
75
76//=======================================================================
77//function : RemoveKnot
78//purpose :
79//=======================================================================
80
81Standard_Boolean BSplCLib::RemoveKnot
82(const Standard_Integer Index,
83 const Standard_Integer Mult,
84 const Standard_Integer Degree,
85 const Standard_Boolean Periodic,
86 const Array1OfPoints& Poles,
0e14656b 87 const TColStd_Array1OfReal* Weights,
7fd59977 88 const TColStd_Array1OfReal& Knots,
89 const TColStd_Array1OfInteger& Mults,
90 Array1OfPoints& NewPoles,
0e14656b 91 TColStd_Array1OfReal* NewWeights,
7fd59977 92 TColStd_Array1OfReal& NewKnots,
93 TColStd_Array1OfInteger& NewMults,
94 const Standard_Real Tolerance)
95{
0e14656b 96 Standard_Boolean rational = Weights != NULL;
7fd59977 97 Standard_Integer dim;
98 dim = Dimension_gen;
99 if (rational) dim++;
100
101 TColStd_Array1OfReal poles(1,dim*Poles.Length());
102 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
103
0e14656b 104 if (rational) PLib::SetPoles(Poles,*Weights,poles);
7fd59977 105 else PLib::SetPoles(Poles,poles);
106
107 if (!RemoveKnot(Index,Mult,Degree,Periodic,dim,
108 poles,Knots,Mults,newpoles,NewKnots,NewMults,Tolerance))
109 return Standard_False;
110
0e14656b 111 if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
7fd59977 112 else PLib::GetPoles(newpoles,NewPoles);
113 return Standard_True;
114}
115
116//=======================================================================
117//function : InsertKnots
118//purpose : insert an array of knots and multiplicities
119//=======================================================================
120
121void BSplCLib::InsertKnots
122(const Standard_Integer Degree,
123 const Standard_Boolean Periodic,
124 const Array1OfPoints& Poles,
0e14656b 125 const TColStd_Array1OfReal* Weights,
7fd59977 126 const TColStd_Array1OfReal& Knots,
127 const TColStd_Array1OfInteger& Mults,
128 const TColStd_Array1OfReal& AddKnots,
0e14656b 129 const TColStd_Array1OfInteger* AddMults,
7fd59977 130 Array1OfPoints& NewPoles,
0e14656b 131 TColStd_Array1OfReal* NewWeights,
7fd59977 132 TColStd_Array1OfReal& NewKnots,
133 TColStd_Array1OfInteger& NewMults,
134 const Standard_Real Epsilon,
135 const Standard_Boolean Add)
136{
0e14656b 137 Standard_Boolean rational = Weights != NULL;
7fd59977 138 Standard_Integer dim;
139 dim = Dimension_gen;
140 if (rational) dim++;
141
142 TColStd_Array1OfReal poles(1,dim*Poles.Length());
143 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
144
0e14656b 145 if (rational) PLib::SetPoles(Poles,*Weights,poles);
7fd59977 146 else PLib::SetPoles(Poles,poles);
147
148 InsertKnots(Degree,Periodic,dim,poles,Knots,Mults,
149 AddKnots,AddMults,newpoles,NewKnots,NewMults,Epsilon,Add);
150
0e14656b 151 if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
7fd59977 152 else PLib::GetPoles(newpoles,NewPoles);
153}
154
155//=======================================================================
156//function : InsertKnot
157//purpose :
158//=======================================================================
159
160void BSplCLib::InsertKnot(const Standard_Integer ,
161 const Standard_Real U,
162 const Standard_Integer UMult,
163 const Standard_Integer Degree,
164 const Standard_Boolean Periodic,
165 const Array1OfPoints& Poles,
0e14656b 166 const TColStd_Array1OfReal* Weights,
7fd59977 167 const TColStd_Array1OfReal& Knots,
168 const TColStd_Array1OfInteger& Mults,
169 Array1OfPoints& NewPoles,
0e14656b 170 TColStd_Array1OfReal* NewWeights)
7fd59977 171{
172 TColStd_Array1OfReal k(1,1);
173 k(1) = U;
174 TColStd_Array1OfInteger m(1,1);
175 m(1) = UMult;
176 TColStd_Array1OfReal nk(1,Knots.Length()+1);
177 TColStd_Array1OfInteger nm(1,Knots.Length()+1);
178 InsertKnots(Degree,Periodic,Poles,Weights,Knots,Mults,
0e14656b 179 k,&m,NewPoles,NewWeights,nk,nm,Epsilon(U));
7fd59977 180}
181
182//=======================================================================
183//function : RaiseMultiplicity
184//purpose :
185//=======================================================================
186
187void BSplCLib::RaiseMultiplicity(const Standard_Integer KnotIndex,
188 const Standard_Integer Mult,
189 const Standard_Integer Degree,
190 const Standard_Boolean Periodic,
191 const Array1OfPoints& Poles,
0e14656b 192 const TColStd_Array1OfReal* Weights,
7fd59977 193 const TColStd_Array1OfReal& Knots,
194 const TColStd_Array1OfInteger& Mults,
195 Array1OfPoints& NewPoles,
0e14656b 196 TColStd_Array1OfReal* NewWeights)
7fd59977 197{
198 TColStd_Array1OfReal k(1,1);
199 k(1) = Knots(KnotIndex);
200 TColStd_Array1OfInteger m(1,1);
201 m(1) = Mult - Mults(KnotIndex);
202 TColStd_Array1OfReal nk(1,Knots.Length());
203 TColStd_Array1OfInteger nm(1,Knots.Length());
204 InsertKnots(Degree,Periodic,Poles,Weights,Knots,Mults,
0e14656b 205 k,&m,NewPoles,NewWeights,nk,nm,Epsilon(k(1)));
7fd59977 206}
207
208//=======================================================================
209//function : IncreaseDegree
210//purpose :
211//=======================================================================
212
213void BSplCLib::IncreaseDegree
214(const Standard_Integer Degree,
215 const Standard_Integer NewDegree,
216 const Standard_Boolean Periodic,
217 const Array1OfPoints& Poles,
0e14656b 218 const TColStd_Array1OfReal* Weights,
7fd59977 219 const TColStd_Array1OfReal& Knots,
220 const TColStd_Array1OfInteger& Mults,
221 Array1OfPoints& NewPoles,
0e14656b 222 TColStd_Array1OfReal* NewWeights,
7fd59977 223 TColStd_Array1OfReal& NewKnots,
224 TColStd_Array1OfInteger& NewMults)
225{
0e14656b 226 Standard_Boolean rational = Weights != NULL;
7fd59977 227 Standard_Integer dim;
228 dim = Dimension_gen;
229 if (rational) dim++;
230
231 TColStd_Array1OfReal poles(1,dim*Poles.Length());
232 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
233
0e14656b 234 if (rational) PLib::SetPoles(Poles,*Weights,poles);
7fd59977 235 else PLib::SetPoles(Poles,poles);
236
237 IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults,
238 newpoles,NewKnots,NewMults);
239
0e14656b 240 if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
7fd59977 241 else PLib::GetPoles(newpoles,NewPoles);
242}
243
244//=======================================================================
245//function : Unperiodize
246//purpose :
247//=======================================================================
248
249void BSplCLib::Unperiodize
250(const Standard_Integer Degree,
251 const TColStd_Array1OfInteger& Mults,
252 const TColStd_Array1OfReal& Knots,
253 const Array1OfPoints& Poles,
0e14656b 254 const TColStd_Array1OfReal* Weights,
7fd59977 255 TColStd_Array1OfInteger& NewMults,
256 TColStd_Array1OfReal& NewKnots,
257 Array1OfPoints& NewPoles,
0e14656b 258 TColStd_Array1OfReal* NewWeights)
7fd59977 259{
0e14656b 260 Standard_Boolean rational = Weights != NULL;
7fd59977 261 Standard_Integer dim;
262 dim = Dimension_gen;
263 if (rational) dim++;
264
265 TColStd_Array1OfReal poles(1,dim*Poles.Length());
266 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
267
0e14656b 268 if (rational) PLib::SetPoles(Poles,*Weights,poles);
7fd59977 269 else PLib::SetPoles(Poles,poles);
270
271 Unperiodize(Degree, dim, Mults, Knots, poles,
272 NewMults,NewKnots, newpoles);
273
0e14656b 274 if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
7fd59977 275 else PLib::GetPoles(newpoles,NewPoles);
276}
277
278//=======================================================================
279//function : Trimming
280//purpose :
281//=======================================================================
282
283void BSplCLib::Trimming(const Standard_Integer Degree,
284 const Standard_Boolean Periodic,
285 const TColStd_Array1OfReal& Knots,
286 const TColStd_Array1OfInteger& Mults,
287 const Array1OfPoints& Poles,
0e14656b 288 const TColStd_Array1OfReal* Weights,
7fd59977 289 const Standard_Real U1,
290 const Standard_Real U2,
291 TColStd_Array1OfReal& NewKnots,
292 TColStd_Array1OfInteger& NewMults,
293 Array1OfPoints& NewPoles,
0e14656b 294 TColStd_Array1OfReal* NewWeights)
7fd59977 295{
0e14656b 296 Standard_Boolean rational = Weights != NULL;
7fd59977 297 Standard_Integer dim;
298 dim = Dimension_gen;
299 if (rational) dim++;
300
301 TColStd_Array1OfReal poles(1,dim*Poles.Length());
302 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
303
0e14656b 304 if (rational) PLib::SetPoles(Poles,*Weights,poles);
7fd59977 305 else PLib::SetPoles(Poles,poles);
306
307 Trimming(Degree, Periodic, dim, Knots, Mults, poles, U1, U2,
308 NewKnots, NewMults, newpoles);
309
0e14656b 310 if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
7fd59977 311 else PLib::GetPoles(newpoles,NewPoles);
312}
313
314//--------------------------------------------------------------------------
315//ELEMENTARY COMPUTATIONS
316//--------------------------------------------------------------------------
317//
318// All the following methods are methods of low level used in BSplCLib
319// but also in BSplSLib (but not in Geom)
320// At the creation time of this package there is no possibility
321// to declare private methods of package and to declare friend
322// methods of package. It could interesting to declare that BSplSLib
323// is a package friend to BSplCLib because it uses all the basis methods
324// of BSplCLib.
325//--------------------------------------------------------------------------
326
327//=======================================================================
328//function : BuildEval
329//purpose : builds the local array for evaluation
330//=======================================================================
331
332void BSplCLib::BuildEval(const Standard_Integer Degree,
333 const Standard_Integer Index,
334 const Array1OfPoints& Poles,
0e14656b 335 const TColStd_Array1OfReal* Weights,
7fd59977 336 Standard_Real& LP)
337{
338 Standard_Real w, *pole = &LP;
339 Standard_Integer PLower = Poles.Lower();
340 Standard_Integer PUpper = Poles.Upper();
341 Standard_Integer i;
342 Standard_Integer ip = PLower + Index - 1;
0e14656b 343 if (Weights == NULL) {
7fd59977 344 for (i = 0; i <= Degree; i++) {
345 ip++;
346 if (ip > PUpper) ip = PLower;
347 const Point& P = Poles(ip);
348 PointToCoords (pole,P,+0);
349 pole += Dimension_gen;
350 }
351 }
352 else {
353 for (i = 0; i <= Degree; i++) {
354 ip++;
355 if (ip > PUpper) ip = PLower;
356 const Point& P = Poles(ip);
0e14656b 357 pole[Dimension_gen] = w = (*Weights)(ip);
7fd59977 358 PointToCoords (pole, P, * w);
359 pole += Dimension_gen + 1;
360 }
361 }
362}
363
364//=======================================================================
365//function : PrepareEval
366//purpose : stores data for Eval in the local arrays
367// dc.poles and dc.knots
368//=======================================================================
369
370static void PrepareEval
371(Standard_Real& u,
372 Standard_Integer& index,
373 Standard_Integer& dim,
374 Standard_Boolean& rational,
375 const Standard_Integer Degree,
376 const Standard_Boolean Periodic,
377 const Array1OfPoints& Poles,
0e14656b 378 const TColStd_Array1OfReal* Weights,
7fd59977 379 const TColStd_Array1OfReal& Knots,
0e14656b 380 const TColStd_Array1OfInteger* Mults,
7fd59977 381 BSplCLib_DataContainer& dc)
382{
383 // Set the Index
384 BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
385
386 // make the knots
387 BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*dc.knots);
0e14656b 388 if (Mults == NULL)
7fd59977 389 index -= Knots.Lower() + Degree;
390 else
0e14656b 391 index = BSplCLib::PoleIndex(Degree,index,Periodic,*Mults);
7fd59977 392
393 // check truly rational
0e14656b 394 rational = (Weights != NULL);
7fd59977 395 if (rational) {
0e14656b 396 Standard_Integer WLower = Weights->Lower() + index;
397 rational = BSplCLib::IsRational(*Weights, WLower, WLower + Degree);
7fd59977 398 }
399
400 // make the poles
401 if (rational) {
402 dim = Dimension_gen + 1;
403 BSplCLib::BuildEval(Degree, index, Poles, Weights , *dc.poles);
404 }
405 else {
406 dim = Dimension_gen;
407 BSplCLib::BuildEval(Degree, index, Poles, BSplCLib::NoWeights(), *dc.poles);
408 }
409}
410
411//=======================================================================
412//function : D0
413//purpose :
414//=======================================================================
415
416void BSplCLib::D0
417(const Standard_Real U,
418 const Standard_Integer Index,
419 const Standard_Integer Degree,
420 const Standard_Boolean Periodic,
421 const Array1OfPoints& Poles,
0e14656b 422 const TColStd_Array1OfReal* Weights,
7fd59977 423 const TColStd_Array1OfReal& Knots,
0e14656b 424 const TColStd_Array1OfInteger* Mults,
7fd59977 425 Point& P)
426{
427// Standard_Integer k,dim,index = Index;
428 Standard_Integer dim,index = Index;
429 Standard_Real u = U;
430 Standard_Boolean rational;
431 BSplCLib_DataContainer dc(Degree);
432 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
433 BSplCLib::Eval(u,Degree,*dc.knots,dim,*dc.poles);
434
435 if (rational) {
436 Standard_Real w = dc.poles[Dimension_gen];
437 CoordsToPoint (P, dc.poles, / w);
438 }
439 else
440 CoordsToPoint (P, dc.poles, );
441}
442
443//=======================================================================
444//function : D1
445//purpose :
446//=======================================================================
447
448void BSplCLib::D1
449(const Standard_Real U,
450 const Standard_Integer Index,
451 const Standard_Integer Degree,
452 const Standard_Boolean Periodic,
453 const Array1OfPoints& Poles,
0e14656b 454 const TColStd_Array1OfReal* Weights,
7fd59977 455 const TColStd_Array1OfReal& Knots,
0e14656b 456 const TColStd_Array1OfInteger* Mults,
7fd59977 457 Point& P,
458 Vector& V)
459{
460 Standard_Integer dim,index = Index;
461 Standard_Real u = U;
462 Standard_Boolean rational;
463 BSplCLib_DataContainer dc(Degree);
464 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
465 BSplCLib::Bohm(u,Degree,1,*dc.knots,dim,*dc.poles);
466 Standard_Real *result = dc.poles;
467 if (rational) {
468 PLib::RationalDerivative(Degree,1,Dimension_gen,*dc.poles,*dc.ders);
469 result = dc.ders;
470 }
471
472 CoordsToPoint (P, result, );
473 CoordsToPoint (V, result + Dimension_gen, );
474}
475
476//=======================================================================
477//function : D2
478//purpose :
479//=======================================================================
480
481void BSplCLib::D2
482(const Standard_Real U,
483 const Standard_Integer Index,
484 const Standard_Integer Degree,
485 const Standard_Boolean Periodic,
486 const Array1OfPoints& Poles,
0e14656b 487 const TColStd_Array1OfReal* Weights,
7fd59977 488 const TColStd_Array1OfReal& Knots,
0e14656b 489 const TColStd_Array1OfInteger* Mults,
7fd59977 490 Point& P,
491 Vector& V1,
492 Vector& V2)
493{
494 Standard_Integer dim,index = Index;
495 Standard_Real u = U;
496 Standard_Boolean rational;
497 BSplCLib_DataContainer dc(Degree);
498 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
499 BSplCLib::Bohm(u,Degree,2,*dc.knots,dim,*dc.poles);
500 Standard_Real *result = dc.poles;
501 if (rational) {
502 PLib::RationalDerivative(Degree,2,Dimension_gen,*dc.poles,*dc.ders);
503 result = dc.ders;
504 }
505
506 CoordsToPoint (P, result, );
507 CoordsToPoint (V1, result + Dimension_gen, );
508 if (!rational && (Degree < 2))
509 NullifyPoint (V2);
510 else
511 CoordsToPoint (V2, result + 2 * Dimension_gen, );
512}
513
514//=======================================================================
515//function : D3
516//purpose :
517//=======================================================================
518
519void BSplCLib::D3
520(const Standard_Real U,
521 const Standard_Integer Index,
522 const Standard_Integer Degree,
523 const Standard_Boolean Periodic,
524 const Array1OfPoints& Poles,
0e14656b 525 const TColStd_Array1OfReal* Weights,
7fd59977 526 const TColStd_Array1OfReal& Knots,
0e14656b 527 const TColStd_Array1OfInteger* Mults,
7fd59977 528 Point& P,
529 Vector& V1,
530 Vector& V2,
531 Vector& V3)
532{
533 Standard_Integer dim,index = Index;
534 Standard_Real u = U;
535 Standard_Boolean rational;
536 BSplCLib_DataContainer dc(Degree);
537 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
538 BSplCLib::Bohm(u,Degree,3,*dc.knots,dim,*dc.poles);
539 Standard_Real *result = dc.poles;
540 if (rational) {
541 PLib::RationalDerivative(Degree,3,Dimension_gen,*dc.poles,*dc.ders);
542 result = dc.ders;
543 }
544
545 CoordsToPoint (P, result, );
546 CoordsToPoint (V1, result + Dimension_gen, );
547 if (!rational && (Degree < 2))
548 NullifyPoint (V2);
549 else
550 CoordsToPoint (V2, result + 2 * Dimension_gen, );
551 if (!rational && (Degree < 3))
552 NullifyPoint (V3);
553 else
554 CoordsToPoint (V3, result + 3 * Dimension_gen, );
555}
556
557//=======================================================================
558//function : DN
559//purpose :
560//=======================================================================
561
562void BSplCLib::DN
563(const Standard_Real U,
564 const Standard_Integer N,
565 const Standard_Integer Index,
566 const Standard_Integer Degree,
567 const Standard_Boolean Periodic,
568 const Array1OfPoints& Poles,
0e14656b 569 const TColStd_Array1OfReal* Weights,
7fd59977 570 const TColStd_Array1OfReal& Knots,
0e14656b 571 const TColStd_Array1OfInteger* Mults,
7fd59977 572 Vector& VN)
573{
574 Standard_Integer dim,index = Index;
575 Standard_Real u = U;
576 Standard_Boolean rational;
577 BSplCLib_DataContainer dc(Degree);
578 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
579 BSplCLib::Bohm(u,Degree,N,*dc.knots,dim,*dc.poles);
580
581 if (rational) {
582 Standard_Real v[Dimension_gen];
583 PLib::RationalDerivative(Degree,N,Dimension_gen,*dc.poles,v[0],Standard_False);
584 CoordsToPoint (VN, v, );
585 }
586 else {
587 if (N > Degree)
588 NullifyPoint (VN);
589 else {
590 Standard_Real *DN = dc.poles + N * Dimension_gen;
591 CoordsToPoint (VN, DN, );
592 }
593 }
594}
595
596//=======================================================================
597//function : Solves a LU factored Matrix
598//purpose :
599//=======================================================================
600
601Standard_Integer
602BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
603 const Standard_Integer UpperBandWidth,
604 const Standard_Integer LowerBandWidth,
605 Array1OfPoints& PolesArray)
606{
607 Standard_Real *PArray ;
608 PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
609
610 return BSplCLib::SolveBandedSystem(Matrix,
611 UpperBandWidth,
612 LowerBandWidth,
613 Dimension_gen,
614 PArray[0]) ;
615}
616
617//=======================================================================
618//function : Solves a LU factored Matrix
619//purpose : if HomogeneousFlag is 1 then the input and the output
620// will be homogeneous that is no division or multiplication
621// by weigths will happen. On the contrary if HomogeneousFlag
622// is 0 then the poles will be multiplied first by the weights
623// and after interpolation they will be devided by the weights
624//=======================================================================
625
626Standard_Integer
627BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
628 const Standard_Integer UpperBandWidth,
629 const Standard_Integer LowerBandWidth,
630 const Standard_Boolean HomogeneousFlag,
631 Array1OfPoints& PolesArray,
632 TColStd_Array1OfReal& WeightsArray)
633{
634 Standard_Real *PArray,
635 * WArray ;
636 PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
637 WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
638 return BSplCLib::SolveBandedSystem(Matrix,
639 UpperBandWidth,
640 LowerBandWidth,
641 HomogeneousFlag,
642 Dimension_gen,
643 PArray[0],
644 WArray[0]) ;
645}
646
647//=======================================================================
648//function : Evaluates a Bspline function : uses the ExtrapMode
649//purpose : the function is extrapolated using the Taylor expansion
650// of degree ExtrapMode[0] to the left and the Taylor
651// expansion of degree ExtrapMode[1] to the right
652// if the HomogeneousFlag == True than the Poles are supposed
653// to be stored homogeneously and the result will also be homogeneous
654// Valid only if Weights
655//=======================================================================
656void BSplCLib::Eval(const Standard_Real Parameter,
657 const Standard_Boolean PeriodicFlag,
658 const Standard_Boolean HomogeneousFlag,
659 Standard_Integer& ExtrapMode,
660 const Standard_Integer Degree,
661 const TColStd_Array1OfReal& FlatKnots,
662 const Array1OfPoints& PolesArray,
663 const TColStd_Array1OfReal& WeightsArray,
664 Point& aPoint,
665 Standard_Real& aWeight)
666{
667 Standard_Real Inverse,
668 P[Dimension_gen],
669 *PArray,
670 *WArray ;
671 Standard_Integer kk ;
672 PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
673 WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
674 if (HomogeneousFlag) {
675 BSplCLib::Eval(Parameter,
676 PeriodicFlag,
677 0,
678 ExtrapMode,
679 Degree,
680 FlatKnots,
681 Dimension_gen,
682 PArray[0],
683 P[0]) ;
684 BSplCLib::Eval(Parameter,
685 PeriodicFlag,
686 0,
687 ExtrapMode,
688 Degree,
689 FlatKnots,
690 1,
691 WArray[0],
692 aWeight) ;
693 }
694 else {
695 BSplCLib::Eval(Parameter,
696 PeriodicFlag,
697 0,
698 ExtrapMode,
699 Degree,
700 FlatKnots,
701 Dimension_gen,
702 PArray[0],
703 WArray[0],
704 P[0],
705 aWeight) ;
706 Inverse = 1.0e0 / aWeight ;
707
708 for (kk = 0 ; kk < Dimension_gen ; kk++) {
709 P[kk] *= Inverse ;
710 }
711 }
712
713 for (kk = 0 ; kk < Dimension_gen ; kk++)
714 aPoint.SetCoord(kk + 1, P[kk]) ;
715}
716
717//=======================================================================
718//function : CacheD0
719//purpose : Evaluates the polynomial cache of the Bspline Curve
720//
721//=======================================================================
722void BSplCLib::CacheD0(const Standard_Real Parameter,
723 const Standard_Integer Degree,
724 const Standard_Real CacheParameter,
725 const Standard_Real SpanLenght,
726 const Array1OfPoints& PolesArray,
0e14656b 727 const TColStd_Array1OfReal* WeightsArray,
7fd59977 728 Point& aPoint)
729{
730 //
731 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
732 // form
733 // the SpanLenght is the normalizing factor so that the polynomial is between
734 // 0 and 1
735 Standard_Real NewParameter, Inverse;
736 Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
737 Standard_Real * myPoint = (Standard_Real *) &aPoint ;
738 NewParameter = (Parameter - CacheParameter) / SpanLenght ;
739 PLib::NoDerivativeEvalPolynomial(NewParameter,
740 Degree,
741 Dimension_gen,
742 Degree * Dimension_gen,
743 PArray[0],
744 myPoint[0]) ;
0e14656b 745 if (WeightsArray != NULL) {
746 const TColStd_Array1OfReal& refWeights = *WeightsArray;
7fd59977 747 Standard_Real *
0e14656b 748 WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
7fd59977 749 PLib::NoDerivativeEvalPolynomial(NewParameter,
750 Degree,
751 1,
752 Degree,
753 WArray[0],
754 Inverse) ;
755
756 Inverse = 1.0e0 / Inverse;
757 ModifyCoords (myPoint, *= Inverse);
758 }
759}
760
761//=======================================================================
762//function : CacheD1
763//purpose : Evaluates the polynomial cache of the Bspline Curve
764//
765//=======================================================================
766void BSplCLib::CacheD1(const Standard_Real Parameter,
767 const Standard_Integer Degree,
768 const Standard_Real CacheParameter,
769 const Standard_Real SpanLenght,
770 const Array1OfPoints& PolesArray,
0e14656b 771 const TColStd_Array1OfReal* WeightsArray,
7fd59977 772 Point& aPoint,
773 Vector& aVector)
774{
775 //
776 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
777 // form
778 // the SpanLenght is the normalizing factor so that the polynomial is between
779 // 0 and 1
780 Standard_Real LocalPDerivatives[Dimension_gen << 1] ;
781// Standard_Real LocalWDerivatives[2], NewParameter, Inverse ;
782 Standard_Real LocalWDerivatives[2], NewParameter ;
783
784 Standard_Real *
785 PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
786 Standard_Real *
787 myPoint = (Standard_Real *) &aPoint ;
788 Standard_Real *
789 myVector = (Standard_Real *) &aVector ;
790 NewParameter = (Parameter - CacheParameter) / SpanLenght ;
791 PLib::EvalPolynomial(NewParameter,
792 1,
793 Degree,
794 Dimension_gen,
795 PArray[0],
796 LocalPDerivatives[0]) ;
797 //
798 // unormalize derivatives since those are computed normalized
799 //
800
801 ModifyCoords (LocalPDerivatives + Dimension_gen, /= SpanLenght);
802
0e14656b 803 if (WeightsArray != NULL) {
804 const TColStd_Array1OfReal& refWeights = *WeightsArray;
7fd59977 805 Standard_Real *
0e14656b 806 WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
7fd59977 807 PLib::EvalPolynomial(NewParameter,
808 1,
809 Degree,
810 1,
811 WArray[0],
812 LocalWDerivatives[0]) ;
813 //
814 // unormalize the result since the polynomial stored in the cache
815 // is normalized between 0 and 1
816 //
817 LocalWDerivatives[1] /= SpanLenght ;
818
819 PLib::RationalDerivatives(1,
820 Dimension_gen,
821 LocalPDerivatives[0],
822 LocalWDerivatives[0],
823 LocalPDerivatives[0]) ;
824 }
825
826 CopyCoords (myPoint, LocalPDerivatives);
827 CopyCoords (myVector, LocalPDerivatives + Dimension_gen);
828}
829
830//=======================================================================
831//function : CacheD2
832//purpose : Evaluates the polynomial cache of the Bspline Curve
833//
834//=======================================================================
835void BSplCLib::CacheD2(const Standard_Real Parameter,
836 const Standard_Integer Degree,
837 const Standard_Real CacheParameter,
838 const Standard_Real SpanLenght,
839 const Array1OfPoints& PolesArray,
0e14656b 840 const TColStd_Array1OfReal* WeightsArray,
7fd59977 841 Point& aPoint,
842 Vector& aVector1,
843 Vector& aVector2)
844{
845 //
846 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
847 // form
848 // the SpanLenght is the normalizing factor so that the polynomial is between
849 // 0 and 1
850 Standard_Integer ii,Index,EndIndex;
851 Standard_Real LocalPDerivatives[(Dimension_gen << 1) + Dimension_gen] ;
852// Standard_Real LocalWDerivatives[3], NewParameter, Factor, Inverse ;
853 Standard_Real LocalWDerivatives[3], NewParameter, Factor ;
854 Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
855 Standard_Real * myPoint = (Standard_Real *) &aPoint ;
856 Standard_Real * myVector1 = (Standard_Real *) &aVector1 ;
857 Standard_Real * myVector2 = (Standard_Real *) &aVector2 ;
858 NewParameter = (Parameter - CacheParameter) / SpanLenght ;
859 PLib::EvalPolynomial(NewParameter,
860 2,
861 Degree,
862 Dimension_gen,
863 PArray[0],
864 LocalPDerivatives[0]) ;
865 //
866 // unormalize derivatives since those are computed normalized
867 //
868 Factor = 1.0e0/ SpanLenght ;
869 Index = Dimension_gen ;
870 EndIndex = Min (2, Degree) ;
871
872 for (ii = 1 ; ii <= EndIndex ; ii++) {
873 ModifyCoords (LocalPDerivatives + Index, *= Factor);
874 Factor /= SpanLenght ;
875 Index += Dimension_gen;
876 }
877
878 Index = (Degree + 1) * Dimension_gen;
879 for (ii = Degree ; ii < 2 ; ii++) {
880 NullifyCoords (LocalPDerivatives + Index);
881 Index += Dimension_gen;
882 }
883
0e14656b 884 if (WeightsArray != NULL) {
885 const TColStd_Array1OfReal& refWeights = *WeightsArray;
7fd59977 886 Standard_Real *
0e14656b 887 WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
7fd59977 888
889 PLib::EvalPolynomial(NewParameter,
890 2,
891 Degree,
892 1,
893 WArray[0],
894 LocalWDerivatives[0]) ;
895
896 for (ii = Degree + 1 ; ii <= 2 ; ii++) {
897 LocalWDerivatives[ii] = 0.0e0 ;
898 }
899 //
900 // unormalize the result since the polynomial stored in the cache
901 // is normalized between 0 and 1
902 //
903 Factor = 1.0e0 / SpanLenght ;
904
905 for (ii = 1 ; ii <= EndIndex ; ii++) {
906 LocalWDerivatives[ii] *= Factor ;
907 Factor /= SpanLenght ;
908 }
909 PLib::RationalDerivatives(2,
910 Dimension_gen,
911 LocalPDerivatives[0],
912 LocalWDerivatives[0],
913 LocalPDerivatives[0]) ;
914 }
915
916 CopyCoords (myPoint, LocalPDerivatives);
917 CopyCoords (myVector1, LocalPDerivatives + Dimension_gen);
918 CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2);
919}
920
921//=======================================================================
922//function : CacheD3
923//purpose : Evaluates the polynomial cache of the Bspline Curve
924//
925//=======================================================================
926void BSplCLib::CacheD3(const Standard_Real Parameter,
927 const Standard_Integer Degree,
928 const Standard_Real CacheParameter,
929 const Standard_Real SpanLenght,
930 const Array1OfPoints& PolesArray,
0e14656b 931 const TColStd_Array1OfReal* WeightsArray,
7fd59977 932 Point& aPoint,
933 Vector& aVector1,
934 Vector& aVector2,
935 Vector& aVector3)
936{
937 //
938 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
939 // form
940 // the SpanLenght is the normalizing factor so that the polynomial is between
941 // 0 and 1
942 Standard_Integer ii, Index, EndIndex;
943 Standard_Real LocalPDerivatives[Dimension_gen << 2] ;
944// Standard_Real LocalWDerivatives[4], Factor, NewParameter, Inverse ;
945 Standard_Real LocalWDerivatives[4], Factor, NewParameter ;
946 Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
947 Standard_Real * myPoint = (Standard_Real *) &aPoint ;
948 Standard_Real * myVector1 = (Standard_Real *) &aVector1 ;
949 Standard_Real * myVector2 = (Standard_Real *) &aVector2 ;
950 Standard_Real * myVector3 = (Standard_Real *) &aVector3 ;
951 NewParameter = (Parameter - CacheParameter) / SpanLenght ;
952 PLib::EvalPolynomial(NewParameter,
953 3,
954 Degree,
955 Dimension_gen,
956 PArray[0],
957 LocalPDerivatives[0]) ;
958
959 Index = (Degree + 1) * Dimension_gen;
960 for (ii = Degree ; ii < 3 ; ii++) {
961 NullifyCoords (LocalPDerivatives + Index);
962 Index += Dimension_gen;
963 }
964
965 //
966 // unormalize derivatives since those are computed normalized
967 //
968 Factor = 1.0e0 / SpanLenght ;
969 Index = Dimension_gen ;
970 EndIndex = Min(3,Degree) ;
971
972 for (ii = 1 ; ii <= EndIndex ; ii++) {
973 ModifyCoords (LocalPDerivatives + Index, *= Factor);
974 Factor /= SpanLenght;
975 Index += Dimension_gen;
976 }
977
0e14656b 978 if (WeightsArray != NULL) {
979 const TColStd_Array1OfReal& refWeights = *WeightsArray;
7fd59977 980 Standard_Real *
0e14656b 981 WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
7fd59977 982
983 PLib::EvalPolynomial(NewParameter,
984 3,
985 Degree,
986 1,
987 WArray[0],
988 LocalWDerivatives[0]) ;
989 //
990 // unormalize the result since the polynomial stored in the cache
991 // is normalized between 0 and 1
992 //
993 Factor = 1.0e0 / SpanLenght ;
994
995 for (ii = 1 ; ii <= EndIndex ; ii++) {
996 LocalWDerivatives[ii] *= Factor ;
997 Factor /= SpanLenght ;
998 }
999
1000 for (ii = (Degree + 1) ; ii <= 3 ; ii++) {
1001 LocalWDerivatives[ii] = 0.0e0 ;
1002 }
1003 PLib::RationalDerivatives(3,
1004 Dimension_gen,
1005 LocalPDerivatives[0],
1006 LocalWDerivatives[0],
1007 LocalPDerivatives[0]) ;
1008 }
1009
1010 CopyCoords (myPoint, LocalPDerivatives);
1011 CopyCoords (myVector1, LocalPDerivatives + Dimension_gen);
1012 CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2);
1013 CopyCoords (myVector3, LocalPDerivatives + Dimension_gen * 3);
1014}
1015
1016//=======================================================================
1017//function : BuildCache
1018//purpose : Stores theTaylor expansion normalized between 0,1 in the
1019// Cache : in case of a rational function the Poles are
1020// stored in homogeneous form
1021//=======================================================================
1022
1023void BSplCLib::BuildCache
1024(const Standard_Real U,
1025 const Standard_Real SpanDomain,
1026 const Standard_Boolean Periodic,
1027 const Standard_Integer Degree,
1028 const TColStd_Array1OfReal& FlatKnots,
1029 const Array1OfPoints& Poles,
0e14656b 1030 const TColStd_Array1OfReal* Weights,
7fd59977 1031 Array1OfPoints& CachePoles,
0e14656b 1032 TColStd_Array1OfReal* CacheWeights)
7fd59977 1033{
1034 Standard_Integer ii,
1035 Dimension,
1036 LocalIndex,
1037 index = 0 ;
1038 Standard_Real u = U,
1039 LocalValue ;
1040 Standard_Boolean rational;
1041
1042 BSplCLib_DataContainer dc(Degree);
1043 PrepareEval(u,
1044 index,
1045 Dimension,
1046 rational,
1047 Degree,
1048 Periodic,
1049 Poles,
1050 Weights,
1051 FlatKnots,
1052 (BSplCLib::NoMults()),
1053 dc);
1054 //
1055 // Watch out : PrepareEval checks if locally the Bspline is polynomial
1056 // therefore rational flag could be set to False even if there are
1057 // Weights. The Dimension is set accordingly.
1058 //
1059
1060 BSplCLib::Bohm(u,Degree,Degree,*dc.knots,Dimension,*dc.poles);
1061
1062 LocalValue = 1.0e0 ;
1063 LocalIndex = 0 ;
1064
1065 if (rational) {
1066
1067 for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1068 CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue);
1069 LocalIndex += Dimension_gen + 1;
1070 LocalValue *= SpanDomain / (Standard_Real) ii ;
1071 }
1072
1073 LocalIndex = Dimension_gen;
1074 LocalValue = 1.0e0 ;
1075 for (ii = 1 ; ii <= Degree + 1 ; ii++) {
0e14656b 1076 (*CacheWeights)(ii) = dc.poles[LocalIndex] * LocalValue ;
7fd59977 1077 LocalIndex += Dimension_gen + 1;
1078 LocalValue *= SpanDomain / (Standard_Real) ii ;
1079 }
1080 }
1081 else {
1082
1083 for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1084 CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue);
1085 LocalIndex += Dimension_gen;
1086 LocalValue *= SpanDomain / (Standard_Real) ii ;
1087 }
1088
0e14656b 1089 if (Weights != NULL) {
7fd59977 1090 for (ii = 1 ; ii <= Degree + 1 ; ii++)
0e14656b 1091 (*CacheWeights)(ii) = 0.0e0 ;
1092 (*CacheWeights)(1) = 1.0e0 ;
7fd59977 1093 }
1094 }
1095}
1096
94f71cad 1097void BSplCLib::BuildCache(const Standard_Real theParameter,
1098 const Standard_Real theSpanDomain,
1099 const Standard_Boolean thePeriodicFlag,
1100 const Standard_Integer theDegree,
a061150b 1101 const Standard_Integer theSpanIndex,
94f71cad 1102 const TColStd_Array1OfReal& theFlatKnots,
1103 const Array1OfPoints& thePoles,
0e14656b 1104 const TColStd_Array1OfReal* theWeights,
94f71cad 1105 TColStd_Array2OfReal& theCacheArray)
1106{
1107 Standard_Real aParam = theParameter;
a061150b 1108 Standard_Integer anIndex = theSpanIndex;
94f71cad 1109 Standard_Integer aDimension;
1110 Standard_Boolean isRational;
1111
1112 BSplCLib_DataContainer dc(theDegree);
1113 PrepareEval(aParam,
1114 anIndex,
1115 aDimension,
1116 isRational,
1117 theDegree,
1118 thePeriodicFlag,
1119 thePoles,
1120 theWeights,
1121 theFlatKnots,
1122 (BSplCLib::NoMults()),
1123 dc);
1124 //
1125 // Watch out : PrepareEval checks if locally the Bspline is polynomial
1126 // therefore rational flag could be set to False even if there are
1127 // Weights. The Dimension is set accordingly.
1128 //
1129
1130 Standard_Integer aCacheShift = // helps to store weights when PrepareEval did not found that the curve is locally polynomial
0e14656b 1131 (theWeights != NULL && !isRational) ? aDimension + 1 : aDimension;
94f71cad 1132
1133 BSplCLib::Bohm(aParam, theDegree, theDegree, *dc.knots, aDimension, *dc.poles);
1134
1135 Standard_Real aCoeff = 1.0;
1136 Standard_Real* aCache = (Standard_Real *) &(theCacheArray(theCacheArray.LowerRow(), theCacheArray.LowerCol()));
1137 Standard_Real* aPolyCoeffs = dc.poles;
1138
1139 for (Standard_Integer i = 0; i <= theDegree; i++)
1140 {
1141 for (Standard_Integer j = 0; j< aDimension; j++)
1142 aCache[j] = aPolyCoeffs[j] * aCoeff;
1143 aCoeff *= theSpanDomain / (i+1);
1144 aPolyCoeffs += aDimension;
1145 aCache += aDimension;
1146 if (aCacheShift > aDimension)
1147 {
1148 aCache[0] = 0.0;
1149 aCache++;
1150 }
1151 }
1152
1153 if (aCacheShift > aDimension)
1154 theCacheArray.SetValue(theCacheArray.LowerRow(), theCacheArray.LowerCol() + aCacheShift - 1, 1.0);
1155}
1156
7fd59977 1157//=======================================================================
1158//function : Interpolate
1159//purpose :
1160//=======================================================================
1161
1162void BSplCLib::Interpolate(const Standard_Integer Degree,
1163 const TColStd_Array1OfReal& FlatKnots,
1164 const TColStd_Array1OfReal& Parameters,
1165 const TColStd_Array1OfInteger& ContactOrderArray,
1166 Array1OfPoints& Poles,
1167 Standard_Integer& InversionProblem)
1168
1169{
1170 Standard_Real *PArray ;
1171
1172 PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1173
1174 BSplCLib::Interpolate(Degree,
1175 FlatKnots,
1176 Parameters,
1177 ContactOrderArray,
1178 Dimension_gen,
1179 PArray[0],
1180 InversionProblem) ;
1181}
1182
1183//=======================================================================
1184//function : Interpolate
1185//purpose :
1186//=======================================================================
1187
1188void BSplCLib::Interpolate(const Standard_Integer Degree,
1189 const TColStd_Array1OfReal& FlatKnots,
1190 const TColStd_Array1OfReal& Parameters,
1191 const TColStd_Array1OfInteger& ContactOrderArray,
1192 Array1OfPoints& Poles,
1193 TColStd_Array1OfReal& Weights,
1194 Standard_Integer& InversionProblem)
1195{
1196 Standard_Real *PArray,
1197 * WArray ;
1198 PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1199 WArray = (Standard_Real *) &Weights(Weights.Lower()) ;
1200 BSplCLib::Interpolate(Degree,
1201 FlatKnots,
1202 Parameters,
1203 ContactOrderArray,
1204 Dimension_gen,
1205 PArray[0],
1206 WArray[0],
1207 InversionProblem) ;
1208}
1209
1210//=======================================================================
1211//function : MovePoint
1212//purpose : Find the new poles which allows an old point (with a
1213// given u as parameter) to reach a new position
1214//=======================================================================
1215
1216void BSplCLib::MovePoint (const Standard_Real U,
1217 const Vector& Displ,
1218 const Standard_Integer Index1,
1219 const Standard_Integer Index2,
1220 const Standard_Integer Degree,
7fd59977 1221 const Array1OfPoints& Poles,
aff73fd5 1222 const TColStd_Array1OfReal* Weights,
7fd59977 1223 const TColStd_Array1OfReal& FlatKnots,
1224 Standard_Integer& FirstIndex,
1225 Standard_Integer& LastIndex,
1226 Array1OfPoints& NewPoles)
1227{
1228 // calculate the BSplineBasis in the parameter U
1229 Standard_Integer FirstNonZeroBsplineIndex;
1230 math_Matrix BSplineBasis(1, 1,
1231 1, Degree+1);
1232 Standard_Integer ErrorCode =
6143f12f 1233 BSplCLib::EvalBsplineBasis(0,
7fd59977 1234 Degree+1,
1235 FlatKnots,
1236 U,
1237 FirstNonZeroBsplineIndex,
1238 BSplineBasis);
1239 if (ErrorCode != 0) {
1240 FirstIndex = 0;
1241 LastIndex = 0;
1242
1243 for (Standard_Integer i=Poles.Lower(); i<=Poles.Upper(); i++) {
1244 NewPoles(i) = Poles(i);
1245 }
1246 return;
1247 }
1248
1249 // find the span which is predominant for parameter U
1250 FirstIndex = FirstNonZeroBsplineIndex;
1251 LastIndex = FirstNonZeroBsplineIndex + Degree ;
1252 if (FirstIndex < Index1) FirstIndex = Index1;
1253 if (LastIndex > Index2) LastIndex = Index2;
1254
1255 Standard_Real maxValue = 0.0;
1256 Standard_Integer i, kk1=0, kk2, ii;
1257
1258 for (i = FirstIndex-FirstNonZeroBsplineIndex+1;
1259 i <= LastIndex-FirstNonZeroBsplineIndex+1; i++) {
1260 if (BSplineBasis(1,i) > maxValue) {
1261 kk1 = i + FirstNonZeroBsplineIndex - 1;
1262 maxValue = BSplineBasis(1,i);
1263 }
1264 }
1265
1266 // find a kk2 if symetriy
1267 kk2 = kk1;
1268 i = kk1 - FirstNonZeroBsplineIndex + 2;
1269 if ((kk1+1) <= LastIndex) {
1270 if (Abs(BSplineBasis(1, kk1-FirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
1271 kk2 = kk1+1;
1272 }
1273 }
1274
1275 // compute the vector of displacement
1276 Standard_Real D1 = 0.0;
1277 Standard_Real D2 = 0.0;
1278 Standard_Real hN, Coef, Dval;
1279
1280 for (i = 1; i <= Degree+1; i++) {
1281 ii = i + FirstNonZeroBsplineIndex - 1;
aff73fd5 1282 if (Weights != NULL) {
1283 hN = Weights->Value(ii)*BSplineBasis(1, i);
7fd59977 1284 D2 += hN;
1285 }
1286 else {
1287 hN = BSplineBasis(1, i);
1288 }
1289 if (ii >= FirstIndex && ii <= LastIndex) {
1290 if (ii < kk1) {
1291 Dval = kk1-ii;
1292 }
1293 else if (ii > kk2) {
1294 Dval = ii - kk2;
1295 }
1296 else {
1297 Dval = 0.0;
1298 }
1299 D1 += 1./(Dval + 1.) * hN;
1300 }
1301 }
1302
aff73fd5 1303 if (Weights != NULL) {
7fd59977 1304 Coef = D2/D1;
1305 }
1306 else {
1307 Coef = 1./D1;
1308 }
1309
1310 // compute the new poles
1311
1312 for (i=Poles.Lower(); i<=Poles.Upper(); i++) {
1313 if (i >= FirstIndex && i <= LastIndex) {
1314 if (i < kk1) {
1315 Dval = kk1-i;
1316 }
1317 else if (i > kk2) {
1318 Dval = i - kk2;
1319 }
1320 else {
1321 Dval = 0.0;
1322 }
1323 NewPoles(i) = Poles(i).Translated((Coef/(Dval+1.))*Displ);
1324 }
1325 else {
1326 NewPoles(i) = Poles(i);
1327 }
1328 }
1329}
1330
1331//=======================================================================
1332//function : MovePoint
1333//purpose : Find the new poles which allows an old point (with a
1334// given u as parameter) to reach a new position
1335//=======================================================================
1336
1337//=======================================================================
1338//function : MovePointAndTangent
1339//purpose :
1340//=======================================================================
1341
1342void BSplCLib::MovePointAndTangent (const Standard_Real U,
1343 const Vector& Delta,
1344 const Vector& DeltaDerivatives,
1345 const Standard_Real Tolerance,
1346 const Standard_Integer Degree,
7fd59977 1347 const Standard_Integer StartingCondition,
1348 const Standard_Integer EndingCondition,
1349 const Array1OfPoints& Poles,
aff73fd5 1350 const TColStd_Array1OfReal* Weights,
7fd59977 1351 const TColStd_Array1OfReal& FlatKnots,
1352 Array1OfPoints& NewPoles,
1353 Standard_Integer & ErrorStatus)
1354{
1355 Standard_Real *delta_array,
1356 *delta_derivative_array,
1357 *poles_array,
1358 *new_poles_array ;
1359
1360 Standard_Integer num_poles ;
1361 num_poles = Poles.Length() ;
1362
1363 if (NewPoles.Length() != num_poles) {
9775fa61 1364 throw Standard_ConstructionError();
7fd59977 1365 }
1366 delta_array = (Standard_Real *) &Delta ;
1367 delta_derivative_array = (Standard_Real *)&DeltaDerivatives ;
1368 poles_array = (Standard_Real *) &Poles(Poles.Lower()) ;
1369
1370 new_poles_array = (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1371 MovePointAndTangent(U,
1372 Dimension_gen,
1373 delta_array[0],
1374 delta_derivative_array[0],
1375 Tolerance,
1376 Degree,
7fd59977 1377 StartingCondition,
1378 EndingCondition,
1379 poles_array[0],
1380 Weights,
1381 FlatKnots,
1382 new_poles_array[0],
1383 ErrorStatus) ;
1384}
1385
1386//=======================================================================
1387//function : Resolution
1388//purpose :
1389//=======================================================================
1390
1391void BSplCLib::Resolution(const Array1OfPoints& Poles,
0e14656b 1392 const TColStd_Array1OfReal* Weights,
7fd59977 1393 const Standard_Integer NumPoles,
1394 const TColStd_Array1OfReal& FlatKnots,
1395 const Standard_Integer Degree,
1396 const Standard_Real Tolerance3D,
1397 Standard_Real& UTolerance)
1398{
1399 Standard_Real *PolesArray ;
1400 PolesArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1401 BSplCLib::Resolution(PolesArray[0],
1402 Dimension_gen,
1403 NumPoles,
1404 Weights,
1405 FlatKnots,
1406 Degree,
1407 Tolerance3D,
1408 UTolerance) ;
1409}
1410
1411//=======================================================================
1412//function : FunctionMultiply
1413//purpose :
1414//=======================================================================
1415
1416void BSplCLib::FunctionMultiply
1417(const BSplCLib_EvaluatorFunction & FunctionPtr,
1418 const Standard_Integer BSplineDegree,
1419 const TColStd_Array1OfReal & BSplineFlatKnots,
1420 const Array1OfPoints & Poles,
1421 const TColStd_Array1OfReal & FlatKnots,
1422 const Standard_Integer NewDegree,
1423 Array1OfPoints & NewPoles,
9fd2d2c3 1424 Standard_Integer & theStatus)
7fd59977 1425{
1426 Standard_Integer num_bspline_poles =
1427 BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1428 Standard_Integer num_new_poles =
1429 FlatKnots.Length() - NewDegree - 1 ;
1430
1431 if (Poles.Length() != num_bspline_poles ||
1432 NewPoles.Length() != num_new_poles) {
9775fa61 1433 throw Standard_ConstructionError();
7fd59977 1434 }
1435 Standard_Real * array_of_poles =
1436 (Standard_Real *) &Poles(Poles.Lower()) ;
1437 Standard_Real * array_of_new_poles =
1438 (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1439 BSplCLib::FunctionMultiply(FunctionPtr,
1440 BSplineDegree,
1441 BSplineFlatKnots,
1442 Dimension_gen,
1443 array_of_poles[0],
1444 FlatKnots,
1445 NewDegree,
1446 array_of_new_poles[0],
9fd2d2c3 1447 theStatus);
7fd59977 1448}
1449
1450//=======================================================================
1451//function : FunctionReparameterise
1452//purpose :
1453//=======================================================================
1454
1455void BSplCLib::FunctionReparameterise
1456(const BSplCLib_EvaluatorFunction & FunctionPtr,
1457 const Standard_Integer BSplineDegree,
1458 const TColStd_Array1OfReal & BSplineFlatKnots,
1459 const Array1OfPoints & Poles,
1460 const TColStd_Array1OfReal & FlatKnots,
1461 const Standard_Integer NewDegree,
1462 Array1OfPoints & NewPoles,
9fd2d2c3 1463 Standard_Integer & theStatus)
7fd59977 1464{
1465 Standard_Integer num_bspline_poles =
1466 BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1467 Standard_Integer num_new_poles =
1468 FlatKnots.Length() - NewDegree - 1 ;
1469
1470 if (Poles.Length() != num_bspline_poles ||
1471 NewPoles.Length() != num_new_poles) {
9775fa61 1472 throw Standard_ConstructionError();
7fd59977 1473 }
1474 Standard_Real * array_of_poles =
1475 (Standard_Real *) &Poles(Poles.Lower()) ;
1476 Standard_Real * array_of_new_poles =
1477 (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1478 BSplCLib::FunctionReparameterise(FunctionPtr,
1479 BSplineDegree,
1480 BSplineFlatKnots,
1481 Dimension_gen,
1482 array_of_poles[0],
1483 FlatKnots,
1484 NewDegree,
1485 array_of_new_poles[0],
9fd2d2c3 1486 theStatus);
7fd59977 1487}