0023909: FPE in BSplSLib::RationalDerivative
[occt.git] / src / BSplSLib / BSplSLib.cxx
1 // Created on: 1991-08-26
2 // Created by: JCV
3 // Copyright (c) 1991-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
5 //
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
10 //
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13 //
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
20
21
22 // Modifed RLE Aug 93 - Complete rewrite
23 // xab  21-Mar-95  implemented cache mecanism
24 // pmn  25-09-96   Interpolation
25 // jct  25-09-96 : Correction de l'alloc de LocalArray dans RationalDerivative.
26 // pmn  07-10-96 : Correction de DN dans le cas rationnal.
27 // pmn  06-02-97 : Correction des poids dans RationalDerivative. (PRO700)
28
29 #include <BSplSLib.ixx>
30 #include <PLib.hxx>
31 #include <PLib_LocalArray.hxx>
32 #include <BSplCLib.hxx>
33 #include <TColgp_Array2OfXYZ.hxx>
34 #include <TColgp_Array1OfXYZ.hxx>
35 #include <TColStd_HArray1OfInteger.hxx>
36 #include <Standard_NotImplemented.hxx>
37 #include <Standard_ConstructionError.hxx>
38 #include <math_Matrix.hxx>
39
40 // for null derivatives
41 static Standard_Real BSplSLib_zero[3] = {0.0, 0.0, 0.0};
42 #ifdef WNT
43 #define M_SQRT2 1.4142135623730950488016887 
44 #endif
45
46 //=======================================================================
47 //struct : BSplCLib_DataContainer 
48 //purpose: Auxiliary structure providing buffers for poles and knots used in
49 //         evaluation of bspline (allocated in the stack)
50 //=======================================================================
51
52 struct BSplSLib_DataContainer
53 {
54   BSplSLib_DataContainer (Standard_Integer UDegree, Standard_Integer VDegree)
55   {
56     Standard_OutOfRange_Raise_if (UDegree > BSplCLib::MaxDegree() ||
57         VDegree > BSplCLib::MaxDegree() || BSplCLib::MaxDegree() > 25,
58         "BSplSLib: bspline degree is greater than maximum supported");
59   }
60
61   Standard_Real poles[4*(25+1)*(25+1)];
62   Standard_Real knots1[2*25];
63   Standard_Real knots2[2*25];
64   Standard_Real ders[48];
65 };
66
67 typedef PLib_LocalArray BSplSLib_LocalArray;
68
69 //**************************************************************************
70 //                     Evaluation methods
71 //**************************************************************************
72
73 //=======================================================================
74 //function : RationalDerivative
75 //purpose  : computes the rational derivatives when whe have the
76 //           the derivatives of the homogeneous numerator and the
77 //           the derivatives of the denominator 
78 //=======================================================================
79
80 void  BSplSLib::RationalDerivative(const Standard_Integer UDeg,
81                                    const Standard_Integer VDeg,
82                                    const Standard_Integer N, 
83                                    const Standard_Integer M, 
84                                    Standard_Real& HDerivatives,
85                                    Standard_Real& RDerivatives,
86                                    const Standard_Boolean All)
87 {
88   //
89   //  if All is True all derivatives are computed. if Not only
90   //  the requested N, M is computed 
91   //
92   //                                           Numerator(u,v) 
93   //   let f(u,v) be a rational function  = ------------------
94   //                                         Denominator(u,v)
95   //
96   //
97   //   Let (N,M) the order of the derivatives we want : then since
98   //   we have :
99   //
100   //         Numerator = f * Denominator 
101   //
102   //   we derive :
103   //
104   //   (N,M)         1           (         (N M)                  (p q)            (N -p M-q) )
105   //  f       =  ------------    (  Numerator   -  SUM SUM  a   * f    * Denominator          )
106   //                       (0,0) (                 p<N q<M   p q                              )
107   //             Denominator
108   //
109   //   with :
110   //
111   //              ( N )  ( M )
112   //    a      =  (   )  (   )
113   //     p q      ( p )  ( q )
114   //
115   //  
116   //   HDerivatives is an array where derivatives are stored in the following form
117   //   Numerator is assumee to have 3 functions that is a vector of dimension
118   //   3
119   // 
120   //             (0,0)           (0,0)                       (0, DegV)           (0, DegV) 
121   //     Numerator     Denominator      ...         Numerator          Denominator
122   //
123   //             (1,0)           (1,0)                       (1, DegV)           (1, DegV) 
124   //     Numerator     Denominator      ...         Numerator          Denominator
125   //
126   //         ...........................................................
127   //
128   //
129   //            (DegU,0)        (DegU,0)                  (DegU, DegV)           (DegU, DegV) 
130   //     Numerator     Denominator      ...         Numerator          Denominator
131   //
132   //
133   Standard_Integer ii,jj,pp,qq,index,index1,index2;
134   Standard_Integer M1,M3,M4,N1,iiM1,iiM3,jjM1,ppM1,ppM3;
135   Standard_Integer MinN,MinN1,MinM,MinM1;
136   Standard_Integer index_u,index_u1,index_v,index_v1,index_w;
137
138   M1 = M + 1;
139   N1 = N + 1;
140   ii = N1 * M1;
141   M3 = (M1 << 1) + M1;
142   M4 = (VDeg + 1) << 2;
143   
144   BSplSLib_LocalArray StoreDerivatives (All ? 0 : ii * 3);
145   Standard_Real *RArray = (All ? &RDerivatives : (Standard_Real*)StoreDerivatives);
146   BSplSLib_LocalArray StoreW (ii);  
147   Standard_Real *HomogeneousArray = &HDerivatives; 
148   Standard_Real denominator,Pii,Pip,Pjq;
149   
150   denominator = 1.0e0 / HomogeneousArray[3];
151   index_u  = 0;
152   index_u1 = 0;
153   if (UDeg < N) MinN = UDeg;
154   else          MinN = N;
155   if (VDeg < M) MinM = VDeg;
156   else          MinM = M;
157   MinN1 = MinN + 1;
158   MinM1 = MinM + 1;
159   iiM1 = - M1;
160
161   for (ii = 0 ; ii < MinN1 ; ii++) {
162     iiM1    += M1;
163     index_v  = index_u;
164     index_v1 = index_u1;
165     index_w  = iiM1;
166     
167     for (jj = 0 ; jj < MinM1 ; jj++) {
168       RArray[index_v++] = HomogeneousArray[index_v1++];
169       RArray[index_v++] = HomogeneousArray[index_v1++];
170       RArray[index_v++] = HomogeneousArray[index_v1++];
171       StoreW[index_w++] = HomogeneousArray[index_v1++];
172     }
173
174     for (jj = MinM1 ; jj < M1 ; jj++) {
175       RArray[index_v++] = 0.;
176       RArray[index_v++] = 0.;
177       RArray[index_v++] = 0.;
178       StoreW[index_w++] = 0.;
179     }
180     index_u1 += M4;
181     index_u  += M3;
182   }
183   index_v = MinN1 * M3;
184   index_w = MinN1 * M1;
185   
186   for (ii = MinN1 ; ii < N1 ; ii++) {
187     
188     for (jj = 0 ; jj < M1 ; jj++) {  
189       RArray[index_v++] = 0.0e0;
190       RArray[index_v++] = 0.0e0;
191       RArray[index_v++] = 0.0e0;
192       StoreW[index_w++] = 0.0e0;
193     }
194   } 
195
196   // ---------------  Calculation ----------------
197
198   iiM1 = - M1;
199   iiM3 = - M3;
200
201   for (ii = 0 ; ii <= N  ; ii++) {
202     iiM1  += M1;
203     iiM3  += M3;
204     index1 = iiM3 - 3;
205     jjM1   = iiM1;
206     
207     for (jj = 0 ; jj <= M ; jj++) {
208       jjM1 ++;
209       ppM1    = - M1;
210       ppM3    = - M3;
211       index1 += 3;
212
213       for (pp = 0 ; pp < ii ; pp++) {
214         ppM1  += M1;
215         ppM3  += M3;
216         index  = ppM3;
217         index2 = jjM1 - ppM1;
218         Pip    = PLib::Bin(ii,pp);
219         
220         for (qq = 0 ; qq <= jj ; qq++) {
221           index2--;
222           Pjq    = Pip * PLib::Bin(jj,qq) * StoreW[index2]; 
223           RArray[index1] -= Pjq * RArray[index]; index++; index1++;
224           RArray[index1] -= Pjq * RArray[index]; index++; index1++;
225           RArray[index1] -= Pjq * RArray[index]; index++;
226           index1 -= 2;
227         }
228       }
229       index  = iiM3;
230       index2 = jj + 1;
231       Pii = PLib::Bin(ii,ii);
232
233       for (qq = 0 ; qq < jj ; qq++) {
234         index2--;
235         Pjq = Pii * PLib::Bin(jj,qq) * StoreW[index2]; 
236         RArray[index1] -= Pjq * RArray[index]; index++; index1++;
237         RArray[index1] -= Pjq * RArray[index]; index++; index1++;
238         RArray[index1] -= Pjq * RArray[index]; index++;
239         index1 -= 2;
240       }
241       RArray[index1] *= denominator; index1++;
242       RArray[index1] *= denominator; index1++;
243       RArray[index1] *= denominator;
244       index1 -= 2;
245     }
246   }
247   if (!All) {
248     RArray = &RDerivatives;
249     index = N * M1 + M;
250     index = (index << 1) + index;
251     RArray[0] = StoreDerivatives[index]; index++;
252     RArray[1] = StoreDerivatives[index]; index++;
253     RArray[2] = StoreDerivatives[index];
254   }
255 }
256
257 //=======================================================================
258 //function : PrepareEval
259 //purpose  : 
260 //=======================================================================
261
262 //
263 // PrepareEval :
264 //
265 // Prepare all data for computing points :
266 //  local arrays of knots
267 //  local array  of poles (multiplied by the weights if rational)
268 //
269 //  The first direction to compute (smaller degree) is returned 
270 //  and the poles are stored according to this direction.
271
272 static Standard_Boolean  PrepareEval
273 (const Standard_Real            U,
274  const Standard_Real            V,
275  const Standard_Integer         Uindex,
276  const Standard_Integer         Vindex,
277  const Standard_Integer         UDegree,
278  const Standard_Integer         VDegree,
279  const Standard_Boolean         URat,
280  const Standard_Boolean         VRat,
281  const Standard_Boolean         UPer,
282  const Standard_Boolean         VPer,
283  const TColgp_Array2OfPnt&      Poles,
284  const TColStd_Array2OfReal&    Weights,
285  const TColStd_Array1OfReal&    UKnots,
286  const TColStd_Array1OfReal&    VKnots,
287  const TColStd_Array1OfInteger& UMults,
288  const TColStd_Array1OfInteger& VMults,
289  Standard_Real& u1,     // first  parameter to use
290  Standard_Real& u2,     // second parameter to use
291  Standard_Integer& d1,  // first degree
292  Standard_Integer& d2,  // second degree
293  Standard_Boolean& rational,
294  BSplSLib_DataContainer& dc)
295 {
296   rational = URat || VRat;
297   Standard_Integer uindex = Uindex;
298   Standard_Integer vindex = Vindex;
299   Standard_Integer UKLower = UKnots.Lower();
300   Standard_Integer UKUpper = UKnots.Upper();
301   Standard_Integer VKLower = VKnots.Lower();
302   Standard_Integer VKUpper = VKnots.Upper();
303   if (UDegree <= VDegree) {
304     // compute the indices
305     if (uindex < UKLower || uindex > UKUpper)
306       BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u1);
307     else u1 = U;
308     if (vindex < VKLower || vindex > VKUpper)
309       BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u2);
310     else u2 = V;
311     // get the knots
312     d1 = UDegree;
313     d2 = VDegree;
314     BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots1);
315     BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots2);
316     if (&UMults == NULL) uindex -= UKLower + UDegree;
317     else                 uindex  = BSplCLib::PoleIndex
318       (UDegree,uindex,UPer,UMults);
319     if (&VMults == NULL) vindex -= VKLower + VDegree;
320     else                 vindex  = BSplCLib::PoleIndex
321       (VDegree,vindex,VPer,VMults);
322     // get the poles
323 //    Standard_Integer i,j,k,ip,jp;
324     Standard_Integer i,j,ip,jp;
325     Standard_Real w, *pole = dc.poles;
326     d1 = UDegree;
327     d2 = VDegree;
328     Standard_Integer PLowerRow = Poles.LowerRow();
329     Standard_Integer PUpperRow = Poles.UpperRow();
330     Standard_Integer PLowerCol = Poles.LowerCol();
331     Standard_Integer PUpperCol = Poles.UpperCol();
332     if (rational) { // verify if locally non rational
333       rational = Standard_False;
334       ip = PLowerRow + uindex;
335       jp = PLowerCol + vindex;
336       w  = Weights.Value(ip,jp);
337       Standard_Real eps = Epsilon(w);
338       Standard_Real dw;
339       
340       for (i = 0; i <= UDegree && !rational; i++) {
341         jp = PLowerCol + vindex;
342         
343         for (j = 0; j <= VDegree && !rational; j++) {
344           dw = Weights.Value(ip,jp) - w;
345           if (dw < 0) dw = - dw;
346           rational = dw > eps;
347           jp++;
348           if (jp > PUpperCol) jp = PLowerCol;
349         }
350         ip++;
351         if (ip > PUpperRow) ip = PLowerRow;
352       }
353     }
354     // copy the poles
355     ip = PLowerRow + uindex;
356     if (rational) {
357     
358       for (i = 0; i <= d1; i++) {
359         jp = PLowerCol + vindex;
360         
361         for (j = 0; j <= d2; j++) {
362           const gp_Pnt& P = Poles  .Value(ip,jp);
363           pole[3] = w     = Weights.Value(ip,jp);
364           pole[0] = P.X() * w;
365           pole[1] = P.Y() * w;
366           pole[2] = P.Z() * w;
367           pole   += 4;
368           jp++;
369           if (jp > PUpperCol) jp = PLowerCol;
370         }
371         ip++;
372         if (ip > PUpperRow) ip = PLowerRow;
373       }
374     }
375     else {
376     
377       for (i = 0; i <= d1; i++) {
378         jp = PLowerCol + vindex;
379         
380         for (j = 0; j <= d2; j++) {
381           const gp_Pnt& P = Poles.Value(ip,jp);
382           pole[0] = P.X();
383           pole[1] = P.Y();
384           pole[2] = P.Z();
385           pole   += 3;
386           jp++;
387           if (jp > PUpperCol) jp = PLowerCol;
388         }
389         ip++;
390         if (ip > PUpperRow) ip = PLowerRow;
391       }
392     }
393     return Standard_True;
394   }
395   else {
396     // compute the indices
397     if (uindex < UKLower || uindex > UKUpper)
398       BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u2);
399     else u2 = U;
400     if (vindex < VKLower || vindex > VKUpper)
401       BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u1);
402     else u1 = V;
403     // get the knots
404     d2 = UDegree;
405     d1 = VDegree;
406     BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots2);
407     BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots1);
408     if (&UMults == NULL) uindex -= UKLower + UDegree;
409     else                 uindex  = BSplCLib::PoleIndex
410       (UDegree,uindex,UPer,UMults);
411     if (&VMults == NULL) vindex -= VKLower + VDegree;
412     else                 vindex  = BSplCLib::PoleIndex
413       (VDegree,vindex,VPer,VMults);
414     // get the poles
415 //    Standard_Integer i,j,k,ip,jp;
416     Standard_Integer i,j,ip,jp;
417     Standard_Real w, *pole = dc.poles;
418     d1 = VDegree;
419     d2 = UDegree;
420     Standard_Integer PLowerRow = Poles.LowerRow();
421     Standard_Integer PUpperRow = Poles.UpperRow();
422     Standard_Integer PLowerCol = Poles.LowerCol();
423     Standard_Integer PUpperCol = Poles.UpperCol();
424     if (rational) { // verify if locally non rational
425       rational = Standard_False;
426       ip = PLowerRow + uindex;
427       jp = PLowerCol + vindex;
428       w  = Weights.Value(ip,jp);
429       Standard_Real eps = Epsilon(w);
430       Standard_Real dw;
431       
432       for (i = 0; i <= UDegree && !rational; i++) {
433         jp = PLowerCol + vindex;
434         
435         for (j = 0; j <= VDegree && !rational; j++) {
436           dw = Weights.Value(ip,jp) - w;
437           if (dw < 0) dw = - dw;
438           rational = dw > eps;
439           jp++;
440           if (jp > PUpperCol) jp = PLowerCol;
441         }
442         ip++;
443         if (ip > PUpperRow) ip = PLowerRow;
444       }
445     }
446     // copy the poles
447     jp = PLowerCol + vindex;
448     if (rational) {
449     
450       for (i = 0; i <= d1; i++) {
451         ip = PLowerRow + uindex;
452         
453         for (j = 0; j <= d2; j++) {
454           const gp_Pnt& P = Poles  .Value(ip,jp);
455           pole[3] = w     = Weights.Value(ip,jp);
456           pole[0] = P.X() * w;
457           pole[1] = P.Y() * w;
458           pole[2] = P.Z() * w;
459           pole   += 4;
460           ip++;
461           if (ip > PUpperRow) ip = PLowerRow;
462         }
463         jp++;
464         if (jp > PUpperCol) jp = PLowerCol;
465       }
466     }
467     else {
468     
469       for (i = 0; i <= d1; i++) {
470         ip = PLowerRow + uindex;
471         
472         for (j = 0; j <= d2; j++) {
473           const gp_Pnt& P = Poles.Value(ip,jp);
474           pole[0] = P.X();
475           pole[1] = P.Y();
476           pole[2] = P.Z();
477           pole   += 3;
478           ip++;
479           if (ip > PUpperRow) ip = PLowerRow;
480         }
481         jp++;
482         if (jp > PUpperCol) jp = PLowerCol;
483       }
484     }
485     return Standard_False;
486   }
487 }
488
489 //=======================================================================
490 //function : D0
491 //purpose  : 
492 //=======================================================================
493
494 void  BSplSLib::D0
495 (const Standard_Real U, 
496  const Standard_Real V, 
497  const Standard_Integer UIndex, 
498  const Standard_Integer VIndex,
499  const TColgp_Array2OfPnt& Poles,
500  const TColStd_Array2OfReal& Weights,
501  const TColStd_Array1OfReal& UKnots,
502  const TColStd_Array1OfReal& VKnots,
503  const TColStd_Array1OfInteger& UMults,
504  const TColStd_Array1OfInteger& VMults,
505  const Standard_Integer UDegree,
506  const Standard_Integer VDegree,
507  const Standard_Boolean URat,
508  const Standard_Boolean VRat,
509  const Standard_Boolean UPer,
510  const Standard_Boolean VPer,
511  gp_Pnt& P)
512 {
513 //  Standard_Integer k ;
514   Standard_Real W ;
515   HomogeneousD0(U,
516                 V, 
517                 UIndex, 
518                 VIndex,
519                 Poles,
520                 Weights,
521                 UKnots,
522                 VKnots,
523                 UMults,
524                 VMults,
525                 UDegree,
526                 VDegree,
527                 URat,
528                 VRat,
529                 UPer,
530                 VPer,
531                 W,
532                 P) ;
533   P.SetX(P.X() / W);
534   P.SetY(P.Y() / W);
535   P.SetZ(P.Z() / W);
536 }
537
538 //=======================================================================
539 //function : D0
540 //purpose  : 
541 //=======================================================================
542
543 void  BSplSLib::HomogeneousD0
544 (const Standard_Real U, 
545  const Standard_Real V, 
546  const Standard_Integer UIndex, 
547  const Standard_Integer VIndex,
548  const TColgp_Array2OfPnt& Poles,
549  const TColStd_Array2OfReal& Weights,
550  const TColStd_Array1OfReal& UKnots,
551  const TColStd_Array1OfReal& VKnots,
552  const TColStd_Array1OfInteger& UMults,
553  const TColStd_Array1OfInteger& VMults,
554  const Standard_Integer UDegree,
555  const Standard_Integer VDegree,
556  const Standard_Boolean URat,
557  const Standard_Boolean VRat,
558  const Standard_Boolean UPer,
559  const Standard_Boolean VPer,
560  Standard_Real & W,
561  gp_Pnt& P)
562 {
563   Standard_Boolean rational;
564 //  Standard_Integer k,dim;
565   Standard_Integer dim;
566   Standard_Real u1,u2;
567   Standard_Integer d1,d2;
568   W = 1.0e0 ;
569   
570   BSplSLib_DataContainer dc (UDegree, VDegree);
571   PrepareEval(U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
572               Poles,Weights,UKnots,VKnots,UMults,VMults,
573               u1,u2,d1,d2,rational,dc);
574   if (rational) {
575     dim = 4;
576     BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
577     BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
578     W = dc.poles[3];
579     P.SetX(dc.poles[0]);
580     P.SetY(dc.poles[1]);
581     P.SetZ(dc.poles[2]);
582   }
583   else {
584     dim = 3;
585     BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
586     BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
587     P.SetX(dc.poles[0]);
588     P.SetY(dc.poles[1]);
589     P.SetZ(dc.poles[2]);
590   }
591 }
592
593 //=======================================================================
594 //function : D1
595 //purpose  : 
596 //=======================================================================
597
598 void  BSplSLib::D1
599 (const Standard_Real U,
600  const Standard_Real V,
601  const Standard_Integer UIndex,
602  const Standard_Integer VIndex,
603  const TColgp_Array2OfPnt& Poles,
604  const TColStd_Array2OfReal& Weights,
605  const TColStd_Array1OfReal& UKnots,
606  const TColStd_Array1OfReal& VKnots,
607  const TColStd_Array1OfInteger& UMults,
608  const TColStd_Array1OfInteger& VMults,
609  const Standard_Integer UDegree,
610  const Standard_Integer VDegree,
611  const Standard_Boolean URat,
612  const Standard_Boolean VRat,
613  const Standard_Boolean UPer,
614  const Standard_Boolean VPer,
615  gp_Pnt& P,
616  gp_Vec& Vu,
617  gp_Vec& Vv)
618 {
619   Standard_Boolean rational;
620 //  Standard_Integer k,dim,dim2;
621   Standard_Integer dim,dim2;
622   Standard_Real u1,u2;
623   Standard_Integer d1,d2;
624   Standard_Real *result, *resVu, *resVv;
625   BSplSLib_DataContainer dc (UDegree, VDegree);
626   if (PrepareEval
627     (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
628      Poles,Weights,UKnots,VKnots,UMults,VMults,
629      u1,u2,d1,d2,rational,dc)) {
630     if (rational) {
631       dim  = 4;
632       dim2 = (d2 + 1) << 2;
633       BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
634       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
635       BSplCLib::Eval(u2,d2,  *dc.knots2,dim ,*(dc.poles + dim2));
636       BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
637       result = dc.ders;
638       resVu  = result + 6;
639       resVv  = result + 3;
640     }
641     else {
642       dim  = 3;
643       dim2 = d2 + 1;
644       dim2 = (dim2 << 1) + dim2;
645       BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
646       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
647       BSplCLib::Eval(u2,d2,  *dc.knots2,dim ,*(dc.poles + dim2));
648       result = dc.poles;
649       resVu  = result + dim2;
650       resVv  = result + 3;
651     }
652   }
653   else {
654     if (rational) {
655       dim  = 4;
656       dim2 = (d2 + 1) << 2;
657       BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
658       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
659       BSplCLib::Eval(u2,d2,  *dc.knots2,dim ,*(dc.poles + dim2));
660       BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
661       result = dc.ders;
662       resVu  = result + 3;
663       resVv  = result + 6;
664     }
665     else {
666       dim  = 3;
667       dim2 = d2 + 1;
668       dim2 = (dim2 << 1) + dim2;
669       BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
670       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
671       BSplCLib::Eval(u2,d2  ,*dc.knots2,dim ,*(dc.poles + dim2));
672       result = dc.poles;
673       resVu  = result + 3;
674       resVv  = result + dim2;
675     }
676   }
677   
678   P .SetX(result[0]);
679   Vu.SetX(resVu [0]);
680   Vv.SetX(resVv [0]);
681   
682   P .SetY(result[1]);
683   Vu.SetY(resVu [1]);
684   Vv.SetY(resVv [1]);
685   
686   P .SetZ(result[2]);
687   Vu.SetZ(resVu [2]);
688   Vv.SetZ(resVv [2]);
689 }
690
691 //=======================================================================
692 //function : D1
693 //purpose  : 
694 //=======================================================================
695
696 void  BSplSLib::HomogeneousD1
697 (const Standard_Real U,
698  const Standard_Real V,
699  const Standard_Integer UIndex,
700  const Standard_Integer VIndex,
701  const TColgp_Array2OfPnt& Poles,
702  const TColStd_Array2OfReal& Weights,
703  const TColStd_Array1OfReal& UKnots,
704  const TColStd_Array1OfReal& VKnots,
705  const TColStd_Array1OfInteger& UMults,
706  const TColStd_Array1OfInteger& VMults,
707  const Standard_Integer UDegree,
708  const Standard_Integer VDegree,
709  const Standard_Boolean URat,
710  const Standard_Boolean VRat,
711  const Standard_Boolean UPer,
712  const Standard_Boolean VPer,
713  gp_Pnt& N,
714  gp_Vec& Nu,
715  gp_Vec& Nv,
716  Standard_Real& D,
717  Standard_Real& Du,
718  Standard_Real& Dv)
719 {
720   Standard_Boolean rational;
721 //  Standard_Integer k,dim;
722   Standard_Integer dim;
723   Standard_Real u1,u2;
724   Standard_Integer d1,d2;
725   
726   D = 1.0e0 ;
727   Du = 0.0e0 ;
728   Dv = 0.0e0 ;
729   BSplSLib_DataContainer dc (UDegree, VDegree);
730   Standard_Boolean ufirst = PrepareEval
731     (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
732      Poles,Weights,UKnots,VKnots,UMults,VMults,
733      u1,u2,d1,d2,rational,dc);
734   dim  = rational ? 4 : 3;
735   
736   BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim * (d2 + 1),*dc.poles);
737   BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim,*dc.poles);
738   BSplCLib::Eval(u2,d2,*dc.knots2,dim,*(dc.poles+dim*(d2+1)));
739   
740   Standard_Real *result, *resVu, *resVv;
741   result = dc.poles;
742   resVu  = result + (ufirst ? dim*(d2+1) : dim);
743   resVv  = result + (ufirst ? dim : dim*(d2+1));
744   
745   N .SetX(result[0]);
746   Nu.SetX(resVu [0]);
747   Nv.SetX(resVv [0]);
748   
749   N .SetY(result[1]);
750   Nu.SetY(resVu [1]);
751   Nv.SetY(resVv [1]);
752   
753   N .SetZ(result[2]);
754   Nu.SetZ(resVu [2]);
755   Nv.SetZ(resVv [2]);
756   
757   if (rational) {
758     D  = result[3];
759     Du = resVu [3];
760     Dv = resVv [3];
761   }
762 }
763
764 //=======================================================================
765 //function : D2
766 //purpose  : 
767 //=======================================================================
768
769 void  BSplSLib::D2
770 (const Standard_Real U,
771  const Standard_Real V,
772  const Standard_Integer UIndex,
773  const Standard_Integer VIndex,
774  const TColgp_Array2OfPnt& Poles,
775  const TColStd_Array2OfReal& Weights,
776  const TColStd_Array1OfReal& UKnots,
777  const TColStd_Array1OfReal& VKnots,
778  const TColStd_Array1OfInteger& UMults,
779  const TColStd_Array1OfInteger& VMults,
780  const Standard_Integer UDegree,
781  const Standard_Integer VDegree,
782  const Standard_Boolean URat,
783  const Standard_Boolean VRat,
784  const Standard_Boolean UPer,
785  const Standard_Boolean VPer,
786  gp_Pnt& P,
787  gp_Vec& Vu,
788  gp_Vec& Vv,
789  gp_Vec& Vuu,
790  gp_Vec& Vvv,
791  gp_Vec& Vuv)
792 {
793   Standard_Boolean rational;
794 //  Standard_Integer k,dim,dim2;
795   Standard_Integer dim,dim2;
796   Standard_Real u1,u2;
797   Standard_Integer d1,d2;
798   Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv;
799   BSplSLib_DataContainer dc (UDegree, VDegree);
800   if (PrepareEval
801       (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
802        Poles,Weights,UKnots,VKnots,UMults,VMults,
803        u1,u2,d1,d2,rational,dc)) {
804     if (rational) {
805       dim = 4;
806       dim2 = (d2 + 1) << 2;
807       BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
808       BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
809       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
810       if (d1 > 1)
811         BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
812       BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
813       result = dc.ders;
814       resVu  = result + 9;
815       resVv  = result + 3;
816       resVuu = result + 18;
817       resVvv = result + 6;
818       resVuv = result + 12;
819     }
820     else {
821       dim = 3;
822       dim2 = d2 + 1;
823       dim2 = (dim2 << 1) + dim2;
824       BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
825       BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
826       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
827       if (d1 > 1)
828         BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
829       result = dc.poles;
830       resVu  = result + dim2;
831       resVv  = result + 3;
832       if (UDegree <= 1) resVuu = BSplSLib_zero;
833       else              resVuu = result + (dim2 << 1);
834       if (VDegree <= 1) resVvv = BSplSLib_zero;
835       else              resVvv = result + 6;
836       resVuv = result + (d2 << 1) + d2 + 6; 
837     }
838   }
839   else {
840     if (rational) {
841       dim = 4;
842       dim2 = (d2 + 1) << 2;
843       BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
844       BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
845       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
846       if (d1 > 1)
847         BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
848       BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
849       result = dc.ders;
850       resVu  = result + 3;
851       resVv  = result + 9;
852       resVuu = result + 6;
853       resVvv = result + 18;
854       resVuv = result + 12;
855     }
856     else {
857       dim = 3;
858       dim2 = d2 + 1;
859       dim2 = (dim2 << 1) + dim2;
860       BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
861       BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
862       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
863       if (d1 > 1)
864         BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
865       result = dc.poles;
866       resVu  = result + 3;
867       resVv  = result + dim2;
868       if (UDegree <= 1) resVuu = BSplSLib_zero;
869       else              resVuu = result + 6;
870       if (VDegree <= 1) resVvv = BSplSLib_zero;
871       else              resVvv = result + (dim2 << 1);
872       resVuv = result + (d2 << 1) + d2 + 6; 
873     }
874   }
875
876   P  .SetX(result[0]);
877   Vu .SetX(resVu [0]);
878   Vv .SetX(resVv [0]);
879   Vuu.SetX(resVuu[0]);
880   Vvv.SetX(resVvv[0]);
881   Vuv.SetX(resVuv[0]);
882   
883   P  .SetY(result[1]);
884   Vu .SetY(resVu [1]);
885   Vv .SetY(resVv [1]);
886   Vuu.SetY(resVuu[1]);
887   Vvv.SetY(resVvv[1]);
888   Vuv.SetY(resVuv[1]);
889   
890   P  .SetZ(result[2]);
891   Vu .SetZ(resVu [2]);
892   Vv .SetZ(resVv [2]);
893   Vuu.SetZ(resVuu[2]);
894   Vvv.SetZ(resVvv[2]);
895   Vuv.SetZ(resVuv[2]);
896 }
897
898 //=======================================================================
899 //function : D3
900 //purpose  : 
901 //=======================================================================
902
903 void  BSplSLib::D3
904 (const Standard_Real U,
905  const Standard_Real V,
906  const Standard_Integer UIndex,
907  const Standard_Integer VIndex,
908  const TColgp_Array2OfPnt& Poles,
909  const TColStd_Array2OfReal& Weights,
910  const TColStd_Array1OfReal& UKnots,
911  const TColStd_Array1OfReal& VKnots,
912  const TColStd_Array1OfInteger& UMults,
913  const TColStd_Array1OfInteger& VMults,
914  const Standard_Integer UDegree,
915  const Standard_Integer VDegree,
916  const Standard_Boolean URat,
917  const Standard_Boolean VRat,
918  const Standard_Boolean UPer,
919  const Standard_Boolean VPer,
920  gp_Pnt& P,
921  gp_Vec& Vu,
922  gp_Vec& Vv,
923  gp_Vec& Vuu,
924  gp_Vec& Vvv,
925  gp_Vec& Vuv,
926  gp_Vec& Vuuu,
927  gp_Vec& Vvvv,
928  gp_Vec& Vuuv,
929  gp_Vec& Vuvv)
930 {
931   Standard_Boolean rational;
932 //  Standard_Integer k,dim,dim2;
933   Standard_Integer dim,dim2;
934   Standard_Real u1,u2;
935   Standard_Integer d1,d2;
936   Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv,
937   *resVuuu, *resVvvv, *resVuuv, *resVuvv;
938   BSplSLib_DataContainer dc (UDegree, VDegree);
939   if (PrepareEval
940     (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
941      Poles,Weights,UKnots,VKnots,UMults,VMults,
942      u1,u2,d1,d2,rational,dc)) {
943     if (rational) {
944       dim = 4;
945       dim2 = (d2 + 1) << 2;
946       BSplCLib::Bohm  (u1,d1,3,*dc.knots1,dim2,*dc.poles);
947       BSplCLib::Bohm  (u2,d2,3,*dc.knots2,dim ,*dc.poles);
948       BSplCLib::Bohm  (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
949       if (d1 > 1)
950         BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
951       if (d1 > 2)
952         BSplCLib::Eval(u2,d2  ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
953       BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
954       result  = dc.ders;
955       resVu   = result + 12;
956       resVv   = result + 3;
957       resVuu  = result + 24;
958       resVvv  = result + 6;
959       resVuv  = result + 15;
960       resVuuu = result + 36;
961       resVvvv = result + 9;
962       resVuuv = result + 27;
963       resVuvv = result + 18;
964     }
965     else {
966       dim = 3;
967       dim2 = (d2 + 1);
968       dim2 = (dim2 << 1) + dim2;
969       BSplCLib::Bohm  (u1,d1,3,*dc.knots1,dim2,*dc.poles);
970       BSplCLib::Bohm  (u2,d2,3,*dc.knots2,dim ,*dc.poles);
971       BSplCLib::Bohm  (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
972       if (d1 > 1)
973         BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
974       if (d1 > 2)
975         BSplCLib::Eval(u2,d2  ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
976       result = dc.poles;
977       resVu  = result + dim2;
978       resVv  = result + 3;
979       if (UDegree <= 1) {
980         resVuu  = BSplSLib_zero;
981         resVuuv = BSplSLib_zero;
982       }
983       else {
984         resVuu  = result + (dim2 << 1);
985         resVuuv = result + (dim2 << 1) + 3;
986       }
987       if (VDegree <= 1) {
988         resVvv  = BSplSLib_zero;
989         resVuvv = BSplSLib_zero;
990       }
991       else {
992         resVvv  = result + 6;
993         resVuvv = result + dim2 + 6;
994       }
995       resVuv = result + (d2 << 1) + d2 + 6;
996       if (UDegree <= 2) resVuuu = BSplSLib_zero;
997       else              resVuuu = result + (dim2 << 1) + dim2;
998       if (VDegree <= 2) resVvvv = BSplSLib_zero;
999       else              resVvvv = result + 9;
1000     }
1001   }
1002   else {
1003     if (rational) {
1004       dim = 4;
1005       dim2 = (d2 + 1) << 2;
1006       BSplCLib::Bohm  (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1007       BSplCLib::Bohm  (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1008       BSplCLib::Bohm  (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1009       if (d1 > 1)
1010         BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1011       if (d1 > 2)
1012         BSplCLib::Eval(u2,d2  ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1013       BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
1014       result  = dc.ders;
1015       resVu   = result + 3;
1016       resVv   = result + 12;
1017       resVuu  = result + 6;
1018       resVvv  = result + 24;
1019       resVuv  = result + 15;
1020       resVuuu = result + 9;
1021       resVvvv = result + 36;
1022       resVuuv = result + 18;
1023       resVuvv = result + 27;
1024     }
1025     else {
1026       dim = 3;
1027       dim2 = (d2 + 1);
1028       dim2 = (dim2 << 1) + dim2;
1029       BSplCLib::Bohm  (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1030       BSplCLib::Bohm  (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1031       BSplCLib::Bohm  (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1032       if (d1 > 1)
1033         BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1034       if (d1 > 2)
1035         BSplCLib::Eval(u2,d2  ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1036       result = dc.poles;
1037       resVu  = result + 3;
1038       resVv  = result + dim2;
1039       if (UDegree <= 1) {
1040         resVuu  = BSplSLib_zero;
1041         resVuuv = BSplSLib_zero;
1042       }
1043       else {
1044         resVuu  = result + 6;
1045         resVuuv = result + dim2 + 6;
1046       }
1047       if (VDegree <= 1) {
1048         resVvv  = BSplSLib_zero;
1049         resVuvv = BSplSLib_zero;
1050       }
1051       else {
1052         resVvv  = result + (dim2 << 1);
1053         resVuvv = result + (dim2 << 1) + 3;
1054       }
1055       resVuv = result + (d2 << 1) + d2 + 6;
1056       if (UDegree <= 2) resVuuu = BSplSLib_zero;
1057       else              resVuuu = result + 9;
1058       if (VDegree <= 2) resVvvv = BSplSLib_zero;
1059       else              resVvvv = result + (dim2 << 1) + dim2;
1060     }
1061   }
1062   
1063   P   .SetX(result [0]);
1064   Vu  .SetX(resVu  [0]);
1065   Vv  .SetX(resVv  [0]);
1066   Vuu .SetX(resVuu [0]);
1067   Vvv .SetX(resVvv [0]);
1068   Vuv .SetX(resVuv [0]);
1069   Vuuu.SetX(resVuuu[0]);
1070   Vvvv.SetX(resVvvv[0]);
1071   Vuuv.SetX(resVuuv[0]);
1072   Vuvv.SetX(resVuvv[0]);
1073   
1074   P   .SetY(result [1]);
1075   Vu  .SetY(resVu  [1]);
1076   Vv  .SetY(resVv  [1]);
1077   Vuu .SetY(resVuu [1]);
1078   Vvv .SetY(resVvv [1]);
1079   Vuv .SetY(resVuv [1]);
1080   Vuuu.SetY(resVuuu[1]);
1081   Vvvv.SetY(resVvvv[1]);
1082   Vuuv.SetY(resVuuv[1]);
1083   Vuvv.SetY(resVuvv[1]);
1084   
1085   P   .SetZ(result [2]);
1086   Vu  .SetZ(resVu  [2]);
1087   Vv  .SetZ(resVv  [2]);
1088   Vuu .SetZ(resVuu [2]);
1089   Vvv .SetZ(resVvv [2]);
1090   Vuv .SetZ(resVuv [2]);
1091   Vuuu.SetZ(resVuuu[2]);
1092   Vvvv.SetZ(resVvvv[2]);
1093   Vuuv.SetZ(resVuuv[2]);
1094   Vuvv.SetZ(resVuvv[2]);
1095 }
1096
1097 //=======================================================================
1098 //function : DN
1099 //purpose  : 
1100 //=======================================================================
1101
1102 void  BSplSLib::DN
1103 (const Standard_Real U,
1104  const Standard_Real V,
1105  const Standard_Integer Nu,
1106  const Standard_Integer Nv,
1107  const Standard_Integer UIndex,
1108  const Standard_Integer VIndex,
1109  const TColgp_Array2OfPnt& Poles,
1110  const TColStd_Array2OfReal& Weights,
1111  const TColStd_Array1OfReal& UKnots,
1112  const TColStd_Array1OfReal& VKnots,
1113  const TColStd_Array1OfInteger& UMults,
1114  const TColStd_Array1OfInteger& VMults,
1115  const Standard_Integer UDegree,
1116  const Standard_Integer VDegree,
1117  const Standard_Boolean URat,
1118  const Standard_Boolean VRat,
1119  const Standard_Boolean UPer,
1120  const Standard_Boolean VPer,
1121  gp_Vec& Vn)
1122 {
1123   Standard_Boolean rational;
1124   Standard_Integer k,dim;
1125   Standard_Real u1,u2;
1126   Standard_Integer d1,d2;
1127   
1128   BSplSLib_DataContainer dc (UDegree, VDegree);
1129   Standard_Boolean ufirst = PrepareEval
1130     (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
1131      Poles,Weights,UKnots,VKnots,UMults,VMults,
1132      u1,u2,d1,d2,rational,dc);
1133   dim  = rational ? 4 : 3;
1134   
1135   if (!rational) {
1136     if ((Nu > UDegree) || (Nv > VDegree)) {
1137       Vn.SetX(0.);
1138       Vn.SetY(0.);
1139       Vn.SetZ(0.);
1140       return;
1141     }
1142   }
1143   
1144   Standard_Integer n1 = ufirst ? Nu : Nv;
1145   Standard_Integer n2 = ufirst ? Nv : Nu;
1146   
1147   BSplCLib::Bohm(u1,d1,n1,*dc.knots1,dim * (d2 + 1),*dc.poles);
1148
1149   for (k = 0; k <= Min(n1,d1); k++) 
1150     BSplCLib::Bohm(u2,d2,n2,*dc.knots2,dim,*(dc.poles+k*dim*(d2+1)));
1151   
1152   Standard_Real *result;
1153   if (rational) {
1154     BSplSLib::RationalDerivative(d1,d2,n1,n2,*dc.poles,*dc.ders,Standard_False);
1155     result = dc.ders; // because Standard_False ci-dessus.
1156     
1157   }
1158   else {
1159     result = dc.poles + (n1 * (d2+1) + n2) * dim ;
1160   }
1161   
1162   Vn.SetX(result[0]);
1163   Vn.SetY(result[1]);
1164   Vn.SetZ(result[2]);
1165 }
1166
1167 //
1168 // Surface modifications
1169 // 
1170 // a surface is processed as a curve of curves.
1171 // i.e : a curve of parameter u where the current point is the set of poles 
1172 //       of the iso.
1173 //
1174
1175 //=======================================================================
1176 //function : Iso
1177 //purpose  : 
1178 //=======================================================================
1179
1180 void  BSplSLib::Iso(const Standard_Real            Param, 
1181                     const Standard_Boolean         IsU,
1182                     const TColgp_Array2OfPnt&      Poles,
1183                     const TColStd_Array2OfReal&    Weights,
1184                     const TColStd_Array1OfReal&    Knots,
1185                     const TColStd_Array1OfInteger& Mults,
1186                     const Standard_Integer         Degree,
1187                     const Standard_Boolean         Periodic,
1188                     TColgp_Array1OfPnt&            CPoles,
1189                     TColStd_Array1OfReal&          CWeights)
1190 {
1191   Standard_Integer index = 0;
1192   Standard_Real    u = Param;
1193   Standard_Boolean rational = &Weights != NULL;
1194   Standard_Integer dim = rational ? 4 : 3;
1195   
1196   // compute local knots
1197   
1198   BSplSLib_LocalArray locknots1 (2*Degree);  
1199   BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
1200   BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*locknots1);
1201   if (&Mults == NULL)
1202     index -= Knots.Lower() + Degree;
1203   else
1204     index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults);
1205   
1206   
1207   // copy the local poles
1208   
1209 //  Standard_Integer f1,l1,f2,l2,i,j,k;
1210   Standard_Integer f1,l1,f2,l2,i,j;
1211   
1212   if (IsU) {
1213     f1 = Poles.LowerRow();
1214     l1 = Poles.UpperRow();
1215     f2 = Poles.LowerCol();
1216     l2 = Poles.UpperCol();
1217   }
1218   else {
1219     f1 = Poles.LowerCol();
1220     l1 = Poles.UpperCol();
1221     f2 = Poles.LowerRow();
1222     l2 = Poles.UpperRow();
1223   }
1224   
1225   BSplSLib_LocalArray locpoles ((Degree+1) * (l2-f2+1) * dim);
1226   
1227   Standard_Real w, *pole = locpoles;
1228   index += f1;
1229
1230   for (i = 0; i <= Degree; i++) {
1231
1232     for (j = f2; j <= l2; j++) {
1233       
1234       const gp_Pnt& P  = IsU ? Poles(index,j)   : Poles(j,index);
1235       if (rational) { 
1236         pole[3] = w      = IsU ? Weights(index,j) : Weights(j,index);
1237         pole[0] = P.X() * w;
1238         pole[1] = P.Y() * w;
1239         pole[2] = P.Z() * w;
1240       }
1241       else {
1242         pole[0] = P.X();
1243         pole[1] = P.Y();
1244         pole[2] = P.Z();
1245       }
1246       pole += dim;
1247     }
1248     index++;
1249     if (index > l1) index = f1;
1250   }
1251   
1252   // compute the iso
1253   BSplCLib::Eval(u,Degree,*locknots1,(l2-f2+1)*dim,*locpoles);
1254   
1255   // get the result
1256   pole = locpoles;
1257
1258   for (i = CPoles.Lower(); i <= CPoles.Upper(); i++) {
1259     gp_Pnt& P = CPoles(i);
1260     if (rational) {
1261       CWeights(i) = w = pole[3];
1262       P.SetX( pole[0] / w);
1263       P.SetY( pole[1] / w);
1264       P.SetZ( pole[2] / w);
1265     }
1266     else {
1267       P.SetX( pole[0]);
1268       P.SetY( pole[1]);
1269       P.SetZ( pole[2]);
1270     }
1271     pole += dim;
1272   }
1273   
1274   // if the input is not rational but weights are wanted
1275   if (!rational && (&CWeights != NULL)) {
1276
1277     for (i = CWeights.Lower(); i <= CWeights.Upper(); i++)
1278       CWeights(i) = 1.;
1279   }
1280 }
1281
1282 //=======================================================================
1283 //function : Reverse
1284 //purpose  : 
1285 //=======================================================================
1286
1287 void  BSplSLib::Reverse(      TColgp_Array2OfPnt& Poles,
1288                         const Standard_Integer    Last, 
1289                         const Standard_Boolean    UDirection)
1290 {
1291   Standard_Integer i,j, l = Last;
1292   if ( UDirection) {
1293     l = Poles.LowerRow() + 
1294       (l - Poles.LowerRow())%(Poles.ColLength());
1295     TColgp_Array2OfPnt temp(0, Poles.ColLength()-1,
1296                             Poles.LowerCol(), Poles.UpperCol());
1297     
1298     for (i = Poles.LowerRow(); i <= l; i++) {
1299
1300       for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1301         temp(l-i,j) = Poles(i,j);
1302       }
1303     }
1304
1305     for (i = l+1; i <= Poles.UpperRow(); i++) {
1306
1307       for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1308         temp(l+Poles.ColLength()-i,j) = Poles(i,j);
1309       }
1310     }
1311
1312     for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1313
1314       for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1315         Poles(i,j) = temp (i-Poles.LowerRow(),j);
1316       }
1317     }
1318   }
1319   else {
1320     l = Poles.LowerCol() + 
1321       (l - Poles.LowerCol())%(Poles.RowLength());
1322     TColgp_Array2OfPnt temp(Poles.LowerRow(), Poles.UpperRow(),
1323                             0, Poles.RowLength()-1);
1324
1325     for (j = Poles.LowerCol(); j <= l; j++) {
1326
1327       for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1328         temp(i,l-j) = Poles(i,j);
1329       }
1330     }
1331
1332     for (j = l+1; j <= Poles.UpperCol(); j++) {
1333
1334       for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1335         temp(i,l+Poles.RowLength()-j) = Poles(i,j);
1336       }
1337     }
1338
1339     for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1340
1341       for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1342         Poles(i,j) = temp (i,j-Poles.LowerCol());
1343       }
1344     }
1345   }
1346 }
1347
1348 //=======================================================================
1349 //function : Reverse
1350 //purpose  : 
1351 //=======================================================================
1352
1353 void  BSplSLib::Reverse(      TColStd_Array2OfReal& Weights, 
1354                         const Standard_Integer      Last, 
1355                         const Standard_Boolean      UDirection)
1356 {
1357   Standard_Integer i,j, l = Last;
1358   if ( UDirection) {
1359     l = Weights.LowerRow() + 
1360       (l - Weights.LowerRow())%(Weights.ColLength());
1361     TColStd_Array2OfReal temp(0, Weights.ColLength()-1,
1362                               Weights.LowerCol(), Weights.UpperCol());
1363
1364     for (i = Weights.LowerRow(); i <= l; i++) {
1365
1366       for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1367         temp(l-i,j) = Weights(i,j);
1368       }
1369     }
1370
1371     for (i = l+1; i <= Weights.UpperRow(); i++) {
1372
1373       for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1374         temp(l+Weights.ColLength()-i,j) = Weights(i,j);
1375       }
1376     }
1377
1378     for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1379
1380       for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1381         Weights(i,j) = temp (i-Weights.LowerRow(),j);
1382       }
1383     }
1384   }
1385   else {
1386     l = Weights.LowerCol() + 
1387       (l - Weights.LowerCol())%(Weights.RowLength());
1388     TColStd_Array2OfReal temp(Weights.LowerRow(), Weights.UpperRow(),
1389                               0, Weights.RowLength()-1);
1390
1391     for (j = Weights.LowerCol(); j <= l; j++) {
1392
1393       for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1394         temp(i,l-j) = Weights(i,j);
1395       }
1396     }
1397
1398     for (j = l+1; j <= Weights.UpperCol(); j++) {
1399
1400       for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1401         temp(i,l+Weights.RowLength()-j) = Weights(i,j);
1402       }
1403     }
1404
1405     for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1406
1407       for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1408         Weights(i,j) = temp (i,j-Weights.LowerCol());
1409       }
1410     }
1411   }
1412 }
1413
1414 //=======================================================================
1415 //function : IsRational
1416 //purpose  : 
1417 //=======================================================================
1418
1419 Standard_Boolean BSplSLib::IsRational
1420 (const TColStd_Array2OfReal& Weights, 
1421  const Standard_Integer      I1, 
1422  const Standard_Integer      I2, 
1423  const Standard_Integer      J1, 
1424  const Standard_Integer      J2, 
1425  const Standard_Real         Epsi)
1426 {
1427   Standard_Real eps = (Epsi > 0.0) ? Epsi : Epsilon(Weights(I1,I2));
1428   Standard_Integer i,j;
1429   Standard_Integer fi = Weights.LowerRow(), li = Weights.ColLength();
1430   Standard_Integer fj = Weights.LowerCol(), lj = Weights.RowLength();
1431
1432   for (i = I1 - fi; i < I2 - fi; i++) {
1433
1434     for (j = J1 - fj; j < J2 - fj; j++) {
1435       if (Abs(Weights(fi+i%li,fj+j%lj)-Weights(fi+(i+1)%li,fj+j%lj))>eps)
1436         return Standard_True;
1437     }
1438   }
1439   return Standard_False;
1440 }
1441
1442 //=======================================================================
1443 //function : SetPoles
1444 //purpose  : 
1445 //=======================================================================
1446
1447 void  BSplSLib::SetPoles(const TColgp_Array2OfPnt&   Poles, 
1448                          TColStd_Array1OfReal& FP,
1449                          const Standard_Boolean      UDirection)
1450 {
1451   Standard_Integer i, j, l = FP.Lower();
1452   Standard_Integer PLowerRow = Poles.LowerRow();
1453   Standard_Integer PUpperRow = Poles.UpperRow();
1454   Standard_Integer PLowerCol = Poles.LowerCol();
1455   Standard_Integer PUpperCol = Poles.UpperCol();
1456   if (UDirection) {
1457
1458     for ( i = PLowerRow; i <= PUpperRow; i++) {
1459
1460       for ( j = PLowerCol; j <= PUpperCol; j++) {
1461         const gp_Pnt& P = Poles.Value(i,j);
1462         FP(l) = P.X(); l++;
1463         FP(l) = P.Y(); l++;
1464         FP(l) = P.Z(); l++;
1465       }
1466     }
1467   }
1468   else {
1469
1470     for ( j = PLowerCol; j <= PUpperCol; j++) {
1471
1472       for ( i = PLowerRow; i <= PUpperRow; i++) {
1473         const gp_Pnt& P = Poles.Value(i,j);
1474         FP(l) = P.X(); l++;
1475         FP(l) = P.Y(); l++;
1476         FP(l) = P.Z(); l++;
1477       }
1478     }
1479   }
1480 }
1481
1482 //=======================================================================
1483 //function : SetPoles
1484 //purpose  : 
1485 //=======================================================================
1486
1487 void  BSplSLib::SetPoles(const TColgp_Array2OfPnt&   Poles, 
1488                          const TColStd_Array2OfReal& Weights, 
1489                          TColStd_Array1OfReal& FP,
1490                          const Standard_Boolean      UDirection)
1491 {
1492   Standard_Integer i, j, l = FP.Lower();
1493   Standard_Integer PLowerRow = Poles.LowerRow();
1494   Standard_Integer PUpperRow = Poles.UpperRow();
1495   Standard_Integer PLowerCol = Poles.LowerCol();
1496   Standard_Integer PUpperCol = Poles.UpperCol();
1497   if (UDirection) {
1498
1499     for ( i = PLowerRow; i <= PUpperRow; i++) {
1500
1501       for ( j = PLowerCol; j <= PUpperCol; j++) {
1502         const gp_Pnt& P = Poles  .Value(i,j);
1503         Standard_Real w = Weights.Value(i,j);
1504         FP(l) = P.X() * w; l++;
1505         FP(l) = P.Y() * w; l++;
1506         FP(l) = P.Z() * w; l++;
1507         FP(l) = w; l++;
1508       }
1509     }
1510   }
1511   else {
1512
1513     for ( j = PLowerCol; j <= PUpperCol; j++) {
1514
1515       for ( i = PLowerRow; i <= PUpperRow; i++) {
1516         const gp_Pnt& P = Poles  .Value(i,j);
1517         Standard_Real w = Weights.Value(i,j);
1518         FP(l) = P.X() * w; l++;
1519         FP(l) = P.Y() * w; l++;
1520         FP(l) = P.Z() * w; l++;
1521         FP(l) = w; l++;
1522       }
1523     }
1524   }
1525 }
1526
1527 //=======================================================================
1528 //function : GetPoles
1529 //purpose  : 
1530 //=======================================================================
1531
1532 void  BSplSLib::GetPoles(const TColStd_Array1OfReal& FP, 
1533                          TColgp_Array2OfPnt&   Poles,
1534                          const Standard_Boolean      UDirection)
1535 {
1536   Standard_Integer i, j, l = FP.Lower();
1537   Standard_Integer PLowerRow = Poles.LowerRow();
1538   Standard_Integer PUpperRow = Poles.UpperRow();
1539   Standard_Integer PLowerCol = Poles.LowerCol();
1540   Standard_Integer PUpperCol = Poles.UpperCol();
1541   if (UDirection) {
1542
1543     for ( i = PLowerRow; i <= PUpperRow; i++) {
1544
1545       for ( j = PLowerCol; j <= PUpperCol; j++) {
1546         gp_Pnt& P = Poles.ChangeValue(i,j);
1547         P.SetX(FP(l)); l++;
1548         P.SetY(FP(l)); l++;
1549         P.SetZ(FP(l)); l++;
1550       }
1551     }
1552   }
1553   else {
1554
1555     for ( j = PLowerCol; j <= PUpperCol; j++) {
1556
1557       for ( i = PLowerRow; i <= PUpperRow; i++) {
1558         gp_Pnt& P = Poles.ChangeValue(i,j);
1559         P.SetX(FP(l)); l++;
1560         P.SetY(FP(l)); l++;
1561         P.SetZ(FP(l)); l++;
1562       }
1563     }
1564   }
1565 }
1566
1567 //=======================================================================
1568 //function : GetPoles
1569 //purpose  : 
1570 //=======================================================================
1571
1572 void  BSplSLib::GetPoles(const TColStd_Array1OfReal& FP, 
1573                          TColgp_Array2OfPnt&   Poles, 
1574                          TColStd_Array2OfReal& Weights,
1575                          const Standard_Boolean      UDirection)
1576 {
1577   Standard_Integer i, j, l = FP.Lower();
1578   Standard_Integer PLowerRow = Poles.LowerRow();
1579   Standard_Integer PUpperRow = Poles.UpperRow();
1580   Standard_Integer PLowerCol = Poles.LowerCol();
1581   Standard_Integer PUpperCol = Poles.UpperCol();
1582   if (UDirection) {
1583
1584     for ( i = PLowerRow; i <= PUpperRow; i++) {
1585
1586       for ( j = PLowerCol; j <= PUpperCol; j++) {
1587         Standard_Real w = FP( l + 3);
1588         Weights(i,j) = w;
1589         gp_Pnt& P = Poles.ChangeValue(i,j);
1590         P.SetX(FP(l) / w); l++;
1591         P.SetY(FP(l) / w); l++;
1592         P.SetZ(FP(l) / w); l++;
1593         l++;
1594       }
1595     }
1596   }
1597   else {
1598
1599     for ( j = PLowerCol; j <= PUpperCol; j++) {
1600
1601       for ( i = PLowerRow; i <= PUpperRow; i++) {
1602         Standard_Real w = FP( l + 3);
1603         Weights(i,j) = w;
1604         gp_Pnt& P = Poles.ChangeValue(i,j);
1605         P.SetX(FP(l) / w); l++;
1606         P.SetY(FP(l) / w); l++;
1607         P.SetZ(FP(l) / w); l++;
1608         l++;
1609       }
1610     }
1611   }
1612 }
1613
1614 //=======================================================================
1615 //function : InsertKnots
1616 //purpose  : 
1617 //=======================================================================
1618
1619 void  BSplSLib::InsertKnots(const Standard_Boolean         UDirection, 
1620                             const Standard_Integer         Degree, 
1621                             const Standard_Boolean         Periodic, 
1622                             const TColgp_Array2OfPnt&      Poles, 
1623                             const TColStd_Array2OfReal&    Weights, 
1624                             const TColStd_Array1OfReal&    Knots, 
1625                             const TColStd_Array1OfInteger& Mults, 
1626                             const TColStd_Array1OfReal&    AddKnots, 
1627                             const TColStd_Array1OfInteger& AddMults, 
1628                             TColgp_Array2OfPnt&      NewPoles,
1629                             TColStd_Array2OfReal&    NewWeights,
1630                             TColStd_Array1OfReal&    NewKnots,
1631                             TColStd_Array1OfInteger& NewMults, 
1632                             const Standard_Real            Epsilon, 
1633                             const Standard_Boolean         Add )
1634 {
1635   Standard_Boolean rational = &Weights != NULL;
1636   Standard_Integer dim = 3;
1637   if (rational) dim++;
1638   
1639   TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1640   TColStd_Array1OfReal 
1641     newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1642   
1643   if (rational) SetPoles(Poles,Weights,poles,UDirection);
1644   else          SetPoles(Poles,poles,UDirection);
1645   
1646   if (UDirection) {
1647     dim *= Poles.RowLength();
1648   }
1649   else {
1650     dim *= Poles.ColLength();
1651   }
1652   BSplCLib::InsertKnots(Degree,Periodic,dim,poles,Knots,Mults,
1653                         AddKnots,AddMults,newpoles,NewKnots,NewMults,
1654                         Epsilon,Add);
1655   
1656   if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1657   else          GetPoles(newpoles,NewPoles,UDirection);
1658 }
1659
1660 //=======================================================================
1661 //function : RemoveKnot
1662 //purpose  : 
1663 //=======================================================================
1664
1665 Standard_Boolean  BSplSLib::RemoveKnot
1666 (const Standard_Boolean         UDirection, 
1667  const Standard_Integer         Index,
1668  const Standard_Integer         Mult,
1669  const Standard_Integer         Degree, 
1670  const Standard_Boolean         Periodic, 
1671  const TColgp_Array2OfPnt&      Poles, 
1672  const TColStd_Array2OfReal&    Weights, 
1673  const TColStd_Array1OfReal&    Knots, 
1674  const TColStd_Array1OfInteger& Mults,
1675  TColgp_Array2OfPnt&      NewPoles,
1676  TColStd_Array2OfReal&    NewWeights,
1677  TColStd_Array1OfReal&    NewKnots,
1678  TColStd_Array1OfInteger& NewMults,
1679  const Standard_Real            Tolerance)
1680 {
1681   Standard_Boolean rational = &Weights != NULL;
1682   Standard_Integer dim = 3;
1683   if (rational) dim++;
1684   
1685   TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1686   TColStd_Array1OfReal 
1687     newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1688   
1689   if (rational) SetPoles(Poles,Weights,poles,UDirection);
1690   else          SetPoles(Poles,poles,UDirection);
1691   
1692   if (UDirection) {
1693     dim *= Poles.RowLength();
1694   }
1695   else {
1696     dim *= Poles.ColLength();
1697   }
1698   
1699   if ( !BSplCLib::RemoveKnot(Index,Mult,Degree,Periodic,dim,
1700                              poles,Knots,Mults,newpoles,NewKnots,NewMults,
1701                              Tolerance))
1702     return Standard_False;
1703   
1704   if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1705   else          GetPoles(newpoles,NewPoles,UDirection);
1706   return Standard_True;
1707 }
1708
1709 //=======================================================================
1710 //function : IncreaseDegree
1711 //purpose  : 
1712 //=======================================================================
1713
1714 void  BSplSLib::IncreaseDegree
1715 (const Standard_Boolean         UDirection,
1716  const Standard_Integer         Degree, 
1717  const Standard_Integer         NewDegree, 
1718  const Standard_Boolean         Periodic, 
1719  const TColgp_Array2OfPnt&      Poles, 
1720  const TColStd_Array2OfReal&    Weights,
1721  const TColStd_Array1OfReal&    Knots,
1722  const TColStd_Array1OfInteger& Mults, 
1723  TColgp_Array2OfPnt&      NewPoles, 
1724  TColStd_Array2OfReal&    NewWeights, 
1725  TColStd_Array1OfReal&    NewKnots, 
1726  TColStd_Array1OfInteger& NewMults)
1727 {
1728   Standard_Boolean rational = &Weights != NULL;
1729   Standard_Integer dim = 3;
1730   if (rational) dim++;
1731   
1732   TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1733   TColStd_Array1OfReal 
1734     newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1735   
1736   if (rational) SetPoles(Poles,Weights,poles,UDirection);
1737   else          SetPoles(Poles,poles,UDirection);
1738   
1739   if (UDirection) {
1740     dim *= Poles.RowLength();
1741   }
1742   else {
1743     dim *= Poles.ColLength();
1744   }
1745   
1746   BSplCLib::IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults,
1747                            newpoles,NewKnots,NewMults);
1748   
1749   if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1750   else          GetPoles(newpoles,NewPoles,UDirection);
1751 }
1752
1753 //=======================================================================
1754 //function : Unperiodize
1755 //purpose  : 
1756 //=======================================================================
1757
1758 void  BSplSLib::Unperiodize
1759 (const Standard_Boolean         UDirection,
1760  const Standard_Integer         Degree,
1761  const TColStd_Array1OfInteger& Mults,
1762  const TColStd_Array1OfReal&    Knots,
1763  const TColgp_Array2OfPnt&      Poles, 
1764  const TColStd_Array2OfReal&    Weights,
1765  TColStd_Array1OfInteger& NewMults,
1766  TColStd_Array1OfReal&    NewKnots,
1767  TColgp_Array2OfPnt&      NewPoles,
1768  TColStd_Array2OfReal&    NewWeights)
1769 {
1770   Standard_Boolean rational = &Weights != NULL;
1771   Standard_Integer dim = 3;
1772   if (rational) dim++;
1773   
1774   TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1775   TColStd_Array1OfReal 
1776     newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1777   
1778   if (rational) SetPoles(Poles,Weights,poles,UDirection);
1779   else          SetPoles(Poles,poles,UDirection);
1780   
1781   if (UDirection) {
1782     dim *= Poles.RowLength();
1783   }
1784   else {
1785     dim *= Poles.ColLength();
1786   }
1787   
1788   BSplCLib::Unperiodize(Degree, dim, Mults, Knots, poles, 
1789                         NewMults, NewKnots, newpoles);
1790   
1791   if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1792   else          GetPoles(newpoles,NewPoles,UDirection);
1793 }
1794
1795 //=======================================================================
1796 //function : BuildCache
1797 //purpose  : Stores theTaylor expansion normalized between 0,1 in the
1798 //           Cache : in case of  a rational function the Poles are
1799 //           stored in homogeneous form 
1800 //=======================================================================
1801
1802 void BSplSLib::BuildCache
1803 (const Standard_Real            U,   
1804  const Standard_Real            V,
1805  const Standard_Real            USpanDomain,
1806  const Standard_Real            VSpanDomain,  
1807  const Standard_Boolean         UPeriodic,
1808  const Standard_Boolean         VPeriodic,
1809  const Standard_Integer         UDegree,
1810  const Standard_Integer         VDegree,
1811  const Standard_Integer         UIndex,
1812  const Standard_Integer         VIndex,
1813  const TColStd_Array1OfReal&    UFlatKnots,   
1814  const TColStd_Array1OfReal&    VFlatKnots,   
1815  const TColgp_Array2OfPnt&      Poles,  
1816  const TColStd_Array2OfReal&    Weights,
1817  TColgp_Array2OfPnt&            CachePoles,
1818  TColStd_Array2OfReal&          CacheWeights)
1819 {  
1820   Standard_Boolean rational,rational_u,rational_v,flag_u_or_v;                  
1821   Standard_Integer kk,d1,d1p1,d2,d2p1,ii,jj,iii,jjj,Index;
1822   Standard_Real u1,min_degree_domain,max_degree_domain,f,factor[2],u2;
1823   if (&Weights != NULL) 
1824     rational_u = rational_v = Standard_True;
1825   else
1826     rational_u = rational_v = Standard_False;
1827   BSplSLib_DataContainer dc (UDegree, VDegree);
1828   flag_u_or_v =
1829     PrepareEval  (U,
1830                   V,
1831                   UIndex,
1832                   VIndex,
1833                   UDegree,
1834                   VDegree,
1835                   rational_u,
1836                   rational_v,
1837                   UPeriodic,
1838                   VPeriodic,
1839                   Poles,
1840                   Weights,
1841                   UFlatKnots,
1842                   VFlatKnots,
1843                   (BSplCLib::NoMults()),
1844                   (BSplCLib::NoMults()),
1845                   u1,
1846                   u2,
1847                   d1,
1848                   d2,
1849                   rational,
1850           dc);
1851   d1p1 = d1 + 1;
1852   d2p1 = d2 + 1;
1853   if (rational) {
1854     BSplCLib::Bohm(u1,d1,d1,*dc.knots1,4 * d2p1,*dc.poles);
1855     
1856     for (kk = 0; kk <= d1 ; kk++) 
1857       BSplCLib::Bohm(u2,d2,d2,*dc.knots2,4,*(dc.poles + kk * 4 * d2p1));
1858     if (flag_u_or_v) {
1859       min_degree_domain = USpanDomain ;
1860       max_degree_domain = VSpanDomain ;
1861     }
1862     else {
1863       min_degree_domain = VSpanDomain ;
1864       max_degree_domain = USpanDomain ;
1865     }
1866     factor[0] = 1.0e0 ;
1867     
1868     for (ii = 0 ; ii <= d2 ; ii++) {
1869       iii = ii + 1;
1870       factor[1] = 1.0e0 ;
1871       
1872       for (jj = 0 ; jj <= d1 ; jj++) {
1873         jjj = jj + 1;
1874         Index = jj * d2p1 + ii ;
1875         Index = Index << 2;
1876         gp_Pnt& P = CachePoles(iii,jjj);
1877         f = factor[0] * factor[1];
1878         P.SetX( f * dc.poles[Index]); Index++;
1879         P.SetY( f * dc.poles[Index]); Index++;
1880         P.SetZ( f * dc.poles[Index]); Index++;
1881         CacheWeights(iii ,jjj) = f * dc.poles[Index] ;
1882         factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
1883       }
1884       factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
1885     }
1886   }
1887   else {
1888     BSplCLib::Bohm(u1,d1,d1,*dc.knots1,3 * d2p1,*dc.poles);
1889     
1890     for (kk = 0; kk <= d1 ; kk++) 
1891       BSplCLib::Bohm(u2,d2,d2,*dc.knots2,3,*(dc.poles + kk * 3 * d2p1));
1892     if (flag_u_or_v) {
1893       min_degree_domain = USpanDomain ;
1894       max_degree_domain = VSpanDomain ;
1895     }
1896     else {
1897       min_degree_domain = VSpanDomain ;
1898       max_degree_domain = USpanDomain ;
1899     }
1900     factor[0] = 1.0e0 ;
1901     
1902     for (ii = 0 ; ii <= d2 ; ii++) {
1903       iii = ii + 1;
1904       factor[1] = 1.0e0 ;
1905       
1906       for (jj = 0 ; jj <= d1 ; jj++) {
1907         jjj = jj + 1;
1908         Index = jj * d2p1 + ii ;
1909         Index = (Index << 1) + Index;
1910         gp_Pnt& P = CachePoles(iii,jjj);
1911         f = factor[0] * factor[1];
1912         P.SetX( f * dc.poles[Index]); Index++;
1913         P.SetY( f * dc.poles[Index]); Index++;
1914         P.SetZ( f * dc.poles[Index]);
1915         factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
1916       }
1917       factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
1918     }
1919     if (&Weights != NULL) {
1920       //
1921       // means that PrepareEval did found out that the surface was 
1922       // locally polynomial but since the surface is constructed
1923       // with some weights we need to set the weight polynomial to constant
1924       // 
1925       
1926       for (ii = 1 ; ii <= d2p1 ; ii++) {
1927         
1928         for (jj = 1 ; jj <= d1p1 ; jj++) {
1929           CacheWeights(ii,jj) = 0.0e0 ;
1930         }
1931       }
1932       CacheWeights(1,1) = 1.0e0 ;
1933     }
1934   }
1935 }
1936
1937 //=======================================================================
1938 //function : CacheD0
1939 //purpose  : Evaluates the polynomial cache of the Bspline Curve
1940 //           
1941 //=======================================================================
1942
1943 void  BSplSLib::CacheD0(const Standard_Real                  UParameter,
1944                         const Standard_Real                  VParameter,
1945                         const  Standard_Integer              UDegree,
1946                         const  Standard_Integer              VDegree,
1947                         const  Standard_Real                 UCacheParameter,
1948                         const  Standard_Real                 VCacheParameter,
1949                         const  Standard_Real                 USpanLenght,
1950                         const  Standard_Real                 VSpanLenght,
1951                         const  TColgp_Array2OfPnt&           PolesArray,
1952                         const  TColStd_Array2OfReal&         WeightsArray,
1953                         gp_Pnt&                              aPoint)
1954 {
1955   //
1956   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
1957   // form
1958   // the SpanLenght     is the normalizing factor so that the polynomial is between
1959   // 0 and 1 
1960   Standard_Integer 
1961 //    ii,
1962   dimension,
1963   min_degree,
1964   max_degree  ;
1965   
1966   Standard_Real 
1967     new_parameter[2] ,
1968   inverse ;
1969   
1970   Standard_Real * 
1971     PArray = (Standard_Real *) 
1972       &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
1973   Standard_Real *
1974     myPoint = (Standard_Real *) &aPoint  ;
1975   if (UDegree <= VDegree) {
1976     min_degree = UDegree ;
1977     max_degree = VDegree ;
1978     new_parameter[1] = (UParameter - UCacheParameter) / USpanLenght ;
1979     new_parameter[0] = (VParameter - VCacheParameter) / VSpanLenght ; 
1980     dimension = 3 * (UDegree + 1) ;
1981   }
1982   else {
1983     min_degree = VDegree ;
1984     max_degree = UDegree ;
1985     new_parameter[0] = (UParameter - UCacheParameter) / USpanLenght ;
1986     new_parameter[1] = (VParameter - VCacheParameter) / VSpanLenght ; 
1987     dimension = 3 * (VDegree + 1) ;
1988   }
1989   BSplSLib_LocalArray locpoles(dimension);
1990   
1991   PLib::NoDerivativeEvalPolynomial(new_parameter[0],
1992                        max_degree,
1993                        dimension,
1994                        max_degree*dimension,
1995                        PArray[0],
1996                        locpoles[0]) ;
1997   
1998   PLib::NoDerivativeEvalPolynomial(new_parameter[1],
1999                        min_degree,
2000                        3,
2001                        (min_degree << 1) + min_degree,
2002                        locpoles[0],
2003                        myPoint[0]) ;
2004   if (&WeightsArray != NULL) {
2005     dimension = min_degree + 1 ;
2006     Standard_Real *
2007       WArray = (Standard_Real *) 
2008         &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2009     PLib::NoDerivativeEvalPolynomial(new_parameter[0],
2010                          max_degree,
2011                          dimension,
2012                          max_degree*dimension,
2013                          WArray[0],
2014                          locpoles[0]) ;
2015     
2016     PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2017                          min_degree,
2018                          1,
2019                          min_degree,
2020                          locpoles[0],
2021                          inverse) ;
2022     inverse = 1.0e0/ inverse ;
2023     
2024     myPoint[0] *= inverse ;
2025     myPoint[1] *= inverse ;
2026     myPoint[2] *= inverse ;
2027   }
2028 }
2029
2030 //=======================================================================
2031 //function : CacheD1
2032 //purpose  : Evaluates the polynomial cache of the Bspline Curve
2033 //           
2034 //=======================================================================
2035
2036 void  BSplSLib::CacheD1(const Standard_Real                  UParameter,
2037                         const Standard_Real                  VParameter,
2038                         const  Standard_Integer              UDegree,
2039                         const  Standard_Integer              VDegree,
2040                         const  Standard_Real                 UCacheParameter,
2041                         const  Standard_Real                 VCacheParameter,
2042                         const  Standard_Real                 USpanLenght,
2043                         const  Standard_Real                 VSpanLenght,
2044                         const  TColgp_Array2OfPnt&           PolesArray,
2045                         const  TColStd_Array2OfReal&         WeightsArray,
2046                         gp_Pnt&                              aPoint,
2047                         gp_Vec&                              aVecU,
2048                         gp_Vec&                              aVecV)
2049 {
2050   //
2051   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2052   // form
2053   // the SpanLenght     is the normalizing factor so that the polynomial is between
2054   // 0 and 1 
2055   Standard_Integer 
2056 //    ii,
2057 //  jj,
2058 //  kk,
2059   dimension,
2060   min_degree,
2061   max_degree  ;
2062   
2063   Standard_Real
2064     inverse_min,
2065   inverse_max, 
2066   new_parameter[2]  ;
2067   
2068   Standard_Real * 
2069     PArray = (Standard_Real *) 
2070       &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2071   Standard_Real local_poles_array[2][2][3],
2072   local_poles_and_weights_array[2][2][4],
2073   local_weights_array[2][2]  ;
2074   Standard_Real * my_vec_min,
2075   * my_vec_max,
2076   * my_point ;
2077   my_point = (Standard_Real *) &aPoint  ;
2078   //
2079   // initialize in case of rational evaluation
2080   // because RationalDerivative will use all
2081   // the coefficients
2082   //
2083   //
2084   if (&WeightsArray != NULL) {
2085
2086     local_poles_array            [0][0][0] = 0.0e0 ;
2087     local_poles_array            [0][0][1] = 0.0e0 ;
2088     local_poles_array            [0][0][2] = 0.0e0 ;
2089     local_weights_array          [0][0]    = 0.0e0 ;
2090     local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2091     local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2092     local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2093     local_poles_and_weights_array[0][0][3] = 0.0e0 ;
2094     
2095     local_poles_array            [0][1][0] = 0.0e0 ;
2096     local_poles_array            [0][1][1] = 0.0e0 ;
2097     local_poles_array            [0][1][2] = 0.0e0 ;
2098     local_weights_array          [0][1]    = 0.0e0 ;
2099     local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2100     local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2101     local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2102     local_poles_and_weights_array[0][1][3] = 0.0e0 ;
2103
2104     local_poles_array            [1][0][0] = 0.0e0 ;
2105     local_poles_array            [1][0][1] = 0.0e0 ;
2106     local_poles_array            [1][0][2] = 0.0e0 ;
2107     local_weights_array          [1][0]    = 0.0e0 ;
2108     local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2109     local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2110     local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2111     local_poles_and_weights_array[1][0][3] = 0.0e0 ;
2112     
2113     local_poles_array            [1][1][0] = 0.0e0 ;
2114     local_poles_array            [1][1][1] = 0.0e0 ;
2115     local_poles_array            [1][1][2] = 0.0e0 ;
2116     local_weights_array          [1][1]    = 0.0e0 ;
2117     local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2118     local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2119     local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2120     local_poles_and_weights_array[1][1][3] = 0.0e0 ;
2121   }
2122
2123   if (UDegree <= VDegree) {
2124     min_degree = UDegree ;
2125     max_degree = VDegree ;
2126     inverse_min = 1.0e0/USpanLenght ;
2127     inverse_max = 1.0e0/VSpanLenght ;
2128     new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ; 
2129     new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2130     
2131     dimension = 3 * (UDegree + 1) ;
2132     my_vec_min = (Standard_Real *) &aVecU ;
2133     my_vec_max = (Standard_Real *) &aVecV ;
2134   }
2135   else {
2136     min_degree = VDegree ;
2137     max_degree = UDegree ;
2138     inverse_min = 1.0e0/VSpanLenght ;
2139     inverse_max = 1.0e0/USpanLenght ;
2140     new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2141     new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ; 
2142     dimension = 3 * (VDegree + 1) ;
2143     my_vec_min = (Standard_Real *) &aVecV ;
2144     my_vec_max = (Standard_Real *) &aVecU ;
2145   }
2146
2147   BSplSLib_LocalArray locpoles (2 * dimension);
2148   
2149   PLib::EvalPolynomial(new_parameter[0],
2150                        1,
2151                        max_degree,
2152                        dimension,
2153                        PArray[0],
2154                        locpoles[0]) ;
2155   
2156   PLib::EvalPolynomial(new_parameter[1],
2157                        1,
2158                        min_degree,
2159                        3,
2160                        locpoles[0],
2161                        local_poles_array[0][0][0]) ;
2162   PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2163                        min_degree,
2164                        3,
2165                        (min_degree << 1) + min_degree,
2166                        locpoles[dimension],
2167                        local_poles_array[1][0][0]) ;
2168   
2169   if (&WeightsArray != NULL) {
2170     dimension = min_degree + 1 ;
2171     Standard_Real *
2172       WArray = (Standard_Real *) 
2173         &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2174     PLib::EvalPolynomial(new_parameter[0],
2175                          1,
2176                          max_degree,
2177                          dimension,
2178                          WArray[0],
2179                          locpoles[0]) ;
2180     
2181     PLib::EvalPolynomial(new_parameter[1],
2182                          1,
2183                          min_degree,
2184                          1,
2185                          locpoles[0],
2186                          local_weights_array[0][0]) ;
2187     PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2188                          min_degree,
2189                          1,
2190                          min_degree,
2191                          locpoles[dimension],
2192                          local_weights_array[1][0]) ;
2193     
2194     local_poles_and_weights_array[0][0][0] = local_poles_array  [0][0][0] ;
2195     local_poles_and_weights_array[0][0][1] = local_poles_array  [0][0][1] ;
2196     local_poles_and_weights_array[0][0][2] = local_poles_array  [0][0][2] ;
2197     local_poles_and_weights_array[0][0][3] = local_weights_array[0][0]    ;
2198     
2199     local_poles_and_weights_array[0][1][0] = local_poles_array  [0][1][0] ;
2200     local_poles_and_weights_array[0][1][1] = local_poles_array  [0][1][1] ;
2201     local_poles_and_weights_array[0][1][2] = local_poles_array  [0][1][2] ;
2202     local_poles_and_weights_array[0][1][3] = local_weights_array[0][1]    ;
2203     
2204     local_poles_and_weights_array[1][0][0] = local_poles_array  [1][0][0] ;
2205     local_poles_and_weights_array[1][0][1] = local_poles_array  [1][0][1] ;
2206     local_poles_and_weights_array[1][0][2] = local_poles_array  [1][0][2] ;
2207     local_poles_and_weights_array[1][0][3] = local_weights_array[1][0]    ;
2208     
2209     local_poles_and_weights_array[1][1][0] = local_poles_array  [1][1][0] ;
2210     local_poles_and_weights_array[1][1][1] = local_poles_array  [1][1][1] ;
2211     local_poles_and_weights_array[1][1][2] = local_poles_array  [1][1][2] ;
2212     local_poles_and_weights_array[1][1][3] = local_weights_array[1][1]    ;
2213
2214     BSplSLib::RationalDerivative(1,
2215                                  1,
2216                                  1,
2217                                  1,
2218                                  local_poles_and_weights_array[0][0][0],
2219                                  local_poles_array[0][0][0]) ;
2220   }
2221   
2222   my_point  [0] = local_poles_array              [0][0][0] ;
2223   my_vec_min[0] = inverse_min * local_poles_array[0][1][0] ;
2224   my_vec_max[0] = inverse_max * local_poles_array[1][0][0] ;
2225   
2226   my_point  [1] = local_poles_array              [0][0][1] ;
2227   my_vec_min[1] = inverse_min * local_poles_array[0][1][1] ;
2228   my_vec_max[1] = inverse_max * local_poles_array[1][0][1] ;
2229   
2230   my_point  [2] = local_poles_array              [0][0][2] ;
2231   my_vec_min[2] = inverse_min * local_poles_array[0][1][2] ;
2232   my_vec_max[2] = inverse_max * local_poles_array[1][0][2] ;
2233 }
2234
2235 //=======================================================================
2236 //function : CacheD2
2237 //purpose  : Evaluates the polynomial cache of the Bspline Curve
2238 //           
2239 //=======================================================================
2240
2241 void  BSplSLib::CacheD2(const Standard_Real                  UParameter,
2242                         const Standard_Real                  VParameter,
2243                         const  Standard_Integer              UDegree,
2244                         const  Standard_Integer              VDegree,
2245                         const  Standard_Real                 UCacheParameter,
2246                         const  Standard_Real                 VCacheParameter,
2247                         const  Standard_Real                 USpanLenght,
2248                         const  Standard_Real                 VSpanLenght,
2249                         const  TColgp_Array2OfPnt&           PolesArray,
2250                         const  TColStd_Array2OfReal&         WeightsArray,
2251                         gp_Pnt&                              aPoint,
2252                         gp_Vec&                              aVecU,
2253                         gp_Vec&                              aVecV,
2254                         gp_Vec&                              aVecUU,
2255                         gp_Vec&                              aVecUV,
2256                         gp_Vec&                              aVecVV)
2257 {
2258   //
2259   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2260   // form
2261   // the SpanLenght     is the normalizing factor so that the polynomial is between
2262   // 0 and 1 
2263   Standard_Integer 
2264     ii,
2265 //  jj,
2266   kk,
2267   index,
2268   dimension,
2269   min_degree,
2270   max_degree  ;
2271   
2272   Standard_Real
2273     inverse_min,
2274   inverse_max, 
2275   new_parameter[2]  ;
2276
2277   Standard_Real * 
2278     PArray = (Standard_Real *) 
2279       &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2280   Standard_Real local_poles_array[3][3][3],
2281   local_poles_and_weights_array[3][3][4],
2282   local_weights_array[3][3]  ;
2283   Standard_Real * my_vec_min,
2284   * my_vec_max,
2285   * my_vec_min_min,
2286   * my_vec_max_max,
2287   * my_vec_min_max,
2288   * my_point ;
2289   my_point = (Standard_Real *) &aPoint  ;
2290   
2291   //
2292   // initialize in case the min and max degree are less than 2
2293   //
2294   local_poles_array[0][0][0] = 0.0e0 ;
2295   local_poles_array[0][0][1] = 0.0e0 ;
2296   local_poles_array[0][0][2] = 0.0e0 ;
2297   local_poles_array[0][1][0] = 0.0e0 ;
2298   local_poles_array[0][1][1] = 0.0e0 ;
2299   local_poles_array[0][1][2] = 0.0e0 ;
2300   local_poles_array[0][2][0] = 0.0e0 ;
2301   local_poles_array[0][2][1] = 0.0e0 ;
2302   local_poles_array[0][2][2] = 0.0e0 ;
2303   
2304   local_poles_array[1][0][0] = 0.0e0 ;
2305   local_poles_array[1][0][1] = 0.0e0 ;
2306   local_poles_array[1][0][2] = 0.0e0 ;
2307   local_poles_array[1][1][0] = 0.0e0 ;
2308   local_poles_array[1][1][1] = 0.0e0 ;
2309   local_poles_array[1][1][2] = 0.0e0 ;
2310   local_poles_array[1][2][0] = 0.0e0 ;
2311   local_poles_array[1][2][1] = 0.0e0 ;
2312   local_poles_array[1][2][2] = 0.0e0 ;
2313   
2314   local_poles_array[2][0][0] = 0.0e0 ;
2315   local_poles_array[2][0][1] = 0.0e0 ;
2316   local_poles_array[2][0][2] = 0.0e0 ;
2317   local_poles_array[2][1][0] = 0.0e0 ;
2318   local_poles_array[2][1][1] = 0.0e0 ;
2319   local_poles_array[2][1][2] = 0.0e0 ;
2320   local_poles_array[2][2][0] = 0.0e0 ;
2321   local_poles_array[2][2][1] = 0.0e0 ;
2322   local_poles_array[2][2][2] = 0.0e0 ;
2323   //
2324   // initialize in case of rational evaluation
2325   // because RationalDerivative will use all
2326   // the coefficients
2327   //
2328   //
2329   if (&WeightsArray != NULL) {
2330     
2331     local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2332     local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2333     local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2334     local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2335     local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2336     local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2337     local_poles_and_weights_array[0][2][0] = 0.0e0 ;
2338     local_poles_and_weights_array[0][2][1] = 0.0e0 ;
2339     local_poles_and_weights_array[0][2][2] = 0.0e0 ;
2340     
2341     local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2342     local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2343     local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2344     local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2345     local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2346     local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2347     local_poles_and_weights_array[1][2][0] = 0.0e0 ;
2348     local_poles_and_weights_array[1][2][1] = 0.0e0 ;
2349     local_poles_and_weights_array[1][2][2] = 0.0e0 ;
2350     
2351     local_poles_and_weights_array[2][0][0] = 0.0e0 ;
2352     local_poles_and_weights_array[2][0][1] = 0.0e0 ;
2353     local_poles_and_weights_array[2][0][2] = 0.0e0 ;
2354     local_poles_and_weights_array[2][1][0] = 0.0e0 ;
2355     local_poles_and_weights_array[2][1][1] = 0.0e0 ;
2356     local_poles_and_weights_array[2][1][2] = 0.0e0 ;
2357     local_poles_and_weights_array[2][2][0] = 0.0e0 ;
2358     local_poles_and_weights_array[2][2][1] = 0.0e0 ;
2359     local_poles_and_weights_array[2][2][2] = 0.0e0 ;
2360     
2361     local_poles_and_weights_array[0][0][3] =
2362       local_weights_array[0][0] = 0.0e0 ;
2363     local_poles_and_weights_array[0][1][3] =
2364       local_weights_array[0][1] = 0.0e0 ;
2365     local_poles_and_weights_array[0][2][3] =
2366       local_weights_array[0][2] = 0.0e0 ;
2367     local_poles_and_weights_array[1][0][3] =
2368       local_weights_array[1][0] = 0.0e0 ;
2369     local_poles_and_weights_array[1][1][3] =
2370       local_weights_array[1][1] = 0.0e0 ;
2371     local_poles_and_weights_array[1][2][3] =
2372       local_weights_array[1][2] = 0.0e0 ;
2373     local_poles_and_weights_array[2][0][3] =
2374       local_weights_array[2][0] = 0.0e0 ;
2375     local_poles_and_weights_array[2][1][3] =
2376       local_weights_array[2][1] = 0.0e0 ;
2377     local_poles_and_weights_array[2][2][3] =
2378       local_weights_array[2][2] = 0.0e0 ;
2379   }
2380
2381   if (UDegree <= VDegree) {
2382     min_degree = UDegree ;
2383     max_degree = VDegree ;
2384     inverse_min = 1.0e0/USpanLenght ;
2385     inverse_max = 1.0e0/VSpanLenght ;
2386     new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ; 
2387     new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2388     
2389     dimension = 3 * (UDegree + 1) ;
2390     my_vec_min     = (Standard_Real *) &aVecU ;
2391     my_vec_max     = (Standard_Real *) &aVecV ;
2392     my_vec_min_min = (Standard_Real *) &aVecUU ;
2393     my_vec_min_max = (Standard_Real *) &aVecUV ;
2394     my_vec_max_max = (Standard_Real *) &aVecVV ;
2395   }
2396   else {
2397     min_degree = VDegree ;
2398     max_degree = UDegree ;
2399     inverse_min = 1.0e0/VSpanLenght ;
2400     inverse_max = 1.0e0/USpanLenght ;
2401     new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2402     new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ; 
2403     dimension = 3 * (VDegree + 1) ;
2404     my_vec_min     = (Standard_Real *) &aVecV ;
2405     my_vec_max     = (Standard_Real *) &aVecU ;
2406     my_vec_min_min = (Standard_Real *) &aVecVV ;
2407     my_vec_min_max = (Standard_Real *) &aVecUV ;
2408     my_vec_max_max = (Standard_Real *) &aVecUU ;
2409   }
2410
2411   BSplSLib_LocalArray locpoles (3 * dimension);
2412   
2413   //
2414   // initialize in case min or max degree are less than 2
2415   //
2416   Standard_Integer MinIndMax = 2;
2417   if ( max_degree < 2) MinIndMax = max_degree;
2418   Standard_Integer MinIndMin = 2;
2419   if ( min_degree < 2) MinIndMin = min_degree;
2420   
2421   index =  MinIndMax * dimension ;
2422
2423   for (ii = MinIndMax ; ii <  3 ; ii++) {
2424     
2425     for (kk = 0 ; kk < dimension ; kk++) {
2426       locpoles[index] = 0.0e0 ;
2427       index += 1 ;
2428     }
2429   }
2430   
2431   PLib::EvalPolynomial(new_parameter[0],
2432                        MinIndMax,
2433                        max_degree,
2434                        dimension,
2435                        PArray[0],
2436                        locpoles[0]) ;
2437   
2438   PLib::EvalPolynomial(new_parameter[1],
2439                        MinIndMin,
2440                        min_degree,
2441                        3,
2442                        locpoles[0],
2443                        local_poles_array[0][0][0]) ;
2444   PLib::EvalPolynomial(new_parameter[1],
2445                        1,
2446                        min_degree,
2447                        3,
2448                        locpoles[dimension],
2449                        local_poles_array[1][0][0]) ;
2450   
2451   PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2452                        min_degree,
2453                        3,
2454                        (min_degree << 1) + min_degree,
2455                        locpoles[dimension + dimension],
2456                        local_poles_array[2][0][0]) ;
2457   
2458   if (&WeightsArray != NULL) {
2459     dimension = min_degree + 1 ;
2460     Standard_Real *
2461       WArray = (Standard_Real *) 
2462         &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2463     PLib::EvalPolynomial(new_parameter[0],
2464                          MinIndMax,
2465                          max_degree,
2466                          dimension,
2467                          WArray[0],
2468                          locpoles[0]) ;
2469     
2470     PLib::EvalPolynomial(new_parameter[1],
2471                          MinIndMin,
2472                          min_degree,
2473                          1,
2474                          locpoles[0],
2475                          local_weights_array[0][0]) ;
2476     PLib::EvalPolynomial(new_parameter[1],
2477                          1,
2478                          min_degree,
2479                          1,
2480                          locpoles[dimension],
2481                          local_weights_array[1][0]) ;
2482     PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2483                          min_degree,
2484                          1,
2485                          min_degree,
2486                          locpoles[dimension + dimension],
2487                          local_weights_array[2][0]) ;
2488     
2489     
2490     local_poles_and_weights_array[0][0][0] = local_poles_array[0][0][0];
2491     local_poles_and_weights_array[0][0][1] = local_poles_array[0][0][1];
2492     local_poles_and_weights_array[0][0][2] = local_poles_array[0][0][2];
2493     local_poles_and_weights_array[0][1][0] = local_poles_array[0][1][0];
2494     local_poles_and_weights_array[0][1][1] = local_poles_array[0][1][1];
2495     local_poles_and_weights_array[0][1][2] = local_poles_array[0][1][2];
2496     local_poles_and_weights_array[0][2][0] = local_poles_array[0][2][0];
2497     local_poles_and_weights_array[0][2][1] = local_poles_array[0][2][1];
2498     local_poles_and_weights_array[0][2][2] = local_poles_array[0][2][2];
2499     
2500     local_poles_and_weights_array[1][0][0] = local_poles_array[1][0][0];
2501     local_poles_and_weights_array[1][0][1] = local_poles_array[1][0][1];
2502     local_poles_and_weights_array[1][0][2] = local_poles_array[1][0][2];
2503     local_poles_and_weights_array[1][1][0] = local_poles_array[1][1][0];
2504     local_poles_and_weights_array[1][1][1] = local_poles_array[1][1][1];
2505     local_poles_and_weights_array[1][1][2] = local_poles_array[1][1][2];
2506     local_poles_and_weights_array[1][2][0] = local_poles_array[1][2][0];
2507     local_poles_and_weights_array[1][2][1] = local_poles_array[1][2][1];
2508     local_poles_and_weights_array[1][2][2] = local_poles_array[1][2][2];
2509     
2510     local_poles_and_weights_array[2][0][0] = local_poles_array[2][0][0];
2511     local_poles_and_weights_array[2][0][1] = local_poles_array[2][0][1];
2512     local_poles_and_weights_array[2][0][2] = local_poles_array[2][0][2];
2513     local_poles_and_weights_array[2][1][0] = local_poles_array[2][1][0];
2514     local_poles_and_weights_array[2][1][1] = local_poles_array[2][1][1];
2515     local_poles_and_weights_array[2][1][2] = local_poles_array[2][1][2];
2516     local_poles_and_weights_array[2][2][0] = local_poles_array[2][2][0];
2517     local_poles_and_weights_array[2][2][1] = local_poles_array[2][2][1];
2518     local_poles_and_weights_array[2][2][2] = local_poles_array[2][2][2];
2519     
2520     
2521     local_poles_and_weights_array[0][0][3] = local_weights_array[0][0];
2522     local_poles_and_weights_array[0][1][3] = local_weights_array[0][1];
2523     local_poles_and_weights_array[0][2][3] = local_weights_array[0][2];
2524     local_poles_and_weights_array[1][0][3] = local_weights_array[1][0];
2525     local_poles_and_weights_array[1][1][3] = local_weights_array[1][1];
2526     local_poles_and_weights_array[1][2][3] = local_weights_array[1][2];
2527     local_poles_and_weights_array[2][0][3] = local_weights_array[2][0];
2528     local_poles_and_weights_array[2][1][3] = local_weights_array[2][1];
2529     local_poles_and_weights_array[2][2][3] = local_weights_array[2][2];
2530     
2531     BSplSLib::RationalDerivative(2,
2532                                  2,
2533                                  2,
2534                                  2,
2535                                  local_poles_and_weights_array[0][0][0],
2536                                  local_poles_array[0][0][0]) ;
2537   }
2538   
2539
2540   Standard_Real minmin = inverse_min * inverse_min;
2541   Standard_Real minmax = inverse_min * inverse_max;
2542   Standard_Real maxmax = inverse_max * inverse_max;
2543   
2544   my_point      [0] = local_poles_array              [0][0][0] ;
2545   my_vec_min    [0] = inverse_min * local_poles_array[0][1][0] ;
2546   my_vec_max    [0] = inverse_max * local_poles_array[1][0][0] ;
2547   my_vec_min_min[0] = minmin * local_poles_array     [0][2][0] ;
2548   my_vec_min_max[0] = minmax * local_poles_array     [1][1][0] ;
2549   my_vec_max_max[0] = maxmax * local_poles_array     [2][0][0] ;
2550
2551   my_point      [1] = local_poles_array              [0][0][1] ;
2552   my_vec_min    [1] = inverse_min * local_poles_array[0][1][1] ;
2553   my_vec_max    [1] = inverse_max * local_poles_array[1][0][1] ;
2554   my_vec_min_min[1] = minmin * local_poles_array     [0][2][1] ;
2555   my_vec_min_max[1] = minmax * local_poles_array     [1][1][1] ;
2556   my_vec_max_max[1] = maxmax * local_poles_array     [2][0][1] ;
2557
2558   my_point      [2] = local_poles_array              [0][0][2] ;
2559   my_vec_min    [2] = inverse_min * local_poles_array[0][1][2] ;
2560   my_vec_max    [2] = inverse_max * local_poles_array[1][0][2] ;
2561   my_vec_min_min[2] = minmin * local_poles_array     [0][2][2] ;
2562   my_vec_min_max[2] = minmax * local_poles_array     [1][1][2] ;
2563   my_vec_max_max[2] = maxmax * local_poles_array     [2][0][2] ;
2564 }
2565
2566 //=======================================================================
2567 //function : MovePoint
2568 //purpose  : Find the new poles which allows  an old point (with a
2569 //           given  u and v as parameters) to reach a new position
2570 //=======================================================================
2571
2572 void BSplSLib::MovePoint (const Standard_Real            U, 
2573                           const Standard_Real            V,
2574                           const gp_Vec&                  Displ,
2575                           const Standard_Integer         UIndex1,
2576                           const Standard_Integer         UIndex2,
2577                           const Standard_Integer         VIndex1,
2578                           const Standard_Integer         VIndex2,
2579                           const Standard_Integer         UDegree,
2580                           const Standard_Integer         VDegree,
2581                           const Standard_Boolean         Rational,
2582                           const TColgp_Array2OfPnt&      Poles,  
2583                           const TColStd_Array2OfReal&    Weights,
2584                           const TColStd_Array1OfReal&    UFlatKnots,
2585                           const TColStd_Array1OfReal&    VFlatKnots,
2586                           Standard_Integer&              UFirstIndex,
2587                           Standard_Integer&              ULastIndex,
2588                           Standard_Integer&              VFirstIndex,
2589                           Standard_Integer&              VLastIndex,
2590                           TColgp_Array2OfPnt&            NewPoles)
2591 {
2592   // calculate the UBSplineBasis in the parameter U
2593   Standard_Integer UFirstNonZeroBsplineIndex;
2594   math_Matrix UBSplineBasis(1, 1,
2595                             1, UDegree+1);
2596   Standard_Integer ErrorCod1 =  BSplCLib::EvalBsplineBasis(1,
2597                                                            0,
2598                                                            UDegree+1,
2599                                                            UFlatKnots,
2600                                                            U,
2601                                                            UFirstNonZeroBsplineIndex,
2602                                                            UBSplineBasis);  
2603   // calculate the VBSplineBasis in the parameter V
2604   Standard_Integer VFirstNonZeroBsplineIndex;
2605   math_Matrix VBSplineBasis(1, 1,
2606                             1, VDegree+1);
2607   Standard_Integer ErrorCod2 =  BSplCLib::EvalBsplineBasis(1,
2608                                                            0,
2609                                                            VDegree+1,
2610                                                            VFlatKnots,
2611                                                            V,
2612                                                            VFirstNonZeroBsplineIndex,
2613                                                            VBSplineBasis);  
2614   if (ErrorCod1 || ErrorCod2) {
2615     UFirstIndex = 0;
2616     ULastIndex = 0;
2617     VFirstIndex = 0;
2618     VLastIndex = 0;
2619     return;
2620   }
2621   
2622   // find the span which is predominant for parameter U
2623   UFirstIndex = UFirstNonZeroBsplineIndex;
2624   ULastIndex = UFirstNonZeroBsplineIndex + UDegree ;
2625   if (UFirstIndex < UIndex1) UFirstIndex = UIndex1;
2626   if (ULastIndex > UIndex2) ULastIndex = UIndex2;
2627
2628   Standard_Real maxValue = 0.0;
2629   Standard_Integer i, ukk1=0, ukk2;
2630
2631   for (i = UFirstIndex-UFirstNonZeroBsplineIndex+1; i <= ULastIndex-UFirstNonZeroBsplineIndex+1; i++) {
2632     if (UBSplineBasis(1,i) > maxValue) {
2633       ukk1 = i + UFirstNonZeroBsplineIndex - 1;
2634       maxValue = UBSplineBasis(1,i);
2635     }
2636   }
2637
2638   // find a ukk2 if symetriy
2639   ukk2 = ukk1;
2640   i = ukk1 - UFirstNonZeroBsplineIndex + 2;
2641   if ((ukk1+1) <= ULastIndex) {
2642     if (Abs(UBSplineBasis(1, ukk1-UFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2643       ukk2 = ukk1+1;
2644     }
2645   }
2646
2647   // find the span which is predominant for parameter V
2648   VFirstIndex = VFirstNonZeroBsplineIndex;
2649   VLastIndex = VFirstNonZeroBsplineIndex + VDegree ;
2650
2651   if (VFirstIndex < VIndex1) VFirstIndex = VIndex1;
2652   if (VLastIndex > VIndex2) VLastIndex = VIndex2;
2653
2654   maxValue = 0.0;
2655   Standard_Integer j, vkk1=0, vkk2;
2656
2657   for (j = VFirstIndex-VFirstNonZeroBsplineIndex+1; j <= VLastIndex-VFirstNonZeroBsplineIndex+1; j++) {
2658     if (VBSplineBasis(1,j) > maxValue) {
2659       vkk1 = j + VFirstNonZeroBsplineIndex - 1;
2660       maxValue = VBSplineBasis(1,j);
2661     }
2662   }
2663
2664   // find a vkk2 if symetriy
2665   vkk2 = vkk1;
2666   j = vkk1 - VFirstNonZeroBsplineIndex + 2;
2667   if ((vkk1+1) <= VLastIndex) {
2668     if (Abs(VBSplineBasis(1, vkk1-VFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2669       vkk2 = vkk1+1;
2670     }
2671   }
2672
2673   // compute the vector of displacement
2674   Standard_Real D1 = 0.0;
2675   Standard_Real D2 = 0.0;
2676   Standard_Real hN, Coef, DvalU, DvalV;
2677
2678   Standard_Integer ii, jj;
2679
2680   for (i = 1; i <= UDegree+1; i++) {
2681     ii = i + UFirstNonZeroBsplineIndex - 1;
2682     if (ii < ukk1) {
2683       DvalU = ukk1-ii;
2684     }
2685     else if (ii > ukk2) {
2686       DvalU = ii - ukk2;
2687     }
2688     else {
2689       DvalU = 0.0;
2690     }
2691
2692     for (j = 1; j <= VDegree+1; j++) {
2693       jj = j + VFirstNonZeroBsplineIndex - 1;
2694       if (Rational) {
2695         hN = Weights(ii, jj)*UBSplineBasis(1, i)*VBSplineBasis(1,j);
2696         D2 += hN;
2697       }
2698       else {
2699         hN = UBSplineBasis(1, i)*VBSplineBasis(1,j);
2700       }
2701       if (ii >= UFirstIndex && ii <= ULastIndex && jj >= VFirstIndex && jj <= VLastIndex) {
2702         if (jj < vkk1) {
2703           DvalV = vkk1-jj;
2704         }
2705         else if (jj > vkk2) {
2706           DvalV = jj - vkk2;
2707         }
2708         else {
2709           DvalV = 0.0;
2710         }
2711         D1 += 1./(DvalU + DvalV + 1.) * hN;
2712       }
2713     }
2714   }
2715   
2716   if (Rational) {
2717     Coef = D2/D1;
2718   }
2719   else {
2720     Coef = 1./D1;
2721   }
2722
2723   // compute the new poles
2724
2725   for (i=Poles.LowerRow(); i<=Poles.UpperRow(); i++) {
2726     if (i < ukk1) {
2727       DvalU = ukk1-i;
2728     }
2729     else if (i > ukk2) {
2730       DvalU = i - ukk2;
2731     }
2732     else {
2733       DvalU = 0.0;
2734     }
2735
2736     for (j=Poles.LowerCol(); j<=Poles.UpperCol(); j++) {
2737       if (i >= UFirstIndex && i <= ULastIndex && j >= VFirstIndex && j <= VLastIndex) {
2738         if (j < vkk1) {
2739           DvalV = vkk1-j;
2740         }
2741         else if (j > vkk2) {
2742           DvalV = j - vkk2;
2743         }
2744         else {
2745           DvalV = 0.0;
2746         }
2747         NewPoles(i,j) = Poles(i,j).Translated((Coef/(DvalU + DvalV + 1.))*Displ);
2748       }
2749       else {
2750         NewPoles(i,j) = Poles(i,j);
2751       }
2752     }
2753   }
2754 }
2755
2756 //=======================================================================
2757 // function : Resolution
2758 // purpose  : this computes an estimate for the maximum of the 
2759 // partial derivatives both in U and in V
2760 //
2761 //
2762 // The calculation resembles at the calculation of curves with 
2763 // additional index for the control point. Let Si,j be the
2764 // control points for ls surface  and  Di,j  the weights.  
2765 // The checking of upper bounds for the partial derivatives 
2766 // will be omitted and Su is the next upper bound in the polynomial case :
2767 //
2768 //
2769 //
2770 //                        |  Si,j - Si-1,j  |
2771 //          d *   Max     |  -------------  |
2772 //                i = 2,n |     ti+d - ti   |
2773 //                i=1.m
2774 //
2775 //
2776 // and in the rational case :
2777 //
2778 //
2779 //
2780 //                         Di,j * (Si,j - Sk,j) - Di-1,j * (Si-1,j - Sk,j)
2781 //   Max   Max       d  *  -----------------------------------------------
2782 // k=1,n  i dans Rj                   ti+d  - ti
2783 // j=1,m
2784 //  ----------------------------------------------------------------------
2785 //
2786 //               Min    Di,j
2787 //              i=1,n
2788 //              j=1,m
2789 //
2790 //
2791 //
2792 // with Rj = {j-d, ....,  j+d+d+1}.
2793 //
2794 //
2795 //=======================================================================
2796
2797 void BSplSLib::Resolution(const TColgp_Array2OfPnt&      Poles,
2798                           const TColStd_Array2OfReal&    Weights,
2799                           const TColStd_Array1OfReal&    UKnots,
2800                           const TColStd_Array1OfReal&    VKnots,
2801                           const TColStd_Array1OfInteger& UMults,
2802                           const TColStd_Array1OfInteger& VMults,
2803                           const Standard_Integer         UDegree,
2804                           const Standard_Integer         VDegree,
2805                           const Standard_Boolean         URational,
2806                           const Standard_Boolean         VRational,
2807                           const Standard_Boolean         UPeriodic,
2808                           const Standard_Boolean         VPeriodic,
2809                           const Standard_Real            Tolerance3D,
2810                           Standard_Real&                 UTolerance,
2811                           Standard_Real&                 VTolerance)
2812 {
2813   Standard_Real Wij,Wmj,Wji,Wjm;
2814   Standard_Real Xij,Xmj,Xji,Xjm,Xpq,Xqp;
2815   Standard_Real Yij,Ymj,Yji,Yjm,Ypq,Yqp;
2816   Standard_Real Zij,Zmj,Zji,Zjm,Zpq,Zqp;
2817   Standard_Real factor,value,min,min_weights=0,inverse,max_derivative[2];
2818
2819   max_derivative[0] = max_derivative[1] = 0.0e0 ; 
2820   
2821   Standard_Integer PRowLength, PColLength;
2822   Standard_Integer ii,jj,pp,qq,ii_index,jj_index,pp_index,qq_index;
2823   Standard_Integer ii_minus,upper[2],lower[2],poles_length[2];
2824   Standard_Integer num_poles[2],num_flat_knots[2];
2825   
2826   num_flat_knots[0] = 
2827     BSplCLib::KnotSequenceLength(UMults,
2828                                  UDegree,
2829                                  UPeriodic) ;
2830   num_flat_knots[1] = 
2831     BSplCLib::KnotSequenceLength(VMults,
2832                                  VDegree,
2833                                  VPeriodic) ;
2834   TColStd_Array1OfReal  flat_knots_in_u(1,num_flat_knots[0]) ;
2835   TColStd_Array1OfReal  flat_knots_in_v(1,num_flat_knots[1]) ;
2836   BSplCLib::KnotSequence(UKnots,
2837                          UMults,
2838                          UDegree,
2839                          UPeriodic,
2840                          flat_knots_in_u) ;
2841   BSplCLib::KnotSequence(VKnots,
2842                          VMults,
2843                          VDegree,
2844                          VPeriodic,
2845                          flat_knots_in_v) ;
2846   PRowLength = Poles.RowLength();
2847   PColLength = Poles.ColLength();
2848   if (URational || VRational) {
2849     Standard_Integer Wsize = PRowLength * PColLength;
2850     const Standard_Real * WG = &Weights(Weights.LowerRow(),Weights.LowerCol());
2851     min_weights = WG[0];
2852     
2853     for (ii = 1 ; ii < Wsize ; ii++) {
2854       min = WG[ii];
2855       if (min_weights > min) min_weights = min;
2856     }
2857   }
2858   Standard_Integer UD1 = UDegree + 1;
2859   Standard_Integer VD1 = VDegree + 1;
2860   num_poles[0] = num_flat_knots[0] - UD1;
2861   num_poles[1] = num_flat_knots[1] - VD1;
2862   poles_length[0] = PColLength;
2863   poles_length[1] = PRowLength;
2864   if(URational) {
2865     Standard_Integer UD2 = UDegree << 1;
2866     Standard_Integer VD2 = VDegree << 1;
2867
2868     for (ii = 2 ; ii <= num_poles[0] ; ii++) {
2869       ii_index = (ii - 1) % poles_length[0] + 1 ;
2870       ii_minus = (ii - 2) % poles_length[0] + 1 ;
2871       inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
2872       inverse = 1.0e0 / inverse ;
2873       lower[0] = ii - UD1;
2874       if (lower[0] < 1) lower[0] = 1;
2875       upper[0] = ii + UD2 + 1;
2876       if (upper[0] > num_poles[0]) upper[0] = num_poles[0];
2877
2878       for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
2879         jj_index = (jj - 1) % poles_length[1] + 1 ;
2880         lower[1] = jj - VD1;
2881         if (lower[1] < 1) lower[1] = 1;
2882         upper[1] = jj + VD2 + 1;
2883         if (upper[1] > num_poles[1]) upper[1] = num_poles[1];
2884         const gp_Pnt& Pij = Poles  .Value(ii_index,jj_index);
2885         Wij               = Weights.Value(ii_index,jj_index);
2886         const gp_Pnt& Pmj = Poles  .Value(ii_minus,jj_index);
2887         Wmj               = Weights.Value(ii_minus,jj_index);
2888         Xij = Pij.X();
2889         Yij = Pij.Y();
2890         Zij = Pij.Z();
2891         Xmj = Pmj.X();
2892         Ymj = Pmj.Y();
2893         Zmj = Pmj.Z();
2894         
2895         for (pp = lower[0] ; pp <= upper[0] ; pp++) {
2896           pp_index = (pp - 1) % poles_length[0] + 1 ;
2897
2898           for (qq = lower[1] ; qq <= upper[1] ; qq++) {
2899             value = 0.0e0 ;
2900             qq_index = (qq - 1) % poles_length[1] + 1 ;
2901             const gp_Pnt& Ppq = Poles.Value(pp_index,qq_index);
2902             Xpq = Ppq.X();
2903             Ypq = Ppq.Y();
2904             Zpq = Ppq.Z();
2905             factor  = (Xpq - Xij) * Wij;
2906             factor -= (Xpq - Xmj) * Wmj;
2907             if (factor < 0) factor = - factor;
2908             value  += factor ;
2909             factor  = (Ypq - Yij) * Wij;
2910             factor -= (Ypq - Ymj) * Wmj;
2911             if (factor < 0) factor = - factor;
2912             value  += factor ;
2913             factor  = (Zpq - Zij) * Wij;
2914             factor -= (Zpq - Zmj) * Wmj;
2915             if (factor < 0) factor = - factor;
2916             value  += factor ;
2917             value *= inverse ;
2918             if (max_derivative[0] < value) max_derivative[0] = value ;
2919           }
2920         }
2921       }
2922     }
2923     max_derivative[0] /= min_weights ;
2924   }
2925   else {
2926
2927     for (ii = 2 ; ii <= num_poles[0] ; ii++) {
2928       ii_index = (ii - 1) % poles_length[0] + 1 ;
2929       ii_minus = (ii - 2) % poles_length[0] + 1 ;
2930       inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
2931       inverse = 1.0e0 / inverse ;
2932
2933       for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
2934         jj_index = (jj - 1) % poles_length[1] + 1 ;
2935         value = 0.0e0 ;
2936         const gp_Pnt& Pij = Poles.Value(ii_index,jj_index);
2937         const gp_Pnt& Pmj = Poles.Value(ii_minus,jj_index);
2938         factor = Pij.X() - Pmj.X();
2939         if (factor < 0) factor = - factor;
2940         value += factor;
2941         factor = Pij.Y() - Pmj.Y();
2942         if (factor < 0) factor = - factor;
2943         value += factor;
2944         factor = Pij.Z() - Pmj.Z();
2945         if (factor < 0) factor = - factor;
2946         value += factor;
2947         value *= inverse ;
2948         if (max_derivative[0] < value) max_derivative[0] = value ;
2949       }
2950     }
2951   }
2952   max_derivative[0] *= UDegree ;
2953   if(VRational) {
2954     Standard_Integer UD2 = UDegree << 1;
2955     Standard_Integer VD2 = VDegree << 1;
2956
2957     for (ii = 2 ; ii <= num_poles[1] ; ii++) {
2958       ii_index = (ii - 1) % poles_length[1] + 1 ;
2959       ii_minus = (ii - 2) % poles_length[1] + 1 ;
2960       inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
2961       inverse = 1.0e0 / inverse ;
2962       lower[0] = ii - VD1;
2963       if (lower[0] < 1) lower[0] = 1;
2964       upper[0] = ii + VD2 + 1;
2965       if (upper[0] > num_poles[1]) upper[0] = num_poles[1];
2966
2967       for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
2968         jj_index = (jj - 1) % poles_length[0] + 1 ;
2969         lower[1] = jj - UD1;
2970         if (lower[1] < 1) lower[1] = 1;
2971         upper[1] = jj + UD2 + 1;
2972         if (upper[1] > num_poles[0]) upper[1] = num_poles[0];
2973         const gp_Pnt& Pji = Poles  .Value(jj_index,ii_index);
2974         Wji               = Weights.Value(jj_index,ii_index);
2975         const gp_Pnt& Pjm = Poles  .Value(jj_index,ii_minus);
2976         Wjm               = Weights.Value(jj_index,ii_minus);
2977         Xji = Pji.X();
2978         Yji = Pji.Y();
2979         Zji = Pji.Z();
2980         Xjm = Pjm.X();
2981         Yjm = Pjm.Y();
2982         Zjm = Pjm.Z();
2983         
2984         for (pp = lower[1] ; pp <= upper[1] ; pp++) {
2985           pp_index = (pp - 1) % poles_length[1] + 1 ;
2986
2987           for (qq = lower[0] ; qq <= upper[0] ; qq++) {
2988             value = 0.0e0 ;
2989             qq_index = (qq - 1) % poles_length[0] + 1 ;
2990             const gp_Pnt& Pqp = Poles.Value(qq_index,pp_index);
2991             Xqp = Pqp.X();
2992             Yqp = Pqp.Y();
2993             Zqp = Pqp.Z();
2994             factor  = (Xqp - Xji) * Wji;
2995             factor -= (Xqp - Xjm) * Wjm;
2996             if (factor < 0) factor = - factor;
2997             value += factor ;
2998             factor  = (Yqp - Yji) * Wji;
2999             factor -= (Yqp - Yjm) * Wjm;
3000             if (factor < 0) factor = - factor;
3001             value += factor ;
3002             factor  = (Zqp - Zji) * Wji;
3003             factor -= (Zqp - Zjm) * Wjm;
3004             if (factor < 0) factor = - factor;
3005             value += factor ;
3006             value *= inverse ;
3007             if (max_derivative[1] < value) max_derivative[1] = value ;
3008           }
3009         }
3010       }
3011     }
3012     max_derivative[1] /= min_weights ;
3013   }
3014   else {
3015
3016     for (ii = 2 ; ii <= num_poles[1] ; ii++) {
3017       ii_index = (ii - 1) % poles_length[1] + 1 ;
3018       ii_minus = (ii - 2) % poles_length[1] + 1 ;
3019       inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
3020       inverse = 1.0e0 / inverse ;
3021
3022       for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
3023         jj_index = (jj - 1) % poles_length[0] + 1 ;
3024         value = 0.0e0 ;
3025         const gp_Pnt& Pji = Poles.Value(jj_index,ii_index);
3026         const gp_Pnt& Pjm = Poles.Value(jj_index,ii_minus);
3027         factor = Pji.X() - Pjm.X() ;
3028         if (factor < 0) factor = - factor;
3029         value += factor;
3030         factor = Pji.Y() - Pjm.Y() ;
3031         if (factor < 0) factor = - factor;
3032         value += factor;
3033         factor = Pji.Z() - Pjm.Z() ;
3034         if (factor < 0) factor = - factor;
3035         value += factor;
3036         value *= inverse ;
3037         if (max_derivative[1] < value) max_derivative[1] = value ;
3038       }
3039     }
3040   }
3041   max_derivative[1] *= VDegree ;
3042   max_derivative[0] *= M_SQRT2 ;
3043   max_derivative[1] *= M_SQRT2 ;
3044   if(max_derivative[0] && max_derivative[1]) { 
3045     UTolerance = Tolerance3D / max_derivative[0] ;
3046     VTolerance = Tolerance3D / max_derivative[1] ;
3047   }
3048   else { 
3049     UTolerance=VTolerance=0.0;
3050 #ifdef DEB
3051     cout<<"ElSLib.cxx : maxderivative = 0.0 "<<endl;
3052 #endif
3053   }
3054 }
3055
3056 //=======================================================================
3057 //function : Interpolate
3058 //purpose  : 
3059 //=======================================================================
3060
3061 void BSplSLib::Interpolate(const Standard_Integer UDegree,
3062                            const Standard_Integer VDegree, 
3063                            const TColStd_Array1OfReal& UFlatKnots,
3064                            const TColStd_Array1OfReal& VFlatKnots,
3065                            const TColStd_Array1OfReal& UParameters,
3066                            const TColStd_Array1OfReal& VParameters,
3067                            TColgp_Array2OfPnt& Poles,
3068                            TColStd_Array2OfReal& Weights, 
3069                            Standard_Integer& InversionProblem)
3070 {
3071   Standard_Integer ii, jj, ll, kk, dimension;
3072   Standard_Integer ULength = UParameters.Length();
3073   Standard_Integer VLength = VParameters.Length();
3074   Standard_Real * poles_array;
3075   
3076   // extraction of iso u
3077   dimension = 4*ULength;
3078   TColStd_Array2OfReal Points(1, VLength, 
3079                               1, dimension);
3080   
3081   Handle(TColStd_HArray1OfInteger) ContactOrder = 
3082     new (TColStd_HArray1OfInteger)(1, VLength);
3083   ContactOrder->Init(0);
3084   
3085   for (ii=1; ii <= VLength; ii++) {
3086
3087     for (jj=1, ll=1; jj<= ULength; jj++, ll+=4) {
3088       Points(ii,ll)   = Poles(jj, ii).X();
3089       Points(ii,ll+1) = Poles(jj, ii).Y();
3090       Points(ii,ll+2) = Poles(jj, ii).Z();
3091       Points(ii,ll+3) = Weights(jj,ii)  ;
3092     }
3093   }
3094
3095   // interpolation of iso u
3096   poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3097   BSplCLib::Interpolate(VDegree,
3098                         VFlatKnots,
3099                         VParameters,
3100                         ContactOrder->Array1(),
3101                         dimension,
3102                         poles_array[0],
3103                         InversionProblem) ;
3104   if (InversionProblem != 0) return;
3105
3106   // extraction of iso v
3107
3108   dimension = VLength*4;
3109   TColStd_Array2OfReal IsoPoles(1, ULength, 
3110                                 1, dimension);
3111   
3112   ContactOrder =  new (TColStd_HArray1OfInteger)(1, ULength);
3113   ContactOrder->Init(0);
3114   poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3115
3116   for (ii=1, kk=1; ii <= ULength; ii++, kk+=4) {
3117
3118     for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3119       IsoPoles (ii,ll)   = Points(jj, kk);
3120       IsoPoles (ii,ll+1) = Points(jj, kk+1);
3121       IsoPoles (ii,ll+2) = Points(jj, kk+2);
3122       IsoPoles (ii,ll+3) = Points(jj, kk+3);
3123     }
3124   }
3125   // interpolation of iso v
3126   BSplCLib::Interpolate(UDegree,
3127                         UFlatKnots,
3128                         UParameters,
3129                         ContactOrder->Array1(),
3130                         dimension,
3131                         poles_array[0],
3132                         InversionProblem);
3133
3134   // return results
3135
3136   for (ii=1; ii <= ULength; ii++) {
3137
3138     for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3139       gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3140       Poles.SetValue(ii, jj, Pnt);
3141       Weights.SetValue(ii,jj,IsoPoles(ii,ll+3)) ;
3142     }
3143   }
3144 }
3145
3146 //=======================================================================
3147 //function : Interpolate
3148 //purpose  : 
3149 //=======================================================================
3150
3151 void BSplSLib::Interpolate(const Standard_Integer UDegree,
3152                            const Standard_Integer VDegree, 
3153                            const TColStd_Array1OfReal& UFlatKnots,
3154                            const TColStd_Array1OfReal& VFlatKnots,
3155                            const TColStd_Array1OfReal& UParameters,
3156                            const TColStd_Array1OfReal& VParameters,
3157                            TColgp_Array2OfPnt& Poles,
3158                            Standard_Integer& InversionProblem)
3159 {
3160   Standard_Integer ii, jj, ll, kk, dimension;
3161   Standard_Integer ULength = UParameters.Length();
3162   Standard_Integer VLength = VParameters.Length();
3163   Standard_Real * poles_array;
3164   
3165   // extraction of iso u
3166   dimension = 3*ULength;
3167   TColStd_Array2OfReal Points(1, VLength, 
3168                               1, dimension);
3169   
3170   Handle(TColStd_HArray1OfInteger) ContactOrder = 
3171     new (TColStd_HArray1OfInteger)(1, VLength);
3172   ContactOrder->Init(0);
3173   
3174   for (ii=1; ii <= VLength; ii++) {
3175     
3176     for (jj=1, ll=1; jj<= ULength; jj++, ll+=3) {
3177       Points(ii,ll)   = Poles(jj, ii).X();
3178       Points(ii,ll+1) = Poles(jj, ii).Y();
3179       Points(ii,ll+2) = Poles(jj, ii).Z();
3180     }
3181   }
3182   
3183   // interpolation of iso u
3184   poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3185   BSplCLib::Interpolate(VDegree,
3186                         VFlatKnots,
3187                         VParameters,
3188                         ContactOrder->Array1(),
3189                         dimension,
3190                         poles_array[0],
3191                         InversionProblem) ;
3192   if (InversionProblem != 0) return;
3193   
3194   // extraction of iso v
3195   
3196   dimension = VLength*3;
3197   TColStd_Array2OfReal IsoPoles(1, ULength, 
3198                                 1, dimension);
3199   
3200   ContactOrder =  new (TColStd_HArray1OfInteger)(1, ULength);
3201   ContactOrder->Init(0);
3202   poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3203   
3204   for (ii=1, kk=1; ii <= ULength; ii++, kk+=3) {
3205
3206     for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3207       IsoPoles (ii,ll)   = Points(jj, kk);
3208       IsoPoles (ii,ll+1) = Points(jj, kk+1);
3209       IsoPoles (ii,ll+2) = Points(jj, kk+2);
3210     }
3211   }
3212   // interpolation of iso v
3213   BSplCLib::Interpolate(UDegree,
3214                         UFlatKnots,
3215                         UParameters,
3216                         ContactOrder->Array1(),
3217                         dimension,
3218                         poles_array[0],
3219                         InversionProblem);
3220   
3221   // return results
3222
3223   for (ii=1; ii <= ULength; ii++) {
3224
3225     for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3226       gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3227       Poles.SetValue(ii, jj, Pnt);
3228     }
3229   }
3230 }
3231
3232 //=======================================================================
3233 //function : FunctionMultiply
3234 //purpose  : 
3235 //=======================================================================
3236
3237 void BSplSLib::FunctionMultiply
3238 (const BSplSLib_EvaluatorFunction&           Function,
3239  const Standard_Integer                      UBSplineDegree,
3240  const Standard_Integer                      VBSplineDegree,
3241  const TColStd_Array1OfReal&                 UBSplineKnots,
3242  const TColStd_Array1OfReal&                 VBSplineKnots,
3243  const TColStd_Array1OfInteger &             UMults,
3244  const TColStd_Array1OfInteger &             VMults,
3245  const TColgp_Array2OfPnt&                   Poles,
3246  const TColStd_Array2OfReal&                 Weights,
3247  const TColStd_Array1OfReal&                 UFlatKnots,
3248  const TColStd_Array1OfReal&                 VFlatKnots,
3249  const Standard_Integer                      UNewDegree,
3250  const Standard_Integer                      VNewDegree,
3251  TColgp_Array2OfPnt&                         NewNumerator,
3252  TColStd_Array2OfReal&                       NewDenominator,
3253  Standard_Integer&                           Status) 
3254 {
3255   Standard_Integer num_uparameters,
3256 //  ii,jj,kk,
3257   ii,jj,
3258   error_code,
3259   num_vparameters ;
3260   Standard_Real    result ;
3261   
3262   num_uparameters = UFlatKnots.Length() - UNewDegree - 1 ;
3263   num_vparameters = VFlatKnots.Length() - VNewDegree - 1 ;
3264   TColStd_Array1OfReal  UParameters(1,num_uparameters) ;
3265   TColStd_Array1OfReal  VParameters(1,num_vparameters) ;
3266   
3267   if ((NewNumerator.ColLength() == num_uparameters) &&
3268       (NewNumerator.RowLength() == num_vparameters) &&
3269       (NewDenominator.ColLength() == num_uparameters) &&
3270       (NewDenominator.RowLength() == num_vparameters)) {
3271     
3272     
3273     BSplCLib::BuildSchoenbergPoints(UNewDegree,
3274                                     UFlatKnots,
3275                                     UParameters) ;
3276     
3277     BSplCLib::BuildSchoenbergPoints(VNewDegree,
3278                                     VFlatKnots,
3279                                     VParameters) ;
3280     
3281     for (ii = 1 ; ii <= num_uparameters ; ii++) {
3282
3283       for (jj = 1 ; jj <= num_vparameters ; jj++) {
3284         HomogeneousD0(UParameters(ii),
3285                       VParameters(jj),
3286                       0,
3287                       0,
3288                       Poles,
3289                       Weights,
3290                       UBSplineKnots,
3291                       VBSplineKnots,
3292                       UMults,
3293                       VMults,
3294                       UBSplineDegree,
3295                       VBSplineDegree,
3296                       Standard_True,
3297                       Standard_True,
3298                       Standard_False,
3299                       Standard_False,
3300                       NewDenominator(ii,jj),
3301                       NewNumerator(ii,jj)) ;
3302         
3303         Function.Evaluate (0,
3304                  UParameters(ii),
3305                  VParameters(jj),
3306                  result,
3307                  error_code) ;
3308         if (error_code) {
3309           Standard_ConstructionError::Raise();
3310         }
3311         gp_Pnt& P = NewNumerator(ii,jj);
3312         P.SetX(P.X() * result);
3313         P.SetY(P.Y() * result);
3314         P.SetZ(P.Z() * result);
3315         NewDenominator(ii,jj) *= result ;
3316       }
3317     }
3318     Interpolate(UNewDegree,
3319                 VNewDegree, 
3320                 UFlatKnots,
3321                 VFlatKnots,
3322                 UParameters,
3323                 VParameters,
3324                 NewNumerator,
3325                 NewDenominator, 
3326                 Status) ;
3327   }
3328   else {
3329     Standard_ConstructionError::Raise();
3330   }
3331 }