0027275: Unused formal parameter in BSplCLib::EvalBsplineBasis
[occt.git] / src / BSplCLib / BSplCLib_2.cxx
1 // Created on: 1995-01-17
2 // Created by: Laurent BOURESCHE
3 // Copyright (c) 1995-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 // xab : modified 15-Mar-95 : added BuildBSplMatrix,FactorBandedMatrix,
18 //                            EvalBsplineBasis,
19 //                            EvalPolynomial : Horners method
20
21 #include <Standard_Stream.hxx>
22
23 #include <BSplCLib.hxx>
24 #include <gp_Mat2d.hxx>
25 #include <PLib.hxx>
26 #include <Standard_ConstructionError.hxx>
27 #include <TColStd_Array1OfReal.hxx>
28 #include <TColStd_Array1OfInteger.hxx>
29 #include <TColStd_HArray1OfReal.hxx>
30 #include <TColStd_HArray1OfInteger.hxx>
31
32 #include <math_Matrix.hxx>
33
34 static const Standard_Integer aGlobalMaxDegree = 25;
35
36 //=======================================================================
37 //struct : BSplCLib_DataContainer 
38 //purpose: Auxiliary structure providing buffers for poles and knots used in
39 //         evaluation of bspline (allocated in the stack)
40 //=======================================================================
41
42 struct BSplCLib_DataContainer 
43 {
44   BSplCLib_DataContainer(Standard_Integer Degree)
45   {
46     (void)Degree; // avoid compiler warning
47     Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree() ||
48         BSplCLib::MaxDegree() > aGlobalMaxDegree,
49         "BSplCLib: bspline degree is greater than maximum supported");
50   }
51
52   Standard_Real poles[2*(25+1)];
53   Standard_Real knots[2*25];
54   Standard_Real ders[4];
55 };
56
57 // methods for 1 dimensional BSplines
58
59 //=======================================================================
60 //function : BuildEval
61 //purpose  : builds the local array for evaluation
62 //=======================================================================
63
64 void  BSplCLib::BuildEval(const Standard_Integer         Degree,
65                           const Standard_Integer         Index,
66                           const TColStd_Array1OfReal&    Poles, 
67                           const TColStd_Array1OfReal*    Weights,
68                           Standard_Real&                 LP)
69 {
70   Standard_Integer PLower = Poles.Lower();
71   Standard_Integer PUpper = Poles.Upper();
72   Standard_Integer i;
73   Standard_Integer ip = PLower + Index - 1;
74   Standard_Real w, *pole = &LP;
75   if (Weights == NULL) {
76     
77     for (i = 0; i <= Degree; i++) {
78       ip++;
79       if (ip > PUpper) ip = PLower;
80       pole[0] = Poles(ip);
81       pole += 1;
82     }
83   }
84   else {
85     
86     for (i = 0; i <= Degree; i++) {
87       ip++;
88       if (ip > PUpper) ip = PLower;
89       pole[1] = w = (*Weights)(ip);
90       pole[0] = Poles(ip) * w;
91       pole += 2;
92     }
93   }
94 }
95
96 //=======================================================================
97 //function : PrepareEval
98 //purpose  : stores data for Eval in the local arrays
99 //           dc.poles and dc.knots
100 //=======================================================================
101
102 static void PrepareEval
103 (Standard_Real&                 u,                  
104  Standard_Integer&              index, 
105  Standard_Integer&              dim,
106  Standard_Boolean&              rational,
107  const Standard_Integer         Degree,     
108  const Standard_Boolean         Periodic,
109  const TColStd_Array1OfReal&    Poles,  
110  const TColStd_Array1OfReal*    Weights,
111  const TColStd_Array1OfReal&    Knots,  
112  const TColStd_Array1OfInteger* Mults,
113  BSplCLib_DataContainer&        dc) 
114 {                    
115   // Set the Index
116   BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
117   
118   // make the knots
119   BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*dc.knots);
120   if (Mults == NULL)
121     index -= Knots.Lower() + Degree;
122   else
123     index = BSplCLib::PoleIndex(Degree,index,Periodic,*Mults);
124   
125   // check truly rational
126   rational = (Weights != NULL);
127   if (rational) {
128     Standard_Integer WLower = Weights->Lower() + index;
129     rational = BSplCLib::IsRational(*Weights, WLower, WLower + Degree);
130   }
131   
132   // make the poles
133   if(rational) {
134     dim = 2;
135     BSplCLib::BuildEval(Degree, index, Poles, Weights              , *dc.poles);
136   }
137   else {
138     dim = 1;
139     BSplCLib::BuildEval(Degree, index, Poles, BSplCLib::NoWeights(), *dc.poles);
140   }
141 }
142
143 //=======================================================================
144 //function : D0
145 //purpose  : 
146 //=======================================================================
147
148 void BSplCLib::D0 
149 (const Standard_Real            U,                  
150  const Standard_Integer         Index,          
151  const Standard_Integer         Degree,     
152  const Standard_Boolean         Periodic,
153  const TColStd_Array1OfReal&    Poles,  
154  const TColStd_Array1OfReal*    Weights,
155  const TColStd_Array1OfReal&    Knots,  
156  const TColStd_Array1OfInteger* Mults,  
157  Standard_Real&                 P) 
158 {                    
159   Standard_Integer dim,index = Index;
160   Standard_Real    u = U;
161   Standard_Boolean rational;
162   BSplCLib_DataContainer dc(Degree);
163   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
164   BSplCLib::Eval(u,Degree,*dc.knots,dim,*dc.poles);
165   if (rational) P = dc.poles[0] / dc.poles[1];
166   else          P = dc.poles[0];
167 }
168
169 //=======================================================================
170 //function : D1
171 //purpose  : 
172 //=======================================================================
173
174 void BSplCLib::D1
175 (const Standard_Real            U,                  
176  const Standard_Integer         Index,          
177  const Standard_Integer         Degree,     
178  const Standard_Boolean         Periodic,
179  const TColStd_Array1OfReal&    Poles,  
180  const TColStd_Array1OfReal*    Weights,
181  const TColStd_Array1OfReal&    Knots,  
182  const TColStd_Array1OfInteger* Mults,  
183  Standard_Real&                 P,
184  Standard_Real&                 V) 
185 {                    
186   Standard_Integer dim,index = Index;
187   Standard_Real    u = U;
188   Standard_Boolean rational;
189   BSplCLib_DataContainer dc(Degree);
190   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
191   BSplCLib::Bohm(u,Degree,1,*dc.knots,dim,*dc.poles);
192   Standard_Real *result = dc.poles;
193   if (rational) {
194     PLib::RationalDerivative(Degree,1,1,*dc.poles,*dc.ders);
195     result = dc.ders;
196   }
197   P = result[0];
198   V = result[1];
199 }
200
201 //=======================================================================
202 //function : D2
203 //purpose  : 
204 //=======================================================================
205
206 void BSplCLib::D2
207 (const Standard_Real            U,                  
208  const Standard_Integer         Index,          
209  const Standard_Integer         Degree,     
210  const Standard_Boolean         Periodic,
211  const TColStd_Array1OfReal&    Poles,  
212  const TColStd_Array1OfReal*    Weights,
213  const TColStd_Array1OfReal&    Knots,  
214  const TColStd_Array1OfInteger* Mults,  
215  Standard_Real&                 P,
216  Standard_Real&                 V1,
217  Standard_Real&                 V2) 
218 {                    
219   Standard_Integer dim,index = Index;
220   Standard_Real    u = U;
221   Standard_Boolean rational;
222   BSplCLib_DataContainer dc(Degree);
223   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
224   BSplCLib::Bohm(u,Degree,2,*dc.knots,dim,*dc.poles);
225   Standard_Real *result = dc.poles;
226   if (rational) {
227     PLib::RationalDerivative(Degree,2,1,*dc.poles,*dc.ders);
228     result = dc.ders;
229   }
230   P  = result[0];
231   V1 = result[1];
232   if (!rational && (Degree < 2)) V2 = 0.;
233   else                           V2 = result[2];
234 }
235
236 //=======================================================================
237 //function : D3
238 //purpose  : 
239 //=======================================================================
240
241 void BSplCLib::D3
242 (const Standard_Real            U,                  
243  const Standard_Integer         Index,          
244  const Standard_Integer         Degree,     
245  const Standard_Boolean         Periodic,
246  const TColStd_Array1OfReal&    Poles,  
247  const TColStd_Array1OfReal*    Weights,
248  const TColStd_Array1OfReal&    Knots,  
249  const TColStd_Array1OfInteger* Mults,  
250  Standard_Real&                 P,
251  Standard_Real&                 V1,
252  Standard_Real&                 V2,
253  Standard_Real&                 V3) 
254 {                    
255   Standard_Integer dim,index = Index;
256   Standard_Real    u = U;
257   Standard_Boolean rational;
258   BSplCLib_DataContainer dc(Degree);
259   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
260   BSplCLib::Bohm(u,Degree,3,*dc.knots,dim,*dc.poles);
261   Standard_Real *result = dc.poles;
262   if (rational) {
263     PLib::RationalDerivative(Degree,3,1,*dc.poles,*dc.ders);
264     result = dc.ders;
265   }
266   P  = result[0];
267   V1 = result[1];
268   if (!rational && (Degree < 2)) V2 = 0.;
269   else                           V2 = result[2];
270   if (!rational && (Degree < 3)) V3 = 0.;
271   else                           V3 = result[3];
272 }
273
274 //=======================================================================
275 //function : DN
276 //purpose  : 
277 //=======================================================================
278
279 void BSplCLib::DN
280 (const Standard_Real            U,                  
281  const Standard_Integer         N,          
282  const Standard_Integer         Index,          
283  const Standard_Integer         Degree,     
284  const Standard_Boolean         Periodic,
285  const TColStd_Array1OfReal&    Poles,  
286  const TColStd_Array1OfReal*    Weights,
287  const TColStd_Array1OfReal&    Knots,  
288  const TColStd_Array1OfInteger* Mults,  
289  Standard_Real&                 VN) 
290 {                    
291   Standard_Integer dim,index = Index;
292   Standard_Real    u = U;
293   Standard_Boolean rational;
294   BSplCLib_DataContainer dc(Degree);
295   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
296   BSplCLib::Bohm(u,Degree,N,*dc.knots,dim,*dc.poles);
297   if (rational) {
298     Standard_Real v;
299     PLib::RationalDerivative(Degree,N,1,*dc.poles,v,Standard_False);
300     VN = v;
301   }
302   else {
303     if (N > Degree) VN = 0.;
304     else            VN = dc.poles[N];
305   }
306 }
307
308 //=======================================================================
309 //function : Build BSpline Matrix
310 //purpose  : Builds the Bspline Matrix
311 //=======================================================================
312
313 Standard_Integer 
314 BSplCLib::BuildBSpMatrix(const  TColStd_Array1OfReal&     Parameters,
315                          const  TColStd_Array1OfInteger&  ContactOrderArray,
316                          const  TColStd_Array1OfReal&     FlatKnots,
317                          const  Standard_Integer          Degree,   
318                          math_Matrix&                     Matrix,
319                          Standard_Integer&                UpperBandWidth,
320                          Standard_Integer&                LowerBandWidth) 
321 {
322   Standard_Integer ii,
323   jj,
324   Index,
325   ErrorCode,
326   ReturnCode = 0,
327   FirstNonZeroBsplineIndex,
328   BandWidth,
329   MaxOrder = aGlobalMaxDegree+1,
330   Order ;
331   
332   math_Matrix   BSplineBasis(1, MaxOrder,
333                              1, MaxOrder) ;
334   
335   Order = Degree + 1 ;
336   UpperBandWidth = Degree ;
337   LowerBandWidth = Degree ;
338   BandWidth = UpperBandWidth + LowerBandWidth + 1 ;
339   if (Matrix.LowerRow() != Parameters.Lower()  || 
340       Matrix.UpperRow() != Parameters.Upper()  ||
341       Matrix.LowerCol() != 1  || 
342       Matrix.UpperCol() != BandWidth) {
343     ReturnCode = 1;
344     goto FINISH ;
345   }
346   
347   for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
348     ErrorCode =
349       BSplCLib::EvalBsplineBasis(ContactOrderArray(ii),
350                                  Order,
351                                  FlatKnots,
352                                  Parameters(ii),
353                                  
354                                  FirstNonZeroBsplineIndex,
355                                  BSplineBasis) ;
356     if (ErrorCode != 0) {
357       ReturnCode = 2 ;
358       goto FINISH ;
359     }
360     Index = LowerBandWidth + 1 + FirstNonZeroBsplineIndex - ii ;
361
362     for (jj = 1 ; jj < Index ; jj++) {
363       Matrix.Value(ii,jj) = 0.0e0 ;
364     }
365
366     for (jj = 1 ; jj <= Order ; jj++) {
367       Matrix.Value(ii,Index) = BSplineBasis(ContactOrderArray(ii) + 1, jj) ;
368       Index += 1 ;
369     }
370     
371     for (jj = Index ; jj <= BandWidth ; jj++) {
372       Matrix.Value(ii,jj) = 0.0e0 ;
373     }
374   } 
375   FINISH : ;
376   return (ReturnCode) ;
377 }
378
379 //=======================================================================
380 //function : Makes LU decompositiomn without Pivoting
381 //purpose  : Builds the Bspline Matrix
382 //=======================================================================
383
384 Standard_Integer 
385 BSplCLib::FactorBandedMatrix(math_Matrix&   Matrix,
386                              const Standard_Integer UpperBandWidth,
387                              const Standard_Integer LowerBandWidth,
388                              Standard_Integer&  PivotIndexProblem) 
389 {
390   Standard_Integer ii,
391   jj,
392   kk,
393   Index,
394   MinIndex,
395   MaxIndex,
396   ReturnCode = 0,
397   BandWidth = UpperBandWidth + LowerBandWidth + 1 ;
398   
399   Standard_Real Inverse ;
400   PivotIndexProblem = 0 ;
401
402   for (ii = Matrix.LowerRow() + 1 ; ii <= Matrix.UpperRow() ; ii++) {
403     MinIndex = ( LowerBandWidth - ii + 2 >= 1 ? LowerBandWidth - ii + 2 : 1) ;
404
405     for (jj = MinIndex ; jj <= LowerBandWidth  ; jj++) {
406       Index = ii - LowerBandWidth + jj - 1 ; 
407       Inverse = Matrix(Index ,LowerBandWidth + 1) ;
408       if (Abs(Inverse) > RealSmall()) {
409         Inverse = -1.0e0/Inverse ;
410       }
411       else {
412         ReturnCode = 1 ;
413         PivotIndexProblem = Index ;
414         goto FINISH ;
415       }
416       Matrix(ii,jj) = Matrix(ii,jj) * Inverse ;
417       MaxIndex = BandWidth + Index - ii ;
418
419       for (kk = jj + 1 ; kk <= MaxIndex ; kk++) {
420         Matrix(ii,kk) += Matrix(ii,jj) * Matrix(Index, kk + ii - Index) ;
421       }
422     }                       
423   }
424   FINISH :
425     return (ReturnCode) ;
426 }
427
428 //=======================================================================
429 //function : Build BSpline Matrix
430 //purpose  : Builds the Bspline Matrix
431 //=======================================================================
432
433 Standard_Integer 
434 BSplCLib::EvalBsplineBasis
435 (const  Standard_Integer              DerivativeRequest,
436  const  Standard_Integer              Order,
437  const  TColStd_Array1OfReal&         FlatKnots,
438  const  Standard_Real                 Parameter,
439  Standard_Integer&             FirstNonZeroBsplineIndex,
440  math_Matrix&                  BsplineBasis,
441  Standard_Boolean              isPeriodic)
442 {
443   // the matrix must have at least DerivativeRequest + 1
444   //   row and Order columns
445   // the result are stored in the following way in
446   // the Bspline matrix 
447   // Let i be the FirstNonZeroBsplineIndex and 
448   // t be the parameter value, k the order of the 
449   // knot vector, r the DerivativeRequest :
450   //   
451   //   B (t)   B (t)                     B (t)
452   //    i       i+1                       i+k-1
453   //   
454   //    (1)     (1)                       (1) 
455   //   B (t)   B (t)                     B (t)
456   //    i       i+1                       i+k-1
457   //  
458   //
459   //
460   //
461   //    (r)     (r)                       (r) 
462   //   B (t)   B (t)                     B (t)
463   //    i       i+1                       i+k-1
464   //
465   Standard_Integer  
466     ReturnCode,
467   ii,
468   pp,
469   qq,
470   ss,
471   NumPoles,
472   LocalRequest ;
473 //  ,Index ;
474   
475   Standard_Real NewParameter,
476   Inverse,
477   Factor,
478   LocalInverse,
479   Saved ;
480 // , *FlatKnotsArray ;
481   
482   ReturnCode = 0 ;
483   FirstNonZeroBsplineIndex = 0 ;
484   LocalRequest = DerivativeRequest ;
485   if (DerivativeRequest >= Order) {
486     LocalRequest = Order - 1 ;
487   }
488   
489   if (BsplineBasis.LowerCol() != 1    ||
490       BsplineBasis.UpperCol() < Order ||
491       BsplineBasis.LowerRow() != 1    ||
492       BsplineBasis.UpperRow() <= LocalRequest) {
493     ReturnCode = 1;
494     goto FINISH ;
495   }
496   NumPoles = FlatKnots.Upper() - FlatKnots.Lower() + 1 - Order ;
497   BSplCLib::LocateParameter(Order - 1, 
498                             FlatKnots,
499                             Parameter,
500                             isPeriodic, 
501                             Order, 
502                             NumPoles+1, 
503                             ii,
504                             NewParameter) ;
505   
506   FirstNonZeroBsplineIndex = ii - Order + 1 ;
507   
508   BsplineBasis(1,1) = 1.0e0 ;
509   LocalRequest = DerivativeRequest ;
510   if (DerivativeRequest >= Order) {
511     LocalRequest = Order - 1 ;
512   }
513
514   for (qq = 2 ; qq <= Order - LocalRequest ; qq++) {
515     BsplineBasis(1,qq) = 0.0e0 ;
516     
517     for (pp = 1 ; pp <= qq - 1 ; pp++) {
518       //
519       // this should be always invertible if ii is correctly computed 
520       //
521       Factor = (Parameter - FlatKnots(ii - qq + pp + 1)) 
522         / (FlatKnots(ii + pp)   - FlatKnots(ii - qq + pp + 1)) ; 
523       Saved = Factor *    BsplineBasis(1,pp) ;
524       BsplineBasis(1,pp) *= (1.0e0 - Factor) ;
525       BsplineBasis(1,pp) += BsplineBasis(1,qq) ;
526       BsplineBasis(1,qq) = Saved ;
527     }
528   }
529   
530   for (qq = Order - LocalRequest + 1 ; qq <= Order ; qq++) {
531     
532     for (pp = 1 ; pp <= qq - 1 ; pp++) {
533       BsplineBasis(Order - qq + 2,pp) = BsplineBasis(1,pp) ;
534     }
535     BsplineBasis(1,qq) = 0.0e0 ;
536
537     for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) {
538       BsplineBasis(Order - ss + 2,qq) = 0.0e0 ;
539     }
540
541     for (pp = 1 ; pp <= qq - 1 ; pp++) {
542       Inverse = 1.0e0 / (FlatKnots(ii + pp)  - FlatKnots(ii - qq + pp + 1)) ;
543       Factor  =  (Parameter - FlatKnots(ii - qq + pp + 1)) * Inverse ;
544       Saved = Factor *                 BsplineBasis(1,pp) ;
545       BsplineBasis(1,pp) *= (1.0e0 - Factor) ;
546       BsplineBasis(1,pp) += BsplineBasis(1,qq) ;
547       BsplineBasis(1,qq) = Saved ;
548       LocalInverse = (Standard_Real) (qq - 1) * Inverse ;
549
550       for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) {
551         Saved = LocalInverse * BsplineBasis(Order - ss + 2, pp) ;
552         BsplineBasis(Order - ss + 2, pp) *= - LocalInverse  ;
553         BsplineBasis(Order - ss + 2, pp) +=   BsplineBasis(Order - ss + 2,qq) ;
554         BsplineBasis(Order - ss + 2,qq) = Saved ;
555       }
556     }
557   }
558   FINISH :
559     return (ReturnCode) ;
560 }
561
562 //=======================================================================
563 //function : MovePointAndTangent
564 //purpose  : 
565 //=======================================================================
566
567 void BSplCLib::MovePointAndTangent(const Standard_Real    U,
568                                    const Standard_Integer ArrayDimension,
569                                    Standard_Real    &Delta,
570                                    Standard_Real    &DeltaDerivatives,
571                                    const Standard_Real    Tolerance,
572                                    const Standard_Integer Degree,
573                                    const Standard_Boolean Rational,
574                                    const Standard_Integer StartingCondition,
575                                    const Standard_Integer EndingCondition,
576                                    Standard_Real&         Poles,
577                                    const TColStd_Array1OfReal&   Weights,
578                                    const TColStd_Array1OfReal&   FlatKnots,
579                                    Standard_Real&        NewPoles,
580                                    Standard_Integer&     ErrorStatus) 
581 {
582   Standard_Integer num_poles,
583   num_knots,
584   ii,
585   jj,
586   conditions,
587   start_num_poles,
588   end_num_poles,
589   index,
590   start_index,
591   end_index,
592   other_index,
593   type,
594   order ;
595   
596   Standard_Real    new_parameter,
597   value,
598   divide,
599   end_value,
600   start_value,
601   *poles_array,
602   *new_poles_array,
603   *delta_array,
604   *derivatives_array,
605   *weights_array ;
606   
607   ErrorStatus = 0 ;
608   weights_array = NULL ;
609   if (Rational) {
610     weights_array = (Standard_Real *) &Weights(Weights.Lower()) ;
611   }
612   
613   poles_array = &Poles ;
614   new_poles_array = &NewPoles ;
615   delta_array = &Delta ;
616   derivatives_array = &DeltaDerivatives ;
617   order = Degree + 1 ;
618   num_knots = FlatKnots.Length() ;
619   num_poles = num_knots - order ;
620   conditions = StartingCondition + EndingCondition + 4 ;
621   //
622   // check validity of input data
623   //
624   if (StartingCondition >= -1 &&  
625       StartingCondition <= Degree &&
626       EndingCondition >= -1 &&
627       EndingCondition <= Degree &&
628       conditions <= num_poles) {
629     //
630     // check the parameter is within bounds 
631     //
632     start_index = FlatKnots.Lower() + Degree ;
633     end_index = FlatKnots.Upper() - Degree ;
634     //
635     //  check if there is enough room to move the poles
636     //
637     conditions = 1 ;
638     if (StartingCondition == -1) {
639       conditions = conditions && (FlatKnots(start_index) <=  U) ;
640     }
641     else {
642       conditions = conditions && (FlatKnots(start_index) + Tolerance < U) ;
643     }
644     if (EndingCondition == -1) {
645       conditions = conditions && (FlatKnots(end_index) >= U) ;
646     }
647     else {
648       conditions = conditions && (FlatKnots(end_index) - Tolerance > U) ; 
649     }
650     
651     if (conditions) {
652       //
653       // build 2 auxialiary functions 
654       // 
655       TColStd_Array1OfReal schoenberg_points(1,num_poles) ;
656       TColStd_Array1OfReal first_function  (1,num_poles) ;
657       TColStd_Array1OfReal second_function (1,num_poles) ;
658       
659       BuildSchoenbergPoints(Degree,
660                             FlatKnots,
661                             schoenberg_points) ;
662       start_index = StartingCondition + 2 ;
663       end_index = num_poles - EndingCondition - 1 ;
664       LocateParameter(schoenberg_points,
665                       U,
666                       Standard_False,
667                       start_index,
668                       end_index,
669                       index,
670                       new_parameter, 
671                       0, 1) ;
672       
673       if (index == start_index) {
674         other_index = index + 1 ;
675       } 
676       else if (index == end_index) {
677         other_index = index -1  ;
678       }
679       else if (U - FlatKnots(index) < FlatKnots(index + 1) - U ) {
680         other_index = index - 1 ;
681       }
682       else {
683         other_index = index + 1 ;
684       }
685       type = 3 ;
686       
687       start_num_poles = StartingCondition + 2 ;
688       end_num_poles   = num_poles - EndingCondition - 1 ;
689       if (start_num_poles == 1) {
690         start_value = schoenberg_points(num_poles) - schoenberg_points(1) ;
691         start_value = schoenberg_points(1) - start_value ;
692       }
693       else {
694         start_value = schoenberg_points(start_num_poles - 1) ;
695       }
696       if (end_num_poles == num_poles) {
697         end_value = schoenberg_points(num_poles) - schoenberg_points(1) ;
698         end_value = schoenberg_points(num_poles) + end_value ;
699       }
700       else {
701         end_value = schoenberg_points(end_num_poles + 1) ;
702       }
703       
704       for (ii = 1 ; ii < start_num_poles ; ii++) {
705         first_function(ii) = 0.e0 ;
706         second_function(ii) = 0.0e0 ;
707       }
708
709       for (ii = end_num_poles + 1 ; ii <= num_poles ; ii++) {
710         first_function(ii) = 0.e0 ;
711         second_function(ii) = 0.0e0 ;
712       }
713       divide = 1.0e0 / (schoenberg_points(index) - start_value) ;
714
715       for (ii = start_num_poles  ; ii <= index ; ii++) {
716         value = schoenberg_points(ii) - start_value ;
717         value *= divide ;
718         first_function(ii) = 1.0e0 ;
719         
720         for (jj = 0 ; jj < type ; jj++) {
721           first_function(ii) *= value ;
722         }
723       }
724       divide = 1.0e0 /(end_value - schoenberg_points(index)) ;
725
726       for (ii = index ; ii <= end_num_poles ; ii++) {
727         value = end_value - schoenberg_points(ii) ;
728         value *= divide ;
729         first_function(ii) = 1.0e0 ;
730
731         for (jj = 0 ; jj < type ; jj++) {
732           first_function(ii) *= value ;
733         }
734       }
735       
736       divide = 1.0e0 / (schoenberg_points(other_index) - start_value) ;
737
738       for (ii = start_num_poles  ; ii <= other_index ; ii++) {
739         value = schoenberg_points(ii) - start_value ;
740         value *= divide ;
741         second_function(ii) = 1.0e0 ;
742
743         for (jj = 0 ; jj < type ; jj++) {
744           second_function(ii) *= value ;
745         }
746       }
747       divide = 1.0e0/( end_value - schoenberg_points(other_index)) ;
748
749       for (ii = other_index ; ii <= end_num_poles ; ii++) {
750         value = end_value - schoenberg_points(ii) ;
751         value *= divide ;
752         second_function(ii) = 1.0e0 ;
753
754         for (jj = 0 ; jj < type ; jj++) {
755           second_function(ii) *= value ;
756         }
757       }
758
759       //
760       //  compute the point and derivatives of both functions
761       //    
762       Standard_Real results[2][2],
763       weights_results[2][2];
764       Standard_Integer extrap_mode[2],
765       derivative_request = 1,
766       dimension = 1 ;
767       Standard_Boolean periodic_flag = Standard_False ;
768       
769       extrap_mode[0] = Degree ;
770       extrap_mode[1] = Degree ;
771       if (Rational) {
772         //
773         // evaluate in homogenised form
774         //
775         Eval(U,
776              periodic_flag,
777              derivative_request,
778              extrap_mode[0],
779              Degree,
780              FlatKnots,
781              dimension,
782              first_function(1),
783              weights_array[0],
784              results[0][0],
785              weights_results[0][0]) ;
786         
787         Eval(U,
788              periodic_flag,
789              derivative_request,
790              extrap_mode[0],
791              Degree,
792              FlatKnots,
793              dimension,
794              second_function(1),
795              weights_array[0],
796              results[1][0],
797              weights_results[1][0]) ;
798         //
799         //  compute the rational derivatives values
800         //       
801
802         for (ii = 0 ; ii < 2 ; ii++) {
803           PLib::RationalDerivatives(1,
804                                     1,
805                                     results[ii][0],
806                                     weights_results[ii][0],
807                                     results[ii][0]) ;
808         }
809       }
810       else {
811         Eval(U,
812              Standard_False,
813              1,
814              extrap_mode[0],
815              Degree,
816              FlatKnots,
817              1,
818              first_function(1),
819              results[0][0]) ;
820         
821         Eval(U,
822              Standard_False,
823              1,
824              extrap_mode[0],
825              Degree,
826              FlatKnots,
827              1,
828              second_function(1),
829              results[1][0]) ;
830       }
831       gp_Mat2d  a_matrix ;
832
833       for (ii = 0 ; ii < 2 ; ii++) {
834
835         for (jj = 0 ; jj < 2 ; jj++) {
836           a_matrix.SetValue(ii+1,jj+1,results[ii][jj]) ;
837         }
838       }
839       a_matrix.Invert() ;
840       TColStd_Array1OfReal the_a_vector(0,ArrayDimension-1) ;
841       TColStd_Array1OfReal the_b_vector(0,ArrayDimension-1) ;
842
843       for( ii = 0 ; ii < ArrayDimension ; ii++) {
844         the_a_vector(ii) = 
845           a_matrix.Value(1,1) * delta_array[ii] +
846             a_matrix.Value(2,1) * derivatives_array[ii] ;
847         the_b_vector(ii) =
848           a_matrix.Value(1,2) * delta_array[ii] +
849             a_matrix.Value(2,2) * derivatives_array[ii] ;
850       }
851       index = 0 ;
852       
853       for (ii = 0 ; ii < num_poles ; ii++) {
854         
855         for (jj = 0 ; jj < ArrayDimension ; jj++) {
856           new_poles_array[index] = poles_array[index] ; 
857           new_poles_array[index] += 
858             first_function(ii+1) * the_a_vector(jj) ;
859           new_poles_array[index] += 
860             second_function(ii+1) * the_b_vector(jj) ;
861           index += 1 ;
862         }
863       }
864     }
865     else {
866       ErrorStatus = 1 ;
867     }
868   }
869   else {
870     ErrorStatus = 2 ;
871   }
872 }
873
874 //=======================================================================
875 //function : FunctionMultiply
876 //purpose  : 
877 //=======================================================================
878
879 void BSplCLib::FunctionMultiply
880 (const BSplCLib_EvaluatorFunction & FunctionPtr,
881  const Standard_Integer             BSplineDegree,
882  const TColStd_Array1OfReal &       BSplineFlatKnots,
883  const Standard_Integer             PolesDimension,
884  Standard_Real &                    Poles,
885  const TColStd_Array1OfReal &       FlatKnots,
886  const Standard_Integer             NewDegree,
887  Standard_Real &                    NewPoles,
888  Standard_Integer &                 Status) 
889 {
890   Standard_Integer ii,
891   jj,
892   index ;
893   Standard_Integer extrap_mode[2],
894   error_code,
895   num_new_poles,
896   derivative_request = 0 ;
897   Standard_Boolean  periodic_flag = Standard_False ;
898   Standard_Real  result,
899   start_end[2],
900   *array_of_poles,
901   *array_of_new_poles ;
902   
903   array_of_poles = (Standard_Real *) &NewPoles ;
904   extrap_mode[0] = 
905     extrap_mode[1] = BSplineDegree ;
906   num_new_poles =
907     FlatKnots.Length() - NewDegree - 1 ;
908   start_end[0] = FlatKnots(NewDegree+1) ;
909   start_end[1] = FlatKnots(num_new_poles+1) ;
910   TColStd_Array1OfReal  parameters(1,num_new_poles) ;
911   TColStd_Array1OfInteger contact_order_array(1,num_new_poles) ;
912   TColStd_Array1OfReal  new_poles_array(1, num_new_poles * PolesDimension) ;
913   
914   array_of_new_poles = 
915     (Standard_Real *) &new_poles_array(1) ;
916   BuildSchoenbergPoints(NewDegree,
917                         FlatKnots,
918                         parameters) ;
919   //
920   // on recadre sur les bornes
921   // 
922   if (parameters(1) < start_end[0]) {
923     parameters(1) = start_end[0] ;
924   }
925   if (parameters(num_new_poles) > start_end[1]) {
926     parameters(num_new_poles) = start_end[1] ;
927   }
928   index = 0 ; 
929
930   for (ii = 1 ; ii <= num_new_poles ; ii++) {
931     contact_order_array(ii) = 0 ;
932     FunctionPtr.Evaluate (contact_order_array(ii),
933                    start_end,
934                    parameters(ii),
935                    result,
936                    error_code);
937     if (error_code) {
938       Status = 1 ;
939       goto FINISH ;
940     }     
941     
942     Eval(parameters(ii),
943          periodic_flag,
944          derivative_request,
945          extrap_mode[0],
946          BSplineDegree,
947          BSplineFlatKnots,
948          PolesDimension,
949          Poles,
950          array_of_new_poles[index]) ;
951     
952     for (jj = 0 ; jj < PolesDimension ; jj++) {
953       array_of_new_poles[index] *= result ;
954       index += 1 ;
955     }
956   }
957   Interpolate(NewDegree,
958               FlatKnots,
959               parameters,
960               contact_order_array,
961               PolesDimension,
962               array_of_new_poles[0],
963               Status) ;
964   
965   for (ii = 0 ; ii < num_new_poles * PolesDimension ; ii++) {
966     array_of_poles[ii] = array_of_new_poles[ii] ;
967     
968   }
969   FINISH :
970     ;
971 }
972
973 //=======================================================================
974 // function : FunctionMultiply
975 //purpose  : 
976 //=======================================================================
977
978 void BSplCLib::FunctionReparameterise
979 (const BSplCLib_EvaluatorFunction & FunctionPtr,
980  const Standard_Integer             BSplineDegree,
981  const TColStd_Array1OfReal &       BSplineFlatKnots,
982  const Standard_Integer             PolesDimension,
983  Standard_Real &                    Poles,
984  const TColStd_Array1OfReal &       FlatKnots,
985  const Standard_Integer             NewDegree,
986  Standard_Real &                    NewPoles,
987  Standard_Integer &                 Status) 
988 {
989   Standard_Integer ii,
990 //  jj,
991   index ;
992   Standard_Integer extrap_mode[2],
993   error_code,
994   num_new_poles,
995   derivative_request = 0 ;
996   Standard_Boolean  periodic_flag = Standard_False ;
997   Standard_Real  result,
998   start_end[2],
999   *array_of_poles,
1000   *array_of_new_poles ;
1001   
1002   array_of_poles = (Standard_Real *) &NewPoles ;
1003   extrap_mode[0] = 
1004     extrap_mode[1] = BSplineDegree ;
1005   num_new_poles =
1006     FlatKnots.Length() - NewDegree - 1 ;
1007   start_end[0] = FlatKnots(NewDegree+1) ;
1008   start_end[1] = FlatKnots(num_new_poles+1) ;
1009   TColStd_Array1OfReal  parameters(1,num_new_poles) ;
1010   TColStd_Array1OfInteger contact_order_array(1,num_new_poles) ;
1011   TColStd_Array1OfReal  new_poles_array(1, num_new_poles * PolesDimension) ;
1012   
1013   array_of_new_poles = 
1014     (Standard_Real *) &new_poles_array(1) ;
1015   BuildSchoenbergPoints(NewDegree,
1016                         FlatKnots,
1017                         parameters) ;
1018   index = 0 ; 
1019
1020   for (ii = 1 ; ii <= num_new_poles ; ii++) {
1021     contact_order_array(ii) = 0 ;
1022     FunctionPtr.Evaluate (contact_order_array(ii),
1023                    start_end,
1024                    parameters(ii),
1025                    result,
1026                    error_code);
1027     if (error_code) {
1028       Status = 1 ;
1029       goto FINISH ;
1030     }     
1031     
1032     Eval(result,
1033          periodic_flag,
1034          derivative_request,
1035          extrap_mode[0],
1036          BSplineDegree,
1037          BSplineFlatKnots,
1038          PolesDimension,
1039          Poles,
1040          array_of_new_poles[index]) ;
1041     index += PolesDimension ;
1042   }
1043   Interpolate(NewDegree,
1044               FlatKnots,
1045               parameters,
1046               contact_order_array,
1047               PolesDimension,
1048               array_of_new_poles[0],
1049               Status) ;
1050   
1051   for (ii = 0 ; ii < num_new_poles * PolesDimension ; ii++) {
1052     array_of_poles[ii] = array_of_new_poles[ii] ;
1053     
1054   }
1055   FINISH :
1056     ;
1057 }
1058
1059 //=======================================================================
1060 //function : FunctionMultiply
1061 //purpose  : 
1062 //=======================================================================
1063
1064 void BSplCLib::FunctionMultiply
1065 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1066  const Standard_Integer              BSplineDegree,
1067  const TColStd_Array1OfReal &        BSplineFlatKnots,
1068  const TColStd_Array1OfReal &        Poles,
1069  const TColStd_Array1OfReal &        FlatKnots,
1070  const Standard_Integer              NewDegree,
1071  TColStd_Array1OfReal &              NewPoles,
1072  Standard_Integer &                  Status) 
1073
1074   Standard_Integer num_bspline_poles =
1075     BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1076   Standard_Integer num_new_poles =
1077     FlatKnots.Length() - NewDegree - 1 ;
1078   
1079   if (Poles.Length() != num_bspline_poles ||
1080       NewPoles.Length() != num_new_poles) {
1081     Standard_ConstructionError::Raise();  
1082   }
1083   Standard_Real  * array_of_poles =
1084     (Standard_Real *) &Poles(Poles.Lower()) ;
1085   Standard_Real  * array_of_new_poles =
1086     (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1087   BSplCLib::FunctionMultiply(FunctionPtr,
1088                              BSplineDegree,
1089                              BSplineFlatKnots,
1090                              1,
1091                              array_of_poles[0],
1092                              FlatKnots,
1093                              NewDegree,
1094                              array_of_new_poles[0],
1095                              Status) ;
1096 }
1097
1098 //=======================================================================
1099 //function : FunctionReparameterise
1100 //purpose  : 
1101 //=======================================================================
1102
1103 void BSplCLib::FunctionReparameterise
1104 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1105  const Standard_Integer              BSplineDegree,
1106  const TColStd_Array1OfReal &        BSplineFlatKnots,
1107  const TColStd_Array1OfReal &        Poles,
1108  const TColStd_Array1OfReal &        FlatKnots,
1109  const Standard_Integer              NewDegree,
1110  TColStd_Array1OfReal &              NewPoles,
1111  Standard_Integer &                  Status) 
1112
1113   Standard_Integer num_bspline_poles =
1114     BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1115   Standard_Integer num_new_poles =
1116     FlatKnots.Length() - NewDegree - 1 ;
1117   
1118   if (Poles.Length() != num_bspline_poles ||
1119       NewPoles.Length() != num_new_poles) {
1120     Standard_ConstructionError::Raise();  
1121   }
1122   Standard_Real  * array_of_poles =
1123     (Standard_Real *) &Poles(Poles.Lower()) ;
1124   Standard_Real  * array_of_new_poles =
1125     (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1126   BSplCLib::FunctionReparameterise(
1127                                    FunctionPtr,
1128                                    BSplineDegree,
1129                                    BSplineFlatKnots,
1130                                    1,
1131                                    array_of_poles[0],
1132                                    FlatKnots,
1133                                    NewDegree,
1134                                    array_of_new_poles[0],
1135                                    Status) ;
1136 }
1137
1138 //=======================================================================
1139 //function : MergeBSplineKnots
1140 //purpose  : 
1141 //=======================================================================
1142 void BSplCLib::MergeBSplineKnots
1143 (const Standard_Real                 Tolerance,
1144  const Standard_Real                 StartValue,
1145  const Standard_Real                 EndValue,
1146  const Standard_Integer              Degree1,
1147  const TColStd_Array1OfReal          & Knots1,
1148  const TColStd_Array1OfInteger       & Mults1,
1149  const Standard_Integer              Degree2,
1150  const TColStd_Array1OfReal          & Knots2,
1151  const TColStd_Array1OfInteger       & Mults2,
1152  Standard_Integer &                  NumPoles,
1153  Handle(TColStd_HArray1OfReal) &     NewKnots,
1154  Handle(TColStd_HArray1OfInteger) &  NewMults) 
1155 {
1156   Standard_Integer ii,
1157   jj,
1158   continuity,
1159   set_mults_flag,
1160   degree,
1161   index,
1162   num_knots ;
1163   if (StartValue < EndValue - Tolerance) {
1164     TColStd_Array1OfReal knots1(1,Knots1.Length()) ;
1165     TColStd_Array1OfReal knots2(1,Knots2.Length()) ;
1166     degree = Degree1 + Degree2 ;
1167     index = 1 ;
1168
1169     for (ii = Knots1.Lower() ; ii <= Knots1.Upper() ; ii++) {
1170       knots1(index) = Knots1(ii) ;
1171       index += 1 ;
1172     }
1173     index = 1 ;
1174
1175     for (ii = Knots2.Lower() ; ii <= Knots2.Upper() ; ii++) {
1176       knots2(index) = Knots2(ii) ;
1177       index += 1 ;
1178     }
1179     BSplCLib::Reparametrize(StartValue,
1180                             EndValue,
1181                             knots1) ;
1182     
1183     BSplCLib::Reparametrize(StartValue,
1184                             EndValue,
1185                             knots2) ;
1186     num_knots = 0 ;
1187     jj =  1 ;
1188
1189     for (ii = 1 ; ii <= knots1.Length() ; ii++) {
1190
1191       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) - Tolerance) {
1192         jj += 1 ;
1193         num_knots += 1 ;
1194       }
1195
1196       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) + Tolerance) {
1197         jj += 1 ;
1198       }
1199       num_knots += 1 ;
1200     }
1201     NewKnots = 
1202       new TColStd_HArray1OfReal(1,num_knots) ;
1203     NewMults =
1204       new TColStd_HArray1OfInteger(1,num_knots) ;
1205     num_knots = 1 ;
1206     jj = 1 ;
1207
1208     for (ii = 1 ; ii <= knots1.Length() ; ii++) {
1209
1210       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) - Tolerance) {
1211         NewKnots->ChangeArray1()(num_knots) = knots2(jj) ;
1212         NewMults->ChangeArray1()(num_knots) = Mults2(jj) + Degree1 ;
1213         jj += 1 ;
1214         num_knots += 1 ;
1215       }
1216       set_mults_flag = 0 ;
1217
1218       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) + Tolerance) {
1219         continuity = Min(Degree1 - Mults1(ii), Degree2 - Mults2(jj)) ;
1220         set_mults_flag = 1 ;
1221         NewMults->ChangeArray1()(num_knots) = degree - continuity ;
1222         jj += 1 ;
1223       }
1224
1225       NewKnots->ChangeArray1()(num_knots) = knots1(ii) ;
1226       if (! set_mults_flag) {
1227         NewMults->ChangeArray1()(num_knots) = Mults1(ii) + Degree2 ;
1228       }
1229       num_knots += 1 ;
1230     }
1231     num_knots -= 1 ;
1232     NewMults->ChangeArray1()(1) = degree + 1 ;
1233     NewMults->ChangeArray1()(num_knots) = degree + 1 ;
1234     index = 0 ;
1235
1236     for (ii = 1 ; ii <= num_knots ; ii++) {
1237       index += NewMults->Value(ii) ;
1238     }
1239     NumPoles = index  - degree - 1  ;
1240   }
1241 }
1242