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 divided 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) {
746 const TColStd_Array1OfReal& refWeights = *WeightsArray;
748 WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
749 PLib::NoDerivativeEvalPolynomial(NewParameter,
756 Inverse = 1.0e0 / Inverse;
757 ModifyCoords (myPoint, *= Inverse);
761 //=======================================================================
763 //purpose : Evaluates the polynomial cache of the Bspline Curve
765 //=======================================================================
766 void 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,
771 const TColStd_Array1OfReal* WeightsArray,
776 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
778 // the SpanLenght is the normalizing factor so that the polynomial is between
780 Standard_Real LocalPDerivatives[Dimension_gen << 1] ;
781 // Standard_Real LocalWDerivatives[2], NewParameter, Inverse ;
782 Standard_Real LocalWDerivatives[2], NewParameter ;
785 PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
787 myPoint = (Standard_Real *) &aPoint ;
789 myVector = (Standard_Real *) &aVector ;
790 NewParameter = (Parameter - CacheParameter) / SpanLenght ;
791 PLib::EvalPolynomial(NewParameter,
796 LocalPDerivatives[0]) ;
798 // unormalize derivatives since those are computed normalized
801 ModifyCoords (LocalPDerivatives + Dimension_gen, /= SpanLenght);
803 if (WeightsArray != NULL) {
804 const TColStd_Array1OfReal& refWeights = *WeightsArray;
806 WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
807 PLib::EvalPolynomial(NewParameter,
812 LocalWDerivatives[0]) ;
814 // unormalize the result since the polynomial stored in the cache
815 // is normalized between 0 and 1
817 LocalWDerivatives[1] /= SpanLenght ;
819 PLib::RationalDerivatives(1,
821 LocalPDerivatives[0],
822 LocalWDerivatives[0],
823 LocalPDerivatives[0]) ;
826 CopyCoords (myPoint, LocalPDerivatives);
827 CopyCoords (myVector, LocalPDerivatives + Dimension_gen);
830 //=======================================================================
832 //purpose : Evaluates the polynomial cache of the Bspline Curve
834 //=======================================================================
835 void 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,
840 const TColStd_Array1OfReal* WeightsArray,
846 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
848 // the SpanLenght is the normalizing factor so that the polynomial is between
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,
864 LocalPDerivatives[0]) ;
866 // unormalize derivatives since those are computed normalized
868 Factor = 1.0e0/ SpanLenght ;
869 Index = Dimension_gen ;
870 EndIndex = Min (2, Degree) ;
872 for (ii = 1 ; ii <= EndIndex ; ii++) {
873 ModifyCoords (LocalPDerivatives + Index, *= Factor);
874 Factor /= SpanLenght ;
875 Index += Dimension_gen;
878 Index = (Degree + 1) * Dimension_gen;
879 for (ii = Degree ; ii < 2 ; ii++) {
880 NullifyCoords (LocalPDerivatives + Index);
881 Index += Dimension_gen;
884 if (WeightsArray != NULL) {
885 const TColStd_Array1OfReal& refWeights = *WeightsArray;
887 WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
889 PLib::EvalPolynomial(NewParameter,
894 LocalWDerivatives[0]) ;
896 for (ii = Degree + 1 ; ii <= 2 ; ii++) {
897 LocalWDerivatives[ii] = 0.0e0 ;
900 // unormalize the result since the polynomial stored in the cache
901 // is normalized between 0 and 1
903 Factor = 1.0e0 / SpanLenght ;
905 for (ii = 1 ; ii <= EndIndex ; ii++) {
906 LocalWDerivatives[ii] *= Factor ;
907 Factor /= SpanLenght ;
909 PLib::RationalDerivatives(2,
911 LocalPDerivatives[0],
912 LocalWDerivatives[0],
913 LocalPDerivatives[0]) ;
916 CopyCoords (myPoint, LocalPDerivatives);
917 CopyCoords (myVector1, LocalPDerivatives + Dimension_gen);
918 CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2);
921 //=======================================================================
923 //purpose : Evaluates the polynomial cache of the Bspline Curve
925 //=======================================================================
926 void 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,
931 const TColStd_Array1OfReal* WeightsArray,
938 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
940 // the SpanLenght is the normalizing factor so that the polynomial is between
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,
957 LocalPDerivatives[0]) ;
959 Index = (Degree + 1) * Dimension_gen;
960 for (ii = Degree ; ii < 3 ; ii++) {
961 NullifyCoords (LocalPDerivatives + Index);
962 Index += Dimension_gen;
966 // unormalize derivatives since those are computed normalized
968 Factor = 1.0e0 / SpanLenght ;
969 Index = Dimension_gen ;
970 EndIndex = Min(3,Degree) ;
972 for (ii = 1 ; ii <= EndIndex ; ii++) {
973 ModifyCoords (LocalPDerivatives + Index, *= Factor);
974 Factor /= SpanLenght;
975 Index += Dimension_gen;
978 if (WeightsArray != NULL) {
979 const TColStd_Array1OfReal& refWeights = *WeightsArray;
981 WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
983 PLib::EvalPolynomial(NewParameter,
988 LocalWDerivatives[0]) ;
990 // unormalize the result since the polynomial stored in the cache
991 // is normalized between 0 and 1
993 Factor = 1.0e0 / SpanLenght ;
995 for (ii = 1 ; ii <= EndIndex ; ii++) {
996 LocalWDerivatives[ii] *= Factor ;
997 Factor /= SpanLenght ;
1000 for (ii = (Degree + 1) ; ii <= 3 ; ii++) {
1001 LocalWDerivatives[ii] = 0.0e0 ;
1003 PLib::RationalDerivatives(3,
1005 LocalPDerivatives[0],
1006 LocalWDerivatives[0],
1007 LocalPDerivatives[0]) ;
1010 CopyCoords (myPoint, LocalPDerivatives);
1011 CopyCoords (myVector1, LocalPDerivatives + Dimension_gen);
1012 CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2);
1013 CopyCoords (myVector3, LocalPDerivatives + Dimension_gen * 3);
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 //=======================================================================
1023 void 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,
1030 const TColStd_Array1OfReal* Weights,
1031 Array1OfPoints& CachePoles,
1032 TColStd_Array1OfReal* CacheWeights)
1034 Standard_Integer ii,
1038 Standard_Real u = U,
1040 Standard_Boolean rational;
1042 BSplCLib_DataContainer dc(Degree);
1052 (BSplCLib::NoMults()),
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.
1060 BSplCLib::Bohm(u,Degree,Degree,*dc.knots,Dimension,*dc.poles);
1062 LocalValue = 1.0e0 ;
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 ;
1073 LocalIndex = Dimension_gen;
1074 LocalValue = 1.0e0 ;
1075 for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1076 (*CacheWeights)(ii) = dc.poles[LocalIndex] * LocalValue ;
1077 LocalIndex += Dimension_gen + 1;
1078 LocalValue *= SpanDomain / (Standard_Real) ii ;
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 ;
1089 if (Weights != NULL) {
1090 for (ii = 1 ; ii <= Degree + 1 ; ii++)
1091 (*CacheWeights)(ii) = 0.0e0 ;
1092 (*CacheWeights)(1) = 1.0e0 ;
1097 void BSplCLib::BuildCache(const Standard_Real theParameter,
1098 const Standard_Real theSpanDomain,
1099 const Standard_Boolean thePeriodicFlag,
1100 const Standard_Integer theDegree,
1101 const Standard_Integer theSpanIndex,
1102 const TColStd_Array1OfReal& theFlatKnots,
1103 const Array1OfPoints& thePoles,
1104 const TColStd_Array1OfReal* theWeights,
1105 TColStd_Array2OfReal& theCacheArray)
1107 Standard_Real aParam = theParameter;
1108 Standard_Integer anIndex = theSpanIndex;
1109 Standard_Integer aDimension;
1110 Standard_Boolean isRational;
1112 BSplCLib_DataContainer dc(theDegree);
1122 (BSplCLib::NoMults()),
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.
1130 Standard_Integer aCacheShift = // helps to store weights when PrepareEval did not found that the curve is locally polynomial
1131 (theWeights != NULL && !isRational) ? aDimension + 1 : aDimension;
1133 BSplCLib::Bohm(aParam, theDegree, theDegree, *dc.knots, aDimension, *dc.poles);
1135 Standard_Real aCoeff = 1.0;
1136 Standard_Real* aCache = (Standard_Real *) &(theCacheArray(theCacheArray.LowerRow(), theCacheArray.LowerCol()));
1137 Standard_Real* aPolyCoeffs = dc.poles;
1139 for (Standard_Integer i = 0; i <= theDegree; i++)
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)
1153 if (aCacheShift > aDimension)
1154 theCacheArray.SetValue(theCacheArray.LowerRow(), theCacheArray.LowerCol() + aCacheShift - 1, 1.0);
1157 //=======================================================================
1158 //function : Interpolate
1160 //=======================================================================
1162 void 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)
1170 Standard_Real *PArray ;
1172 PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1174 BSplCLib::Interpolate(Degree,
1183 //=======================================================================
1184 //function : Interpolate
1186 //=======================================================================
1188 void 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)
1196 Standard_Real *PArray,
1198 PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1199 WArray = (Standard_Real *) &Weights(Weights.Lower()) ;
1200 BSplCLib::Interpolate(Degree,
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 //=======================================================================
1216 void 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,
1221 const Array1OfPoints& Poles,
1222 const TColStd_Array1OfReal* Weights,
1223 const TColStd_Array1OfReal& FlatKnots,
1224 Standard_Integer& FirstIndex,
1225 Standard_Integer& LastIndex,
1226 Array1OfPoints& NewPoles)
1228 // calculate the BSplineBasis in the parameter U
1229 Standard_Integer FirstNonZeroBsplineIndex;
1230 math_Matrix BSplineBasis(1, 1,
1232 Standard_Integer ErrorCode =
1233 BSplCLib::EvalBsplineBasis(0,
1237 FirstNonZeroBsplineIndex,
1239 if (ErrorCode != 0) {
1243 for (Standard_Integer i=Poles.Lower(); i<=Poles.Upper(); i++) {
1244 NewPoles(i) = Poles(i);
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;
1255 Standard_Real maxValue = 0.0;
1256 Standard_Integer i, kk1=0, kk2, ii;
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);
1266 // find a kk2 if symetriy
1268 i = kk1 - FirstNonZeroBsplineIndex + 2;
1269 if ((kk1+1) <= LastIndex) {
1270 if (Abs(BSplineBasis(1, kk1-FirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
1275 // compute the vector of displacement
1276 Standard_Real D1 = 0.0;
1277 Standard_Real D2 = 0.0;
1278 Standard_Real hN, Coef, Dval;
1280 for (i = 1; i <= Degree+1; i++) {
1281 ii = i + FirstNonZeroBsplineIndex - 1;
1282 if (Weights != NULL) {
1283 hN = Weights->Value(ii)*BSplineBasis(1, i);
1287 hN = BSplineBasis(1, i);
1289 if (ii >= FirstIndex && ii <= LastIndex) {
1293 else if (ii > kk2) {
1299 D1 += 1./(Dval + 1.) * hN;
1303 if (Weights != NULL) {
1310 // compute the new poles
1312 for (i=Poles.Lower(); i<=Poles.Upper(); i++) {
1313 if (i >= FirstIndex && i <= LastIndex) {
1323 NewPoles(i) = Poles(i).Translated((Coef/(Dval+1.))*Displ);
1326 NewPoles(i) = Poles(i);
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 //=======================================================================
1337 //=======================================================================
1338 //function : MovePointAndTangent
1340 //=======================================================================
1342 void BSplCLib::MovePointAndTangent (const Standard_Real U,
1343 const Vector& Delta,
1344 const Vector& DeltaDerivatives,
1345 const Standard_Real Tolerance,
1346 const Standard_Integer Degree,
1347 const Standard_Integer StartingCondition,
1348 const Standard_Integer EndingCondition,
1349 const Array1OfPoints& Poles,
1350 const TColStd_Array1OfReal* Weights,
1351 const TColStd_Array1OfReal& FlatKnots,
1352 Array1OfPoints& NewPoles,
1353 Standard_Integer & ErrorStatus)
1355 Standard_Real *delta_array,
1356 *delta_derivative_array,
1360 Standard_Integer num_poles ;
1361 num_poles = Poles.Length() ;
1363 if (NewPoles.Length() != num_poles) {
1364 throw Standard_ConstructionError();
1366 delta_array = (Standard_Real *) &Delta ;
1367 delta_derivative_array = (Standard_Real *)&DeltaDerivatives ;
1368 poles_array = (Standard_Real *) &Poles(Poles.Lower()) ;
1370 new_poles_array = (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1371 MovePointAndTangent(U,
1374 delta_derivative_array[0],
1386 //=======================================================================
1387 //function : Resolution
1389 //=======================================================================
1391 void BSplCLib::Resolution(const Array1OfPoints& Poles,
1392 const TColStd_Array1OfReal* Weights,
1393 const Standard_Integer NumPoles,
1394 const TColStd_Array1OfReal& FlatKnots,
1395 const Standard_Integer Degree,
1396 const Standard_Real Tolerance3D,
1397 Standard_Real& UTolerance)
1399 Standard_Real *PolesArray ;
1400 PolesArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1401 BSplCLib::Resolution(PolesArray[0],
1411 //=======================================================================
1412 //function : FunctionMultiply
1414 //=======================================================================
1416 void 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,
1424 Standard_Integer & theStatus)
1426 Standard_Integer num_bspline_poles =
1427 BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1428 Standard_Integer num_new_poles =
1429 FlatKnots.Length() - NewDegree - 1 ;
1431 if (Poles.Length() != num_bspline_poles ||
1432 NewPoles.Length() != num_new_poles) {
1433 throw Standard_ConstructionError();
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,
1446 array_of_new_poles[0],
1450 //=======================================================================
1451 //function : FunctionReparameterise
1453 //=======================================================================
1455 void 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,
1463 Standard_Integer & theStatus)
1465 Standard_Integer num_bspline_poles =
1466 BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1467 Standard_Integer num_new_poles =
1468 FlatKnots.Length() - NewDegree - 1 ;
1470 if (Poles.Length() != num_bspline_poles ||
1471 NewPoles.Length() != num_new_poles) {
1472 throw Standard_ConstructionError();
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,
1485 array_of_new_poles[0],