0028838: Configuration - undefine macros coming from X11 headers in place of collision
[occt.git] / src / BSplCLib / BSplCLib_CurveComputation.gxx
1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15 // xab : modified 15-Mar-94 added EvalBSplineMatrix, FactorBandedMatrix, SolveBandedSystem
16 //                          Eval, BuildCache, cacheD0, cacheD1, cacheD2, cacheD3, LocalMatrix,
17 //                          EvalPolynomial, RationalDerivatives.
18 // xab : 22-Mar-94 : fixed rational problem in BuildCache
19 // xab : 12-Mar-96 : added MovePointAndTangent
20 #include <TColStd_Array1OfInteger.hxx>
21 #include <TColStd_Array1OfReal.hxx>
22 #include <TColStd_Array2OfReal.hxx>
23 #include <gp_Vec2d.hxx>
24 #include <Standard_ConstructionError.hxx>
25 #include <PLib.hxx>
26 #include <math_Matrix.hxx>
27
28 //=======================================================================
29 //struct : BSplCLib_DataContainer 
30 //purpose: Auxiliary structure providing buffers for poles and knots used in
31 //         evaluation of bspline (allocated in the stack)
32 //=======================================================================
33
34 struct BSplCLib_DataContainer 
35 {
36   BSplCLib_DataContainer(Standard_Integer Degree) 
37   {
38     (void)Degree; // avoid compiler warning
39     Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree() ||
40         BSplCLib::MaxDegree() > 25,
41         "BSplCLib: bspline degree is greater than maximum supported");
42   }
43
44   Standard_Real poles[(25+1)*(Dimension_gen+1)];
45   Standard_Real knots[2*25];
46   Standard_Real ders[Dimension_gen*4];
47 };
48
49 //=======================================================================
50 //function : Reverse
51 //purpose  : 
52 //=======================================================================
53
54 void  BSplCLib::Reverse(Array1OfPoints& Poles,
55                         const Standard_Integer L)
56 {
57   Standard_Integer i, l = L;
58   l = Poles.Lower() + (l-Poles.Lower()) % (Poles.Upper()-Poles.Lower()+1);
59   
60   Array1OfPoints temp(0,Poles.Length()-1);
61   
62   for (i = Poles.Lower(); i <= l; i++)
63     temp(l-i) = Poles(i);
64   
65   for (i = l+1; i <= Poles.Upper(); i++)
66     temp(l-Poles.Lower()+Poles.Upper()-i+1) = Poles(i);
67   
68   for (i = Poles.Lower(); i <= Poles.Upper(); i++)
69     Poles(i) = temp(i-Poles.Lower());
70 }
71
72 //
73 // CURVES MODIFICATIONS
74 //
75
76 //=======================================================================
77 //function : RemoveKnot
78 //purpose  : 
79 //=======================================================================
80
81 Standard_Boolean BSplCLib::RemoveKnot 
82 (const Standard_Integer         Index,       
83  const Standard_Integer         Mult,        
84  const Standard_Integer         Degree,  
85  const Standard_Boolean         Periodic,
86  const Array1OfPoints&          Poles,
87  const TColStd_Array1OfReal*    Weights,
88  const TColStd_Array1OfReal&    Knots,  
89  const TColStd_Array1OfInteger& Mults,
90  Array1OfPoints&                NewPoles,   
91  TColStd_Array1OfReal*          NewWeights,
92  TColStd_Array1OfReal&          NewKnots,  
93  TColStd_Array1OfInteger&       NewMults,
94  const Standard_Real            Tolerance) 
95 {
96   Standard_Boolean rational = Weights  != NULL;
97   Standard_Integer dim;
98   dim = Dimension_gen;
99   if (rational) dim++;
100   
101   TColStd_Array1OfReal poles(1,dim*Poles.Length());
102   TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
103   
104   if (rational) PLib::SetPoles(Poles,*Weights,poles);
105   else          PLib::SetPoles(Poles,poles);
106   
107   if (!RemoveKnot(Index,Mult,Degree,Periodic,dim,
108                   poles,Knots,Mults,newpoles,NewKnots,NewMults,Tolerance))
109     return Standard_False;
110
111   if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
112   else          PLib::GetPoles(newpoles,NewPoles);
113   return Standard_True;
114 }
115
116 //=======================================================================
117 //function : InsertKnots
118 //purpose  : insert an array of knots and multiplicities
119 //=======================================================================
120
121 void BSplCLib::InsertKnots
122 (const Standard_Integer         Degree, 
123  const Standard_Boolean         Periodic,
124  const Array1OfPoints&          Poles,  
125  const TColStd_Array1OfReal*    Weights,  
126  const TColStd_Array1OfReal&    Knots,    
127  const TColStd_Array1OfInteger& Mults, 
128  const TColStd_Array1OfReal&    AddKnots,    
129  const TColStd_Array1OfInteger* AddMults, 
130  Array1OfPoints&                NewPoles,     
131  TColStd_Array1OfReal*          NewWeights,
132  TColStd_Array1OfReal&          NewKnots,    
133  TColStd_Array1OfInteger&       NewMults, 
134  const Standard_Real            Epsilon,
135  const Standard_Boolean         Add) 
136 {
137   Standard_Boolean rational = Weights  != NULL;
138   Standard_Integer dim;
139   dim = Dimension_gen;
140   if (rational) dim++;
141   
142   TColStd_Array1OfReal poles(1,dim*Poles.Length());
143   TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
144   
145   if (rational) PLib::SetPoles(Poles,*Weights,poles);
146   else          PLib::SetPoles(Poles,poles);
147   
148   InsertKnots(Degree,Periodic,dim,poles,Knots,Mults,
149               AddKnots,AddMults,newpoles,NewKnots,NewMults,Epsilon,Add);
150   
151   if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
152   else          PLib::GetPoles(newpoles,NewPoles);
153 }
154
155 //=======================================================================
156 //function : InsertKnot
157 //purpose  : 
158 //=======================================================================
159
160 void  BSplCLib::InsertKnot(const Standard_Integer , 
161                            const Standard_Real U,
162                            const Standard_Integer UMult, 
163                            const Standard_Integer Degree, 
164                            const Standard_Boolean Periodic, 
165                            const Array1OfPoints& Poles, 
166                            const TColStd_Array1OfReal* Weights, 
167                            const TColStd_Array1OfReal& Knots, 
168                            const TColStd_Array1OfInteger& Mults, 
169                            Array1OfPoints& NewPoles, 
170                            TColStd_Array1OfReal* NewWeights)
171 {
172   TColStd_Array1OfReal k(1,1);
173   k(1) = U;
174   TColStd_Array1OfInteger m(1,1);
175   m(1) = UMult;
176   TColStd_Array1OfReal    nk(1,Knots.Length()+1);
177   TColStd_Array1OfInteger nm(1,Knots.Length()+1);
178   InsertKnots(Degree,Periodic,Poles,Weights,Knots,Mults,
179               k,&m,NewPoles,NewWeights,nk,nm,Epsilon(U));
180 }
181
182 //=======================================================================
183 //function : RaiseMultiplicity
184 //purpose  : 
185 //=======================================================================
186
187 void  BSplCLib::RaiseMultiplicity(const Standard_Integer KnotIndex,
188                                   const Standard_Integer Mult,
189                                   const Standard_Integer Degree, 
190                                   const Standard_Boolean Periodic,
191                                   const Array1OfPoints& Poles,
192                                   const TColStd_Array1OfReal* Weights, 
193                                   const TColStd_Array1OfReal& Knots, 
194                                   const TColStd_Array1OfInteger& Mults, 
195                                   Array1OfPoints& NewPoles,
196                                   TColStd_Array1OfReal* NewWeights)
197 {
198   TColStd_Array1OfReal k(1,1);
199   k(1) = Knots(KnotIndex);
200   TColStd_Array1OfInteger m(1,1);
201   m(1) = Mult - Mults(KnotIndex);
202   TColStd_Array1OfReal    nk(1,Knots.Length());
203   TColStd_Array1OfInteger nm(1,Knots.Length());
204   InsertKnots(Degree,Periodic,Poles,Weights,Knots,Mults,
205               k,&m,NewPoles,NewWeights,nk,nm,Epsilon(k(1)));
206 }
207
208 //=======================================================================
209 //function : IncreaseDegree
210 //purpose  : 
211 //=======================================================================
212
213 void BSplCLib::IncreaseDegree 
214 (const Standard_Integer          Degree,
215  const Standard_Integer          NewDegree,
216  const Standard_Boolean          Periodic,
217  const Array1OfPoints&           Poles,
218  const TColStd_Array1OfReal*     Weights,
219  const TColStd_Array1OfReal&     Knots,
220  const TColStd_Array1OfInteger&  Mults,
221  Array1OfPoints&                 NewPoles,
222  TColStd_Array1OfReal*           NewWeights,
223  TColStd_Array1OfReal&           NewKnots,
224  TColStd_Array1OfInteger&        NewMults) 
225 {
226   Standard_Boolean rational = Weights  != NULL;
227   Standard_Integer dim;
228   dim = Dimension_gen;
229   if (rational) dim++;
230   
231   TColStd_Array1OfReal poles(1,dim*Poles.Length());
232   TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
233   
234   if (rational) PLib::SetPoles(Poles,*Weights,poles);
235   else          PLib::SetPoles(Poles,poles);
236   
237   IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults,
238                  newpoles,NewKnots,NewMults);
239   
240   if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
241   else          PLib::GetPoles(newpoles,NewPoles);
242 }
243
244 //=======================================================================
245 //function : Unperiodize
246 //purpose  : 
247 //=======================================================================
248
249 void  BSplCLib::Unperiodize
250 (const Standard_Integer         Degree, 
251  const TColStd_Array1OfInteger& Mults,
252  const TColStd_Array1OfReal&    Knots,
253  const Array1OfPoints&          Poles,
254  const TColStd_Array1OfReal*    Weights,
255  TColStd_Array1OfInteger& NewMults,
256  TColStd_Array1OfReal&    NewKnots,
257  Array1OfPoints&          NewPoles,
258  TColStd_Array1OfReal*    NewWeights)
259 {
260   Standard_Boolean rational = Weights  != NULL;
261   Standard_Integer dim;
262   dim = Dimension_gen;
263   if (rational) dim++;
264   
265   TColStd_Array1OfReal poles(1,dim*Poles.Length());
266   TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
267   
268   if (rational) PLib::SetPoles(Poles,*Weights,poles);
269   else          PLib::SetPoles(Poles,poles);
270   
271   Unperiodize(Degree, dim, Mults, Knots, poles,
272               NewMults,NewKnots, newpoles);
273   
274   if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
275   else          PLib::GetPoles(newpoles,NewPoles);
276 }
277
278 //=======================================================================
279 //function : Trimming
280 //purpose  : 
281 //=======================================================================
282
283 void BSplCLib::Trimming(const Standard_Integer         Degree, 
284                         const Standard_Boolean         Periodic,
285                         const TColStd_Array1OfReal&    Knots, 
286                         const TColStd_Array1OfInteger& Mults,
287                         const Array1OfPoints&          Poles,
288                         const TColStd_Array1OfReal*    Weights, 
289                         const Standard_Real            U1,
290                         const Standard_Real            U2, 
291                         TColStd_Array1OfReal&    NewKnots,
292                         TColStd_Array1OfInteger& NewMults,
293                         Array1OfPoints&          NewPoles,
294                         TColStd_Array1OfReal*    NewWeights)
295 {
296   Standard_Boolean rational = Weights  != NULL;
297   Standard_Integer dim;
298   dim = Dimension_gen;
299   if (rational) dim++;
300   
301   TColStd_Array1OfReal poles(1,dim*Poles.Length());
302   TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
303   
304   if (rational) PLib::SetPoles(Poles,*Weights,poles);
305   else          PLib::SetPoles(Poles,poles);
306   
307   Trimming(Degree, Periodic, dim, Knots, Mults, poles, U1, U2,
308            NewKnots, NewMults, newpoles);
309   
310   if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
311   else          PLib::GetPoles(newpoles,NewPoles);
312 }
313
314 //--------------------------------------------------------------------------
315 //ELEMENTARY COMPUTATIONS
316 //--------------------------------------------------------------------------
317 //
318 // All the following methods are methods of low level used in BSplCLib 
319 // but also in BSplSLib (but not in Geom)
320 // At the creation time of this package there is no possibility 
321 // to declare private methods of package and to declare friend
322 // methods of package.  It could interesting to declare that BSplSLib
323 // is a package friend to BSplCLib because it uses all the basis methods
324 // of BSplCLib.
325 //--------------------------------------------------------------------------
326
327 //=======================================================================
328 //function : BuildEval
329 //purpose  : builds the local array for evaluation
330 //=======================================================================
331
332 void  BSplCLib::BuildEval(const Standard_Integer         Degree,
333                           const Standard_Integer         Index,
334                           const Array1OfPoints&          Poles, 
335                           const TColStd_Array1OfReal*    Weights,
336                           Standard_Real&                 LP)
337 {
338   Standard_Real w, *pole = &LP;
339   Standard_Integer PLower = Poles.Lower();
340   Standard_Integer PUpper = Poles.Upper();
341   Standard_Integer i;
342   Standard_Integer ip = PLower + Index - 1;
343   if (Weights == NULL) {
344     for (i = 0; i <= Degree; i++) {
345       ip++;
346       if (ip > PUpper) ip = PLower;
347       const Point& P = Poles(ip);
348       PointToCoords (pole,P,+0);
349       pole += Dimension_gen;
350     }
351   }
352   else {
353     for (i = 0; i <= Degree; i++) {
354       ip++;
355       if (ip > PUpper) ip = PLower;
356       const Point& P = Poles(ip);
357       pole[Dimension_gen] = w = (*Weights)(ip);
358       PointToCoords (pole, P, * w);
359       pole += Dimension_gen + 1;
360     }
361   }
362 }
363
364 //=======================================================================
365 //function : PrepareEval
366 //purpose  : stores data for Eval in the local arrays
367 //           dc.poles and dc.knots
368 //=======================================================================
369
370 static void PrepareEval
371 (Standard_Real&                 u,                  
372  Standard_Integer&              index, 
373  Standard_Integer&              dim,
374  Standard_Boolean&              rational,
375  const Standard_Integer         Degree,     
376  const Standard_Boolean         Periodic,
377  const Array1OfPoints&          Poles,  
378  const TColStd_Array1OfReal*    Weights,
379  const TColStd_Array1OfReal&    Knots,  
380  const TColStd_Array1OfInteger* Mults,
381  BSplCLib_DataContainer& dc) 
382 {
383   // Set the Index
384   BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
385
386   // make the knots
387   BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*dc.knots);
388   if (Mults == NULL)
389     index -= Knots.Lower() + Degree;
390   else
391     index = BSplCLib::PoleIndex(Degree,index,Periodic,*Mults);
392   
393   // check truly rational
394   rational = (Weights != NULL);
395   if (rational) {
396     Standard_Integer WLower = Weights->Lower() + index;
397     rational = BSplCLib::IsRational(*Weights, WLower, WLower + Degree);
398   }
399   
400   // make the poles
401   if (rational) {
402     dim = Dimension_gen + 1;
403     BSplCLib::BuildEval(Degree, index, Poles, Weights              , *dc.poles);
404   }
405   else {
406     dim = Dimension_gen; 
407     BSplCLib::BuildEval(Degree, index, Poles, BSplCLib::NoWeights(), *dc.poles);
408   }
409 }
410
411 //=======================================================================
412 //function : D0
413 //purpose  : 
414 //=======================================================================
415
416 void BSplCLib::D0 
417 (const Standard_Real            U,                  
418  const Standard_Integer         Index,          
419  const Standard_Integer         Degree,     
420  const Standard_Boolean         Periodic,
421  const Array1OfPoints&          Poles,  
422  const TColStd_Array1OfReal*    Weights,
423  const TColStd_Array1OfReal&    Knots,  
424  const TColStd_Array1OfInteger* Mults,  
425  Point&                         P) 
426 {                    
427 //  Standard_Integer k,dim,index = Index;
428   Standard_Integer dim,index = Index;
429   Standard_Real    u = U;
430   Standard_Boolean rational;
431   BSplCLib_DataContainer dc(Degree);
432   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
433   BSplCLib::Eval(u,Degree,*dc.knots,dim,*dc.poles);
434
435   if (rational) {
436     Standard_Real w = dc.poles[Dimension_gen];
437     CoordsToPoint (P, dc.poles, / w);
438   }
439   else
440     CoordsToPoint (P, dc.poles, );
441 }
442
443 //=======================================================================
444 //function : D1
445 //purpose  : 
446 //=======================================================================
447
448 void BSplCLib::D1
449 (const Standard_Real            U,                  
450  const Standard_Integer         Index,          
451  const Standard_Integer         Degree,     
452  const Standard_Boolean         Periodic,
453  const Array1OfPoints&          Poles,  
454  const TColStd_Array1OfReal*    Weights,
455  const TColStd_Array1OfReal&    Knots,  
456  const TColStd_Array1OfInteger* Mults,  
457  Point&                         P,
458  Vector&                        V) 
459 {                    
460   Standard_Integer dim,index = Index;
461   Standard_Real    u = U;
462   Standard_Boolean rational;
463   BSplCLib_DataContainer dc(Degree);
464   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
465   BSplCLib::Bohm(u,Degree,1,*dc.knots,dim,*dc.poles);
466   Standard_Real *result = dc.poles;
467   if (rational) {
468     PLib::RationalDerivative(Degree,1,Dimension_gen,*dc.poles,*dc.ders);
469     result = dc.ders;
470   }
471   
472   CoordsToPoint (P, result, );
473   CoordsToPoint (V, result + Dimension_gen, );
474 }
475
476 //=======================================================================
477 //function : D2
478 //purpose  : 
479 //=======================================================================
480
481 void BSplCLib::D2
482 (const Standard_Real            U,                  
483  const Standard_Integer         Index,          
484  const Standard_Integer         Degree,     
485  const Standard_Boolean         Periodic,
486  const Array1OfPoints&          Poles,  
487  const TColStd_Array1OfReal*    Weights,
488  const TColStd_Array1OfReal&    Knots,  
489  const TColStd_Array1OfInteger* Mults,  
490  Point&                         P,
491  Vector&                        V1,
492  Vector&                        V2) 
493 {                    
494   Standard_Integer dim,index = Index;
495   Standard_Real    u = U;
496   Standard_Boolean rational;
497   BSplCLib_DataContainer dc(Degree);
498   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
499   BSplCLib::Bohm(u,Degree,2,*dc.knots,dim,*dc.poles);
500   Standard_Real *result = dc.poles;
501   if (rational) {
502     PLib::RationalDerivative(Degree,2,Dimension_gen,*dc.poles,*dc.ders);
503     result = dc.ders;
504   }
505   
506   CoordsToPoint (P,  result, );
507   CoordsToPoint (V1, result + Dimension_gen, );
508   if (!rational && (Degree < 2)) 
509     NullifyPoint (V2);
510   else
511     CoordsToPoint (V2, result + 2 * Dimension_gen, );  
512 }
513
514 //=======================================================================
515 //function : D3
516 //purpose  : 
517 //=======================================================================
518
519 void BSplCLib::D3
520 (const Standard_Real            U,                  
521  const Standard_Integer         Index,          
522  const Standard_Integer         Degree,     
523  const Standard_Boolean         Periodic,
524  const Array1OfPoints&          Poles,  
525  const TColStd_Array1OfReal*    Weights,
526  const TColStd_Array1OfReal&    Knots,  
527  const TColStd_Array1OfInteger* Mults,  
528  Point&                         P,
529  Vector&                        V1,
530  Vector&                        V2,
531  Vector&                        V3) 
532 {                    
533   Standard_Integer dim,index = Index;
534   Standard_Real    u = U;
535   Standard_Boolean rational;
536   BSplCLib_DataContainer dc(Degree);
537   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
538   BSplCLib::Bohm(u,Degree,3,*dc.knots,dim,*dc.poles);
539   Standard_Real *result = dc.poles;
540   if (rational) {
541     PLib::RationalDerivative(Degree,3,Dimension_gen,*dc.poles,*dc.ders);
542     result = dc.ders;
543   }
544
545   CoordsToPoint (P,  result, );
546   CoordsToPoint (V1, result + Dimension_gen, );
547   if (!rational && (Degree < 2)) 
548     NullifyPoint (V2);
549   else
550     CoordsToPoint (V2, result + 2 * Dimension_gen, );  
551   if (!rational && (Degree < 3)) 
552     NullifyPoint (V3);
553   else
554     CoordsToPoint (V3, result + 3 * Dimension_gen, );  
555 }
556
557 //=======================================================================
558 //function : DN
559 //purpose  : 
560 //=======================================================================
561
562 void BSplCLib::DN
563 (const Standard_Real            U,                  
564  const Standard_Integer         N,          
565  const Standard_Integer         Index,          
566  const Standard_Integer         Degree,     
567  const Standard_Boolean         Periodic,
568  const Array1OfPoints&          Poles,  
569  const TColStd_Array1OfReal*    Weights,
570  const TColStd_Array1OfReal&    Knots,  
571  const TColStd_Array1OfInteger* Mults,  
572  Vector&                        VN) 
573 {                    
574   Standard_Integer dim,index = Index;
575   Standard_Real    u = U;
576   Standard_Boolean rational;
577   BSplCLib_DataContainer dc(Degree);
578   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
579   BSplCLib::Bohm(u,Degree,N,*dc.knots,dim,*dc.poles);
580
581   if (rational) {
582     Standard_Real v[Dimension_gen];
583     PLib::RationalDerivative(Degree,N,Dimension_gen,*dc.poles,v[0],Standard_False);
584     CoordsToPoint (VN, v, );
585   }
586   else {
587     if (N > Degree) 
588       NullifyPoint (VN);
589     else {
590       Standard_Real *DN = dc.poles + N * Dimension_gen;
591       CoordsToPoint (VN, DN, );
592     }
593   }
594 }
595
596 //=======================================================================
597 //function : Solves a LU factored Matrix 
598 //purpose  : 
599 //=======================================================================
600
601 Standard_Integer 
602 BSplCLib::SolveBandedSystem(const math_Matrix&  Matrix,
603                             const Standard_Integer UpperBandWidth,
604                             const Standard_Integer LowerBandWidth,
605                             Array1OfPoints&   PolesArray) 
606 {
607   Standard_Real *PArray ;
608   PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
609   
610   return BSplCLib::SolveBandedSystem(Matrix,
611                                      UpperBandWidth,
612                                      LowerBandWidth,
613                                      Dimension_gen,
614                                      PArray[0]) ;
615 }
616
617 //=======================================================================
618 //function : Solves a LU factored Matrix 
619 //purpose  : if HomogeneousFlag is 1 then the input and the output
620 //           will be homogeneous that is no division or multiplication
621 //           by weigths will happen. On the contrary if HomogeneousFlag
622 //           is 0 then the poles will be multiplied first by the weights
623 //           and after interpolation they will be devided by the weights
624 //=======================================================================
625
626 Standard_Integer 
627 BSplCLib::SolveBandedSystem(const math_Matrix&  Matrix,
628                             const Standard_Integer UpperBandWidth,
629                             const Standard_Integer LowerBandWidth,
630                             const Standard_Boolean HomogeneousFlag,
631                             Array1OfPoints&        PolesArray,
632                             TColStd_Array1OfReal&  WeightsArray) 
633 {
634   Standard_Real *PArray,
635   * WArray ;
636   PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
637   WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
638   return BSplCLib::SolveBandedSystem(Matrix,
639                                      UpperBandWidth,
640                                      LowerBandWidth,
641                                      HomogeneousFlag,
642                                      Dimension_gen,
643                                      PArray[0],
644                                      WArray[0]) ;
645 }
646
647 //=======================================================================
648 //function : Evaluates a Bspline function : uses the ExtrapMode 
649 //purpose  : the function is extrapolated using the Taylor expansion
650 //           of degree ExtrapMode[0] to the left and the Taylor
651 //           expansion of degree ExtrapMode[1] to the right 
652 //   if the HomogeneousFlag == True than the Poles are supposed
653 //   to be stored homogeneously and the result will also be homogeneous
654 //   Valid only if Weights 
655 //=======================================================================
656 void  BSplCLib::Eval(const Standard_Real                   Parameter,
657                      const Standard_Boolean                PeriodicFlag,
658                      const Standard_Boolean                HomogeneousFlag,
659                      Standard_Integer&                     ExtrapMode,
660                      const Standard_Integer                Degree,
661                      const  TColStd_Array1OfReal&          FlatKnots, 
662                      const  Array1OfPoints&                PolesArray,
663                      const  TColStd_Array1OfReal&          WeightsArray,
664                      Point&                                aPoint,
665                      Standard_Real&                        aWeight)
666 {
667   Standard_Real  Inverse,
668   P[Dimension_gen],
669   *PArray,
670   *WArray ;
671   Standard_Integer kk ;
672   PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
673   WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
674   if (HomogeneousFlag) {
675     BSplCLib::Eval(Parameter,
676                    PeriodicFlag,
677                    0,
678                    ExtrapMode,
679                    Degree,
680                    FlatKnots,
681                    Dimension_gen,
682                    PArray[0],
683                    P[0]) ;
684     BSplCLib::Eval(Parameter,
685                    PeriodicFlag,
686                    0,
687                    ExtrapMode,
688                    Degree,
689                    FlatKnots,
690                    1,
691                    WArray[0],
692                    aWeight) ;
693   }
694   else {
695     BSplCLib::Eval(Parameter,
696                    PeriodicFlag,
697                    0,
698                    ExtrapMode,
699                    Degree,
700                    FlatKnots,
701                    Dimension_gen,
702                    PArray[0],
703                    WArray[0],
704                    P[0],
705                    aWeight) ;
706     Inverse = 1.0e0 / aWeight ;
707
708     for (kk = 0 ; kk < Dimension_gen ; kk++) {
709       P[kk] *= Inverse ;
710     } 
711   }
712   
713   for (kk = 0 ; kk < Dimension_gen ; kk++) 
714     aPoint.SetCoord(kk + 1, P[kk]) ;
715 }
716
717 //=======================================================================
718 //function : CacheD0
719 //purpose  : Evaluates the polynomial cache of the Bspline Curve
720 //           
721 //=======================================================================
722 void  BSplCLib::CacheD0(const Standard_Real                  Parameter,
723                         const  Standard_Integer               Degree,
724                         const  Standard_Real                  CacheParameter,
725                         const  Standard_Real                  SpanLenght,
726                         const  Array1OfPoints&                PolesArray,
727                         const  TColStd_Array1OfReal*          WeightsArray,
728                         Point&                                aPoint)
729 {
730   //
731   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
732   // form
733   // the SpanLenght     is the normalizing factor so that the polynomial is between
734   // 0 and 1 
735   Standard_Real NewParameter, Inverse;
736   Standard_Real * PArray  = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
737   Standard_Real * myPoint = (Standard_Real *) &aPoint  ;
738   NewParameter = (Parameter - CacheParameter) / SpanLenght ; 
739   PLib::NoDerivativeEvalPolynomial(NewParameter,
740                        Degree,
741                        Dimension_gen,
742                        Degree * Dimension_gen,
743                        PArray[0],
744                        myPoint[0]) ;
745   if (WeightsArray != NULL) {
746     const TColStd_Array1OfReal& refWeights = *WeightsArray;
747     Standard_Real *
748       WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
749     PLib::NoDerivativeEvalPolynomial(NewParameter,
750                          Degree,
751                          1,
752                          Degree,
753                          WArray[0],
754                          Inverse) ;
755     
756     Inverse = 1.0e0 / Inverse;
757     ModifyCoords (myPoint, *= Inverse);
758   }
759 }
760
761 //=======================================================================
762 //function : CacheD1
763 //purpose  : Evaluates the polynomial cache of the Bspline Curve
764 //           
765 //=======================================================================
766 void  BSplCLib::CacheD1(const Standard_Real                  Parameter,
767                         const Standard_Integer                Degree,
768                         const  Standard_Real                  CacheParameter,
769                         const  Standard_Real                  SpanLenght,
770                         const  Array1OfPoints&                PolesArray,
771                         const  TColStd_Array1OfReal*          WeightsArray,
772                         Point&                                aPoint,
773                         Vector&                               aVector) 
774 {
775   //
776   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
777   // form
778   // the SpanLenght     is the normalizing factor so that the polynomial is between
779   // 0 and 1 
780   Standard_Real LocalPDerivatives[Dimension_gen << 1] ;
781 //  Standard_Real LocalWDerivatives[2], NewParameter, Inverse ;
782   Standard_Real LocalWDerivatives[2], NewParameter ;
783   
784   Standard_Real * 
785     PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
786   Standard_Real *
787     myPoint = (Standard_Real *) &aPoint  ;
788   Standard_Real *
789     myVector = (Standard_Real *) &aVector ;
790   NewParameter = (Parameter - CacheParameter) / SpanLenght ; 
791   PLib::EvalPolynomial(NewParameter,
792                        1,
793                        Degree,
794                        Dimension_gen,
795                        PArray[0],
796                        LocalPDerivatives[0]) ;
797   // 
798   // unormalize derivatives since those are computed normalized
799   //
800
801   ModifyCoords (LocalPDerivatives + Dimension_gen, /= SpanLenght);
802   
803   if (WeightsArray != NULL) {
804     const TColStd_Array1OfReal& refWeights = *WeightsArray;
805     Standard_Real *
806       WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
807     PLib::EvalPolynomial(NewParameter,
808                          1,
809                          Degree,
810                          1,
811                          WArray[0],
812                          LocalWDerivatives[0]) ;
813     // 
814     // unormalize the result since the polynomial stored in the cache
815     // is normalized between 0 and 1
816     //
817     LocalWDerivatives[1] /= SpanLenght ;
818     
819     PLib::RationalDerivatives(1,
820                               Dimension_gen,
821                               LocalPDerivatives[0],
822                               LocalWDerivatives[0],
823                               LocalPDerivatives[0]) ;
824   }
825   
826   CopyCoords (myPoint,  LocalPDerivatives);
827   CopyCoords (myVector, LocalPDerivatives + Dimension_gen);
828 }
829
830 //=======================================================================
831 //function : CacheD2
832 //purpose  : Evaluates the polynomial cache of the Bspline Curve
833 //           
834 //=======================================================================
835 void  BSplCLib::CacheD2(const Standard_Real                  Parameter,
836                         const Standard_Integer                Degree,
837                         const  Standard_Real                  CacheParameter,
838                         const  Standard_Real                  SpanLenght,
839                         const  Array1OfPoints&                PolesArray,
840                         const  TColStd_Array1OfReal*          WeightsArray,
841                         Point&                                aPoint,
842                         Vector&                               aVector1, 
843                         Vector&                               aVector2) 
844 {
845   //
846   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
847   // form
848   // the SpanLenght     is the normalizing factor so that the polynomial is between
849   // 0 and 1 
850   Standard_Integer ii,Index,EndIndex;
851   Standard_Real LocalPDerivatives[(Dimension_gen << 1) + Dimension_gen] ;
852 //  Standard_Real LocalWDerivatives[3], NewParameter, Factor, Inverse ;
853   Standard_Real LocalWDerivatives[3], NewParameter, Factor ;
854   Standard_Real * PArray    = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
855   Standard_Real * myPoint   = (Standard_Real *) &aPoint  ;
856   Standard_Real * myVector1 = (Standard_Real *) &aVector1 ;
857   Standard_Real * myVector2 = (Standard_Real *) &aVector2 ;
858   NewParameter = (Parameter - CacheParameter) / SpanLenght ; 
859   PLib::EvalPolynomial(NewParameter,
860                        2,
861                        Degree,
862                        Dimension_gen,
863                        PArray[0],
864                        LocalPDerivatives[0]) ;
865   // 
866   // unormalize derivatives since those are computed normalized
867   //
868   Factor = 1.0e0/ SpanLenght ;
869   Index = Dimension_gen ;
870   EndIndex = Min (2, Degree) ;
871
872   for (ii = 1 ; ii <= EndIndex ; ii++) {
873     ModifyCoords (LocalPDerivatives + Index, *= Factor);
874     Factor /= SpanLenght ;
875     Index  += Dimension_gen;
876   }
877
878   Index = (Degree + 1) * Dimension_gen;
879   for (ii = Degree ; ii < 2 ; ii++) {
880     NullifyCoords (LocalPDerivatives + Index);
881     Index += Dimension_gen;
882   }
883
884   if (WeightsArray != NULL) {
885     const TColStd_Array1OfReal& refWeights = *WeightsArray;
886     Standard_Real *
887       WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
888     
889     PLib::EvalPolynomial(NewParameter,
890                          2,
891                          Degree,
892                          1,
893                          WArray[0],
894                          LocalWDerivatives[0]) ;
895
896     for (ii = Degree + 1  ; ii <= 2 ; ii++) {
897       LocalWDerivatives[ii] = 0.0e0 ;
898     }
899     // 
900     // unormalize the result since the polynomial stored in the cache
901     // is normalized between 0 and 1
902     //
903     Factor = 1.0e0 / SpanLenght ;
904
905     for (ii = 1 ; ii <= EndIndex ; ii++) { 
906       LocalWDerivatives[ii] *= Factor ;
907       Factor /= SpanLenght ;
908     }
909     PLib::RationalDerivatives(2,
910                               Dimension_gen,
911                               LocalPDerivatives[0],
912                               LocalWDerivatives[0],
913                               LocalPDerivatives[0]) ;
914   }
915
916   CopyCoords (myPoint,   LocalPDerivatives);
917   CopyCoords (myVector1, LocalPDerivatives + Dimension_gen);
918   CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2);
919 }
920
921 //=======================================================================
922 //function : CacheD3
923 //purpose  : Evaluates the polynomial cache of the Bspline Curve
924 //           
925 //=======================================================================
926 void  BSplCLib::CacheD3(const Standard_Real                  Parameter,
927                         const Standard_Integer                Degree,
928                         const  Standard_Real                  CacheParameter,
929                         const  Standard_Real                  SpanLenght,
930                         const  Array1OfPoints&                PolesArray,
931                         const  TColStd_Array1OfReal*          WeightsArray,
932                         Point&                                aPoint,
933                         Vector&                               aVector1, 
934                         Vector&                               aVector2,
935                         Vector&                               aVector3) 
936 {
937   //
938   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
939   // form
940   // the SpanLenght     is the normalizing factor so that the polynomial is between
941   // 0 and 1 
942   Standard_Integer ii, Index, EndIndex;
943   Standard_Real LocalPDerivatives[Dimension_gen << 2] ;
944 //  Standard_Real LocalWDerivatives[4], Factor, NewParameter, Inverse ;
945   Standard_Real LocalWDerivatives[4], Factor, NewParameter ;
946   Standard_Real * PArray    = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
947   Standard_Real * myPoint   = (Standard_Real *) &aPoint  ;
948   Standard_Real * myVector1 = (Standard_Real *) &aVector1 ;
949   Standard_Real * myVector2 = (Standard_Real *) &aVector2 ;
950   Standard_Real * myVector3 = (Standard_Real *) &aVector3 ;
951   NewParameter = (Parameter - CacheParameter) / SpanLenght ; 
952   PLib::EvalPolynomial(NewParameter,
953                        3,
954                        Degree,
955                        Dimension_gen,
956                        PArray[0],
957                        LocalPDerivatives[0]) ;
958
959   Index = (Degree + 1) * Dimension_gen;
960   for (ii = Degree ; ii < 3 ; ii++) {
961     NullifyCoords (LocalPDerivatives + Index);
962     Index += Dimension_gen;
963   }
964   
965   // 
966   // unormalize derivatives since those are computed normalized
967   //
968   Factor = 1.0e0 / SpanLenght ;
969   Index = Dimension_gen ;
970   EndIndex = Min(3,Degree) ;
971
972   for (ii = 1 ; ii <= EndIndex ; ii++) { 
973     ModifyCoords (LocalPDerivatives + Index, *= Factor);
974     Factor /= SpanLenght;
975     Index  += Dimension_gen;
976   }
977   
978   if (WeightsArray != NULL) {
979     const TColStd_Array1OfReal& refWeights = *WeightsArray;
980     Standard_Real *
981       WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
982     
983     PLib::EvalPolynomial(NewParameter,
984                          3,
985                          Degree,
986                          1,
987                          WArray[0],
988                          LocalWDerivatives[0]) ;
989     // 
990     // unormalize the result since the polynomial stored in the cache
991     // is normalized between 0 and 1
992     //
993     Factor = 1.0e0 / SpanLenght ;
994
995     for (ii = 1 ; ii <= EndIndex ; ii++) { 
996       LocalWDerivatives[ii] *= Factor ;
997       Factor /= SpanLenght ;
998     }
999
1000     for (ii = (Degree + 1)  ; ii <= 3 ; ii++) {
1001       LocalWDerivatives[ii] = 0.0e0 ;
1002     }
1003     PLib::RationalDerivatives(3,
1004                               Dimension_gen,
1005                               LocalPDerivatives[0],
1006                               LocalWDerivatives[0],
1007                               LocalPDerivatives[0]) ;
1008   }
1009   
1010   CopyCoords (myPoint,   LocalPDerivatives);
1011   CopyCoords (myVector1, LocalPDerivatives + Dimension_gen);
1012   CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2);
1013   CopyCoords (myVector3, LocalPDerivatives + Dimension_gen * 3);
1014 }
1015
1016 //=======================================================================
1017 //function : BuildCache
1018 //purpose  : Stores theTaylor expansion normalized between 0,1 in the
1019 //           Cache : in case of  a rational function the Poles are
1020 //           stored in homogeneous form 
1021 //=======================================================================
1022
1023 void BSplCLib::BuildCache
1024 (const Standard_Real            U,   
1025  const Standard_Real            SpanDomain,    
1026  const Standard_Boolean         Periodic,
1027  const Standard_Integer         Degree,     
1028  const TColStd_Array1OfReal&    FlatKnots,   
1029  const Array1OfPoints&          Poles,  
1030  const TColStd_Array1OfReal*    Weights,
1031  Array1OfPoints&                CachePoles,
1032  TColStd_Array1OfReal*          CacheWeights)
1033 {                    
1034   Standard_Integer ii,
1035   Dimension,
1036   LocalIndex,
1037   index = 0 ;
1038   Standard_Real    u = U,
1039   LocalValue ;
1040   Standard_Boolean rational;
1041   
1042   BSplCLib_DataContainer dc(Degree);
1043   PrepareEval(u,
1044               index,
1045               Dimension,
1046               rational,
1047               Degree,
1048               Periodic,
1049               Poles,
1050               Weights,
1051               FlatKnots,
1052               (BSplCLib::NoMults()),
1053               dc);
1054   // 
1055   // Watch out : PrepareEval checks if locally the Bspline is polynomial
1056   // therefore rational flag could be set to False even if there are 
1057   // Weights. The Dimension is set accordingly.
1058   //
1059   
1060   BSplCLib::Bohm(u,Degree,Degree,*dc.knots,Dimension,*dc.poles);
1061   
1062   LocalValue = 1.0e0 ;
1063   LocalIndex = 0 ;
1064
1065   if (rational) {
1066  
1067     for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1068       CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue);
1069       LocalIndex += Dimension_gen + 1;
1070       LocalValue *= SpanDomain / (Standard_Real) ii ;
1071     }
1072
1073     LocalIndex = Dimension_gen;
1074     LocalValue = 1.0e0 ;      
1075     for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1076       (*CacheWeights)(ii) = dc.poles[LocalIndex] * LocalValue ;
1077       LocalIndex += Dimension_gen + 1;
1078       LocalValue *= SpanDomain / (Standard_Real) ii ;
1079     }
1080   }
1081   else {
1082     
1083     for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1084       CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue);
1085       LocalIndex += Dimension_gen;
1086       LocalValue *= SpanDomain / (Standard_Real) ii ;
1087     }
1088
1089     if (Weights != NULL) { 
1090       for (ii = 1 ; ii <= Degree + 1 ; ii++)
1091         (*CacheWeights)(ii) = 0.0e0 ;
1092       (*CacheWeights)(1) = 1.0e0 ;
1093     }
1094   }
1095 }
1096
1097 void BSplCLib::BuildCache(const Standard_Real          theParameter,
1098                           const Standard_Real          theSpanDomain,
1099                           const Standard_Boolean       thePeriodicFlag,
1100                           const Standard_Integer       theDegree,
1101                           const Standard_Integer       theSpanIndex,
1102                           const TColStd_Array1OfReal&  theFlatKnots,
1103                           const Array1OfPoints&        thePoles,
1104                           const TColStd_Array1OfReal*  theWeights,
1105                                 TColStd_Array2OfReal&  theCacheArray)
1106 {
1107   Standard_Real    aParam = theParameter;
1108   Standard_Integer anIndex = theSpanIndex;
1109   Standard_Integer aDimension;
1110   Standard_Boolean isRational;
1111   
1112   BSplCLib_DataContainer dc(theDegree);
1113   PrepareEval(aParam,
1114               anIndex,
1115               aDimension,
1116               isRational,
1117               theDegree,
1118               thePeriodicFlag,
1119               thePoles,
1120               theWeights,
1121               theFlatKnots,
1122               (BSplCLib::NoMults()),
1123               dc);
1124   //
1125   // Watch out : PrepareEval checks if locally the Bspline is polynomial
1126   // therefore rational flag could be set to False even if there are 
1127   // Weights. The Dimension is set accordingly.
1128   //
1129
1130   Standard_Integer aCacheShift = // helps to store weights when PrepareEval did not found that the curve is locally polynomial
1131     (theWeights != NULL && !isRational) ? aDimension + 1 : aDimension;
1132
1133   BSplCLib::Bohm(aParam, theDegree, theDegree, *dc.knots, aDimension, *dc.poles);
1134
1135   Standard_Real aCoeff = 1.0;
1136   Standard_Real* aCache = (Standard_Real *) &(theCacheArray(theCacheArray.LowerRow(), theCacheArray.LowerCol()));
1137   Standard_Real* aPolyCoeffs = dc.poles;
1138
1139   for (Standard_Integer i = 0; i <= theDegree; i++)
1140   {
1141     for (Standard_Integer j = 0; j< aDimension; j++)
1142       aCache[j] = aPolyCoeffs[j] * aCoeff;
1143     aCoeff *= theSpanDomain / (i+1);
1144     aPolyCoeffs += aDimension;
1145     aCache += aDimension;
1146     if (aCacheShift > aDimension)
1147     {
1148       aCache[0] = 0.0;
1149       aCache++;
1150     }
1151   }
1152
1153   if (aCacheShift > aDimension)
1154     theCacheArray.SetValue(theCacheArray.LowerRow(), theCacheArray.LowerCol() + aCacheShift - 1, 1.0);
1155 }
1156
1157 //=======================================================================
1158 //function : Interpolate
1159 //purpose  : 
1160 //=======================================================================
1161
1162 void  BSplCLib::Interpolate(const Standard_Integer         Degree,
1163                             const TColStd_Array1OfReal&    FlatKnots,
1164                             const TColStd_Array1OfReal&    Parameters,
1165                             const TColStd_Array1OfInteger& ContactOrderArray,
1166                             Array1OfPoints&                Poles,
1167                             Standard_Integer&              InversionProblem) 
1168      
1169 {
1170   Standard_Real *PArray ;
1171   
1172   PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1173   
1174   BSplCLib::Interpolate(Degree,
1175                         FlatKnots,
1176                         Parameters,
1177                         ContactOrderArray,
1178                         Dimension_gen,
1179                         PArray[0],
1180                         InversionProblem) ;
1181 }
1182
1183 //=======================================================================
1184 //function : Interpolate
1185 //purpose  : 
1186 //=======================================================================
1187
1188 void  BSplCLib::Interpolate(const Standard_Integer         Degree,
1189                             const TColStd_Array1OfReal&    FlatKnots,
1190                             const TColStd_Array1OfReal&    Parameters,
1191                             const TColStd_Array1OfInteger& ContactOrderArray,
1192                             Array1OfPoints&                Poles,
1193                             TColStd_Array1OfReal&          Weights,
1194                             Standard_Integer&              InversionProblem) 
1195 {
1196   Standard_Real *PArray,
1197   * WArray ;
1198   PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1199   WArray = (Standard_Real *) &Weights(Weights.Lower()) ;
1200   BSplCLib::Interpolate(Degree,
1201                         FlatKnots,
1202                         Parameters,
1203                         ContactOrderArray,
1204                         Dimension_gen,
1205                         PArray[0],
1206                         WArray[0],
1207                         InversionProblem) ;
1208 }
1209
1210 //=======================================================================
1211 //function : MovePoint
1212 //purpose  : Find the new poles which allows  an old point (with a
1213 //           given  u as parameter) to reach a new position
1214 //=======================================================================
1215
1216 void BSplCLib::MovePoint (const Standard_Real            U, 
1217                           const Vector&                  Displ,
1218                           const Standard_Integer         Index1,
1219                           const Standard_Integer         Index2,
1220                           const Standard_Integer         Degree,
1221                           const Standard_Boolean         Rational,
1222                           const Array1OfPoints&          Poles,  
1223                           const TColStd_Array1OfReal&    Weights,
1224                           const TColStd_Array1OfReal&    FlatKnots,
1225                           Standard_Integer&              FirstIndex,
1226                           Standard_Integer&              LastIndex,
1227                           Array1OfPoints&                NewPoles)
1228 {
1229   // calculate the BSplineBasis in the parameter U
1230   Standard_Integer FirstNonZeroBsplineIndex;
1231   math_Matrix BSplineBasis(1, 1,
1232                            1, Degree+1);
1233   Standard_Integer ErrorCode =
1234     BSplCLib::EvalBsplineBasis(0,
1235                                Degree+1,
1236                                FlatKnots,
1237                                U,
1238                                FirstNonZeroBsplineIndex,
1239                                BSplineBasis);  
1240   if (ErrorCode != 0) {
1241     FirstIndex = 0;
1242     LastIndex = 0;
1243
1244     for (Standard_Integer i=Poles.Lower(); i<=Poles.Upper(); i++) {
1245       NewPoles(i) = Poles(i);
1246     }
1247     return;
1248   }
1249   
1250   // find the span which is predominant for parameter U
1251   FirstIndex = FirstNonZeroBsplineIndex;
1252   LastIndex = FirstNonZeroBsplineIndex + Degree ;
1253   if (FirstIndex < Index1) FirstIndex = Index1;
1254   if (LastIndex > Index2) LastIndex = Index2;
1255   
1256   Standard_Real maxValue = 0.0;
1257   Standard_Integer i, kk1=0, kk2, ii;
1258
1259   for (i = FirstIndex-FirstNonZeroBsplineIndex+1;
1260        i <= LastIndex-FirstNonZeroBsplineIndex+1; i++) {
1261     if (BSplineBasis(1,i) > maxValue) {
1262       kk1 = i + FirstNonZeroBsplineIndex - 1;
1263       maxValue = BSplineBasis(1,i);
1264     }
1265   }
1266   
1267   // find a kk2 if symetriy
1268   kk2 = kk1;
1269   i = kk1 - FirstNonZeroBsplineIndex + 2;
1270   if ((kk1+1) <= LastIndex) {
1271     if (Abs(BSplineBasis(1, kk1-FirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
1272       kk2 = kk1+1;
1273     }
1274   }
1275   
1276   // compute the vector of displacement
1277   Standard_Real D1 = 0.0;
1278   Standard_Real D2 = 0.0;
1279   Standard_Real hN, Coef, Dval;
1280   
1281   for (i = 1; i <= Degree+1; i++) {
1282     ii = i + FirstNonZeroBsplineIndex - 1;
1283     if (Rational) {
1284       hN = Weights(ii)*BSplineBasis(1, i);
1285       D2 += hN;
1286     }
1287     else {
1288       hN = BSplineBasis(1, i);
1289     }
1290     if (ii >= FirstIndex && ii <= LastIndex) {
1291       if (ii < kk1) {
1292         Dval = kk1-ii;
1293       }
1294       else if (ii > kk2) {
1295         Dval = ii - kk2;
1296       }
1297       else {
1298         Dval = 0.0;
1299       }
1300       D1 += 1./(Dval + 1.) * hN;
1301     }
1302   }
1303   
1304   if (Rational) {
1305     Coef = D2/D1;
1306   }
1307   else {
1308     Coef = 1./D1;
1309   }
1310   
1311   // compute the new poles
1312
1313   for (i=Poles.Lower(); i<=Poles.Upper(); i++) {
1314     if (i >= FirstIndex && i <= LastIndex) {
1315       if (i < kk1) {
1316         Dval = kk1-i;
1317       }
1318       else if (i > kk2) {
1319         Dval = i - kk2;
1320       }
1321       else {
1322         Dval = 0.0;
1323       }
1324       NewPoles(i) = Poles(i).Translated((Coef/(Dval+1.))*Displ);
1325     }
1326     else {
1327       NewPoles(i) = Poles(i);
1328     }
1329   }
1330 }
1331
1332 //=======================================================================
1333 //function : MovePoint
1334 //purpose  : Find the new poles which allows  an old point (with a
1335 //           given  u as parameter) to reach a new position
1336 //=======================================================================
1337
1338 //=======================================================================
1339 //function : MovePointAndTangent
1340 //purpose  : 
1341 //=======================================================================
1342
1343 void BSplCLib::MovePointAndTangent (const Standard_Real      U, 
1344                                     const Vector&            Delta,
1345                                     const Vector&            DeltaDerivatives,
1346                                     const Standard_Real      Tolerance,
1347                                     const Standard_Integer   Degree,
1348                                     const Standard_Boolean   Rational,
1349                                     const Standard_Integer   StartingCondition,
1350                                     const Standard_Integer   EndingCondition,
1351                                     const Array1OfPoints&    Poles, 
1352                                     const TColStd_Array1OfReal&    Weights,
1353                                     const TColStd_Array1OfReal&    FlatKnots,
1354                                     Array1OfPoints&                NewPoles,
1355                                     Standard_Integer &             ErrorStatus)
1356 {
1357   Standard_Real *delta_array,
1358   *delta_derivative_array,
1359   *poles_array,
1360   *new_poles_array ;
1361   
1362   Standard_Integer num_poles ;
1363   num_poles = Poles.Length() ;
1364   
1365   if (NewPoles.Length() != num_poles) {
1366     throw Standard_ConstructionError();
1367   }
1368   delta_array = (Standard_Real *)  &Delta ;
1369   delta_derivative_array = (Standard_Real *)&DeltaDerivatives ;
1370   poles_array = (Standard_Real *) &Poles(Poles.Lower()) ;
1371   
1372   new_poles_array = (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1373   MovePointAndTangent(U,
1374                       Dimension_gen,
1375                       delta_array[0],
1376                       delta_derivative_array[0],
1377                       Tolerance,
1378                       Degree,
1379                       Rational,
1380                       StartingCondition,
1381                       EndingCondition,
1382                       poles_array[0],
1383                       Weights,
1384                       FlatKnots,
1385                       new_poles_array[0],
1386                       ErrorStatus) ;
1387 }
1388
1389 //=======================================================================
1390 //function : Resolution
1391 //purpose  : 
1392 //=======================================================================
1393
1394 void BSplCLib::Resolution(const Array1OfPoints&          Poles,  
1395                           const TColStd_Array1OfReal*    Weights,
1396                           const Standard_Integer         NumPoles,
1397                           const TColStd_Array1OfReal&    FlatKnots, 
1398                           const Standard_Integer         Degree,
1399                           const Standard_Real            Tolerance3D,
1400                           Standard_Real&                 UTolerance) 
1401 {
1402   Standard_Real *PolesArray ;
1403   PolesArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1404   BSplCLib::Resolution(PolesArray[0],
1405                        Dimension_gen,
1406                        NumPoles,
1407                        Weights,
1408                        FlatKnots,
1409                        Degree,
1410                        Tolerance3D,
1411                        UTolerance) ;
1412 }
1413
1414 //=======================================================================
1415 //function : FunctionMultiply
1416 //purpose  : 
1417 //=======================================================================
1418
1419 void BSplCLib::FunctionMultiply
1420 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1421  const Standard_Integer              BSplineDegree,
1422  const TColStd_Array1OfReal &        BSplineFlatKnots,
1423  const Array1OfPoints &              Poles,
1424  const TColStd_Array1OfReal &        FlatKnots,
1425  const Standard_Integer              NewDegree,
1426  Array1OfPoints &                   NewPoles,
1427  Standard_Integer &                  theStatus)
1428
1429   Standard_Integer num_bspline_poles =
1430     BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1431   Standard_Integer num_new_poles =
1432     FlatKnots.Length() - NewDegree - 1 ;
1433   
1434   if (Poles.Length() != num_bspline_poles ||
1435       NewPoles.Length() != num_new_poles) {
1436     throw Standard_ConstructionError();
1437   }
1438   Standard_Real  * array_of_poles =
1439     (Standard_Real *) &Poles(Poles.Lower()) ;
1440   Standard_Real  * array_of_new_poles =
1441     (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1442   BSplCLib::FunctionMultiply(FunctionPtr,
1443                              BSplineDegree,
1444                              BSplineFlatKnots,
1445                              Dimension_gen,
1446                              array_of_poles[0],
1447                              FlatKnots,
1448                              NewDegree,
1449                              array_of_new_poles[0],
1450                              theStatus);
1451 }
1452
1453 //=======================================================================
1454 //function : FunctionReparameterise
1455 //purpose  : 
1456 //=======================================================================
1457
1458 void BSplCLib::FunctionReparameterise
1459 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1460  const Standard_Integer              BSplineDegree,
1461  const TColStd_Array1OfReal &        BSplineFlatKnots,
1462  const Array1OfPoints &              Poles,
1463  const TColStd_Array1OfReal &        FlatKnots,
1464  const Standard_Integer              NewDegree,
1465  Array1OfPoints &                   NewPoles,
1466  Standard_Integer &                  theStatus)
1467
1468   Standard_Integer num_bspline_poles =
1469     BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1470   Standard_Integer num_new_poles =
1471     FlatKnots.Length() - NewDegree - 1 ;
1472   
1473   if (Poles.Length() != num_bspline_poles ||
1474       NewPoles.Length() != num_new_poles) {
1475     throw Standard_ConstructionError();
1476   }
1477   Standard_Real  * array_of_poles =
1478     (Standard_Real *) &Poles(Poles.Lower()) ;
1479   Standard_Real  * array_of_new_poles =
1480     (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1481   BSplCLib::FunctionReparameterise(FunctionPtr,
1482                                    BSplineDegree,
1483                                    BSplineFlatKnots,
1484                                    Dimension_gen,
1485                                    array_of_poles[0],
1486                                    FlatKnots,
1487                                    NewDegree,
1488                                    array_of_new_poles[0],
1489                                    theStatus);
1490 }