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