1 // Created on: 1995-01-17
2 // Created by: Laurent BOURESCHE
3 // Copyright (c) 1995-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
21 // xab : modified 15-Mar-95 : added BuildBSplMatrix,FactorBandedMatrix,
23 // EvalPolynomial : Horners method
25 #include <Standard_Stream.hxx>
27 #include <BSplCLib.hxx>
28 #include <gp_Mat2d.hxx>
30 #include <Standard_ConstructionError.hxx>
31 #include <TColStd_Array1OfReal.hxx>
32 #include <TColStd_Array1OfInteger.hxx>
33 #include <TColStd_HArray1OfReal.hxx>
34 #include <TColStd_HArray1OfInteger.hxx>
36 #include <math_Matrix.hxx>
38 //=======================================================================
39 //struct : BSplCLib_DataContainer
40 //purpose: Auxiliary structure providing buffers for poles and knots used in
41 // evaluation of bspline (allocated in the stack)
42 //=======================================================================
44 struct BSplCLib_DataContainer
46 BSplCLib_DataContainer(Standard_Integer Degree)
48 (void)Degree; // avoid compiler warning
49 Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree() ||
50 BSplCLib::MaxDegree() > 25,
51 "BSplCLib: bspline degree is greater than maximum supported");
54 Standard_Real poles[2*(25+1)];
55 Standard_Real knots[2*25];
56 Standard_Real ders[4];
59 // methods for 1 dimensional BSplines
61 //=======================================================================
62 //function : BuildEval
63 //purpose : builds the local array for evaluation
64 //=======================================================================
66 void BSplCLib::BuildEval(const Standard_Integer Degree,
67 const Standard_Integer Index,
68 const TColStd_Array1OfReal& Poles,
69 const TColStd_Array1OfReal& Weights,
72 Standard_Integer PLower = Poles.Lower();
73 Standard_Integer PUpper = Poles.Upper();
75 Standard_Integer ip = PLower + Index - 1;
76 Standard_Real w, *pole = &LP;
77 if (&Weights == NULL) {
79 for (i = 0; i <= Degree; i++) {
81 if (ip > PUpper) ip = PLower;
88 for (i = 0; i <= Degree; i++) {
90 if (ip > PUpper) ip = PLower;
91 pole[1] = w = Weights(ip);
92 pole[0] = Poles(ip) * w;
98 //=======================================================================
99 //function : PrepareEval
100 //purpose : stores data for Eval in the local arrays
101 // dc.poles and dc.knots
102 //=======================================================================
104 static void PrepareEval
106 Standard_Integer& index,
107 Standard_Integer& dim,
108 Standard_Boolean& rational,
109 const Standard_Integer Degree,
110 const Standard_Boolean Periodic,
111 const TColStd_Array1OfReal& Poles,
112 const TColStd_Array1OfReal& Weights,
113 const TColStd_Array1OfReal& Knots,
114 const TColStd_Array1OfInteger& Mults,
115 BSplCLib_DataContainer& dc)
118 BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
121 BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*dc.knots);
123 index -= Knots.Lower() + Degree;
125 index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults);
127 // check truly rational
128 rational = (&Weights != NULL);
130 Standard_Integer WLower = Weights.Lower() + index;
131 rational = BSplCLib::IsRational(Weights, WLower, WLower + Degree);
137 BSplCLib::BuildEval(Degree, index, Poles, Weights , *dc.poles);
141 BSplCLib::BuildEval(Degree, index, Poles, BSplCLib::NoWeights(), *dc.poles);
145 //=======================================================================
148 //=======================================================================
151 (const Standard_Real U,
152 const Standard_Integer Index,
153 const Standard_Integer Degree,
154 const Standard_Boolean Periodic,
155 const TColStd_Array1OfReal& Poles,
156 const TColStd_Array1OfReal& Weights,
157 const TColStd_Array1OfReal& Knots,
158 const TColStd_Array1OfInteger& Mults,
161 Standard_Integer dim,index = Index;
163 Standard_Boolean rational;
164 BSplCLib_DataContainer dc(Degree);
165 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
166 BSplCLib::Eval(u,Degree,*dc.knots,dim,*dc.poles);
167 if (rational) P = dc.poles[0] / dc.poles[1];
168 else P = dc.poles[0];
171 //=======================================================================
174 //=======================================================================
177 (const Standard_Real U,
178 const Standard_Integer Index,
179 const Standard_Integer Degree,
180 const Standard_Boolean Periodic,
181 const TColStd_Array1OfReal& Poles,
182 const TColStd_Array1OfReal& Weights,
183 const TColStd_Array1OfReal& Knots,
184 const TColStd_Array1OfInteger& Mults,
188 Standard_Integer dim,index = Index;
190 Standard_Boolean rational;
191 BSplCLib_DataContainer dc(Degree);
192 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
193 BSplCLib::Bohm(u,Degree,1,*dc.knots,dim,*dc.poles);
194 Standard_Real *result = dc.poles;
196 PLib::RationalDerivative(Degree,1,1,*dc.poles,*dc.ders);
203 //=======================================================================
206 //=======================================================================
209 (const Standard_Real U,
210 const Standard_Integer Index,
211 const Standard_Integer Degree,
212 const Standard_Boolean Periodic,
213 const TColStd_Array1OfReal& Poles,
214 const TColStd_Array1OfReal& Weights,
215 const TColStd_Array1OfReal& Knots,
216 const TColStd_Array1OfInteger& Mults,
221 Standard_Integer dim,index = Index;
223 Standard_Boolean rational;
224 BSplCLib_DataContainer dc(Degree);
225 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
226 BSplCLib::Bohm(u,Degree,2,*dc.knots,dim,*dc.poles);
227 Standard_Real *result = dc.poles;
229 PLib::RationalDerivative(Degree,2,1,*dc.poles,*dc.ders);
234 if (!rational && (Degree < 2)) V2 = 0.;
238 //=======================================================================
241 //=======================================================================
244 (const Standard_Real U,
245 const Standard_Integer Index,
246 const Standard_Integer Degree,
247 const Standard_Boolean Periodic,
248 const TColStd_Array1OfReal& Poles,
249 const TColStd_Array1OfReal& Weights,
250 const TColStd_Array1OfReal& Knots,
251 const TColStd_Array1OfInteger& Mults,
257 Standard_Integer dim,index = Index;
259 Standard_Boolean rational;
260 BSplCLib_DataContainer dc(Degree);
261 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
262 BSplCLib::Bohm(u,Degree,3,*dc.knots,dim,*dc.poles);
263 Standard_Real *result = dc.poles;
265 PLib::RationalDerivative(Degree,3,1,*dc.poles,*dc.ders);
270 if (!rational && (Degree < 2)) V2 = 0.;
272 if (!rational && (Degree < 3)) V3 = 0.;
276 //=======================================================================
279 //=======================================================================
282 (const Standard_Real U,
283 const Standard_Integer N,
284 const Standard_Integer Index,
285 const Standard_Integer Degree,
286 const Standard_Boolean Periodic,
287 const TColStd_Array1OfReal& Poles,
288 const TColStd_Array1OfReal& Weights,
289 const TColStd_Array1OfReal& Knots,
290 const TColStd_Array1OfInteger& Mults,
293 Standard_Integer dim,index = Index;
295 Standard_Boolean rational;
296 BSplCLib_DataContainer dc(Degree);
297 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
298 BSplCLib::Bohm(u,Degree,N,*dc.knots,dim,*dc.poles);
301 PLib::RationalDerivative(Degree,N,1,*dc.poles,v,Standard_False);
305 if (N > Degree) VN = 0.;
306 else VN = dc.poles[N];
310 //=======================================================================
311 //function : Build BSpline Matrix
312 //purpose : Builds the Bspline Matrix
313 //=======================================================================
316 BSplCLib::BuildBSpMatrix(const TColStd_Array1OfReal& Parameters,
317 const TColStd_Array1OfInteger& ContactOrderArray,
318 const TColStd_Array1OfReal& FlatKnots,
319 const Standard_Integer Degree,
321 Standard_Integer& UpperBandWidth,
322 Standard_Integer& LowerBandWidth)
329 FirstNonZeroBsplineIndex,
334 math_Matrix BSplineBasis(1, MaxOrder,
338 UpperBandWidth = Degree ;
339 LowerBandWidth = Degree ;
340 BandWidth = UpperBandWidth + LowerBandWidth + 1 ;
341 if (Matrix.LowerRow() != Parameters.Lower() ||
342 Matrix.UpperRow() != Parameters.Upper() ||
343 Matrix.LowerCol() != 1 ||
344 Matrix.UpperCol() != BandWidth) {
349 for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
351 BSplCLib::EvalBsplineBasis(1,
352 ContactOrderArray(ii),
357 FirstNonZeroBsplineIndex,
359 if (ErrorCode != 0) {
363 Index = LowerBandWidth + 1 + FirstNonZeroBsplineIndex - ii ;
365 for (jj = 1 ; jj < Index ; jj++) {
366 Matrix.Value(ii,jj) = 0.0e0 ;
369 for (jj = 1 ; jj <= Order ; jj++) {
370 Matrix.Value(ii,Index) = BSplineBasis(ContactOrderArray(ii) + 1, jj) ;
374 for (jj = Index ; jj <= BandWidth ; jj++) {
375 Matrix.Value(ii,jj) = 0.0e0 ;
379 return (ReturnCode) ;
382 //=======================================================================
383 //function : Makes LU decompositiomn without Pivoting
384 //purpose : Builds the Bspline Matrix
385 //=======================================================================
388 BSplCLib::FactorBandedMatrix(math_Matrix& Matrix,
389 const Standard_Integer UpperBandWidth,
390 const Standard_Integer LowerBandWidth,
391 Standard_Integer& PivotIndexProblem)
400 BandWidth = UpperBandWidth + LowerBandWidth + 1 ;
402 Standard_Real Inverse ;
403 PivotIndexProblem = 0 ;
405 for (ii = Matrix.LowerRow() + 1 ; ii <= Matrix.UpperRow() ; ii++) {
406 MinIndex = ( LowerBandWidth - ii + 2 >= 1 ? LowerBandWidth - ii + 2 : 1) ;
408 for (jj = MinIndex ; jj <= LowerBandWidth ; jj++) {
409 Index = ii - LowerBandWidth + jj - 1 ;
410 Inverse = Matrix(Index ,LowerBandWidth + 1) ;
411 if (Abs(Inverse) > RealSmall()) {
412 Inverse = -1.0e0/Inverse ;
416 PivotIndexProblem = Index ;
419 Matrix(ii,jj) = Matrix(ii,jj) * Inverse ;
420 MaxIndex = BandWidth + Index - ii ;
422 for (kk = jj + 1 ; kk <= MaxIndex ; kk++) {
423 Matrix(ii,kk) += Matrix(ii,jj) * Matrix(Index, kk + ii - Index) ;
428 return (ReturnCode) ;
431 //=======================================================================
432 //function : Build BSpline Matrix
433 //purpose : Builds the Bspline Matrix
434 //=======================================================================
437 BSplCLib::EvalBsplineBasis
438 //(const Standard_Integer Side, // = 1 rigth side, -1 left side
439 (const Standard_Integer , // = 1 rigth side, -1 left side
440 const Standard_Integer DerivativeRequest,
441 const Standard_Integer Order,
442 const TColStd_Array1OfReal& FlatKnots,
443 const Standard_Real Parameter,
444 Standard_Integer& FirstNonZeroBsplineIndex,
445 math_Matrix& BsplineBasis)
447 // the matrix must have at least DerivativeRequest + 1
448 // row and Order columns
449 // the result are stored in the following way in
450 // the Bspline matrix
451 // Let i be the FirstNonZeroBsplineIndex and
452 // t be the parameter value, k the order of the
453 // knot vector, r the DerivativeRequest :
479 Standard_Real NewParameter,
484 // , *FlatKnotsArray ;
487 FirstNonZeroBsplineIndex = 0 ;
488 LocalRequest = DerivativeRequest ;
489 if (DerivativeRequest >= Order) {
490 LocalRequest = Order - 1 ;
493 if (BsplineBasis.LowerCol() != 1 ||
494 BsplineBasis.UpperCol() < Order ||
495 BsplineBasis.LowerRow() != 1 ||
496 BsplineBasis.UpperRow() <= LocalRequest) {
500 NumPoles = FlatKnots.Upper() - FlatKnots.Lower() + 1 - Order ;
501 BSplCLib::LocateParameter(Order - 1,
510 FirstNonZeroBsplineIndex = ii - Order + 1 ;
512 BsplineBasis(1,1) = 1.0e0 ;
513 LocalRequest = DerivativeRequest ;
514 if (DerivativeRequest >= Order) {
515 LocalRequest = Order - 1 ;
518 for (qq = 2 ; qq <= Order - LocalRequest ; qq++) {
519 BsplineBasis(1,qq) = 0.0e0 ;
521 for (pp = 1 ; pp <= qq - 1 ; pp++) {
523 // this should be always invertible if ii is correctly computed
525 Factor = (Parameter - FlatKnots(ii - qq + pp + 1))
526 / (FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1)) ;
527 Saved = Factor * BsplineBasis(1,pp) ;
528 BsplineBasis(1,pp) *= (1.0e0 - Factor) ;
529 BsplineBasis(1,pp) += BsplineBasis(1,qq) ;
530 BsplineBasis(1,qq) = Saved ;
534 for (qq = Order - LocalRequest + 1 ; qq <= Order ; qq++) {
536 for (pp = 1 ; pp <= qq - 1 ; pp++) {
537 BsplineBasis(Order - qq + 2,pp) = BsplineBasis(1,pp) ;
539 BsplineBasis(1,qq) = 0.0e0 ;
541 for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) {
542 BsplineBasis(Order - ss + 2,qq) = 0.0e0 ;
545 for (pp = 1 ; pp <= qq - 1 ; pp++) {
546 Inverse = 1.0e0 / (FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1)) ;
547 Factor = (Parameter - FlatKnots(ii - qq + pp + 1)) * Inverse ;
548 Saved = Factor * BsplineBasis(1,pp) ;
549 BsplineBasis(1,pp) *= (1.0e0 - Factor) ;
550 BsplineBasis(1,pp) += BsplineBasis(1,qq) ;
551 BsplineBasis(1,qq) = Saved ;
552 LocalInverse = (Standard_Real) (qq - 1) * Inverse ;
554 for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) {
555 Saved = LocalInverse * BsplineBasis(Order - ss + 2, pp) ;
556 BsplineBasis(Order - ss + 2, pp) *= - LocalInverse ;
557 BsplineBasis(Order - ss + 2, pp) += BsplineBasis(Order - ss + 2,qq) ;
558 BsplineBasis(Order - ss + 2,qq) = Saved ;
563 return (ReturnCode) ;
566 //=======================================================================
567 //function : MovePointAndTangent
569 //=======================================================================
571 void BSplCLib::MovePointAndTangent(const Standard_Real U,
572 const Standard_Integer ArrayDimension,
573 Standard_Real &Delta,
574 Standard_Real &DeltaDerivatives,
575 const Standard_Real Tolerance,
576 const Standard_Integer Degree,
577 const Standard_Boolean Rational,
578 const Standard_Integer StartingCondition,
579 const Standard_Integer EndingCondition,
580 Standard_Real& Poles,
581 const TColStd_Array1OfReal& Weights,
582 const TColStd_Array1OfReal& FlatKnots,
583 Standard_Real& NewPoles,
584 Standard_Integer& ErrorStatus)
586 Standard_Integer num_poles,
600 Standard_Real new_parameter,
612 weights_array = NULL ;
614 weights_array = (Standard_Real *) &Weights(Weights.Lower()) ;
617 poles_array = &Poles ;
618 new_poles_array = &NewPoles ;
619 delta_array = &Delta ;
620 derivatives_array = &DeltaDerivatives ;
622 num_knots = FlatKnots.Length() ;
623 num_poles = num_knots - order ;
624 conditions = StartingCondition + EndingCondition + 4 ;
626 // check validity of input data
628 if (StartingCondition >= -1 &&
629 StartingCondition <= Degree &&
630 EndingCondition >= -1 &&
631 EndingCondition <= Degree &&
632 conditions <= num_poles) {
634 // check the parameter is within bounds
636 start_index = FlatKnots.Lower() + Degree ;
637 end_index = FlatKnots.Upper() - Degree ;
639 // check if there is enough room to move the poles
642 if (StartingCondition == -1) {
643 conditions = conditions && (FlatKnots(start_index) <= U) ;
646 conditions = conditions && (FlatKnots(start_index) + Tolerance < U) ;
648 if (EndingCondition == -1) {
649 conditions = conditions && (FlatKnots(end_index) >= U) ;
652 conditions = conditions && (FlatKnots(end_index) - Tolerance > U) ;
657 // build 2 auxialiary functions
659 TColStd_Array1OfReal schoenberg_points(1,num_poles) ;
660 TColStd_Array1OfReal first_function (1,num_poles) ;
661 TColStd_Array1OfReal second_function (1,num_poles) ;
663 BuildSchoenbergPoints(Degree,
666 start_index = StartingCondition + 2 ;
667 end_index = num_poles - EndingCondition - 1 ;
668 LocateParameter(schoenberg_points,
677 if (index == start_index) {
678 other_index = index + 1 ;
680 else if (index == end_index) {
681 other_index = index -1 ;
683 else if (U - FlatKnots(index) < FlatKnots(index + 1) - U ) {
684 other_index = index - 1 ;
687 other_index = index + 1 ;
691 start_num_poles = StartingCondition + 2 ;
692 end_num_poles = num_poles - EndingCondition - 1 ;
693 if (start_num_poles == 1) {
694 start_value = schoenberg_points(num_poles) - schoenberg_points(1) ;
695 start_value = schoenberg_points(1) - start_value ;
698 start_value = schoenberg_points(start_num_poles - 1) ;
700 if (end_num_poles == num_poles) {
701 end_value = schoenberg_points(num_poles) - schoenberg_points(1) ;
702 end_value = schoenberg_points(num_poles) + end_value ;
705 end_value = schoenberg_points(end_num_poles + 1) ;
708 for (ii = 1 ; ii < start_num_poles ; ii++) {
709 first_function(ii) = 0.e0 ;
710 second_function(ii) = 0.0e0 ;
713 for (ii = end_num_poles + 1 ; ii <= num_poles ; ii++) {
714 first_function(ii) = 0.e0 ;
715 second_function(ii) = 0.0e0 ;
717 divide = 1.0e0 / (schoenberg_points(index) - start_value) ;
719 for (ii = start_num_poles ; ii <= index ; ii++) {
720 value = schoenberg_points(ii) - start_value ;
722 first_function(ii) = 1.0e0 ;
724 for (jj = 0 ; jj < type ; jj++) {
725 first_function(ii) *= value ;
728 divide = 1.0e0 /(end_value - schoenberg_points(index)) ;
730 for (ii = index ; ii <= end_num_poles ; ii++) {
731 value = end_value - schoenberg_points(ii) ;
733 first_function(ii) = 1.0e0 ;
735 for (jj = 0 ; jj < type ; jj++) {
736 first_function(ii) *= value ;
740 divide = 1.0e0 / (schoenberg_points(other_index) - start_value) ;
742 for (ii = start_num_poles ; ii <= other_index ; ii++) {
743 value = schoenberg_points(ii) - start_value ;
745 second_function(ii) = 1.0e0 ;
747 for (jj = 0 ; jj < type ; jj++) {
748 second_function(ii) *= value ;
751 divide = 1.0e0/( end_value - schoenberg_points(other_index)) ;
753 for (ii = other_index ; ii <= end_num_poles ; ii++) {
754 value = end_value - schoenberg_points(ii) ;
756 second_function(ii) = 1.0e0 ;
758 for (jj = 0 ; jj < type ; jj++) {
759 second_function(ii) *= value ;
764 // compute the point and derivatives of both functions
766 Standard_Real results[2][2],
767 weights_results[2][2];
768 Standard_Integer extrap_mode[2],
769 derivative_request = 1,
771 Standard_Boolean periodic_flag = Standard_False ;
773 extrap_mode[0] = Degree ;
774 extrap_mode[1] = Degree ;
777 // evaluate in homogenised form
789 weights_results[0][0]) ;
801 weights_results[1][0]) ;
803 // compute the rational derivatives values
806 for (ii = 0 ; ii < 2 ; ii++) {
807 PLib::RationalDerivatives(1,
810 weights_results[ii][0],
837 for (ii = 0 ; ii < 2 ; ii++) {
839 for (jj = 0 ; jj < 2 ; jj++) {
840 a_matrix.SetValue(ii+1,jj+1,results[ii][jj]) ;
844 TColStd_Array1OfReal the_a_vector(0,ArrayDimension-1) ;
845 TColStd_Array1OfReal the_b_vector(0,ArrayDimension-1) ;
847 for( ii = 0 ; ii < ArrayDimension ; ii++) {
849 a_matrix.Value(1,1) * delta_array[ii] +
850 a_matrix.Value(2,1) * derivatives_array[ii] ;
852 a_matrix.Value(1,2) * delta_array[ii] +
853 a_matrix.Value(2,2) * derivatives_array[ii] ;
857 for (ii = 0 ; ii < num_poles ; ii++) {
859 for (jj = 0 ; jj < ArrayDimension ; jj++) {
860 new_poles_array[index] = poles_array[index] ;
861 new_poles_array[index] +=
862 first_function(ii+1) * the_a_vector(jj) ;
863 new_poles_array[index] +=
864 second_function(ii+1) * the_b_vector(jj) ;
878 //=======================================================================
879 //function : FunctionMultiply
881 //=======================================================================
883 void BSplCLib::FunctionMultiply
884 (const BSplCLib_EvaluatorFunction & FunctionPtr,
885 const Standard_Integer BSplineDegree,
886 const TColStd_Array1OfReal & BSplineFlatKnots,
887 const Standard_Integer PolesDimension,
888 Standard_Real & Poles,
889 const TColStd_Array1OfReal & FlatKnots,
890 const Standard_Integer NewDegree,
891 Standard_Real & NewPoles,
892 Standard_Integer & Status)
897 Standard_Integer extrap_mode[2],
900 derivative_request = 0 ;
901 Standard_Boolean periodic_flag = Standard_False ;
902 Standard_Real result,
905 *array_of_new_poles ;
907 array_of_poles = (Standard_Real *) &NewPoles ;
909 extrap_mode[1] = BSplineDegree ;
911 FlatKnots.Length() - NewDegree - 1 ;
912 start_end[0] = FlatKnots(NewDegree+1) ;
913 start_end[1] = FlatKnots(num_new_poles+1) ;
914 TColStd_Array1OfReal parameters(1,num_new_poles) ;
915 TColStd_Array1OfInteger contact_order_array(1,num_new_poles) ;
916 TColStd_Array1OfReal new_poles_array(1, num_new_poles * PolesDimension) ;
919 (Standard_Real *) &new_poles_array(1) ;
920 BuildSchoenbergPoints(NewDegree,
924 // on recadre sur les bornes
926 if (parameters(1) < start_end[0]) {
927 parameters(1) = start_end[0] ;
929 if (parameters(num_new_poles) > start_end[1]) {
930 parameters(num_new_poles) = start_end[1] ;
934 for (ii = 1 ; ii <= num_new_poles ; ii++) {
935 contact_order_array(ii) = 0 ;
936 FunctionPtr.Evaluate (contact_order_array(ii),
954 array_of_new_poles[index]) ;
956 for (jj = 0 ; jj < PolesDimension ; jj++) {
957 array_of_new_poles[index] *= result ;
961 Interpolate(NewDegree,
966 array_of_new_poles[0],
969 for (ii = 0 ; ii < num_new_poles * PolesDimension ; ii++) {
970 array_of_poles[ii] = array_of_new_poles[ii] ;
977 //=======================================================================
978 // function : FunctionMultiply
980 //=======================================================================
982 void BSplCLib::FunctionReparameterise
983 (const BSplCLib_EvaluatorFunction & FunctionPtr,
984 const Standard_Integer BSplineDegree,
985 const TColStd_Array1OfReal & BSplineFlatKnots,
986 const Standard_Integer PolesDimension,
987 Standard_Real & Poles,
988 const TColStd_Array1OfReal & FlatKnots,
989 const Standard_Integer NewDegree,
990 Standard_Real & NewPoles,
991 Standard_Integer & Status)
996 Standard_Integer extrap_mode[2],
999 derivative_request = 0 ;
1000 Standard_Boolean periodic_flag = Standard_False ;
1001 Standard_Real result,
1004 *array_of_new_poles ;
1006 array_of_poles = (Standard_Real *) &NewPoles ;
1008 extrap_mode[1] = BSplineDegree ;
1010 FlatKnots.Length() - NewDegree - 1 ;
1011 start_end[0] = FlatKnots(NewDegree+1) ;
1012 start_end[1] = FlatKnots(num_new_poles+1) ;
1013 TColStd_Array1OfReal parameters(1,num_new_poles) ;
1014 TColStd_Array1OfInteger contact_order_array(1,num_new_poles) ;
1015 TColStd_Array1OfReal new_poles_array(1, num_new_poles * PolesDimension) ;
1017 array_of_new_poles =
1018 (Standard_Real *) &new_poles_array(1) ;
1019 BuildSchoenbergPoints(NewDegree,
1024 for (ii = 1 ; ii <= num_new_poles ; ii++) {
1025 contact_order_array(ii) = 0 ;
1026 FunctionPtr.Evaluate (contact_order_array(ii),
1044 array_of_new_poles[index]) ;
1045 index += PolesDimension ;
1047 Interpolate(NewDegree,
1050 contact_order_array,
1052 array_of_new_poles[0],
1055 for (ii = 0 ; ii < num_new_poles * PolesDimension ; ii++) {
1056 array_of_poles[ii] = array_of_new_poles[ii] ;
1063 //=======================================================================
1064 //function : FunctionMultiply
1066 //=======================================================================
1068 void BSplCLib::FunctionMultiply
1069 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1070 const Standard_Integer BSplineDegree,
1071 const TColStd_Array1OfReal & BSplineFlatKnots,
1072 const TColStd_Array1OfReal & Poles,
1073 const TColStd_Array1OfReal & FlatKnots,
1074 const Standard_Integer NewDegree,
1075 TColStd_Array1OfReal & NewPoles,
1076 Standard_Integer & Status)
1078 Standard_Integer num_bspline_poles =
1079 BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1080 Standard_Integer num_new_poles =
1081 FlatKnots.Length() - NewDegree - 1 ;
1083 if (Poles.Length() != num_bspline_poles ||
1084 NewPoles.Length() != num_new_poles) {
1085 Standard_ConstructionError::Raise();
1087 Standard_Real * array_of_poles =
1088 (Standard_Real *) &Poles(Poles.Lower()) ;
1089 Standard_Real * array_of_new_poles =
1090 (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1091 BSplCLib::FunctionMultiply(FunctionPtr,
1098 array_of_new_poles[0],
1102 //=======================================================================
1103 //function : FunctionReparameterise
1105 //=======================================================================
1107 void BSplCLib::FunctionReparameterise
1108 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1109 const Standard_Integer BSplineDegree,
1110 const TColStd_Array1OfReal & BSplineFlatKnots,
1111 const TColStd_Array1OfReal & Poles,
1112 const TColStd_Array1OfReal & FlatKnots,
1113 const Standard_Integer NewDegree,
1114 TColStd_Array1OfReal & NewPoles,
1115 Standard_Integer & Status)
1117 Standard_Integer num_bspline_poles =
1118 BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1119 Standard_Integer num_new_poles =
1120 FlatKnots.Length() - NewDegree - 1 ;
1122 if (Poles.Length() != num_bspline_poles ||
1123 NewPoles.Length() != num_new_poles) {
1124 Standard_ConstructionError::Raise();
1126 Standard_Real * array_of_poles =
1127 (Standard_Real *) &Poles(Poles.Lower()) ;
1128 Standard_Real * array_of_new_poles =
1129 (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1130 BSplCLib::FunctionReparameterise(
1138 array_of_new_poles[0],
1142 //=======================================================================
1143 //function : MergeBSplineKnots
1145 //=======================================================================
1146 void BSplCLib::MergeBSplineKnots
1147 (const Standard_Real Tolerance,
1148 const Standard_Real StartValue,
1149 const Standard_Real EndValue,
1150 const Standard_Integer Degree1,
1151 const TColStd_Array1OfReal & Knots1,
1152 const TColStd_Array1OfInteger & Mults1,
1153 const Standard_Integer Degree2,
1154 const TColStd_Array1OfReal & Knots2,
1155 const TColStd_Array1OfInteger & Mults2,
1156 Standard_Integer & NumPoles,
1157 Handle(TColStd_HArray1OfReal) & NewKnots,
1158 Handle(TColStd_HArray1OfInteger) & NewMults)
1160 Standard_Integer ii,
1167 if (StartValue < EndValue - Tolerance) {
1168 TColStd_Array1OfReal knots1(1,Knots1.Length()) ;
1169 TColStd_Array1OfReal knots2(1,Knots2.Length()) ;
1170 degree = Degree1 + Degree2 ;
1173 for (ii = Knots1.Lower() ; ii <= Knots1.Upper() ; ii++) {
1174 knots1(index) = Knots1(ii) ;
1179 for (ii = Knots2.Lower() ; ii <= Knots2.Upper() ; ii++) {
1180 knots2(index) = Knots2(ii) ;
1183 BSplCLib::Reparametrize(StartValue,
1187 BSplCLib::Reparametrize(StartValue,
1193 for (ii = 1 ; ii <= knots1.Length() ; ii++) {
1195 while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) - Tolerance) {
1200 while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) + Tolerance) {
1206 new TColStd_HArray1OfReal(1,num_knots) ;
1208 new TColStd_HArray1OfInteger(1,num_knots) ;
1212 for (ii = 1 ; ii <= knots1.Length() ; ii++) {
1214 while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) - Tolerance) {
1215 NewKnots->ChangeArray1()(num_knots) = knots2(jj) ;
1216 NewMults->ChangeArray1()(num_knots) = Mults2(jj) + Degree1 ;
1220 set_mults_flag = 0 ;
1222 while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) + Tolerance) {
1223 continuity = Min(Degree1 - Mults1(ii), Degree2 - Mults2(jj)) ;
1224 set_mults_flag = 1 ;
1225 NewMults->ChangeArray1()(num_knots) = degree - continuity ;
1229 NewKnots->ChangeArray1()(num_knots) = knots1(ii) ;
1230 if (! set_mults_flag) {
1231 NewMults->ChangeArray1()(num_knots) = Mults1(ii) + Degree2 ;
1236 NewMults->ChangeArray1()(1) = degree + 1 ;
1237 NewMults->ChangeArray1()(num_knots) = degree + 1 ;
1240 for (ii = 1 ; ii <= num_knots ; ii++) {
1241 index += NewMults->Value(ii) ;
1243 NumPoles = index - degree - 1 ;