1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
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
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.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
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
20 #include <TColStd_Array1OfInteger.hxx>
21 #include <TColStd_Array1OfReal.hxx>
22 #include <TColStd_Array2OfReal.hxx>
23 #include <gp_Vec2d.hxx>
24 #include <Standard_ConstructionError.hxx>
26 #include <math_Matrix.hxx>
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 //=======================================================================
34 struct BSplCLib_DataContainer
36 BSplCLib_DataContainer(Standard_Integer Degree)
38 (void)Degree; // avoid compiler warning
39 Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree() ||
40 BSplCLib::MaxDegree() > 25,
41 "BSplCLib: bspline degree is greater than maximum supported");
44 Standard_Real poles[(25+1)*(Dimension_gen+1)];
45 Standard_Real knots[2*25];
46 Standard_Real ders[Dimension_gen*4];
49 //=======================================================================
52 //=======================================================================
54 void BSplCLib::Reverse(Array1OfPoints& Poles,
55 const Standard_Integer L)
57 Standard_Integer i, l = L;
58 l = Poles.Lower() + (l-Poles.Lower()) % (Poles.Upper()-Poles.Lower()+1);
60 Array1OfPoints temp(0,Poles.Length()-1);
62 for (i = Poles.Lower(); i <= l; i++)
65 for (i = l+1; i <= Poles.Upper(); i++)
66 temp(l-Poles.Lower()+Poles.Upper()-i+1) = Poles(i);
68 for (i = Poles.Lower(); i <= Poles.Upper(); i++)
69 Poles(i) = temp(i-Poles.Lower());
73 // CURVES MODIFICATIONS
76 //=======================================================================
77 //function : RemoveKnot
79 //=======================================================================
81 Standard_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,
87 const TColStd_Array1OfReal& Weights,
88 const TColStd_Array1OfReal& Knots,
89 const TColStd_Array1OfInteger& Mults,
90 Array1OfPoints& NewPoles,
91 TColStd_Array1OfReal& NewWeights,
92 TColStd_Array1OfReal& NewKnots,
93 TColStd_Array1OfInteger& NewMults,
94 const Standard_Real Tolerance)
96 Standard_Boolean rational = &Weights != NULL;
101 TColStd_Array1OfReal poles(1,dim*Poles.Length());
102 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
104 if (rational) PLib::SetPoles(Poles,Weights,poles);
105 else PLib::SetPoles(Poles,poles);
107 if (!RemoveKnot(Index,Mult,Degree,Periodic,dim,
108 poles,Knots,Mults,newpoles,NewKnots,NewMults,Tolerance))
109 return Standard_False;
111 if (rational) PLib::GetPoles(newpoles,NewPoles,NewWeights);
112 else PLib::GetPoles(newpoles,NewPoles);
113 return Standard_True;
116 //=======================================================================
117 //function : InsertKnots
118 //purpose : insert an array of knots and multiplicities
119 //=======================================================================
121 void BSplCLib::InsertKnots
122 (const Standard_Integer Degree,
123 const Standard_Boolean Periodic,
124 const Array1OfPoints& Poles,
125 const TColStd_Array1OfReal& Weights,
126 const TColStd_Array1OfReal& Knots,
127 const TColStd_Array1OfInteger& Mults,
128 const TColStd_Array1OfReal& AddKnots,
129 const TColStd_Array1OfInteger& AddMults,
130 Array1OfPoints& NewPoles,
131 TColStd_Array1OfReal& NewWeights,
132 TColStd_Array1OfReal& NewKnots,
133 TColStd_Array1OfInteger& NewMults,
134 const Standard_Real Epsilon,
135 const Standard_Boolean Add)
137 Standard_Boolean rational = &Weights != NULL;
138 Standard_Integer dim;
142 TColStd_Array1OfReal poles(1,dim*Poles.Length());
143 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
145 if (rational) PLib::SetPoles(Poles,Weights,poles);
146 else PLib::SetPoles(Poles,poles);
148 InsertKnots(Degree,Periodic,dim,poles,Knots,Mults,
149 AddKnots,AddMults,newpoles,NewKnots,NewMults,Epsilon,Add);
151 if (rational) PLib::GetPoles(newpoles,NewPoles,NewWeights);
152 else PLib::GetPoles(newpoles,NewPoles);
155 //=======================================================================
156 //function : InsertKnot
158 //=======================================================================
160 void 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,
166 const TColStd_Array1OfReal& Weights,
167 const TColStd_Array1OfReal& Knots,
168 const TColStd_Array1OfInteger& Mults,
169 Array1OfPoints& NewPoles,
170 TColStd_Array1OfReal& NewWeights)
172 TColStd_Array1OfReal k(1,1);
174 TColStd_Array1OfInteger m(1,1);
176 TColStd_Array1OfReal nk(1,Knots.Length()+1);
177 TColStd_Array1OfInteger nm(1,Knots.Length()+1);
178 InsertKnots(Degree,Periodic,Poles,Weights,Knots,Mults,
179 k,m,NewPoles,NewWeights,nk,nm,Epsilon(U));
182 //=======================================================================
183 //function : RaiseMultiplicity
185 //=======================================================================
187 void 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,
192 const TColStd_Array1OfReal& Weights,
193 const TColStd_Array1OfReal& Knots,
194 const TColStd_Array1OfInteger& Mults,
195 Array1OfPoints& NewPoles,
196 TColStd_Array1OfReal& NewWeights)
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,
205 k,m,NewPoles,NewWeights,nk,nm,Epsilon(k(1)));
208 //=======================================================================
209 //function : IncreaseDegree
211 //=======================================================================
213 void BSplCLib::IncreaseDegree
214 (const Standard_Integer Degree,
215 const Standard_Integer NewDegree,
216 const Standard_Boolean Periodic,
217 const Array1OfPoints& Poles,
218 const TColStd_Array1OfReal& Weights,
219 const TColStd_Array1OfReal& Knots,
220 const TColStd_Array1OfInteger& Mults,
221 Array1OfPoints& NewPoles,
222 TColStd_Array1OfReal& NewWeights,
223 TColStd_Array1OfReal& NewKnots,
224 TColStd_Array1OfInteger& NewMults)
226 Standard_Boolean rational = &Weights != NULL;
227 Standard_Integer dim;
231 TColStd_Array1OfReal poles(1,dim*Poles.Length());
232 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
234 if (rational) PLib::SetPoles(Poles,Weights,poles);
235 else PLib::SetPoles(Poles,poles);
237 IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults,
238 newpoles,NewKnots,NewMults);
240 if (rational) PLib::GetPoles(newpoles,NewPoles,NewWeights);
241 else PLib::GetPoles(newpoles,NewPoles);
244 //=======================================================================
245 //function : Unperiodize
247 //=======================================================================
249 void BSplCLib::Unperiodize
250 (const Standard_Integer Degree,
251 const TColStd_Array1OfInteger& Mults,
252 const TColStd_Array1OfReal& Knots,
253 const Array1OfPoints& Poles,
254 const TColStd_Array1OfReal& Weights,
255 TColStd_Array1OfInteger& NewMults,
256 TColStd_Array1OfReal& NewKnots,
257 Array1OfPoints& NewPoles,
258 TColStd_Array1OfReal& NewWeights)
260 Standard_Boolean rational = &Weights != NULL;
261 Standard_Integer dim;
265 TColStd_Array1OfReal poles(1,dim*Poles.Length());
266 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
268 if (rational) PLib::SetPoles(Poles,Weights,poles);
269 else PLib::SetPoles(Poles,poles);
271 Unperiodize(Degree, dim, Mults, Knots, poles,
272 NewMults,NewKnots, newpoles);
274 if (rational) PLib::GetPoles(newpoles,NewPoles,NewWeights);
275 else PLib::GetPoles(newpoles,NewPoles);
278 //=======================================================================
279 //function : Trimming
281 //=======================================================================
283 void 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,
288 const TColStd_Array1OfReal& Weights,
289 const Standard_Real U1,
290 const Standard_Real U2,
291 TColStd_Array1OfReal& NewKnots,
292 TColStd_Array1OfInteger& NewMults,
293 Array1OfPoints& NewPoles,
294 TColStd_Array1OfReal& NewWeights)
296 Standard_Boolean rational = &Weights != NULL;
297 Standard_Integer dim;
301 TColStd_Array1OfReal poles(1,dim*Poles.Length());
302 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
304 if (rational) PLib::SetPoles(Poles,Weights,poles);
305 else PLib::SetPoles(Poles,poles);
307 Trimming(Degree, Periodic, dim, Knots, Mults, poles, U1, U2,
308 NewKnots, NewMults, newpoles);
310 if (rational) PLib::GetPoles(newpoles,NewPoles,NewWeights);
311 else PLib::GetPoles(newpoles,NewPoles);
314 //--------------------------------------------------------------------------
315 //ELEMENTARY COMPUTATIONS
316 //--------------------------------------------------------------------------
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
325 //--------------------------------------------------------------------------
327 //=======================================================================
328 //function : BuildEval
329 //purpose : builds the local array for evaluation
330 //=======================================================================
332 void BSplCLib::BuildEval(const Standard_Integer Degree,
333 const Standard_Integer Index,
334 const Array1OfPoints& Poles,
335 const TColStd_Array1OfReal& Weights,
338 Standard_Real w, *pole = &LP;
339 Standard_Integer PLower = Poles.Lower();
340 Standard_Integer PUpper = Poles.Upper();
342 Standard_Integer ip = PLower + Index - 1;
343 if (&Weights == NULL) {
344 for (i = 0; i <= Degree; i++) {
346 if (ip > PUpper) ip = PLower;
347 const Point& P = Poles(ip);
348 PointToCoords (pole,P,+0);
349 pole += Dimension_gen;
353 for (i = 0; i <= Degree; i++) {
355 if (ip > PUpper) ip = PLower;
356 const Point& P = Poles(ip);
357 pole[Dimension_gen] = w = Weights(ip);
358 PointToCoords (pole, P, * w);
359 pole += Dimension_gen + 1;
364 //=======================================================================
365 //function : PrepareEval
366 //purpose : stores data for Eval in the local arrays
367 // dc.poles and dc.knots
368 //=======================================================================
370 static void PrepareEval
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,
378 const TColStd_Array1OfReal& Weights,
379 const TColStd_Array1OfReal& Knots,
380 const TColStd_Array1OfInteger& Mults,
381 BSplCLib_DataContainer& dc)
384 BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
387 BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*dc.knots);
389 index -= Knots.Lower() + Degree;
391 index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults);
393 // check truly rational
394 rational = (&Weights != NULL);
396 Standard_Integer WLower = Weights.Lower() + index;
397 rational = BSplCLib::IsRational(Weights, WLower, WLower + Degree);
402 dim = Dimension_gen + 1;
403 BSplCLib::BuildEval(Degree, index, Poles, Weights , *dc.poles);
407 BSplCLib::BuildEval(Degree, index, Poles, BSplCLib::NoWeights(), *dc.poles);
411 //=======================================================================
414 //=======================================================================
417 (const Standard_Real U,
418 const Standard_Integer Index,
419 const Standard_Integer Degree,
420 const Standard_Boolean Periodic,
421 const Array1OfPoints& Poles,
422 const TColStd_Array1OfReal& Weights,
423 const TColStd_Array1OfReal& Knots,
424 const TColStd_Array1OfInteger& Mults,
427 // Standard_Integer k,dim,index = Index;
428 Standard_Integer dim,index = Index;
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);
436 Standard_Real w = dc.poles[Dimension_gen];
437 CoordsToPoint (P, dc.poles, / w);
440 CoordsToPoint (P, dc.poles, );
443 //=======================================================================
446 //=======================================================================
449 (const Standard_Real U,
450 const Standard_Integer Index,
451 const Standard_Integer Degree,
452 const Standard_Boolean Periodic,
453 const Array1OfPoints& Poles,
454 const TColStd_Array1OfReal& Weights,
455 const TColStd_Array1OfReal& Knots,
456 const TColStd_Array1OfInteger& Mults,
460 Standard_Integer dim,index = Index;
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;
468 PLib::RationalDerivative(Degree,1,Dimension_gen,*dc.poles,*dc.ders);
472 CoordsToPoint (P, result, );
473 CoordsToPoint (V, result + Dimension_gen, );
476 //=======================================================================
479 //=======================================================================
482 (const Standard_Real U,
483 const Standard_Integer Index,
484 const Standard_Integer Degree,
485 const Standard_Boolean Periodic,
486 const Array1OfPoints& Poles,
487 const TColStd_Array1OfReal& Weights,
488 const TColStd_Array1OfReal& Knots,
489 const TColStd_Array1OfInteger& Mults,
494 Standard_Integer dim,index = Index;
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;
502 PLib::RationalDerivative(Degree,2,Dimension_gen,*dc.poles,*dc.ders);
506 CoordsToPoint (P, result, );
507 CoordsToPoint (V1, result + Dimension_gen, );
508 if (!rational && (Degree < 2))
511 CoordsToPoint (V2, result + 2 * Dimension_gen, );
514 //=======================================================================
517 //=======================================================================
520 (const Standard_Real U,
521 const Standard_Integer Index,
522 const Standard_Integer Degree,
523 const Standard_Boolean Periodic,
524 const Array1OfPoints& Poles,
525 const TColStd_Array1OfReal& Weights,
526 const TColStd_Array1OfReal& Knots,
527 const TColStd_Array1OfInteger& Mults,
533 Standard_Integer dim,index = Index;
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;
541 PLib::RationalDerivative(Degree,3,Dimension_gen,*dc.poles,*dc.ders);
545 CoordsToPoint (P, result, );
546 CoordsToPoint (V1, result + Dimension_gen, );
547 if (!rational && (Degree < 2))
550 CoordsToPoint (V2, result + 2 * Dimension_gen, );
551 if (!rational && (Degree < 3))
554 CoordsToPoint (V3, result + 3 * Dimension_gen, );
557 //=======================================================================
560 //=======================================================================
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,
569 const TColStd_Array1OfReal& Weights,
570 const TColStd_Array1OfReal& Knots,
571 const TColStd_Array1OfInteger& Mults,
574 Standard_Integer dim,index = Index;
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);
582 Standard_Real v[Dimension_gen];
583 PLib::RationalDerivative(Degree,N,Dimension_gen,*dc.poles,v[0],Standard_False);
584 CoordsToPoint (VN, v, );
590 Standard_Real *DN = dc.poles + N * Dimension_gen;
591 CoordsToPoint (VN, DN, );
596 //=======================================================================
597 //function : Solves a LU factored Matrix
599 //=======================================================================
602 BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
603 const Standard_Integer UpperBandWidth,
604 const Standard_Integer LowerBandWidth,
605 Array1OfPoints& PolesArray)
607 Standard_Real *PArray ;
608 PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
610 return BSplCLib::SolveBandedSystem(Matrix,
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 //=======================================================================
627 BSplCLib::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)
634 Standard_Real *PArray,
636 PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
637 WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
638 return BSplCLib::SolveBandedSystem(Matrix,
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 //=======================================================================
656 void 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,
665 Standard_Real& aWeight)
667 Standard_Real Inverse,
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,
684 BSplCLib::Eval(Parameter,
695 BSplCLib::Eval(Parameter,
706 Inverse = 1.0e0 / aWeight ;
708 for (kk = 0 ; kk < Dimension_gen ; kk++) {
713 for (kk = 0 ; kk < Dimension_gen ; kk++)
714 aPoint.SetCoord(kk + 1, P[kk]) ;
717 //=======================================================================
719 //purpose : Evaluates the polynomial cache of the Bspline Curve
721 //=======================================================================
722 void 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,
727 const TColStd_Array1OfReal& WeightsArray,
731 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
733 // the SpanLenght is the normalizing factor so that the polynomial is between
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,
742 Degree * Dimension_gen,
745 if (&WeightsArray != NULL) {
747 WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
748 PLib::NoDerivativeEvalPolynomial(NewParameter,
755 Inverse = 1.0e0 / Inverse;
756 ModifyCoords (myPoint, *= Inverse);
760 //=======================================================================
762 //purpose : Evaluates the polynomial cache of the Bspline Curve
764 //=======================================================================
765 void BSplCLib::CacheD1(const Standard_Real Parameter,
766 const Standard_Integer Degree,
767 const Standard_Real CacheParameter,
768 const Standard_Real SpanLenght,
769 const Array1OfPoints& PolesArray,
770 const TColStd_Array1OfReal& WeightsArray,
775 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
777 // the SpanLenght is the normalizing factor so that the polynomial is between
779 Standard_Real LocalPDerivatives[Dimension_gen << 1] ;
780 // Standard_Real LocalWDerivatives[2], NewParameter, Inverse ;
781 Standard_Real LocalWDerivatives[2], NewParameter ;
784 PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
786 myPoint = (Standard_Real *) &aPoint ;
788 myVector = (Standard_Real *) &aVector ;
789 NewParameter = (Parameter - CacheParameter) / SpanLenght ;
790 PLib::EvalPolynomial(NewParameter,
795 LocalPDerivatives[0]) ;
797 // unormalize derivatives since those are computed normalized
800 ModifyCoords (LocalPDerivatives + Dimension_gen, /= SpanLenght);
802 if (&WeightsArray != NULL) {
804 WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
805 PLib::EvalPolynomial(NewParameter,
810 LocalWDerivatives[0]) ;
812 // unormalize the result since the polynomial stored in the cache
813 // is normalized between 0 and 1
815 LocalWDerivatives[1] /= SpanLenght ;
817 PLib::RationalDerivatives(1,
819 LocalPDerivatives[0],
820 LocalWDerivatives[0],
821 LocalPDerivatives[0]) ;
824 CopyCoords (myPoint, LocalPDerivatives);
825 CopyCoords (myVector, LocalPDerivatives + Dimension_gen);
828 //=======================================================================
830 //purpose : Evaluates the polynomial cache of the Bspline Curve
832 //=======================================================================
833 void BSplCLib::CacheD2(const Standard_Real Parameter,
834 const Standard_Integer Degree,
835 const Standard_Real CacheParameter,
836 const Standard_Real SpanLenght,
837 const Array1OfPoints& PolesArray,
838 const TColStd_Array1OfReal& WeightsArray,
844 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
846 // the SpanLenght is the normalizing factor so that the polynomial is between
848 Standard_Integer ii,Index,EndIndex;
849 Standard_Real LocalPDerivatives[(Dimension_gen << 1) + Dimension_gen] ;
850 // Standard_Real LocalWDerivatives[3], NewParameter, Factor, Inverse ;
851 Standard_Real LocalWDerivatives[3], NewParameter, Factor ;
852 Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
853 Standard_Real * myPoint = (Standard_Real *) &aPoint ;
854 Standard_Real * myVector1 = (Standard_Real *) &aVector1 ;
855 Standard_Real * myVector2 = (Standard_Real *) &aVector2 ;
856 NewParameter = (Parameter - CacheParameter) / SpanLenght ;
857 PLib::EvalPolynomial(NewParameter,
862 LocalPDerivatives[0]) ;
864 // unormalize derivatives since those are computed normalized
866 Factor = 1.0e0/ SpanLenght ;
867 Index = Dimension_gen ;
868 EndIndex = Min (2, Degree) ;
870 for (ii = 1 ; ii <= EndIndex ; ii++) {
871 ModifyCoords (LocalPDerivatives + Index, *= Factor);
872 Factor /= SpanLenght ;
873 Index += Dimension_gen;
876 Index = (Degree + 1) * Dimension_gen;
877 for (ii = Degree ; ii < 2 ; ii++) {
878 NullifyCoords (LocalPDerivatives + Index);
879 Index += Dimension_gen;
882 if (&WeightsArray != NULL) {
884 WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
886 PLib::EvalPolynomial(NewParameter,
891 LocalWDerivatives[0]) ;
893 for (ii = Degree + 1 ; ii <= 2 ; ii++) {
894 LocalWDerivatives[ii] = 0.0e0 ;
897 // unormalize the result since the polynomial stored in the cache
898 // is normalized between 0 and 1
900 Factor = 1.0e0 / SpanLenght ;
902 for (ii = 1 ; ii <= EndIndex ; ii++) {
903 LocalWDerivatives[ii] *= Factor ;
904 Factor /= SpanLenght ;
906 PLib::RationalDerivatives(2,
908 LocalPDerivatives[0],
909 LocalWDerivatives[0],
910 LocalPDerivatives[0]) ;
913 CopyCoords (myPoint, LocalPDerivatives);
914 CopyCoords (myVector1, LocalPDerivatives + Dimension_gen);
915 CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2);
918 //=======================================================================
920 //purpose : Evaluates the polynomial cache of the Bspline Curve
922 //=======================================================================
923 void BSplCLib::CacheD3(const Standard_Real Parameter,
924 const Standard_Integer Degree,
925 const Standard_Real CacheParameter,
926 const Standard_Real SpanLenght,
927 const Array1OfPoints& PolesArray,
928 const TColStd_Array1OfReal& WeightsArray,
935 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
937 // the SpanLenght is the normalizing factor so that the polynomial is between
939 Standard_Integer ii, Index, EndIndex;
940 Standard_Real LocalPDerivatives[Dimension_gen << 2] ;
941 // Standard_Real LocalWDerivatives[4], Factor, NewParameter, Inverse ;
942 Standard_Real LocalWDerivatives[4], Factor, NewParameter ;
943 Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
944 Standard_Real * myPoint = (Standard_Real *) &aPoint ;
945 Standard_Real * myVector1 = (Standard_Real *) &aVector1 ;
946 Standard_Real * myVector2 = (Standard_Real *) &aVector2 ;
947 Standard_Real * myVector3 = (Standard_Real *) &aVector3 ;
948 NewParameter = (Parameter - CacheParameter) / SpanLenght ;
949 PLib::EvalPolynomial(NewParameter,
954 LocalPDerivatives[0]) ;
956 Index = (Degree + 1) * Dimension_gen;
957 for (ii = Degree ; ii < 3 ; ii++) {
958 NullifyCoords (LocalPDerivatives + Index);
959 Index += Dimension_gen;
963 // unormalize derivatives since those are computed normalized
965 Factor = 1.0e0 / SpanLenght ;
966 Index = Dimension_gen ;
967 EndIndex = Min(3,Degree) ;
969 for (ii = 1 ; ii <= EndIndex ; ii++) {
970 ModifyCoords (LocalPDerivatives + Index, *= Factor);
971 Factor /= SpanLenght;
972 Index += Dimension_gen;
975 if (&WeightsArray != NULL) {
977 WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
979 PLib::EvalPolynomial(NewParameter,
984 LocalWDerivatives[0]) ;
986 // unormalize the result since the polynomial stored in the cache
987 // is normalized between 0 and 1
989 Factor = 1.0e0 / SpanLenght ;
991 for (ii = 1 ; ii <= EndIndex ; ii++) {
992 LocalWDerivatives[ii] *= Factor ;
993 Factor /= SpanLenght ;
996 for (ii = (Degree + 1) ; ii <= 3 ; ii++) {
997 LocalWDerivatives[ii] = 0.0e0 ;
999 PLib::RationalDerivatives(3,
1001 LocalPDerivatives[0],
1002 LocalWDerivatives[0],
1003 LocalPDerivatives[0]) ;
1006 CopyCoords (myPoint, LocalPDerivatives);
1007 CopyCoords (myVector1, LocalPDerivatives + Dimension_gen);
1008 CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2);
1009 CopyCoords (myVector3, LocalPDerivatives + Dimension_gen * 3);
1012 //=======================================================================
1013 //function : BuildCache
1014 //purpose : Stores theTaylor expansion normalized between 0,1 in the
1015 // Cache : in case of a rational function the Poles are
1016 // stored in homogeneous form
1017 //=======================================================================
1019 void BSplCLib::BuildCache
1020 (const Standard_Real U,
1021 const Standard_Real SpanDomain,
1022 const Standard_Boolean Periodic,
1023 const Standard_Integer Degree,
1024 const TColStd_Array1OfReal& FlatKnots,
1025 const Array1OfPoints& Poles,
1026 const TColStd_Array1OfReal& Weights,
1027 Array1OfPoints& CachePoles,
1028 TColStd_Array1OfReal& CacheWeights)
1030 Standard_Integer ii,
1034 Standard_Real u = U,
1036 Standard_Boolean rational;
1038 BSplCLib_DataContainer dc(Degree);
1048 (BSplCLib::NoMults()),
1051 // Watch out : PrepareEval checks if locally the Bspline is polynomial
1052 // therefore rational flag could be set to False even if there are
1053 // Weights. The Dimension is set accordingly.
1056 BSplCLib::Bohm(u,Degree,Degree,*dc.knots,Dimension,*dc.poles);
1058 LocalValue = 1.0e0 ;
1063 for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1064 CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue);
1065 LocalIndex += Dimension_gen + 1;
1066 LocalValue *= SpanDomain / (Standard_Real) ii ;
1069 LocalIndex = Dimension_gen;
1070 LocalValue = 1.0e0 ;
1071 for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1072 CacheWeights(ii) = dc.poles[LocalIndex] * LocalValue ;
1073 LocalIndex += Dimension_gen + 1;
1074 LocalValue *= SpanDomain / (Standard_Real) ii ;
1079 for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1080 CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue);
1081 LocalIndex += Dimension_gen;
1082 LocalValue *= SpanDomain / (Standard_Real) ii ;
1085 if (&Weights != NULL) {
1086 for (ii = 1 ; ii <= Degree + 1 ; ii++)
1087 CacheWeights(ii) = 0.0e0 ;
1088 CacheWeights(1) = 1.0e0 ;
1093 void BSplCLib::BuildCache(const Standard_Real theParameter,
1094 const Standard_Real theSpanDomain,
1095 const Standard_Boolean thePeriodicFlag,
1096 const Standard_Integer theDegree,
1097 const TColStd_Array1OfReal& theFlatKnots,
1098 const Array1OfPoints& thePoles,
1099 const TColStd_Array1OfReal& theWeights,
1100 TColStd_Array2OfReal& theCacheArray)
1102 Standard_Real aParam = theParameter;
1103 Standard_Integer anIndex = 0;
1104 Standard_Integer aDimension;
1105 Standard_Boolean isRational;
1107 BSplCLib_DataContainer dc(theDegree);
1117 (BSplCLib::NoMults()),
1120 // Watch out : PrepareEval checks if locally the Bspline is polynomial
1121 // therefore rational flag could be set to False even if there are
1122 // Weights. The Dimension is set accordingly.
1125 Standard_Integer aCacheShift = // helps to store weights when PrepareEval did not found that the curve is locally polynomial
1126 (&theWeights != NULL && !isRational) ? aDimension + 1 : aDimension;
1128 BSplCLib::Bohm(aParam, theDegree, theDegree, *dc.knots, aDimension, *dc.poles);
1130 Standard_Real aCoeff = 1.0;
1131 Standard_Real* aCache = (Standard_Real *) &(theCacheArray(theCacheArray.LowerRow(), theCacheArray.LowerCol()));
1132 Standard_Real* aPolyCoeffs = dc.poles;
1134 for (Standard_Integer i = 0; i <= theDegree; i++)
1136 for (Standard_Integer j = 0; j< aDimension; j++)
1137 aCache[j] = aPolyCoeffs[j] * aCoeff;
1138 aCoeff *= theSpanDomain / (i+1);
1139 aPolyCoeffs += aDimension;
1140 aCache += aDimension;
1141 if (aCacheShift > aDimension)
1148 if (aCacheShift > aDimension)
1149 theCacheArray.SetValue(theCacheArray.LowerRow(), theCacheArray.LowerCol() + aCacheShift - 1, 1.0);
1152 //=======================================================================
1153 //function : Interpolate
1155 //=======================================================================
1157 void BSplCLib::Interpolate(const Standard_Integer Degree,
1158 const TColStd_Array1OfReal& FlatKnots,
1159 const TColStd_Array1OfReal& Parameters,
1160 const TColStd_Array1OfInteger& ContactOrderArray,
1161 Array1OfPoints& Poles,
1162 Standard_Integer& InversionProblem)
1165 Standard_Real *PArray ;
1167 PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1169 BSplCLib::Interpolate(Degree,
1178 //=======================================================================
1179 //function : Interpolate
1181 //=======================================================================
1183 void BSplCLib::Interpolate(const Standard_Integer Degree,
1184 const TColStd_Array1OfReal& FlatKnots,
1185 const TColStd_Array1OfReal& Parameters,
1186 const TColStd_Array1OfInteger& ContactOrderArray,
1187 Array1OfPoints& Poles,
1188 TColStd_Array1OfReal& Weights,
1189 Standard_Integer& InversionProblem)
1191 Standard_Real *PArray,
1193 PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1194 WArray = (Standard_Real *) &Weights(Weights.Lower()) ;
1195 BSplCLib::Interpolate(Degree,
1205 //=======================================================================
1206 //function : MovePoint
1207 //purpose : Find the new poles which allows an old point (with a
1208 // given u as parameter) to reach a new position
1209 //=======================================================================
1211 void BSplCLib::MovePoint (const Standard_Real U,
1212 const Vector& Displ,
1213 const Standard_Integer Index1,
1214 const Standard_Integer Index2,
1215 const Standard_Integer Degree,
1216 const Standard_Boolean Rational,
1217 const Array1OfPoints& Poles,
1218 const TColStd_Array1OfReal& Weights,
1219 const TColStd_Array1OfReal& FlatKnots,
1220 Standard_Integer& FirstIndex,
1221 Standard_Integer& LastIndex,
1222 Array1OfPoints& NewPoles)
1224 // calculate the BSplineBasis in the parameter U
1225 Standard_Integer FirstNonZeroBsplineIndex;
1226 math_Matrix BSplineBasis(1, 1,
1228 Standard_Integer ErrorCode =
1229 BSplCLib::EvalBsplineBasis(1,
1234 FirstNonZeroBsplineIndex,
1236 if (ErrorCode != 0) {
1240 for (Standard_Integer i=Poles.Lower(); i<=Poles.Upper(); i++) {
1241 NewPoles(i) = Poles(i);
1246 // find the span which is predominant for parameter U
1247 FirstIndex = FirstNonZeroBsplineIndex;
1248 LastIndex = FirstNonZeroBsplineIndex + Degree ;
1249 if (FirstIndex < Index1) FirstIndex = Index1;
1250 if (LastIndex > Index2) LastIndex = Index2;
1252 Standard_Real maxValue = 0.0;
1253 Standard_Integer i, kk1=0, kk2, ii;
1255 for (i = FirstIndex-FirstNonZeroBsplineIndex+1;
1256 i <= LastIndex-FirstNonZeroBsplineIndex+1; i++) {
1257 if (BSplineBasis(1,i) > maxValue) {
1258 kk1 = i + FirstNonZeroBsplineIndex - 1;
1259 maxValue = BSplineBasis(1,i);
1263 // find a kk2 if symetriy
1265 i = kk1 - FirstNonZeroBsplineIndex + 2;
1266 if ((kk1+1) <= LastIndex) {
1267 if (Abs(BSplineBasis(1, kk1-FirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
1272 // compute the vector of displacement
1273 Standard_Real D1 = 0.0;
1274 Standard_Real D2 = 0.0;
1275 Standard_Real hN, Coef, Dval;
1277 for (i = 1; i <= Degree+1; i++) {
1278 ii = i + FirstNonZeroBsplineIndex - 1;
1280 hN = Weights(ii)*BSplineBasis(1, i);
1284 hN = BSplineBasis(1, i);
1286 if (ii >= FirstIndex && ii <= LastIndex) {
1290 else if (ii > kk2) {
1296 D1 += 1./(Dval + 1.) * hN;
1307 // compute the new poles
1309 for (i=Poles.Lower(); i<=Poles.Upper(); i++) {
1310 if (i >= FirstIndex && i <= LastIndex) {
1320 NewPoles(i) = Poles(i).Translated((Coef/(Dval+1.))*Displ);
1323 NewPoles(i) = Poles(i);
1328 //=======================================================================
1329 //function : MovePoint
1330 //purpose : Find the new poles which allows an old point (with a
1331 // given u as parameter) to reach a new position
1332 //=======================================================================
1334 //=======================================================================
1335 //function : MovePointAndTangent
1337 //=======================================================================
1339 void BSplCLib::MovePointAndTangent (const Standard_Real U,
1340 const Vector& Delta,
1341 const Vector& DeltaDerivatives,
1342 const Standard_Real Tolerance,
1343 const Standard_Integer Degree,
1344 const Standard_Boolean Rational,
1345 const Standard_Integer StartingCondition,
1346 const Standard_Integer EndingCondition,
1347 const Array1OfPoints& Poles,
1348 const TColStd_Array1OfReal& Weights,
1349 const TColStd_Array1OfReal& FlatKnots,
1350 Array1OfPoints& NewPoles,
1351 Standard_Integer & ErrorStatus)
1353 Standard_Real *delta_array,
1354 *delta_derivative_array,
1358 Standard_Integer num_poles ;
1359 num_poles = Poles.Length() ;
1361 if (NewPoles.Length() != num_poles) {
1362 Standard_ConstructionError::Raise();
1364 delta_array = (Standard_Real *) &Delta ;
1365 delta_derivative_array = (Standard_Real *)&DeltaDerivatives ;
1366 poles_array = (Standard_Real *) &Poles(Poles.Lower()) ;
1368 new_poles_array = (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1369 MovePointAndTangent(U,
1372 delta_derivative_array[0],
1385 //=======================================================================
1386 //function : Resolution
1388 //=======================================================================
1390 void BSplCLib::Resolution(const Array1OfPoints& Poles,
1391 const TColStd_Array1OfReal& Weights,
1392 const Standard_Integer NumPoles,
1393 const TColStd_Array1OfReal& FlatKnots,
1394 const Standard_Integer Degree,
1395 const Standard_Real Tolerance3D,
1396 Standard_Real& UTolerance)
1398 Standard_Real *PolesArray ;
1399 PolesArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1400 BSplCLib::Resolution(PolesArray[0],
1410 //=======================================================================
1411 //function : FunctionMultiply
1413 //=======================================================================
1415 void BSplCLib::FunctionMultiply
1416 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1417 const Standard_Integer BSplineDegree,
1418 const TColStd_Array1OfReal & BSplineFlatKnots,
1419 const Array1OfPoints & Poles,
1420 const TColStd_Array1OfReal & FlatKnots,
1421 const Standard_Integer NewDegree,
1422 Array1OfPoints & NewPoles,
1423 Standard_Integer & Status)
1425 Standard_Integer num_bspline_poles =
1426 BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1427 Standard_Integer num_new_poles =
1428 FlatKnots.Length() - NewDegree - 1 ;
1430 if (Poles.Length() != num_bspline_poles ||
1431 NewPoles.Length() != num_new_poles) {
1432 Standard_ConstructionError::Raise();
1434 Standard_Real * array_of_poles =
1435 (Standard_Real *) &Poles(Poles.Lower()) ;
1436 Standard_Real * array_of_new_poles =
1437 (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1438 BSplCLib::FunctionMultiply(FunctionPtr,
1445 array_of_new_poles[0],
1449 //=======================================================================
1450 //function : FunctionReparameterise
1452 //=======================================================================
1454 void BSplCLib::FunctionReparameterise
1455 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1456 const Standard_Integer BSplineDegree,
1457 const TColStd_Array1OfReal & BSplineFlatKnots,
1458 const Array1OfPoints & Poles,
1459 const TColStd_Array1OfReal & FlatKnots,
1460 const Standard_Integer NewDegree,
1461 Array1OfPoints & NewPoles,
1462 Standard_Integer & Status)
1464 Standard_Integer num_bspline_poles =
1465 BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1466 Standard_Integer num_new_poles =
1467 FlatKnots.Length() - NewDegree - 1 ;
1469 if (Poles.Length() != num_bspline_poles ||
1470 NewPoles.Length() != num_new_poles) {
1471 Standard_ConstructionError::Raise();
1473 Standard_Real * array_of_poles =
1474 (Standard_Real *) &Poles(Poles.Lower()) ;
1475 Standard_Real * array_of_new_poles =
1476 (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1477 BSplCLib::FunctionReparameterise(FunctionPtr,
1484 array_of_new_poles[0],