f5c421640f148e7a8ced775706026065abe3df8b
[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 //=======================================================================
35 //struct : BSplCLib_DataContainer 
36 //purpose: Auxiliary structure providing buffers for poles and knots used in
37 //         evaluation of bspline (allocated in the stack)
38 //=======================================================================
39
40 struct BSplCLib_DataContainer 
41 {
42   BSplCLib_DataContainer(Standard_Integer Degree)
43   {
44     (void)Degree; // avoid compiler warning
45     Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree() ||
46         BSplCLib::MaxDegree() > 25,
47         "BSplCLib: bspline degree is greater than maximum supported");
48   }
49
50   Standard_Real poles[2*(25+1)];
51   Standard_Real knots[2*25];
52   Standard_Real ders[4];
53 };
54
55 // methods for 1 dimensional BSplines
56
57 //=======================================================================
58 //function : BuildEval
59 //purpose  : builds the local array for evaluation
60 //=======================================================================
61
62 void  BSplCLib::BuildEval(const Standard_Integer         Degree,
63                           const Standard_Integer         Index,
64                           const TColStd_Array1OfReal&    Poles, 
65                           const TColStd_Array1OfReal&    Weights,
66                           Standard_Real&                 LP)
67 {
68   Standard_Integer PLower = Poles.Lower();
69   Standard_Integer PUpper = Poles.Upper();
70   Standard_Integer i;
71   Standard_Integer ip = PLower + Index - 1;
72   Standard_Real w, *pole = &LP;
73   if (&Weights == NULL) {
74     
75     for (i = 0; i <= Degree; i++) {
76       ip++;
77       if (ip > PUpper) ip = PLower;
78       pole[0] = Poles(ip);
79       pole += 1;
80     }
81   }
82   else {
83     
84     for (i = 0; i <= Degree; i++) {
85       ip++;
86       if (ip > PUpper) ip = PLower;
87       pole[1] = w = Weights(ip);
88       pole[0] = Poles(ip) * w;
89       pole += 2;
90     }
91   }
92 }
93
94 //=======================================================================
95 //function : PrepareEval
96 //purpose  : stores data for Eval in the local arrays
97 //           dc.poles and dc.knots
98 //=======================================================================
99
100 static void PrepareEval
101 (Standard_Real&                 u,                  
102  Standard_Integer&              index, 
103  Standard_Integer&              dim,
104  Standard_Boolean&              rational,
105  const Standard_Integer         Degree,     
106  const Standard_Boolean         Periodic,
107  const TColStd_Array1OfReal&    Poles,  
108  const TColStd_Array1OfReal&    Weights,
109  const TColStd_Array1OfReal&    Knots,  
110  const TColStd_Array1OfInteger& Mults,
111  BSplCLib_DataContainer&        dc) 
112 {                    
113   // Set the Index
114   BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
115   
116   // make the knots
117   BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*dc.knots);
118   if (&Mults == NULL)
119     index -= Knots.Lower() + Degree;
120   else
121     index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults);
122   
123   // check truly rational
124   rational = (&Weights != NULL);
125   if (rational) {
126     Standard_Integer WLower = Weights.Lower() + index;
127     rational = BSplCLib::IsRational(Weights, WLower, WLower + Degree);
128   }
129   
130   // make the poles
131   if(rational) {
132     dim = 2;
133     BSplCLib::BuildEval(Degree, index, Poles, Weights              , *dc.poles);
134   }
135   else {
136     dim = 1;
137     BSplCLib::BuildEval(Degree, index, Poles, BSplCLib::NoWeights(), *dc.poles);
138   }
139 }
140
141 //=======================================================================
142 //function : D0
143 //purpose  : 
144 //=======================================================================
145
146 void BSplCLib::D0 
147 (const Standard_Real            U,                  
148  const Standard_Integer         Index,          
149  const Standard_Integer         Degree,     
150  const Standard_Boolean         Periodic,
151  const TColStd_Array1OfReal&    Poles,  
152  const TColStd_Array1OfReal&    Weights,
153  const TColStd_Array1OfReal&    Knots,  
154  const TColStd_Array1OfInteger& Mults,  
155  Standard_Real&                 P) 
156 {                    
157   Standard_Integer dim,index = Index;
158   Standard_Real    u = U;
159   Standard_Boolean rational;
160   BSplCLib_DataContainer dc(Degree);
161   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
162   BSplCLib::Eval(u,Degree,*dc.knots,dim,*dc.poles);
163   if (rational) P = dc.poles[0] / dc.poles[1];
164   else          P = dc.poles[0];
165 }
166
167 //=======================================================================
168 //function : D1
169 //purpose  : 
170 //=======================================================================
171
172 void BSplCLib::D1
173 (const Standard_Real            U,                  
174  const Standard_Integer         Index,          
175  const Standard_Integer         Degree,     
176  const Standard_Boolean         Periodic,
177  const TColStd_Array1OfReal&    Poles,  
178  const TColStd_Array1OfReal&    Weights,
179  const TColStd_Array1OfReal&    Knots,  
180  const TColStd_Array1OfInteger& Mults,  
181  Standard_Real&                 P,
182  Standard_Real&                 V) 
183 {                    
184   Standard_Integer dim,index = Index;
185   Standard_Real    u = U;
186   Standard_Boolean rational;
187   BSplCLib_DataContainer dc(Degree);
188   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
189   BSplCLib::Bohm(u,Degree,1,*dc.knots,dim,*dc.poles);
190   Standard_Real *result = dc.poles;
191   if (rational) {
192     PLib::RationalDerivative(Degree,1,1,*dc.poles,*dc.ders);
193     result = dc.ders;
194   }
195   P = result[0];
196   V = result[1];
197 }
198
199 //=======================================================================
200 //function : D2
201 //purpose  : 
202 //=======================================================================
203
204 void BSplCLib::D2
205 (const Standard_Real            U,                  
206  const Standard_Integer         Index,          
207  const Standard_Integer         Degree,     
208  const Standard_Boolean         Periodic,
209  const TColStd_Array1OfReal&    Poles,  
210  const TColStd_Array1OfReal&    Weights,
211  const TColStd_Array1OfReal&    Knots,  
212  const TColStd_Array1OfInteger& Mults,  
213  Standard_Real&                 P,
214  Standard_Real&                 V1,
215  Standard_Real&                 V2) 
216 {                    
217   Standard_Integer dim,index = Index;
218   Standard_Real    u = U;
219   Standard_Boolean rational;
220   BSplCLib_DataContainer dc(Degree);
221   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
222   BSplCLib::Bohm(u,Degree,2,*dc.knots,dim,*dc.poles);
223   Standard_Real *result = dc.poles;
224   if (rational) {
225     PLib::RationalDerivative(Degree,2,1,*dc.poles,*dc.ders);
226     result = dc.ders;
227   }
228   P  = result[0];
229   V1 = result[1];
230   if (!rational && (Degree < 2)) V2 = 0.;
231   else                           V2 = result[2];
232 }
233
234 //=======================================================================
235 //function : D3
236 //purpose  : 
237 //=======================================================================
238
239 void BSplCLib::D3
240 (const Standard_Real            U,                  
241  const Standard_Integer         Index,          
242  const Standard_Integer         Degree,     
243  const Standard_Boolean         Periodic,
244  const TColStd_Array1OfReal&    Poles,  
245  const TColStd_Array1OfReal&    Weights,
246  const TColStd_Array1OfReal&    Knots,  
247  const TColStd_Array1OfInteger& Mults,  
248  Standard_Real&                 P,
249  Standard_Real&                 V1,
250  Standard_Real&                 V2,
251  Standard_Real&                 V3) 
252 {                    
253   Standard_Integer dim,index = Index;
254   Standard_Real    u = U;
255   Standard_Boolean rational;
256   BSplCLib_DataContainer dc(Degree);
257   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
258   BSplCLib::Bohm(u,Degree,3,*dc.knots,dim,*dc.poles);
259   Standard_Real *result = dc.poles;
260   if (rational) {
261     PLib::RationalDerivative(Degree,3,1,*dc.poles,*dc.ders);
262     result = dc.ders;
263   }
264   P  = result[0];
265   V1 = result[1];
266   if (!rational && (Degree < 2)) V2 = 0.;
267   else                           V2 = result[2];
268   if (!rational && (Degree < 3)) V3 = 0.;
269   else                           V3 = result[3];
270 }
271
272 //=======================================================================
273 //function : DN
274 //purpose  : 
275 //=======================================================================
276
277 void BSplCLib::DN
278 (const Standard_Real            U,                  
279  const Standard_Integer         N,          
280  const Standard_Integer         Index,          
281  const Standard_Integer         Degree,     
282  const Standard_Boolean         Periodic,
283  const TColStd_Array1OfReal&    Poles,  
284  const TColStd_Array1OfReal&    Weights,
285  const TColStd_Array1OfReal&    Knots,  
286  const TColStd_Array1OfInteger& Mults,  
287  Standard_Real&                 VN) 
288 {                    
289   Standard_Integer dim,index = Index;
290   Standard_Real    u = U;
291   Standard_Boolean rational;
292   BSplCLib_DataContainer dc(Degree);
293   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
294   BSplCLib::Bohm(u,Degree,N,*dc.knots,dim,*dc.poles);
295   if (rational) {
296     Standard_Real v;
297     PLib::RationalDerivative(Degree,N,1,*dc.poles,v,Standard_False);
298     VN = v;
299   }
300   else {
301     if (N > Degree) VN = 0.;
302     else            VN = dc.poles[N];
303   }
304 }
305
306 //=======================================================================
307 //function : Build BSpline Matrix
308 //purpose  : Builds the Bspline Matrix
309 //=======================================================================
310
311 Standard_Integer 
312 BSplCLib::BuildBSpMatrix(const  TColStd_Array1OfReal&     Parameters,
313                          const  TColStd_Array1OfInteger&  ContactOrderArray,
314                          const  TColStd_Array1OfReal&     FlatKnots,
315                          const  Standard_Integer          Degree,   
316                          math_Matrix&                     Matrix,
317                          Standard_Integer&                UpperBandWidth,
318                          Standard_Integer&                LowerBandWidth) 
319 {
320   Standard_Integer ii,
321   jj,
322   Index,
323   ErrorCode,
324   ReturnCode = 0,
325   FirstNonZeroBsplineIndex,
326   BandWidth,
327   MaxOrder = 21,
328   Order ;
329   
330   math_Matrix   BSplineBasis(1, MaxOrder,
331                              1, MaxOrder) ;
332   
333   Order = Degree + 1 ;
334   UpperBandWidth = Degree ;
335   LowerBandWidth = Degree ;
336   BandWidth = UpperBandWidth + LowerBandWidth + 1 ;
337   if (Matrix.LowerRow() != Parameters.Lower()  || 
338       Matrix.UpperRow() != Parameters.Upper()  ||
339       Matrix.LowerCol() != 1  || 
340       Matrix.UpperCol() != BandWidth) {
341     ReturnCode = 1;
342     goto FINISH ;
343   }
344   
345   for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
346     ErrorCode =
347       BSplCLib::EvalBsplineBasis(1,
348                                  ContactOrderArray(ii),
349                                  Order,
350                                  FlatKnots,
351                                  Parameters(ii),
352                                  
353                                  FirstNonZeroBsplineIndex,
354                                  BSplineBasis) ;
355     if (ErrorCode != 0) {
356       ReturnCode = 2 ;
357       goto FINISH ;
358     }
359     Index = LowerBandWidth + 1 + FirstNonZeroBsplineIndex - ii ;
360
361     for (jj = 1 ; jj < Index ; jj++) {
362       Matrix.Value(ii,jj) = 0.0e0 ;
363     }
364
365     for (jj = 1 ; jj <= Order ; jj++) {
366       Matrix.Value(ii,Index) = BSplineBasis(ContactOrderArray(ii) + 1, jj) ;
367       Index += 1 ;
368     }
369     
370     for (jj = Index ; jj <= BandWidth ; jj++) {
371       Matrix.Value(ii,jj) = 0.0e0 ;
372     }
373   } 
374   FINISH : ;
375   return (ReturnCode) ;
376 }
377
378 //=======================================================================
379 //function : Makes LU decompositiomn without Pivoting
380 //purpose  : Builds the Bspline Matrix
381 //=======================================================================
382
383 Standard_Integer 
384 BSplCLib::FactorBandedMatrix(math_Matrix&   Matrix,
385                              const Standard_Integer UpperBandWidth,
386                              const Standard_Integer LowerBandWidth,
387                              Standard_Integer&  PivotIndexProblem) 
388 {
389   Standard_Integer ii,
390   jj,
391   kk,
392   Index,
393   MinIndex,
394   MaxIndex,
395   ReturnCode = 0,
396   BandWidth = UpperBandWidth + LowerBandWidth + 1 ;
397   
398   Standard_Real Inverse ;
399   PivotIndexProblem = 0 ;
400
401   for (ii = Matrix.LowerRow() + 1 ; ii <= Matrix.UpperRow() ; ii++) {
402     MinIndex = ( LowerBandWidth - ii + 2 >= 1 ? LowerBandWidth - ii + 2 : 1) ;
403
404     for (jj = MinIndex ; jj <= LowerBandWidth  ; jj++) {
405       Index = ii - LowerBandWidth + jj - 1 ; 
406       Inverse = Matrix(Index ,LowerBandWidth + 1) ;
407       if (Abs(Inverse) > RealSmall()) {
408         Inverse = -1.0e0/Inverse ;
409       }
410       else {
411         ReturnCode = 1 ;
412         PivotIndexProblem = Index ;
413         goto FINISH ;
414       }
415       Matrix(ii,jj) = Matrix(ii,jj) * Inverse ;
416       MaxIndex = BandWidth + Index - ii ;
417
418       for (kk = jj + 1 ; kk <= MaxIndex ; kk++) {
419         Matrix(ii,kk) += Matrix(ii,jj) * Matrix(Index, kk + ii - Index) ;
420       }
421     }                       
422   }
423   FINISH :
424     return (ReturnCode) ;
425 }
426
427 //=======================================================================
428 //function : Build BSpline Matrix
429 //purpose  : Builds the Bspline Matrix
430 //=======================================================================
431
432 Standard_Integer 
433 BSplCLib::EvalBsplineBasis
434 //(const Standard_Integer              Side, // = 1 rigth side, -1 left side 
435 (const Standard_Integer              , // = 1 rigth side, -1 left side 
436  const  Standard_Integer              DerivativeRequest,
437  const  Standard_Integer              Order,
438  const  TColStd_Array1OfReal&         FlatKnots,
439  const  Standard_Real                 Parameter,
440  Standard_Integer&             FirstNonZeroBsplineIndex,
441  math_Matrix&                  BsplineBasis)
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                             Standard_False, 
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