1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
4 // The content of this file is subject to the Open CASCADE Technology Public
5 // License Version 6.5 (the "License"). You may not use the content of this file
6 // except in compliance with the License. Please obtain a copy of the License
7 // at http://www.opencascade.org and read it completely before using this file.
9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
12 // The Original Code and all software distributed under the License is
13 // distributed on an "AS IS" basis, without warranty of any kind, and the
14 // Initial Developer hereby disclaims all such warranties, including without
15 // limitation, any warranties of merchantability, fitness for a particular
16 // purpose or non-infringement. Please see the License for the specific terms
17 // and conditions governing the rights and limitations under the License.
19 // xab : modified 15-Mar-94 added EvalBSplineMatrix, FactorBandedMatrix, SolveBandedSystem
20 // Eval, BuildCache, cacheD0, cacheD1, cacheD2, cacheD3, LocalMatrix,
21 // EvalPolynomial, RationalDerivatives.
22 // xab : 22-Mar-94 : fixed rational problem in BuildCache
23 // xab : 12-Mar-96 : added MovePointAndTangent
24 #include <TColStd_Array1OfInteger.hxx>
25 #include <TColStd_Array1OfReal.hxx>
26 #include <gp_Vec2d.hxx>
27 #include <Standard_ConstructionError.hxx>
29 #include <math_Matrix.hxx>
31 //=======================================================================
32 //struct : BSplCLib_DataContainer
33 //purpose: Auxiliary structure providing buffers for poles and knots used in
34 // evaluation of bspline (allocated in the stack)
35 //=======================================================================
37 struct BSplCLib_DataContainer
39 BSplCLib_DataContainer(Standard_Integer Degree)
41 (void)Degree; // avoid compiler warning
42 Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree() ||
43 BSplCLib::MaxDegree() > 25,
44 "BSplCLib: bspline degree is greater than maximum supported");
47 Standard_Real poles[(25+1)*(Dimension_gen+1)];
48 Standard_Real knots[2*25];
49 Standard_Real ders[Dimension_gen*4];
52 //=======================================================================
55 //=======================================================================
57 void BSplCLib::Reverse(Array1OfPoints& Poles,
58 const Standard_Integer L)
60 Standard_Integer i, l = L;
61 l = Poles.Lower() + (l-Poles.Lower()) % (Poles.Upper()-Poles.Lower()+1);
63 Array1OfPoints temp(0,Poles.Length()-1);
65 for (i = Poles.Lower(); i <= l; i++)
68 for (i = l+1; i <= Poles.Upper(); i++)
69 temp(l-Poles.Lower()+Poles.Upper()-i+1) = Poles(i);
71 for (i = Poles.Lower(); i <= Poles.Upper(); i++)
72 Poles(i) = temp(i-Poles.Lower());
76 // CURVES MODIFICATIONS
79 //=======================================================================
80 //function : RemoveKnot
82 //=======================================================================
84 Standard_Boolean BSplCLib::RemoveKnot
85 (const Standard_Integer Index,
86 const Standard_Integer Mult,
87 const Standard_Integer Degree,
88 const Standard_Boolean Periodic,
89 const Array1OfPoints& Poles,
90 const TColStd_Array1OfReal& Weights,
91 const TColStd_Array1OfReal& Knots,
92 const TColStd_Array1OfInteger& Mults,
93 Array1OfPoints& NewPoles,
94 TColStd_Array1OfReal& NewWeights,
95 TColStd_Array1OfReal& NewKnots,
96 TColStd_Array1OfInteger& NewMults,
97 const Standard_Real Tolerance)
99 Standard_Boolean rational = &Weights != NULL;
100 Standard_Integer dim;
104 TColStd_Array1OfReal poles(1,dim*Poles.Length());
105 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
107 if (rational) PLib::SetPoles(Poles,Weights,poles);
108 else PLib::SetPoles(Poles,poles);
110 if (!RemoveKnot(Index,Mult,Degree,Periodic,dim,
111 poles,Knots,Mults,newpoles,NewKnots,NewMults,Tolerance))
112 return Standard_False;
114 if (rational) PLib::GetPoles(newpoles,NewPoles,NewWeights);
115 else PLib::GetPoles(newpoles,NewPoles);
116 return Standard_True;
119 //=======================================================================
120 //function : InsertKnots
121 //purpose : insert an array of knots and multiplicities
122 //=======================================================================
124 void BSplCLib::InsertKnots
125 (const Standard_Integer Degree,
126 const Standard_Boolean Periodic,
127 const Array1OfPoints& Poles,
128 const TColStd_Array1OfReal& Weights,
129 const TColStd_Array1OfReal& Knots,
130 const TColStd_Array1OfInteger& Mults,
131 const TColStd_Array1OfReal& AddKnots,
132 const TColStd_Array1OfInteger& AddMults,
133 Array1OfPoints& NewPoles,
134 TColStd_Array1OfReal& NewWeights,
135 TColStd_Array1OfReal& NewKnots,
136 TColStd_Array1OfInteger& NewMults,
137 const Standard_Real Epsilon,
138 const Standard_Boolean Add)
140 Standard_Boolean rational = &Weights != NULL;
141 Standard_Integer dim;
145 TColStd_Array1OfReal poles(1,dim*Poles.Length());
146 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
148 if (rational) PLib::SetPoles(Poles,Weights,poles);
149 else PLib::SetPoles(Poles,poles);
151 InsertKnots(Degree,Periodic,dim,poles,Knots,Mults,
152 AddKnots,AddMults,newpoles,NewKnots,NewMults,Epsilon,Add);
154 if (rational) PLib::GetPoles(newpoles,NewPoles,NewWeights);
155 else PLib::GetPoles(newpoles,NewPoles);
158 //=======================================================================
159 //function : InsertKnot
161 //=======================================================================
163 void BSplCLib::InsertKnot(const Standard_Integer ,
164 const Standard_Real U,
165 const Standard_Integer UMult,
166 const Standard_Integer Degree,
167 const Standard_Boolean Periodic,
168 const Array1OfPoints& Poles,
169 const TColStd_Array1OfReal& Weights,
170 const TColStd_Array1OfReal& Knots,
171 const TColStd_Array1OfInteger& Mults,
172 Array1OfPoints& NewPoles,
173 TColStd_Array1OfReal& NewWeights)
175 TColStd_Array1OfReal k(1,1);
177 TColStd_Array1OfInteger m(1,1);
179 TColStd_Array1OfReal nk(1,Knots.Length()+1);
180 TColStd_Array1OfInteger nm(1,Knots.Length()+1);
181 InsertKnots(Degree,Periodic,Poles,Weights,Knots,Mults,
182 k,m,NewPoles,NewWeights,nk,nm,Epsilon(U));
185 //=======================================================================
186 //function : RaiseMultiplicity
188 //=======================================================================
190 void BSplCLib::RaiseMultiplicity(const Standard_Integer KnotIndex,
191 const Standard_Integer Mult,
192 const Standard_Integer Degree,
193 const Standard_Boolean Periodic,
194 const Array1OfPoints& Poles,
195 const TColStd_Array1OfReal& Weights,
196 const TColStd_Array1OfReal& Knots,
197 const TColStd_Array1OfInteger& Mults,
198 Array1OfPoints& NewPoles,
199 TColStd_Array1OfReal& NewWeights)
201 TColStd_Array1OfReal k(1,1);
202 k(1) = Knots(KnotIndex);
203 TColStd_Array1OfInteger m(1,1);
204 m(1) = Mult - Mults(KnotIndex);
205 TColStd_Array1OfReal nk(1,Knots.Length());
206 TColStd_Array1OfInteger nm(1,Knots.Length());
207 InsertKnots(Degree,Periodic,Poles,Weights,Knots,Mults,
208 k,m,NewPoles,NewWeights,nk,nm,Epsilon(k(1)));
211 //=======================================================================
212 //function : IncreaseDegree
214 //=======================================================================
216 void BSplCLib::IncreaseDegree
217 (const Standard_Integer Degree,
218 const Standard_Integer NewDegree,
219 const Standard_Boolean Periodic,
220 const Array1OfPoints& Poles,
221 const TColStd_Array1OfReal& Weights,
222 const TColStd_Array1OfReal& Knots,
223 const TColStd_Array1OfInteger& Mults,
224 Array1OfPoints& NewPoles,
225 TColStd_Array1OfReal& NewWeights,
226 TColStd_Array1OfReal& NewKnots,
227 TColStd_Array1OfInteger& NewMults)
229 Standard_Boolean rational = &Weights != NULL;
230 Standard_Integer dim;
234 TColStd_Array1OfReal poles(1,dim*Poles.Length());
235 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
237 if (rational) PLib::SetPoles(Poles,Weights,poles);
238 else PLib::SetPoles(Poles,poles);
240 IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults,
241 newpoles,NewKnots,NewMults);
243 if (rational) PLib::GetPoles(newpoles,NewPoles,NewWeights);
244 else PLib::GetPoles(newpoles,NewPoles);
247 //=======================================================================
248 //function : Unperiodize
250 //=======================================================================
252 void BSplCLib::Unperiodize
253 (const Standard_Integer Degree,
254 const TColStd_Array1OfInteger& Mults,
255 const TColStd_Array1OfReal& Knots,
256 const Array1OfPoints& Poles,
257 const TColStd_Array1OfReal& Weights,
258 TColStd_Array1OfInteger& NewMults,
259 TColStd_Array1OfReal& NewKnots,
260 Array1OfPoints& NewPoles,
261 TColStd_Array1OfReal& NewWeights)
263 Standard_Boolean rational = &Weights != NULL;
264 Standard_Integer dim;
268 TColStd_Array1OfReal poles(1,dim*Poles.Length());
269 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
271 if (rational) PLib::SetPoles(Poles,Weights,poles);
272 else PLib::SetPoles(Poles,poles);
274 Unperiodize(Degree, dim, Mults, Knots, poles,
275 NewMults,NewKnots, newpoles);
277 if (rational) PLib::GetPoles(newpoles,NewPoles,NewWeights);
278 else PLib::GetPoles(newpoles,NewPoles);
281 //=======================================================================
282 //function : Trimming
284 //=======================================================================
286 void BSplCLib::Trimming(const Standard_Integer Degree,
287 const Standard_Boolean Periodic,
288 const TColStd_Array1OfReal& Knots,
289 const TColStd_Array1OfInteger& Mults,
290 const Array1OfPoints& Poles,
291 const TColStd_Array1OfReal& Weights,
292 const Standard_Real U1,
293 const Standard_Real U2,
294 TColStd_Array1OfReal& NewKnots,
295 TColStd_Array1OfInteger& NewMults,
296 Array1OfPoints& NewPoles,
297 TColStd_Array1OfReal& NewWeights)
299 Standard_Boolean rational = &Weights != NULL;
300 Standard_Integer dim;
304 TColStd_Array1OfReal poles(1,dim*Poles.Length());
305 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
307 if (rational) PLib::SetPoles(Poles,Weights,poles);
308 else PLib::SetPoles(Poles,poles);
310 Trimming(Degree, Periodic, dim, Knots, Mults, poles, U1, U2,
311 NewKnots, NewMults, newpoles);
313 if (rational) PLib::GetPoles(newpoles,NewPoles,NewWeights);
314 else PLib::GetPoles(newpoles,NewPoles);
317 //--------------------------------------------------------------------------
318 //ELEMENTARY COMPUTATIONS
319 //--------------------------------------------------------------------------
321 // All the following methods are methods of low level used in BSplCLib
322 // but also in BSplSLib (but not in Geom)
323 // At the creation time of this package there is no possibility
324 // to declare private methods of package and to declare friend
325 // methods of package. It could interesting to declare that BSplSLib
326 // is a package friend to BSplCLib because it uses all the basis methods
328 //--------------------------------------------------------------------------
330 //=======================================================================
331 //function : BuildEval
332 //purpose : builds the local array for evaluation
333 //=======================================================================
335 void BSplCLib::BuildEval(const Standard_Integer Degree,
336 const Standard_Integer Index,
337 const Array1OfPoints& Poles,
338 const TColStd_Array1OfReal& Weights,
341 Standard_Real w, *pole = &LP;
342 Standard_Integer PLower = Poles.Lower();
343 Standard_Integer PUpper = Poles.Upper();
345 Standard_Integer ip = PLower + Index - 1;
346 if (&Weights == NULL) {
347 for (i = 0; i <= Degree; i++) {
349 if (ip > PUpper) ip = PLower;
350 const Point& P = Poles(ip);
351 PointToCoords (pole,P,+0);
352 pole += Dimension_gen;
356 for (i = 0; i <= Degree; i++) {
358 if (ip > PUpper) ip = PLower;
359 const Point& P = Poles(ip);
360 pole[Dimension_gen] = w = Weights(ip);
361 PointToCoords (pole, P, * w);
362 pole += Dimension_gen + 1;
367 //=======================================================================
368 //function : PrepareEval
369 //purpose : stores data for Eval in the local arrays
370 // dc.poles and dc.knots
371 //=======================================================================
373 static void PrepareEval
375 Standard_Integer& index,
376 Standard_Integer& dim,
377 Standard_Boolean& rational,
378 const Standard_Integer Degree,
379 const Standard_Boolean Periodic,
380 const Array1OfPoints& Poles,
381 const TColStd_Array1OfReal& Weights,
382 const TColStd_Array1OfReal& Knots,
383 const TColStd_Array1OfInteger& Mults,
384 BSplCLib_DataContainer& dc)
387 BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
390 BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*dc.knots);
392 index -= Knots.Lower() + Degree;
394 index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults);
396 // check truly rational
397 rational = (&Weights != NULL);
399 Standard_Integer WLower = Weights.Lower() + index;
400 rational = BSplCLib::IsRational(Weights, WLower, WLower + Degree);
405 dim = Dimension_gen + 1;
406 BSplCLib::BuildEval(Degree, index, Poles, Weights , *dc.poles);
410 BSplCLib::BuildEval(Degree, index, Poles, BSplCLib::NoWeights(), *dc.poles);
414 //=======================================================================
417 //=======================================================================
420 (const Standard_Real U,
421 const Standard_Integer Index,
422 const Standard_Integer Degree,
423 const Standard_Boolean Periodic,
424 const Array1OfPoints& Poles,
425 const TColStd_Array1OfReal& Weights,
426 const TColStd_Array1OfReal& Knots,
427 const TColStd_Array1OfInteger& Mults,
430 // Standard_Integer k,dim,index = Index;
431 Standard_Integer dim,index = Index;
433 Standard_Boolean rational;
434 BSplCLib_DataContainer dc(Degree);
435 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
436 BSplCLib::Eval(u,Degree,*dc.knots,dim,*dc.poles);
439 Standard_Real w = dc.poles[Dimension_gen];
440 CoordsToPoint (P, dc.poles, / w);
443 CoordsToPoint (P, dc.poles, );
446 //=======================================================================
449 //=======================================================================
452 (const Standard_Real U,
453 const Standard_Integer Index,
454 const Standard_Integer Degree,
455 const Standard_Boolean Periodic,
456 const Array1OfPoints& Poles,
457 const TColStd_Array1OfReal& Weights,
458 const TColStd_Array1OfReal& Knots,
459 const TColStd_Array1OfInteger& Mults,
463 Standard_Integer dim,index = Index;
465 Standard_Boolean rational;
466 BSplCLib_DataContainer dc(Degree);
467 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
468 BSplCLib::Bohm(u,Degree,1,*dc.knots,dim,*dc.poles);
469 Standard_Real *result = dc.poles;
471 PLib::RationalDerivative(Degree,1,Dimension_gen,*dc.poles,*dc.ders);
475 CoordsToPoint (P, result, );
476 CoordsToPoint (V, result + Dimension_gen, );
479 //=======================================================================
482 //=======================================================================
485 (const Standard_Real U,
486 const Standard_Integer Index,
487 const Standard_Integer Degree,
488 const Standard_Boolean Periodic,
489 const Array1OfPoints& Poles,
490 const TColStd_Array1OfReal& Weights,
491 const TColStd_Array1OfReal& Knots,
492 const TColStd_Array1OfInteger& Mults,
497 Standard_Integer dim,index = Index;
499 Standard_Boolean rational;
500 BSplCLib_DataContainer dc(Degree);
501 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
502 BSplCLib::Bohm(u,Degree,2,*dc.knots,dim,*dc.poles);
503 Standard_Real *result = dc.poles;
505 PLib::RationalDerivative(Degree,2,Dimension_gen,*dc.poles,*dc.ders);
509 CoordsToPoint (P, result, );
510 CoordsToPoint (V1, result + Dimension_gen, );
511 if (!rational && (Degree < 2))
514 CoordsToPoint (V2, result + 2 * Dimension_gen, );
517 //=======================================================================
520 //=======================================================================
523 (const Standard_Real U,
524 const Standard_Integer Index,
525 const Standard_Integer Degree,
526 const Standard_Boolean Periodic,
527 const Array1OfPoints& Poles,
528 const TColStd_Array1OfReal& Weights,
529 const TColStd_Array1OfReal& Knots,
530 const TColStd_Array1OfInteger& Mults,
536 Standard_Integer dim,index = Index;
538 Standard_Boolean rational;
539 BSplCLib_DataContainer dc(Degree);
540 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
541 BSplCLib::Bohm(u,Degree,3,*dc.knots,dim,*dc.poles);
542 Standard_Real *result = dc.poles;
544 PLib::RationalDerivative(Degree,3,Dimension_gen,*dc.poles,*dc.ders);
548 CoordsToPoint (P, result, );
549 CoordsToPoint (V1, result + Dimension_gen, );
550 if (!rational && (Degree < 2))
553 CoordsToPoint (V2, result + 2 * Dimension_gen, );
554 if (!rational && (Degree < 3))
557 CoordsToPoint (V3, result + 3 * Dimension_gen, );
560 //=======================================================================
563 //=======================================================================
566 (const Standard_Real U,
567 const Standard_Integer N,
568 const Standard_Integer Index,
569 const Standard_Integer Degree,
570 const Standard_Boolean Periodic,
571 const Array1OfPoints& Poles,
572 const TColStd_Array1OfReal& Weights,
573 const TColStd_Array1OfReal& Knots,
574 const TColStd_Array1OfInteger& Mults,
577 Standard_Integer dim,index = Index;
579 Standard_Boolean rational;
580 BSplCLib_DataContainer dc(Degree);
581 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
582 BSplCLib::Bohm(u,Degree,N,*dc.knots,dim,*dc.poles);
585 Standard_Real v[Dimension_gen];
586 PLib::RationalDerivative(Degree,N,Dimension_gen,*dc.poles,v[0],Standard_False);
587 CoordsToPoint (VN, v, );
593 Standard_Real *DN = dc.poles + N * Dimension_gen;
594 CoordsToPoint (VN, DN, );
599 //=======================================================================
600 //function : Solves a LU factored Matrix
602 //=======================================================================
605 BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
606 const Standard_Integer UpperBandWidth,
607 const Standard_Integer LowerBandWidth,
608 Array1OfPoints& PolesArray)
610 Standard_Real *PArray ;
611 PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
613 return BSplCLib::SolveBandedSystem(Matrix,
620 //=======================================================================
621 //function : Solves a LU factored Matrix
622 //purpose : if HomogeneousFlag is 1 then the input and the output
623 // will be homogeneous that is no division or multiplication
624 // by weigths will happen. On the contrary if HomogeneousFlag
625 // is 0 then the poles will be multiplied first by the weights
626 // and after interpolation they will be devided by the weights
627 //=======================================================================
630 BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
631 const Standard_Integer UpperBandWidth,
632 const Standard_Integer LowerBandWidth,
633 const Standard_Boolean HomogeneousFlag,
634 Array1OfPoints& PolesArray,
635 TColStd_Array1OfReal& WeightsArray)
637 Standard_Real *PArray,
639 PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
640 WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
641 return BSplCLib::SolveBandedSystem(Matrix,
650 //=======================================================================
651 //function : Evaluates a Bspline function : uses the ExtrapMode
652 //purpose : the function is extrapolated using the Taylor expansion
653 // of degree ExtrapMode[0] to the left and the Taylor
654 // expansion of degree ExtrapMode[1] to the right
655 // if the HomogeneousFlag == True than the Poles are supposed
656 // to be stored homogeneously and the result will also be homogeneous
657 // Valid only if Weights
658 //=======================================================================
659 void BSplCLib::Eval(const Standard_Real Parameter,
660 const Standard_Boolean PeriodicFlag,
661 const Standard_Boolean HomogeneousFlag,
662 Standard_Integer& ExtrapMode,
663 const Standard_Integer Degree,
664 const TColStd_Array1OfReal& FlatKnots,
665 const Array1OfPoints& PolesArray,
666 const TColStd_Array1OfReal& WeightsArray,
668 Standard_Real& aWeight)
670 Standard_Real Inverse,
674 Standard_Integer kk ;
675 PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
676 WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
677 if (HomogeneousFlag) {
678 BSplCLib::Eval(Parameter,
687 BSplCLib::Eval(Parameter,
698 BSplCLib::Eval(Parameter,
709 Inverse = 1.0e0 / aWeight ;
711 for (kk = 0 ; kk < Dimension_gen ; kk++) {
716 for (kk = 0 ; kk < Dimension_gen ; kk++)
717 aPoint.SetCoord(kk + 1, P[kk]) ;
720 //=======================================================================
722 //purpose : Evaluates the polynomial cache of the Bspline Curve
724 //=======================================================================
725 void BSplCLib::CacheD0(const Standard_Real Parameter,
726 const Standard_Integer Degree,
727 const Standard_Real CacheParameter,
728 const Standard_Real SpanLenght,
729 const Array1OfPoints& PolesArray,
730 const TColStd_Array1OfReal& WeightsArray,
734 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
736 // the SpanLenght is the normalizing factor so that the polynomial is between
738 Standard_Real NewParameter, Inverse;
739 Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
740 Standard_Real * myPoint = (Standard_Real *) &aPoint ;
741 NewParameter = (Parameter - CacheParameter) / SpanLenght ;
742 PLib::NoDerivativeEvalPolynomial(NewParameter,
745 Degree * Dimension_gen,
748 if (&WeightsArray != NULL) {
750 WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
751 PLib::NoDerivativeEvalPolynomial(NewParameter,
758 Inverse = 1.0e0 / Inverse;
759 ModifyCoords (myPoint, *= Inverse);
763 //=======================================================================
765 //purpose : Evaluates the polynomial cache of the Bspline Curve
767 //=======================================================================
768 void BSplCLib::CacheD1(const Standard_Real Parameter,
769 const Standard_Integer Degree,
770 const Standard_Real CacheParameter,
771 const Standard_Real SpanLenght,
772 const Array1OfPoints& PolesArray,
773 const TColStd_Array1OfReal& WeightsArray,
778 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
780 // the SpanLenght is the normalizing factor so that the polynomial is between
782 Standard_Real LocalPDerivatives[Dimension_gen << 1] ;
783 // Standard_Real LocalWDerivatives[2], NewParameter, Inverse ;
784 Standard_Real LocalWDerivatives[2], NewParameter ;
787 PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
789 myPoint = (Standard_Real *) &aPoint ;
791 myVector = (Standard_Real *) &aVector ;
792 NewParameter = (Parameter - CacheParameter) / SpanLenght ;
793 PLib::EvalPolynomial(NewParameter,
798 LocalPDerivatives[0]) ;
800 // unormalize derivatives since those are computed normalized
803 ModifyCoords (LocalPDerivatives + Dimension_gen, /= SpanLenght);
805 if (&WeightsArray != NULL) {
807 WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
808 PLib::EvalPolynomial(NewParameter,
813 LocalWDerivatives[0]) ;
815 // unormalize the result since the polynomial stored in the cache
816 // is normalized between 0 and 1
818 LocalWDerivatives[1] /= SpanLenght ;
820 PLib::RationalDerivatives(1,
822 LocalPDerivatives[0],
823 LocalWDerivatives[0],
824 LocalPDerivatives[0]) ;
827 CopyCoords (myPoint, LocalPDerivatives);
828 CopyCoords (myVector, LocalPDerivatives + Dimension_gen);
831 //=======================================================================
833 //purpose : Evaluates the polynomial cache of the Bspline Curve
835 //=======================================================================
836 void BSplCLib::CacheD2(const Standard_Real Parameter,
837 const Standard_Integer Degree,
838 const Standard_Real CacheParameter,
839 const Standard_Real SpanLenght,
840 const Array1OfPoints& PolesArray,
841 const TColStd_Array1OfReal& WeightsArray,
847 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
849 // the SpanLenght is the normalizing factor so that the polynomial is between
851 Standard_Integer ii,Index,EndIndex;
852 Standard_Real LocalPDerivatives[(Dimension_gen << 1) + Dimension_gen] ;
853 // Standard_Real LocalWDerivatives[3], NewParameter, Factor, Inverse ;
854 Standard_Real LocalWDerivatives[3], NewParameter, Factor ;
855 Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
856 Standard_Real * myPoint = (Standard_Real *) &aPoint ;
857 Standard_Real * myVector1 = (Standard_Real *) &aVector1 ;
858 Standard_Real * myVector2 = (Standard_Real *) &aVector2 ;
859 NewParameter = (Parameter - CacheParameter) / SpanLenght ;
860 PLib::EvalPolynomial(NewParameter,
865 LocalPDerivatives[0]) ;
867 // unormalize derivatives since those are computed normalized
869 Factor = 1.0e0/ SpanLenght ;
870 Index = Dimension_gen ;
871 EndIndex = Min (2, Degree) ;
873 for (ii = 1 ; ii <= EndIndex ; ii++) {
874 ModifyCoords (LocalPDerivatives + Index, *= Factor);
875 Factor /= SpanLenght ;
876 Index += Dimension_gen;
879 Index = (Degree + 1) * Dimension_gen;
880 for (ii = Degree ; ii < 2 ; ii++) {
881 NullifyCoords (LocalPDerivatives + Index);
882 Index += Dimension_gen;
885 if (&WeightsArray != NULL) {
887 WArray = (Standard_Real *) &WeightsArray(WeightsArray.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) {
980 WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
982 PLib::EvalPolynomial(NewParameter,
987 LocalWDerivatives[0]) ;
989 // unormalize the result since the polynomial stored in the cache
990 // is normalized between 0 and 1
992 Factor = 1.0e0 / SpanLenght ;
994 for (ii = 1 ; ii <= EndIndex ; ii++) {
995 LocalWDerivatives[ii] *= Factor ;
996 Factor /= SpanLenght ;
999 for (ii = (Degree + 1) ; ii <= 3 ; ii++) {
1000 LocalWDerivatives[ii] = 0.0e0 ;
1002 PLib::RationalDerivatives(3,
1004 LocalPDerivatives[0],
1005 LocalWDerivatives[0],
1006 LocalPDerivatives[0]) ;
1009 CopyCoords (myPoint, LocalPDerivatives);
1010 CopyCoords (myVector1, LocalPDerivatives + Dimension_gen);
1011 CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2);
1012 CopyCoords (myVector3, LocalPDerivatives + Dimension_gen * 3);
1015 //=======================================================================
1016 //function : BuildCache
1017 //purpose : Stores theTaylor expansion normalized between 0,1 in the
1018 // Cache : in case of a rational function the Poles are
1019 // stored in homogeneous form
1020 //=======================================================================
1022 void BSplCLib::BuildCache
1023 (const Standard_Real U,
1024 const Standard_Real SpanDomain,
1025 const Standard_Boolean Periodic,
1026 const Standard_Integer Degree,
1027 const TColStd_Array1OfReal& FlatKnots,
1028 const Array1OfPoints& Poles,
1029 const TColStd_Array1OfReal& Weights,
1030 Array1OfPoints& CachePoles,
1031 TColStd_Array1OfReal& CacheWeights)
1033 Standard_Integer ii,
1037 Standard_Real u = U,
1039 Standard_Boolean rational;
1041 BSplCLib_DataContainer dc(Degree);
1051 (BSplCLib::NoMults()),
1054 // Watch out : PrepareEval checks if locally the Bspline is polynomial
1055 // therefore rational flag could be set to False even if there are
1056 // Weights. The Dimension is set accordingly.
1059 BSplCLib::Bohm(u,Degree,Degree,*dc.knots,Dimension,*dc.poles);
1061 LocalValue = 1.0e0 ;
1066 for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1067 CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue);
1068 LocalIndex += Dimension_gen + 1;
1069 LocalValue *= SpanDomain / (Standard_Real) ii ;
1072 LocalIndex = Dimension_gen;
1073 LocalValue = 1.0e0 ;
1074 for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1075 CacheWeights(ii) = dc.poles[LocalIndex] * LocalValue ;
1076 LocalIndex += Dimension_gen + 1;
1077 LocalValue *= SpanDomain / (Standard_Real) ii ;
1082 for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1083 CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue);
1084 LocalIndex += Dimension_gen;
1085 LocalValue *= SpanDomain / (Standard_Real) ii ;
1088 if (&Weights != NULL) {
1089 for (ii = 1 ; ii <= Degree + 1 ; ii++)
1090 CacheWeights(ii) = 0.0e0 ;
1091 CacheWeights(1) = 1.0e0 ;
1096 //=======================================================================
1097 //function : Interpolate
1099 //=======================================================================
1101 void BSplCLib::Interpolate(const Standard_Integer Degree,
1102 const TColStd_Array1OfReal& FlatKnots,
1103 const TColStd_Array1OfReal& Parameters,
1104 const TColStd_Array1OfInteger& ContactOrderArray,
1105 Array1OfPoints& Poles,
1106 Standard_Integer& InversionProblem)
1109 Standard_Real *PArray ;
1111 PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1113 BSplCLib::Interpolate(Degree,
1122 //=======================================================================
1123 //function : Interpolate
1125 //=======================================================================
1127 void BSplCLib::Interpolate(const Standard_Integer Degree,
1128 const TColStd_Array1OfReal& FlatKnots,
1129 const TColStd_Array1OfReal& Parameters,
1130 const TColStd_Array1OfInteger& ContactOrderArray,
1131 Array1OfPoints& Poles,
1132 TColStd_Array1OfReal& Weights,
1133 Standard_Integer& InversionProblem)
1135 Standard_Real *PArray,
1137 PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1138 WArray = (Standard_Real *) &Weights(Weights.Lower()) ;
1139 BSplCLib::Interpolate(Degree,
1149 //=======================================================================
1150 //function : MovePoint
1151 //purpose : Find the new poles which allows an old point (with a
1152 // given u as parameter) to reach a new position
1153 //=======================================================================
1155 void BSplCLib::MovePoint (const Standard_Real U,
1156 const Vector& Displ,
1157 const Standard_Integer Index1,
1158 const Standard_Integer Index2,
1159 const Standard_Integer Degree,
1160 const Standard_Boolean Rational,
1161 const Array1OfPoints& Poles,
1162 const TColStd_Array1OfReal& Weights,
1163 const TColStd_Array1OfReal& FlatKnots,
1164 Standard_Integer& FirstIndex,
1165 Standard_Integer& LastIndex,
1166 Array1OfPoints& NewPoles)
1168 // calculate the BSplineBasis in the parameter U
1169 Standard_Integer FirstNonZeroBsplineIndex;
1170 math_Matrix BSplineBasis(1, 1,
1172 Standard_Integer ErrorCode =
1173 BSplCLib::EvalBsplineBasis(1,
1178 FirstNonZeroBsplineIndex,
1180 if (ErrorCode != 0) {
1184 for (Standard_Integer i=Poles.Lower(); i<=Poles.Upper(); i++) {
1185 NewPoles(i) = Poles(i);
1190 // find the span which is predominant for parameter U
1191 FirstIndex = FirstNonZeroBsplineIndex;
1192 LastIndex = FirstNonZeroBsplineIndex + Degree ;
1193 if (FirstIndex < Index1) FirstIndex = Index1;
1194 if (LastIndex > Index2) LastIndex = Index2;
1196 Standard_Real maxValue = 0.0;
1197 Standard_Integer i, kk1=0, kk2, ii;
1199 for (i = FirstIndex-FirstNonZeroBsplineIndex+1;
1200 i <= LastIndex-FirstNonZeroBsplineIndex+1; i++) {
1201 if (BSplineBasis(1,i) > maxValue) {
1202 kk1 = i + FirstNonZeroBsplineIndex - 1;
1203 maxValue = BSplineBasis(1,i);
1207 // find a kk2 if symetriy
1209 i = kk1 - FirstNonZeroBsplineIndex + 2;
1210 if ((kk1+1) <= LastIndex) {
1211 if (Abs(BSplineBasis(1, kk1-FirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
1216 // compute the vector of displacement
1217 Standard_Real D1 = 0.0;
1218 Standard_Real D2 = 0.0;
1219 Standard_Real hN, Coef, Dval;
1221 for (i = 1; i <= Degree+1; i++) {
1222 ii = i + FirstNonZeroBsplineIndex - 1;
1224 hN = Weights(ii)*BSplineBasis(1, i);
1228 hN = BSplineBasis(1, i);
1230 if (ii >= FirstIndex && ii <= LastIndex) {
1234 else if (ii > kk2) {
1240 D1 += 1./(Dval + 1.) * hN;
1251 // compute the new poles
1253 for (i=Poles.Lower(); i<=Poles.Upper(); i++) {
1254 if (i >= FirstIndex && i <= LastIndex) {
1264 NewPoles(i) = Poles(i).Translated((Coef/(Dval+1.))*Displ);
1267 NewPoles(i) = Poles(i);
1272 //=======================================================================
1273 //function : MovePoint
1274 //purpose : Find the new poles which allows an old point (with a
1275 // given u as parameter) to reach a new position
1276 //=======================================================================
1278 //=======================================================================
1279 //function : MovePointAndTangent
1281 //=======================================================================
1283 void BSplCLib::MovePointAndTangent (const Standard_Real U,
1284 const Vector& Delta,
1285 const Vector& DeltaDerivatives,
1286 const Standard_Real Tolerance,
1287 const Standard_Integer Degree,
1288 const Standard_Boolean Rational,
1289 const Standard_Integer StartingCondition,
1290 const Standard_Integer EndingCondition,
1291 const Array1OfPoints& Poles,
1292 const TColStd_Array1OfReal& Weights,
1293 const TColStd_Array1OfReal& FlatKnots,
1294 Array1OfPoints& NewPoles,
1295 Standard_Integer & ErrorStatus)
1297 Standard_Real *delta_array,
1298 *delta_derivative_array,
1302 Standard_Integer num_poles ;
1303 num_poles = Poles.Length() ;
1305 if (NewPoles.Length() != num_poles) {
1306 Standard_ConstructionError::Raise();
1308 delta_array = (Standard_Real *) &Delta ;
1309 delta_derivative_array = (Standard_Real *)&DeltaDerivatives ;
1310 poles_array = (Standard_Real *) &Poles(Poles.Lower()) ;
1312 new_poles_array = (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1313 MovePointAndTangent(U,
1316 delta_derivative_array[0],
1329 //=======================================================================
1330 //function : Resolution
1332 //=======================================================================
1334 void BSplCLib::Resolution(const Array1OfPoints& Poles,
1335 const TColStd_Array1OfReal& Weights,
1336 const Standard_Integer NumPoles,
1337 const TColStd_Array1OfReal& FlatKnots,
1338 const Standard_Integer Degree,
1339 const Standard_Real Tolerance3D,
1340 Standard_Real& UTolerance)
1342 Standard_Real *PolesArray ;
1343 PolesArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1344 BSplCLib::Resolution(PolesArray[0],
1354 //=======================================================================
1355 //function : FunctionMultiply
1357 //=======================================================================
1359 void BSplCLib::FunctionMultiply
1360 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1361 const Standard_Integer BSplineDegree,
1362 const TColStd_Array1OfReal & BSplineFlatKnots,
1363 const Array1OfPoints & Poles,
1364 const TColStd_Array1OfReal & FlatKnots,
1365 const Standard_Integer NewDegree,
1366 Array1OfPoints & NewPoles,
1367 Standard_Integer & Status)
1369 Standard_Integer num_bspline_poles =
1370 BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1371 Standard_Integer num_new_poles =
1372 FlatKnots.Length() - NewDegree - 1 ;
1374 if (Poles.Length() != num_bspline_poles ||
1375 NewPoles.Length() != num_new_poles) {
1376 Standard_ConstructionError::Raise();
1378 Standard_Real * array_of_poles =
1379 (Standard_Real *) &Poles(Poles.Lower()) ;
1380 Standard_Real * array_of_new_poles =
1381 (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1382 BSplCLib::FunctionMultiply(FunctionPtr,
1389 array_of_new_poles[0],
1393 //=======================================================================
1394 //function : FunctionReparameterise
1396 //=======================================================================
1398 void BSplCLib::FunctionReparameterise
1399 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1400 const Standard_Integer BSplineDegree,
1401 const TColStd_Array1OfReal & BSplineFlatKnots,
1402 const Array1OfPoints & Poles,
1403 const TColStd_Array1OfReal & FlatKnots,
1404 const Standard_Integer NewDegree,
1405 Array1OfPoints & NewPoles,
1406 Standard_Integer & Status)
1408 Standard_Integer num_bspline_poles =
1409 BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1410 Standard_Integer num_new_poles =
1411 FlatKnots.Length() - NewDegree - 1 ;
1413 if (Poles.Length() != num_bspline_poles ||
1414 NewPoles.Length() != num_new_poles) {
1415 Standard_ConstructionError::Raise();
1417 Standard_Real * array_of_poles =
1418 (Standard_Real *) &Poles(Poles.Lower()) ;
1419 Standard_Real * array_of_new_poles =
1420 (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1421 BSplCLib::FunctionReparameterise(FunctionPtr,
1428 array_of_new_poles[0],