0026431: Can't cut a sphere from a cylinder
[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(1,
350                                  ContactOrderArray(ii),
351                                  Order,
352                                  FlatKnots,
353                                  Parameters(ii),
354                                  
355                                  FirstNonZeroBsplineIndex,
356                                  BSplineBasis) ;
357     if (ErrorCode != 0) {
358       ReturnCode = 2 ;
359       goto FINISH ;
360     }
361     Index = LowerBandWidth + 1 + FirstNonZeroBsplineIndex - ii ;
362
363     for (jj = 1 ; jj < Index ; jj++) {
364       Matrix.Value(ii,jj) = 0.0e0 ;
365     }
366
367     for (jj = 1 ; jj <= Order ; jj++) {
368       Matrix.Value(ii,Index) = BSplineBasis(ContactOrderArray(ii) + 1, jj) ;
369       Index += 1 ;
370     }
371     
372     for (jj = Index ; jj <= BandWidth ; jj++) {
373       Matrix.Value(ii,jj) = 0.0e0 ;
374     }
375   } 
376   FINISH : ;
377   return (ReturnCode) ;
378 }
379
380 //=======================================================================
381 //function : Makes LU decompositiomn without Pivoting
382 //purpose  : Builds the Bspline Matrix
383 //=======================================================================
384
385 Standard_Integer 
386 BSplCLib::FactorBandedMatrix(math_Matrix&   Matrix,
387                              const Standard_Integer UpperBandWidth,
388                              const Standard_Integer LowerBandWidth,
389                              Standard_Integer&  PivotIndexProblem) 
390 {
391   Standard_Integer ii,
392   jj,
393   kk,
394   Index,
395   MinIndex,
396   MaxIndex,
397   ReturnCode = 0,
398   BandWidth = UpperBandWidth + LowerBandWidth + 1 ;
399   
400   Standard_Real Inverse ;
401   PivotIndexProblem = 0 ;
402
403   for (ii = Matrix.LowerRow() + 1 ; ii <= Matrix.UpperRow() ; ii++) {
404     MinIndex = ( LowerBandWidth - ii + 2 >= 1 ? LowerBandWidth - ii + 2 : 1) ;
405
406     for (jj = MinIndex ; jj <= LowerBandWidth  ; jj++) {
407       Index = ii - LowerBandWidth + jj - 1 ; 
408       Inverse = Matrix(Index ,LowerBandWidth + 1) ;
409       if (Abs(Inverse) > RealSmall()) {
410         Inverse = -1.0e0/Inverse ;
411       }
412       else {
413         ReturnCode = 1 ;
414         PivotIndexProblem = Index ;
415         goto FINISH ;
416       }
417       Matrix(ii,jj) = Matrix(ii,jj) * Inverse ;
418       MaxIndex = BandWidth + Index - ii ;
419
420       for (kk = jj + 1 ; kk <= MaxIndex ; kk++) {
421         Matrix(ii,kk) += Matrix(ii,jj) * Matrix(Index, kk + ii - Index) ;
422       }
423     }                       
424   }
425   FINISH :
426     return (ReturnCode) ;
427 }
428
429 //=======================================================================
430 //function : Build BSpline Matrix
431 //purpose  : Builds the Bspline Matrix
432 //=======================================================================
433
434 Standard_Integer 
435 BSplCLib::EvalBsplineBasis
436 //(const Standard_Integer              Side, // = 1 rigth side, -1 left side 
437 (const Standard_Integer              , // = 1 rigth side, -1 left side 
438  const  Standard_Integer              DerivativeRequest,
439  const  Standard_Integer              Order,
440  const  TColStd_Array1OfReal&         FlatKnots,
441  const  Standard_Real                 Parameter,
442  Standard_Integer&             FirstNonZeroBsplineIndex,
443  math_Matrix&                  BsplineBasis,
444  Standard_Boolean              isPeriodic)
445 {
446   // the matrix must have at least DerivativeRequest + 1
447   //   row and Order columns
448   // the result are stored in the following way in
449   // the Bspline matrix 
450   // Let i be the FirstNonZeroBsplineIndex and 
451   // t be the parameter value, k the order of the 
452   // knot vector, r the DerivativeRequest :
453   //   
454   //   B (t)   B (t)                     B (t)
455   //    i       i+1                       i+k-1
456   //   
457   //    (1)     (1)                       (1) 
458   //   B (t)   B (t)                     B (t)
459   //    i       i+1                       i+k-1
460   //  
461   //
462   //
463   //
464   //    (r)     (r)                       (r) 
465   //   B (t)   B (t)                     B (t)
466   //    i       i+1                       i+k-1
467   //
468   Standard_Integer  
469     ReturnCode,
470   ii,
471   pp,
472   qq,
473   ss,
474   NumPoles,
475   LocalRequest ;
476 //  ,Index ;
477   
478   Standard_Real NewParameter,
479   Inverse,
480   Factor,
481   LocalInverse,
482   Saved ;
483 // , *FlatKnotsArray ;
484   
485   ReturnCode = 0 ;
486   FirstNonZeroBsplineIndex = 0 ;
487   LocalRequest = DerivativeRequest ;
488   if (DerivativeRequest >= Order) {
489     LocalRequest = Order - 1 ;
490   }
491   
492   if (BsplineBasis.LowerCol() != 1    ||
493       BsplineBasis.UpperCol() < Order ||
494       BsplineBasis.LowerRow() != 1    ||
495       BsplineBasis.UpperRow() <= LocalRequest) {
496     ReturnCode = 1;
497     goto FINISH ;
498   }
499   NumPoles = FlatKnots.Upper() - FlatKnots.Lower() + 1 - Order ;
500   BSplCLib::LocateParameter(Order - 1, 
501                             FlatKnots,
502                             Parameter,
503                             isPeriodic, 
504                             Order, 
505                             NumPoles+1, 
506                             ii,
507                             NewParameter) ;
508   
509   FirstNonZeroBsplineIndex = ii - Order + 1 ;
510   
511   BsplineBasis(1,1) = 1.0e0 ;
512   LocalRequest = DerivativeRequest ;
513   if (DerivativeRequest >= Order) {
514     LocalRequest = Order - 1 ;
515   }
516
517   for (qq = 2 ; qq <= Order - LocalRequest ; qq++) {
518     BsplineBasis(1,qq) = 0.0e0 ;
519     
520     for (pp = 1 ; pp <= qq - 1 ; pp++) {
521       //
522       // this should be always invertible if ii is correctly computed 
523       //
524       Factor = (Parameter - FlatKnots(ii - qq + pp + 1)) 
525         / (FlatKnots(ii + pp)   - FlatKnots(ii - qq + pp + 1)) ; 
526       Saved = Factor *    BsplineBasis(1,pp) ;
527       BsplineBasis(1,pp) *= (1.0e0 - Factor) ;
528       BsplineBasis(1,pp) += BsplineBasis(1,qq) ;
529       BsplineBasis(1,qq) = Saved ;
530     }
531   }
532   
533   for (qq = Order - LocalRequest + 1 ; qq <= Order ; qq++) {
534     
535     for (pp = 1 ; pp <= qq - 1 ; pp++) {
536       BsplineBasis(Order - qq + 2,pp) = BsplineBasis(1,pp) ;
537     }
538     BsplineBasis(1,qq) = 0.0e0 ;
539
540     for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) {
541       BsplineBasis(Order - ss + 2,qq) = 0.0e0 ;
542     }
543
544     for (pp = 1 ; pp <= qq - 1 ; pp++) {
545       Inverse = 1.0e0 / (FlatKnots(ii + pp)  - FlatKnots(ii - qq + pp + 1)) ;
546       Factor  =  (Parameter - FlatKnots(ii - qq + pp + 1)) * Inverse ;
547       Saved = Factor *                 BsplineBasis(1,pp) ;
548       BsplineBasis(1,pp) *= (1.0e0 - Factor) ;
549       BsplineBasis(1,pp) += BsplineBasis(1,qq) ;
550       BsplineBasis(1,qq) = Saved ;
551       LocalInverse = (Standard_Real) (qq - 1) * Inverse ;
552
553       for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) {
554         Saved = LocalInverse * BsplineBasis(Order - ss + 2, pp) ;
555         BsplineBasis(Order - ss + 2, pp) *= - LocalInverse  ;
556         BsplineBasis(Order - ss + 2, pp) +=   BsplineBasis(Order - ss + 2,qq) ;
557         BsplineBasis(Order - ss + 2,qq) = Saved ;
558       }
559     }
560   }
561   FINISH :
562     return (ReturnCode) ;
563 }
564
565 //=======================================================================
566 //function : MovePointAndTangent
567 //purpose  : 
568 //=======================================================================
569
570 void BSplCLib::MovePointAndTangent(const Standard_Real    U,
571                                    const Standard_Integer ArrayDimension,
572                                    Standard_Real    &Delta,
573                                    Standard_Real    &DeltaDerivatives,
574                                    const Standard_Real    Tolerance,
575                                    const Standard_Integer Degree,
576                                    const Standard_Boolean Rational,
577                                    const Standard_Integer StartingCondition,
578                                    const Standard_Integer EndingCondition,
579                                    Standard_Real&         Poles,
580                                    const TColStd_Array1OfReal&   Weights,
581                                    const TColStd_Array1OfReal&   FlatKnots,
582                                    Standard_Real&        NewPoles,
583                                    Standard_Integer&     ErrorStatus) 
584 {
585   Standard_Integer num_poles,
586   num_knots,
587   ii,
588   jj,
589   conditions,
590   start_num_poles,
591   end_num_poles,
592   index,
593   start_index,
594   end_index,
595   other_index,
596   type,
597   order ;
598   
599   Standard_Real    new_parameter,
600   value,
601   divide,
602   end_value,
603   start_value,
604   *poles_array,
605   *new_poles_array,
606   *delta_array,
607   *derivatives_array,
608   *weights_array ;
609   
610   ErrorStatus = 0 ;
611   weights_array = NULL ;
612   if (Rational) {
613     weights_array = (Standard_Real *) &Weights(Weights.Lower()) ;
614   }
615   
616   poles_array = &Poles ;
617   new_poles_array = &NewPoles ;
618   delta_array = &Delta ;
619   derivatives_array = &DeltaDerivatives ;
620   order = Degree + 1 ;
621   num_knots = FlatKnots.Length() ;
622   num_poles = num_knots - order ;
623   conditions = StartingCondition + EndingCondition + 4 ;
624   //
625   // check validity of input data
626   //
627   if (StartingCondition >= -1 &&  
628       StartingCondition <= Degree &&
629       EndingCondition >= -1 &&
630       EndingCondition <= Degree &&
631       conditions <= num_poles) {
632     //
633     // check the parameter is within bounds 
634     //
635     start_index = FlatKnots.Lower() + Degree ;
636     end_index = FlatKnots.Upper() - Degree ;
637     //
638     //  check if there is enough room to move the poles
639     //
640     conditions = 1 ;
641     if (StartingCondition == -1) {
642       conditions = conditions && (FlatKnots(start_index) <=  U) ;
643     }
644     else {
645       conditions = conditions && (FlatKnots(start_index) + Tolerance < U) ;
646     }
647     if (EndingCondition == -1) {
648       conditions = conditions && (FlatKnots(end_index) >= U) ;
649     }
650     else {
651       conditions = conditions && (FlatKnots(end_index) - Tolerance > U) ; 
652     }
653     
654     if (conditions) {
655       //
656       // build 2 auxialiary functions 
657       // 
658       TColStd_Array1OfReal schoenberg_points(1,num_poles) ;
659       TColStd_Array1OfReal first_function  (1,num_poles) ;
660       TColStd_Array1OfReal second_function (1,num_poles) ;
661       
662       BuildSchoenbergPoints(Degree,
663                             FlatKnots,
664                             schoenberg_points) ;
665       start_index = StartingCondition + 2 ;
666       end_index = num_poles - EndingCondition - 1 ;
667       LocateParameter(schoenberg_points,
668                       U,
669                       Standard_False,
670                       start_index,
671                       end_index,
672                       index,
673                       new_parameter, 
674                       0, 1) ;
675       
676       if (index == start_index) {
677         other_index = index + 1 ;
678       } 
679       else if (index == end_index) {
680         other_index = index -1  ;
681       }
682       else if (U - FlatKnots(index) < FlatKnots(index + 1) - U ) {
683         other_index = index - 1 ;
684       }
685       else {
686         other_index = index + 1 ;
687       }
688       type = 3 ;
689       
690       start_num_poles = StartingCondition + 2 ;
691       end_num_poles   = num_poles - EndingCondition - 1 ;
692       if (start_num_poles == 1) {
693         start_value = schoenberg_points(num_poles) - schoenberg_points(1) ;
694         start_value = schoenberg_points(1) - start_value ;
695       }
696       else {
697         start_value = schoenberg_points(start_num_poles - 1) ;
698       }
699       if (end_num_poles == num_poles) {
700         end_value = schoenberg_points(num_poles) - schoenberg_points(1) ;
701         end_value = schoenberg_points(num_poles) + end_value ;
702       }
703       else {
704         end_value = schoenberg_points(end_num_poles + 1) ;
705       }
706       
707       for (ii = 1 ; ii < start_num_poles ; ii++) {
708         first_function(ii) = 0.e0 ;
709         second_function(ii) = 0.0e0 ;
710       }
711
712       for (ii = end_num_poles + 1 ; ii <= num_poles ; ii++) {
713         first_function(ii) = 0.e0 ;
714         second_function(ii) = 0.0e0 ;
715       }
716       divide = 1.0e0 / (schoenberg_points(index) - start_value) ;
717
718       for (ii = start_num_poles  ; ii <= index ; ii++) {
719         value = schoenberg_points(ii) - start_value ;
720         value *= divide ;
721         first_function(ii) = 1.0e0 ;
722         
723         for (jj = 0 ; jj < type ; jj++) {
724           first_function(ii) *= value ;
725         }
726       }
727       divide = 1.0e0 /(end_value - schoenberg_points(index)) ;
728
729       for (ii = index ; ii <= end_num_poles ; ii++) {
730         value = end_value - schoenberg_points(ii) ;
731         value *= divide ;
732         first_function(ii) = 1.0e0 ;
733
734         for (jj = 0 ; jj < type ; jj++) {
735           first_function(ii) *= value ;
736         }
737       }
738       
739       divide = 1.0e0 / (schoenberg_points(other_index) - start_value) ;
740
741       for (ii = start_num_poles  ; ii <= other_index ; ii++) {
742         value = schoenberg_points(ii) - start_value ;
743         value *= divide ;
744         second_function(ii) = 1.0e0 ;
745
746         for (jj = 0 ; jj < type ; jj++) {
747           second_function(ii) *= value ;
748         }
749       }
750       divide = 1.0e0/( end_value - schoenberg_points(other_index)) ;
751
752       for (ii = other_index ; ii <= end_num_poles ; ii++) {
753         value = end_value - schoenberg_points(ii) ;
754         value *= divide ;
755         second_function(ii) = 1.0e0 ;
756
757         for (jj = 0 ; jj < type ; jj++) {
758           second_function(ii) *= value ;
759         }
760       }
761
762       //
763       //  compute the point and derivatives of both functions
764       //    
765       Standard_Real results[2][2],
766       weights_results[2][2];
767       Standard_Integer extrap_mode[2],
768       derivative_request = 1,
769       dimension = 1 ;
770       Standard_Boolean periodic_flag = Standard_False ;
771       
772       extrap_mode[0] = Degree ;
773       extrap_mode[1] = Degree ;
774       if (Rational) {
775         //
776         // evaluate in homogenised form
777         //
778         Eval(U,
779              periodic_flag,
780              derivative_request,
781              extrap_mode[0],
782              Degree,
783              FlatKnots,
784              dimension,
785              first_function(1),
786              weights_array[0],
787              results[0][0],
788              weights_results[0][0]) ;
789         
790         Eval(U,
791              periodic_flag,
792              derivative_request,
793              extrap_mode[0],
794              Degree,
795              FlatKnots,
796              dimension,
797              second_function(1),
798              weights_array[0],
799              results[1][0],
800              weights_results[1][0]) ;
801         //
802         //  compute the rational derivatives values
803         //       
804
805         for (ii = 0 ; ii < 2 ; ii++) {
806           PLib::RationalDerivatives(1,
807                                     1,
808                                     results[ii][0],
809                                     weights_results[ii][0],
810                                     results[ii][0]) ;
811         }
812       }
813       else {
814         Eval(U,
815              Standard_False,
816              1,
817              extrap_mode[0],
818              Degree,
819              FlatKnots,
820              1,
821              first_function(1),
822              results[0][0]) ;
823         
824         Eval(U,
825              Standard_False,
826              1,
827              extrap_mode[0],
828              Degree,
829              FlatKnots,
830              1,
831              second_function(1),
832              results[1][0]) ;
833       }
834       gp_Mat2d  a_matrix ;
835
836       for (ii = 0 ; ii < 2 ; ii++) {
837
838         for (jj = 0 ; jj < 2 ; jj++) {
839           a_matrix.SetValue(ii+1,jj+1,results[ii][jj]) ;
840         }
841       }
842       a_matrix.Invert() ;
843       TColStd_Array1OfReal the_a_vector(0,ArrayDimension-1) ;
844       TColStd_Array1OfReal the_b_vector(0,ArrayDimension-1) ;
845
846       for( ii = 0 ; ii < ArrayDimension ; ii++) {
847         the_a_vector(ii) = 
848           a_matrix.Value(1,1) * delta_array[ii] +
849             a_matrix.Value(2,1) * derivatives_array[ii] ;
850         the_b_vector(ii) =
851           a_matrix.Value(1,2) * delta_array[ii] +
852             a_matrix.Value(2,2) * derivatives_array[ii] ;
853       }
854       index = 0 ;
855       
856       for (ii = 0 ; ii < num_poles ; ii++) {
857         
858         for (jj = 0 ; jj < ArrayDimension ; jj++) {
859           new_poles_array[index] = poles_array[index] ; 
860           new_poles_array[index] += 
861             first_function(ii+1) * the_a_vector(jj) ;
862           new_poles_array[index] += 
863             second_function(ii+1) * the_b_vector(jj) ;
864           index += 1 ;
865         }
866       }
867     }
868     else {
869       ErrorStatus = 1 ;
870     }
871   }
872   else {
873     ErrorStatus = 2 ;
874   }
875 }
876
877 //=======================================================================
878 //function : FunctionMultiply
879 //purpose  : 
880 //=======================================================================
881
882 void BSplCLib::FunctionMultiply
883 (const BSplCLib_EvaluatorFunction & FunctionPtr,
884  const Standard_Integer             BSplineDegree,
885  const TColStd_Array1OfReal &       BSplineFlatKnots,
886  const Standard_Integer             PolesDimension,
887  Standard_Real &                    Poles,
888  const TColStd_Array1OfReal &       FlatKnots,
889  const Standard_Integer             NewDegree,
890  Standard_Real &                    NewPoles,
891  Standard_Integer &                 Status) 
892 {
893   Standard_Integer ii,
894   jj,
895   index ;
896   Standard_Integer extrap_mode[2],
897   error_code,
898   num_new_poles,
899   derivative_request = 0 ;
900   Standard_Boolean  periodic_flag = Standard_False ;
901   Standard_Real  result,
902   start_end[2],
903   *array_of_poles,
904   *array_of_new_poles ;
905   
906   array_of_poles = (Standard_Real *) &NewPoles ;
907   extrap_mode[0] = 
908     extrap_mode[1] = BSplineDegree ;
909   num_new_poles =
910     FlatKnots.Length() - NewDegree - 1 ;
911   start_end[0] = FlatKnots(NewDegree+1) ;
912   start_end[1] = FlatKnots(num_new_poles+1) ;
913   TColStd_Array1OfReal  parameters(1,num_new_poles) ;
914   TColStd_Array1OfInteger contact_order_array(1,num_new_poles) ;
915   TColStd_Array1OfReal  new_poles_array(1, num_new_poles * PolesDimension) ;
916   
917   array_of_new_poles = 
918     (Standard_Real *) &new_poles_array(1) ;
919   BuildSchoenbergPoints(NewDegree,
920                         FlatKnots,
921                         parameters) ;
922   //
923   // on recadre sur les bornes
924   // 
925   if (parameters(1) < start_end[0]) {
926     parameters(1) = start_end[0] ;
927   }
928   if (parameters(num_new_poles) > start_end[1]) {
929     parameters(num_new_poles) = start_end[1] ;
930   }
931   index = 0 ; 
932
933   for (ii = 1 ; ii <= num_new_poles ; ii++) {
934     contact_order_array(ii) = 0 ;
935     FunctionPtr.Evaluate (contact_order_array(ii),
936                    start_end,
937                    parameters(ii),
938                    result,
939                    error_code);
940     if (error_code) {
941       Status = 1 ;
942       goto FINISH ;
943     }     
944     
945     Eval(parameters(ii),
946          periodic_flag,
947          derivative_request,
948          extrap_mode[0],
949          BSplineDegree,
950          BSplineFlatKnots,
951          PolesDimension,
952          Poles,
953          array_of_new_poles[index]) ;
954     
955     for (jj = 0 ; jj < PolesDimension ; jj++) {
956       array_of_new_poles[index] *= result ;
957       index += 1 ;
958     }
959   }
960   Interpolate(NewDegree,
961               FlatKnots,
962               parameters,
963               contact_order_array,
964               PolesDimension,
965               array_of_new_poles[0],
966               Status) ;
967   
968   for (ii = 0 ; ii < num_new_poles * PolesDimension ; ii++) {
969     array_of_poles[ii] = array_of_new_poles[ii] ;
970     
971   }
972   FINISH :
973     ;
974 }
975
976 //=======================================================================
977 // function : FunctionMultiply
978 //purpose  : 
979 //=======================================================================
980
981 void BSplCLib::FunctionReparameterise
982 (const BSplCLib_EvaluatorFunction & FunctionPtr,
983  const Standard_Integer             BSplineDegree,
984  const TColStd_Array1OfReal &       BSplineFlatKnots,
985  const Standard_Integer             PolesDimension,
986  Standard_Real &                    Poles,
987  const TColStd_Array1OfReal &       FlatKnots,
988  const Standard_Integer             NewDegree,
989  Standard_Real &                    NewPoles,
990  Standard_Integer &                 Status) 
991 {
992   Standard_Integer ii,
993 //  jj,
994   index ;
995   Standard_Integer extrap_mode[2],
996   error_code,
997   num_new_poles,
998   derivative_request = 0 ;
999   Standard_Boolean  periodic_flag = Standard_False ;
1000   Standard_Real  result,
1001   start_end[2],
1002   *array_of_poles,
1003   *array_of_new_poles ;
1004   
1005   array_of_poles = (Standard_Real *) &NewPoles ;
1006   extrap_mode[0] = 
1007     extrap_mode[1] = BSplineDegree ;
1008   num_new_poles =
1009     FlatKnots.Length() - NewDegree - 1 ;
1010   start_end[0] = FlatKnots(NewDegree+1) ;
1011   start_end[1] = FlatKnots(num_new_poles+1) ;
1012   TColStd_Array1OfReal  parameters(1,num_new_poles) ;
1013   TColStd_Array1OfInteger contact_order_array(1,num_new_poles) ;
1014   TColStd_Array1OfReal  new_poles_array(1, num_new_poles * PolesDimension) ;
1015   
1016   array_of_new_poles = 
1017     (Standard_Real *) &new_poles_array(1) ;
1018   BuildSchoenbergPoints(NewDegree,
1019                         FlatKnots,
1020                         parameters) ;
1021   index = 0 ; 
1022
1023   for (ii = 1 ; ii <= num_new_poles ; ii++) {
1024     contact_order_array(ii) = 0 ;
1025     FunctionPtr.Evaluate (contact_order_array(ii),
1026                    start_end,
1027                    parameters(ii),
1028                    result,
1029                    error_code);
1030     if (error_code) {
1031       Status = 1 ;
1032       goto FINISH ;
1033     }     
1034     
1035     Eval(result,
1036          periodic_flag,
1037          derivative_request,
1038          extrap_mode[0],
1039          BSplineDegree,
1040          BSplineFlatKnots,
1041          PolesDimension,
1042          Poles,
1043          array_of_new_poles[index]) ;
1044     index += PolesDimension ;
1045   }
1046   Interpolate(NewDegree,
1047               FlatKnots,
1048               parameters,
1049               contact_order_array,
1050               PolesDimension,
1051               array_of_new_poles[0],
1052               Status) ;
1053   
1054   for (ii = 0 ; ii < num_new_poles * PolesDimension ; ii++) {
1055     array_of_poles[ii] = array_of_new_poles[ii] ;
1056     
1057   }
1058   FINISH :
1059     ;
1060 }
1061
1062 //=======================================================================
1063 //function : FunctionMultiply
1064 //purpose  : 
1065 //=======================================================================
1066
1067 void BSplCLib::FunctionMultiply
1068 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1069  const Standard_Integer              BSplineDegree,
1070  const TColStd_Array1OfReal &        BSplineFlatKnots,
1071  const TColStd_Array1OfReal &        Poles,
1072  const TColStd_Array1OfReal &        FlatKnots,
1073  const Standard_Integer              NewDegree,
1074  TColStd_Array1OfReal &              NewPoles,
1075  Standard_Integer &                  Status) 
1076
1077   Standard_Integer num_bspline_poles =
1078     BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1079   Standard_Integer num_new_poles =
1080     FlatKnots.Length() - NewDegree - 1 ;
1081   
1082   if (Poles.Length() != num_bspline_poles ||
1083       NewPoles.Length() != num_new_poles) {
1084     Standard_ConstructionError::Raise();  
1085   }
1086   Standard_Real  * array_of_poles =
1087     (Standard_Real *) &Poles(Poles.Lower()) ;
1088   Standard_Real  * array_of_new_poles =
1089     (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1090   BSplCLib::FunctionMultiply(FunctionPtr,
1091                              BSplineDegree,
1092                              BSplineFlatKnots,
1093                              1,
1094                              array_of_poles[0],
1095                              FlatKnots,
1096                              NewDegree,
1097                              array_of_new_poles[0],
1098                              Status) ;
1099 }
1100
1101 //=======================================================================
1102 //function : FunctionReparameterise
1103 //purpose  : 
1104 //=======================================================================
1105
1106 void BSplCLib::FunctionReparameterise
1107 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1108  const Standard_Integer              BSplineDegree,
1109  const TColStd_Array1OfReal &        BSplineFlatKnots,
1110  const TColStd_Array1OfReal &        Poles,
1111  const TColStd_Array1OfReal &        FlatKnots,
1112  const Standard_Integer              NewDegree,
1113  TColStd_Array1OfReal &              NewPoles,
1114  Standard_Integer &                  Status) 
1115
1116   Standard_Integer num_bspline_poles =
1117     BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1118   Standard_Integer num_new_poles =
1119     FlatKnots.Length() - NewDegree - 1 ;
1120   
1121   if (Poles.Length() != num_bspline_poles ||
1122       NewPoles.Length() != num_new_poles) {
1123     Standard_ConstructionError::Raise();  
1124   }
1125   Standard_Real  * array_of_poles =
1126     (Standard_Real *) &Poles(Poles.Lower()) ;
1127   Standard_Real  * array_of_new_poles =
1128     (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1129   BSplCLib::FunctionReparameterise(
1130                                    FunctionPtr,
1131                                    BSplineDegree,
1132                                    BSplineFlatKnots,
1133                                    1,
1134                                    array_of_poles[0],
1135                                    FlatKnots,
1136                                    NewDegree,
1137                                    array_of_new_poles[0],
1138                                    Status) ;
1139 }
1140
1141 //=======================================================================
1142 //function : MergeBSplineKnots
1143 //purpose  : 
1144 //=======================================================================
1145 void BSplCLib::MergeBSplineKnots
1146 (const Standard_Real                 Tolerance,
1147  const Standard_Real                 StartValue,
1148  const Standard_Real                 EndValue,
1149  const Standard_Integer              Degree1,
1150  const TColStd_Array1OfReal          & Knots1,
1151  const TColStd_Array1OfInteger       & Mults1,
1152  const Standard_Integer              Degree2,
1153  const TColStd_Array1OfReal          & Knots2,
1154  const TColStd_Array1OfInteger       & Mults2,
1155  Standard_Integer &                  NumPoles,
1156  Handle(TColStd_HArray1OfReal) &     NewKnots,
1157  Handle(TColStd_HArray1OfInteger) &  NewMults) 
1158 {
1159   Standard_Integer ii,
1160   jj,
1161   continuity,
1162   set_mults_flag,
1163   degree,
1164   index,
1165   num_knots ;
1166   if (StartValue < EndValue - Tolerance) {
1167     TColStd_Array1OfReal knots1(1,Knots1.Length()) ;
1168     TColStd_Array1OfReal knots2(1,Knots2.Length()) ;
1169     degree = Degree1 + Degree2 ;
1170     index = 1 ;
1171
1172     for (ii = Knots1.Lower() ; ii <= Knots1.Upper() ; ii++) {
1173       knots1(index) = Knots1(ii) ;
1174       index += 1 ;
1175     }
1176     index = 1 ;
1177
1178     for (ii = Knots2.Lower() ; ii <= Knots2.Upper() ; ii++) {
1179       knots2(index) = Knots2(ii) ;
1180       index += 1 ;
1181     }
1182     BSplCLib::Reparametrize(StartValue,
1183                             EndValue,
1184                             knots1) ;
1185     
1186     BSplCLib::Reparametrize(StartValue,
1187                             EndValue,
1188                             knots2) ;
1189     num_knots = 0 ;
1190     jj =  1 ;
1191
1192     for (ii = 1 ; ii <= knots1.Length() ; ii++) {
1193
1194       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) - Tolerance) {
1195         jj += 1 ;
1196         num_knots += 1 ;
1197       }
1198
1199       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) + Tolerance) {
1200         jj += 1 ;
1201       }
1202       num_knots += 1 ;
1203     }
1204     NewKnots = 
1205       new TColStd_HArray1OfReal(1,num_knots) ;
1206     NewMults =
1207       new TColStd_HArray1OfInteger(1,num_knots) ;
1208     num_knots = 1 ;
1209     jj = 1 ;
1210
1211     for (ii = 1 ; ii <= knots1.Length() ; ii++) {
1212
1213       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) - Tolerance) {
1214         NewKnots->ChangeArray1()(num_knots) = knots2(jj) ;
1215         NewMults->ChangeArray1()(num_knots) = Mults2(jj) + Degree1 ;
1216         jj += 1 ;
1217         num_knots += 1 ;
1218       }
1219       set_mults_flag = 0 ;
1220
1221       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) + Tolerance) {
1222         continuity = Min(Degree1 - Mults1(ii), Degree2 - Mults2(jj)) ;
1223         set_mults_flag = 1 ;
1224         NewMults->ChangeArray1()(num_knots) = degree - continuity ;
1225         jj += 1 ;
1226       }
1227
1228       NewKnots->ChangeArray1()(num_knots) = knots1(ii) ;
1229       if (! set_mults_flag) {
1230         NewMults->ChangeArray1()(num_knots) = Mults1(ii) + Degree2 ;
1231       }
1232       num_knots += 1 ;
1233     }
1234     num_knots -= 1 ;
1235     NewMults->ChangeArray1()(1) = degree + 1 ;
1236     NewMults->ChangeArray1()(num_knots) = degree + 1 ;
1237     index = 0 ;
1238
1239     for (ii = 1 ; ii <= num_knots ; ii++) {
1240       index += NewMults->Value(ii) ;
1241     }
1242     NumPoles = index  - degree - 1  ;
1243   }
1244 }
1245