90fff069666065bf022c0c59fb20d1200ae8b72a
[occt.git] / src / GeomLib / GeomLib.cxx
1 // Created on: 1993-07-07
2 // Created by: Jean Claude VAUTHIER
3 // Copyright (c) 1993-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 // Version:     
18 //pmn 24/09/96 Ajout du prolongement de courbe.
19 //              jct 15/04/97 Ajout du prolongement de surface.
20 //              jct 24/04/97 simplification ou suppression de calculs
21 //                           inutiles dans ExtendSurfByLength
22 //                           correction de Tbord et Continuity=0 accepte
23 //                           correction du calcul de lambda et appel a
24 //                           TangExtendToConstraint avec lambmin au lieu de 1.
25 //                           correction du passage Sr rat --> BSp nD
26 //              xab 26/06/97 treatement partiel anulation des derivees 
27 //                           partiels du denonimateur des Surfaces BSplines Rationnelles
28 //                           dans le cas de valeurs proportionnelles des denominateurs
29 //                           en umin umax et/ou vmin vmax.
30 //              pmn 4/07/97  Gestion de la continuite dans BuildCurve3d (PRO9097)
31
32 //              xab 10/07/97 on revient en arriere sur l'ajout du 26/06/97
33 //              pmn 26/09/97 Ajout des parametres d'approx dans BuildCurve3d
34 //              xab 29/09/97 on reintegre l'ajout du 26/06/97
35 //              pmn 31/10/97 Ajoute AdjustExtremity
36 //              jct 26/11/98 blindage dans ExtendSurf qd NTgte = 0 (CTS21288)
37 //              jct 19/01/99 traitement de la periodicite dans ExtendSurf
38 // Design:       
39 // Warning:      None   
40 // References:   None   
41 // Language:     C++2.0 
42 // Purpose:     
43
44 // Declarations:        
45
46 #include <GeomLib.ixx>
47
48 #include <Precision.hxx>
49 #include <GeomConvert.hxx>
50 #include <Hermit.hxx>
51 #include <Standard_NotImplemented.hxx>
52 #include <GeomLib_MakeCurvefromApprox.hxx>
53 #include <GeomLib_DenominatorMultiplier.hxx>
54 #include <GeomLib_DenominatorMultiplierPtr.hxx>
55 #include <GeomLib_PolyFunc.hxx>
56 #include <GeomLib_LogSample.hxx>
57
58 #include <AdvApprox_ApproxAFunction.hxx>
59 #include <AdvApprox_PrefAndRec.hxx>
60
61 #include <Adaptor2d_HCurve2d.hxx>
62 #include <Adaptor3d_HCurve.hxx>
63 #include <Adaptor3d_HSurface.hxx>
64 #include <Adaptor3d_CurveOnSurface.hxx>
65 #include <Geom2dAdaptor_Curve.hxx>
66 #include <GeomAdaptor_Surface.hxx>
67 #include <GeomAdaptor_HSurface.hxx>
68 #include <Geom2dAdaptor_HCurve.hxx>
69 #include <Geom2dAdaptor_GHCurve.hxx>
70
71 #include <Geom2d_BSplineCurve.hxx>
72 #include <Geom_BSplineCurve.hxx>
73 #include <Geom2d_BezierCurve.hxx>
74 #include <Geom_BezierCurve.hxx>
75 #include <Geom_RectangularTrimmedSurface.hxx>
76 #include <Geom_Plane.hxx>
77 #include <Geom_Line.hxx>
78 #include <Geom2d_Line.hxx>
79 #include <Geom_Circle.hxx>
80 #include <Geom2d_Circle.hxx>
81 #include <Geom_Ellipse.hxx>
82 #include <Geom2d_Ellipse.hxx>
83 #include <Geom_Parabola.hxx>
84 #include <Geom2d_Parabola.hxx>
85 #include <Geom_Hyperbola.hxx>
86 #include <Geom2d_Hyperbola.hxx>
87 #include <Geom_TrimmedCurve.hxx>
88 #include <Geom2d_TrimmedCurve.hxx>
89 #include <Geom_OffsetCurve.hxx>
90 #include <Geom2d_OffsetCurve.hxx>
91 #include <Geom_BezierSurface.hxx>
92 #include <Geom_BSplineSurface.hxx>
93
94 #include <BSplCLib.hxx>
95 #include <BSplSLib.hxx>
96 #include <PLib.hxx>
97 #include <math_Matrix.hxx>
98 #include <math_Vector.hxx>
99 #include <math_Jacobi.hxx>
100 #include <math.hxx>
101 #include <math_FunctionAllRoots.hxx>
102 #include <math_FunctionSample.hxx>
103
104 #include <TColStd_HArray1OfReal.hxx>
105 #include <TColgp_Array1OfPnt.hxx>
106 #include <TColgp_Array1OfVec.hxx>
107 #include <TColgp_Array2OfPnt.hxx>
108 #include <TColgp_HArray2OfPnt.hxx>
109 #include <TColgp_Array1OfPnt2d.hxx>
110 #include <TColgp_Array1OfXYZ.hxx>
111 #include <TColStd_Array1OfReal.hxx>
112 #include <TColStd_Array2OfReal.hxx>
113 #include <TColStd_HArray2OfReal.hxx>
114 #include <TColStd_Array1OfInteger.hxx>
115
116 #include <gp_TrsfForm.hxx>
117 #include <gp_Lin.hxx>
118 #include <gp_Lin2d.hxx>
119 #include <gp_Circ.hxx>
120 #include <gp_Circ2d.hxx>
121 #include <gp_Elips.hxx>
122 #include <gp_Elips2d.hxx>
123 #include <gp_Hypr.hxx>
124 #include <gp_Hypr2d.hxx>
125 #include <gp_Parab.hxx>
126 #include <gp_Parab2d.hxx>
127 #include <gp_GTrsf2d.hxx>
128 #include <gp_Trsf2d.hxx>
129
130 #include <ElCLib.hxx>
131 #include <Geom2dConvert.hxx>
132 #include <GeomConvert_CompCurveToBSplineCurve.hxx>
133 #include <GeomConvert_ApproxSurface.hxx>
134
135 #include <CSLib.hxx>
136 #include <CSLib_NormalStatus.hxx>
137
138
139 #include <Standard_ConstructionError.hxx>
140
141 //=======================================================================
142 //function : ComputeLambda
143 //purpose  : Calcul le facteur lambda qui minimise la variation de vittesse
144 //           sur une interpolation d'hermite d'ordre (i,0)
145 //=======================================================================
146 static void ComputeLambda(const math_Matrix& Constraint,
147                           const math_Matrix& Hermit,
148                           const Standard_Real Length,
149                           Standard_Real& Lambda )
150 {
151   Standard_Integer size = Hermit.RowNumber();
152   Standard_Integer Continuity = size-2;
153   Standard_Integer ii, jj, ip, pp;
154
155   //Minimization
156   math_Matrix HDer(1, size-1, 1, size);
157   for (jj=1; jj<=size; jj++) {
158     for (ii=1; ii<size;ii++) {
159       HDer(ii, jj) = ii*Hermit(jj, ii+1);
160     }
161   }
162
163   math_Vector V(1, size);
164   math_Vector Vec1(1, Constraint.RowNumber());
165   math_Vector Vec2(1, Constraint.RowNumber());
166   math_Vector Vec3(1, Constraint.RowNumber()); 
167   math_Vector Vec4(1, Constraint.RowNumber());  
168
169   Standard_Real * polynome = &HDer(1,1);
170   Standard_Real * valhder =  &V(1);
171   Vec2 =  Constraint.Col(2);
172   Vec2 /= Length;
173   Standard_Real t,  squared1 = Vec2.Norm2(), GW;
174 //  math_Matrix Vec(1, Constraint.RowNumber(), 1, size-1);
175 //  gp_Vec Vfirst(p0.XYZ()), Vlast(Point.XYZ());
176 //  TColgp_Array1OfVec Der(2, 4);
177 //  Der(2) = d1; Der(3) = d2; Der(4) = d3;
178
179   Standard_Integer GOrdre = 4 + 4*Continuity, 
180                    DDim=Continuity*(Continuity+2);
181   math_Vector GaussP(1, GOrdre), GaussW(1, GOrdre), 
182               pol2(1, 2*Continuity+1), 
183               pol4(1, 4*Continuity+1);
184   math::GaussPoints(GOrdre, GaussP);
185   math::GaussWeights (GOrdre, GaussW);
186   pol4.Init(0.);
187
188   for (ip=1; ip<=GOrdre; ip++) {
189     t = (GaussP(ip)+1.)/2;
190     GW = GaussW(ip);
191     PLib::NoDerivativeEvalPolynomial(t ,  Continuity, Continuity+2, DDim,
192                                      polynome[0], valhder[0]);
193     V /= Length; //Normalisation   
194
195     //                      i
196     // C'(t) = SUM Vi*Lambda 
197     Vec1 = Constraint.Col(1);
198     Vec1 *= V(1);
199     Vec1 += V(size)*Constraint.Col(size);
200     Vec2 = Constraint.Col(2);
201     Vec2 *= V(2);
202     if (Continuity > 1) {
203       Vec3 = Constraint.Col(3);
204       Vec3 *= V(3);
205       if (Continuity > 2) {
206         Vec4 = Constraint.Col(4);
207         Vec4 *= V(4);  
208       }
209     }
210     
211     
212     //   2          2
213     // C'(t) - C'(0)
214
215     pol2(1) = Vec1.Norm2();
216     pol2(2) = 2*(Vec1.Multiplied(Vec2));
217     pol2(3) = Vec2.Norm2() - squared1;
218     if (Continuity>1) { 
219       pol2(3) += 2*(Vec1.Multiplied(Vec3));
220       pol2(4) =  2*(Vec2.Multiplied(Vec3));
221       pol2(5) =  Vec3.Norm2();
222       if (Continuity>2) {
223         pol2(4)+= 2*(Vec1.Multiplied(Vec4));
224         pol2(5)+= 2*(Vec2.Multiplied(Vec4));
225         pol2(6) = 2*(Vec3.Multiplied(Vec4));
226         pol2(7) = Vec4.Norm2();
227       }
228     }
229
230     //                     2      2  2
231     // Integrale de ( C'(t) - C'(0) )
232     for (ii=1; ii<=pol2.Length(); ii++) {
233       pp = ii;
234       for(jj=1; jj<ii; jj++, pp++) {
235         pol4(pp) += 2*GW*pol2(ii)*pol2(jj);
236       }
237       pol4(2*ii-1) += GW*Pow(pol2(ii), 2);
238     }
239   }
240
241   Standard_Real EMin, E;
242   PLib::NoDerivativeEvalPolynomial(Lambda , pol4.Length()-1, 1, 
243                                    pol4.Length()-1,
244                                    pol4(1), EMin); 
245
246   if (EMin > Precision::Confusion()) {
247     // Recheche des extrema de la fonction
248     GeomLib_PolyFunc FF(pol4);
249     GeomLib_LogSample S(Lambda/1000, 50*Lambda, 100);
250     math_FunctionAllRoots Solve(FF, S, Precision::Confusion(), 
251                                 Precision::Confusion()*(Length+1),
252                                 1.e-15);
253     if (Solve.IsDone()) {
254       for (ii=1; ii<=Solve.NbPoints(); ii++) {
255         t = Solve.GetPoint(ii);
256         PLib::NoDerivativeEvalPolynomial(t , pol4.Length()-1, 1, 
257                                          pol4.Length()-1,
258                                          pol4(1), E);
259         if (E < EMin) {
260           Lambda = t;
261           EMin = E;
262         }
263       }
264     }
265   }
266 }
267
268 #include <Extrema_LocateExtPC.hxx>
269 //=======================================================================
270 //function : RemovePointsFromArray
271 //purpose  : 
272 //=======================================================================
273
274 void GeomLib::RemovePointsFromArray(const Standard_Integer NumPoints,
275                                     const TColStd_Array1OfReal& InParameters,
276                                     Handle(TColStd_HArray1OfReal)& OutParameters) 
277 {
278  Standard_Integer ii,
279    jj,
280    add_one_point,
281    loc_num_points,
282    num_points,
283    index ;
284  Standard_Real delta,
285    current_parameter ;
286
287    loc_num_points = Max(0,NumPoints-2) ;
288    delta = InParameters(InParameters.Upper()) - InParameters(InParameters.Lower()) ;
289    delta /= (Standard_Real) (loc_num_points + 1) ;
290    num_points = 1 ;
291    current_parameter = InParameters(InParameters.Lower()) + delta * 0.5e0 ;
292    ii = InParameters.Lower() + 1 ;
293    for (jj = 0 ; ii < InParameters.Upper() && jj < NumPoints ; jj++) {
294      add_one_point = 0 ;
295      while ( ii < InParameters.Upper() && InParameters(ii) < current_parameter) {
296        ii += 1 ;
297        add_one_point = 1 ;
298      }
299      num_points += add_one_point ;
300      current_parameter += delta ;
301    }
302    if (NumPoints <= 2) {
303      num_points = 2 ;
304    }
305    index = 2 ;
306    current_parameter = InParameters(InParameters.Lower()) + delta * 0.5e0 ;
307    OutParameters = 
308      new TColStd_HArray1OfReal(1,num_points) ;
309    OutParameters->ChangeArray1()(1) = InParameters(InParameters.Lower()) ;
310    ii = InParameters.Lower() + 1 ;
311    for (jj = 0 ; ii < InParameters.Upper() && jj < NumPoints ; jj++) {
312      add_one_point = 0 ;
313      while (ii < InParameters.Upper() && InParameters(ii) < current_parameter) {
314        ii += 1 ;
315        add_one_point = 1 ;
316      }
317      if (add_one_point && index <= num_points) {
318        OutParameters->ChangeArray1()(index) = InParameters(ii-1) ;
319        index += 1 ;
320      }
321      current_parameter += delta ;
322    }
323    OutParameters->ChangeArray1()(num_points) = InParameters(InParameters.Upper()) ;
324 }
325 //=======================================================================
326 //function : DensifyArray1OfReal
327 //purpose  : 
328 //=======================================================================
329
330 void GeomLib::DensifyArray1OfReal(const Standard_Integer MinNumPoints,
331                                   const TColStd_Array1OfReal& InParameters,
332                                   Handle(TColStd_HArray1OfReal)& OutParameters) 
333 {
334  Standard_Integer ii,
335    in_order,
336    num_points,
337    num_parameters_to_add,
338    index ;
339  Standard_Real delta,
340    current_parameter ;
341
342  in_order = 1 ;
343  if (MinNumPoints > InParameters.Length()) {
344
345    //
346    // checks the paramaters are in increasing order
347    // 
348    for (ii = InParameters.Lower() ; ii < InParameters.Upper() ; ii++) {
349      if (InParameters(ii) > InParameters(ii+1)) {
350        in_order = 0 ;
351        break ;
352      }
353    }
354    if (in_order) {
355      num_parameters_to_add = MinNumPoints - InParameters.Length()  ;
356      delta = InParameters(InParameters.Upper()) - InParameters(InParameters.Lower()) ;
357      delta /= (Standard_Real) (num_parameters_to_add + 1) ;
358      num_points = MinNumPoints ;
359      OutParameters = 
360        new TColStd_HArray1OfReal(1,num_points) ;
361      index = 1 ;
362      current_parameter = InParameters(InParameters.Lower()) ;
363      OutParameters->ChangeArray1()(index) = current_parameter ;
364      index += 1 ;
365      current_parameter += delta ; 
366      for (ii = InParameters.Lower() + 1 ; index <= num_points && ii <= InParameters.Upper() ; ii++) {
367        while (current_parameter < InParameters(ii) && index <= num_points) {
368          OutParameters->ChangeArray1()(index) = current_parameter ;
369          index += 1 ;
370          current_parameter += delta ;
371        }
372        if (index <= num_points) { 
373          OutParameters->ChangeArray1()(index) = InParameters(ii) ;
374        }
375        index += 1 ;
376      }
377      //
378      // beware of roundoff !
379      //
380      OutParameters->ChangeArray1()(num_points) = InParameters(InParameters.Upper()) ;
381    }
382    else {
383      index = 1 ;
384      num_points = InParameters.Length() ;
385      OutParameters = 
386        new TColStd_HArray1OfReal(1,num_points) ;
387      for (ii = InParameters.Lower()  ; ii <= InParameters.Upper() ; ii++) {
388        OutParameters->ChangeArray1()(index) = InParameters(ii) ;
389        index += 1 ;
390      }
391    }
392  }
393  else {
394    index = 1 ;
395    num_points = InParameters.Length() ;
396    OutParameters = 
397      new TColStd_HArray1OfReal(1,num_points) ;
398    for (ii = InParameters.Lower()  ; ii <= InParameters.Upper() ; ii++) {
399      OutParameters->ChangeArray1()(index) = InParameters(ii) ;
400      index += 1 ;
401    }
402  }
403 }
404
405 //=======================================================================
406 //function : FuseIntervals
407 //purpose  : 
408 //=======================================================================
409 void GeomLib::FuseIntervals(const  TColStd_Array1OfReal& I1,
410                             const  TColStd_Array1OfReal& I2,
411                             TColStd_SequenceOfReal&  Seq,
412                             const Standard_Real  Epspar) 
413 {
414  Standard_Integer ind1=1, ind2=1;
415  Standard_Real    v1, v2;
416 // Initialisations : les IND1 et IND2 pointent sur le 1er element
417 // de chacune des 2 tables a traiter.INDS pointe sur le dernier
418 // element cree de TABSOR
419
420
421 //--- On remplit TABSOR en parcourant TABLE1 et TABLE2 simultanement ---
422 //------------------ en eliminant les occurrences multiples ------------
423
424  while ((ind1<=I1.Upper()) && (ind2<=I2.Upper())) {
425       v1 = I1(ind1);
426       v2 = I2(ind2);
427       if (Abs(v1-v2)<= Epspar) {
428 // Ici les elements de I1 et I2 conviennent .
429          Seq.Append((v1+v2)/2);
430          ind1++;
431          ind2++;
432        }
433       else if (v1 < v2) {
434         // Ici l' element de I1 convient.
435          Seq.Append(v1);
436          ind1++;
437        }
438       else {
439 // Ici l' element de TABLE2 convient.
440          Seq.Append(v2);
441          ind2++;
442        }
443     }
444
445   if (ind1>I1.Upper()) { 
446 //----- Ici I1 est epuise, on complete avec la fin de TABLE2 -------
447
448     for (; ind2<=I2.Upper(); ind2++) {
449       Seq.Append(I2(ind2));
450     }
451   }
452
453   if (ind2>I2.Upper()) { 
454 //----- Ici I2 est epuise, on complete avec la fin de I1 -------
455     for (; ind1<=I1.Upper(); ind1++) {
456       Seq.Append(I1(ind1));
457     }
458   } 
459 }
460
461
462 //=======================================================================
463 //function : EvalMaxParametricDistance
464 //purpose  : 
465 //=======================================================================
466
467 void GeomLib::EvalMaxParametricDistance(const Adaptor3d_Curve& ACurve,
468                                const Adaptor3d_Curve& AReferenceCurve,
469 //                             const Standard_Real  Tolerance,
470                                const Standard_Real  ,
471                                const TColStd_Array1OfReal& Parameters,
472                                Standard_Real& MaxDistance) 
473 {
474   Standard_Integer ii ;
475
476   Standard_Real max_squared = 0.0e0,
477 //    tolerance_squared,
478     local_distance_squared ;
479
480 //  tolerance_squared = Tolerance * Tolerance ;
481   gp_Pnt Point1 ;
482   gp_Pnt Point2 ;
483   for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
484     ACurve.D0(Parameters(ii),
485               Point1) ;
486     AReferenceCurve.D0(Parameters(ii),
487                        Point2) ;
488     local_distance_squared =
489       Point1.SquareDistance (Point2) ;
490     max_squared = Max(max_squared,local_distance_squared) ;
491   }
492   if (max_squared > 0.0e0) {
493     MaxDistance = sqrt(max_squared) ;
494   }
495   else {
496     MaxDistance = 0.0e0 ;
497   }
498   
499 }
500 //=======================================================================
501 //function : EvalMaxDistanceAlongParameter
502 //purpose  : 
503 //=======================================================================
504
505 void GeomLib::EvalMaxDistanceAlongParameter(const Adaptor3d_Curve& ACurve,
506                                const Adaptor3d_Curve& AReferenceCurve,
507                                const Standard_Real  Tolerance,
508                                const TColStd_Array1OfReal& Parameters,
509                                Standard_Real& MaxDistance) 
510 {
511   Standard_Integer ii ;
512   Standard_Real max_squared = 0.0e0,
513     tolerance_squared = Tolerance * Tolerance,
514     other_parameter,
515     para_tolerance,
516     local_distance_squared ;
517   gp_Pnt Point1 ;
518   gp_Pnt Point2 ;
519
520
521
522   para_tolerance = 
523     AReferenceCurve.Resolution(Tolerance) ;
524   other_parameter = Parameters(Parameters.Lower()) ;
525   ACurve.D0(other_parameter,
526             Point1) ;
527   Extrema_LocateExtPC a_projector(Point1,
528                                   AReferenceCurve,
529                                   other_parameter,
530                                   para_tolerance) ;
531   for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
532     ACurve.D0(Parameters(ii),
533               Point1) ;
534     AReferenceCurve.D0(Parameters(ii),
535                        Point2) ;
536     local_distance_squared =
537       Point1.SquareDistance (Point2) ;
538     
539     local_distance_squared =
540       Point1.SquareDistance (Point2) ;
541     
542     
543     if (local_distance_squared > tolerance_squared) {
544       
545       
546       a_projector.Perform(Point1,
547                           other_parameter) ;
548       if (a_projector.IsDone()) {
549         other_parameter =
550           a_projector.Point().Parameter() ;
551         AReferenceCurve.D0(other_parameter,
552                            Point2) ;
553         local_distance_squared =
554           Point1.SquareDistance (Point2) ;
555       }
556       else {
557         local_distance_squared = 0.0e0 ;
558         other_parameter = Parameters(ii) ;
559       }
560     }
561     else {
562       other_parameter = Parameters(ii) ;
563     }
564     
565     
566     max_squared = Max(max_squared,local_distance_squared) ;
567   }
568   if (max_squared > tolerance_squared) {
569     MaxDistance = sqrt(max_squared) ;
570   }
571   else {
572     MaxDistance = Tolerance ;
573   }
574 }
575
576
577
578 // Aliases:     
579
580 // Global data definitions:     
581
582 // Methods :
583
584
585 //=======================================================================
586 //function : To3d
587 //purpose  : 
588 //=======================================================================
589
590 Handle(Geom_Curve) GeomLib::To3d (const gp_Ax2&               Position,
591                                   const Handle(Geom2d_Curve)& Curve2d  ) {
592   Handle(Geom_Curve) Curve3d;
593   Handle(Standard_Type) KindOfCurve = Curve2d->DynamicType();
594
595   if (KindOfCurve == STANDARD_TYPE (Geom2d_TrimmedCurve)) {
596     Handle(Geom2d_TrimmedCurve) Ct =
597       Handle(Geom2d_TrimmedCurve)::DownCast(Curve2d);
598     Standard_Real U1 = Ct->FirstParameter ();
599     Standard_Real U2 = Ct->LastParameter  ();
600     Handle(Geom2d_Curve) CBasis2d = Ct->BasisCurve();
601     Handle(Geom_Curve) CC = GeomLib::To3d(Position, CBasis2d);
602     Curve3d = new Geom_TrimmedCurve (CC, U1, U2);
603   }
604   else if (KindOfCurve == STANDARD_TYPE (Geom2d_OffsetCurve)) {
605     Handle(Geom2d_OffsetCurve) Co =
606       Handle(Geom2d_OffsetCurve)::DownCast(Curve2d);
607     Standard_Real Offset = Co->Offset();
608     Handle(Geom2d_Curve) CBasis2d = Co->BasisCurve();
609     Handle(Geom_Curve) CC = GeomLib::To3d(Position, CBasis2d);
610     Curve3d = new Geom_OffsetCurve (CC, Offset, Position.Direction());
611   }
612   else if (KindOfCurve == STANDARD_TYPE (Geom2d_BezierCurve)) {
613     Handle(Geom2d_BezierCurve) CBez2d = 
614       Handle(Geom2d_BezierCurve)::DownCast (Curve2d);
615     Standard_Integer Nbpoles = CBez2d->NbPoles ();
616     TColgp_Array1OfPnt2d Poles2d (1, Nbpoles);
617     CBez2d->Poles (Poles2d);
618     TColgp_Array1OfPnt Poles3d (1, Nbpoles);
619     for (Standard_Integer i = 1; i <= Nbpoles; i++) {
620       Poles3d (i) = ElCLib::To3d (Position, Poles2d (i));
621     }
622     Handle(Geom_BezierCurve) CBez3d;
623     if (CBez2d->IsRational()) {
624       TColStd_Array1OfReal TheWeights (1, Nbpoles);
625       CBez2d->Weights (TheWeights);
626       CBez3d = new Geom_BezierCurve (Poles3d, TheWeights);
627     }
628     else {
629       CBez3d = new Geom_BezierCurve (Poles3d);
630     }
631     Curve3d = CBez3d;
632   }
633   else if (KindOfCurve == STANDARD_TYPE (Geom2d_BSplineCurve)) {
634     Handle(Geom2d_BSplineCurve) CBSpl2d =
635       Handle(Geom2d_BSplineCurve)::DownCast (Curve2d);
636     Standard_Integer Nbpoles   = CBSpl2d->NbPoles ();
637     Standard_Integer Nbknots   = CBSpl2d->NbKnots ();
638     Standard_Integer TheDegree = CBSpl2d->Degree ();
639     Standard_Boolean IsPeriodic = CBSpl2d->IsPeriodic();
640     TColgp_Array1OfPnt2d Poles2d (1, Nbpoles);
641     CBSpl2d->Poles (Poles2d);
642     TColgp_Array1OfPnt Poles3d (1, Nbpoles);
643     for (Standard_Integer i = 1; i <= Nbpoles; i++) {
644       Poles3d (i) = ElCLib::To3d (Position, Poles2d (i));
645     }
646     TColStd_Array1OfReal    TheKnots (1, Nbknots);
647     TColStd_Array1OfInteger TheMults (1, Nbknots);
648     CBSpl2d->Knots (TheKnots);
649     CBSpl2d->Multiplicities (TheMults);
650     Handle(Geom_BSplineCurve) CBSpl3d;
651     if (CBSpl2d->IsRational()) {
652       TColStd_Array1OfReal TheWeights (1, Nbpoles);
653       CBSpl2d->Weights (TheWeights);
654       CBSpl3d = new Geom_BSplineCurve (Poles3d, TheWeights, TheKnots, TheMults, TheDegree, IsPeriodic);
655     }
656     else {
657       CBSpl3d = new Geom_BSplineCurve (Poles3d, TheKnots, TheMults, TheDegree, IsPeriodic);
658     }
659     Curve3d = CBSpl3d;
660   }
661   else if (KindOfCurve == STANDARD_TYPE (Geom2d_Line)) {
662     Handle(Geom2d_Line) Line2d = Handle(Geom2d_Line)::DownCast (Curve2d);
663     gp_Lin2d L2d = Line2d->Lin2d();
664     gp_Lin   L3d = ElCLib::To3d (Position, L2d);
665     Handle(Geom_Line) GeomL3d = new Geom_Line (L3d);
666     Curve3d = GeomL3d;
667   }
668   else if (KindOfCurve == STANDARD_TYPE (Geom2d_Circle)) {
669     Handle(Geom2d_Circle) Circle2d = 
670       Handle(Geom2d_Circle)::DownCast (Curve2d);
671     gp_Circ2d C2d = Circle2d->Circ2d();
672     gp_Circ   C3d = ElCLib::To3d (Position, C2d);
673     Handle(Geom_Circle) GeomC3d = new Geom_Circle (C3d);
674     Curve3d = GeomC3d;
675   }
676   else if (KindOfCurve == STANDARD_TYPE (Geom2d_Ellipse)) {
677     Handle(Geom2d_Ellipse) Ellipse2d =
678       Handle(Geom2d_Ellipse)::DownCast (Curve2d);
679     gp_Elips2d E2d = Ellipse2d->Elips2d ();
680     gp_Elips   E3d = ElCLib::To3d (Position, E2d);
681     Handle(Geom_Ellipse) GeomE3d = new Geom_Ellipse (E3d);
682     Curve3d = GeomE3d;
683   }
684   else if (KindOfCurve == STANDARD_TYPE (Geom2d_Parabola)) {
685     Handle(Geom2d_Parabola) Parabola2d =
686       Handle(Geom2d_Parabola)::DownCast (Curve2d);
687     gp_Parab2d Prb2d = Parabola2d->Parab2d ();
688     gp_Parab   Prb3d = ElCLib::To3d (Position, Prb2d);
689     Handle(Geom_Parabola) GeomPrb3d = new Geom_Parabola (Prb3d);
690     Curve3d = GeomPrb3d;
691   }
692   else if (KindOfCurve == STANDARD_TYPE (Geom2d_Hyperbola)) {
693     Handle(Geom2d_Hyperbola) Hyperbola2d =
694       Handle(Geom2d_Hyperbola)::DownCast (Curve2d);
695     gp_Hypr2d H2d = Hyperbola2d->Hypr2d ();
696     gp_Hypr   H3d = ElCLib::To3d (Position, H2d);
697     Handle(Geom_Hyperbola) GeomH3d = new Geom_Hyperbola (H3d);
698     Curve3d = GeomH3d;
699   }
700   else {
701     Standard_NotImplemented::Raise();
702   }
703   
704   return Curve3d;
705 }
706
707
708
709 //=======================================================================
710 //function : GTransform
711 //purpose  : 
712 //=======================================================================
713
714 Handle(Geom2d_Curve) GeomLib::GTransform(const Handle(Geom2d_Curve)& Curve, 
715                                          const gp_GTrsf2d&           GTrsf)
716 {
717   gp_TrsfForm Form = GTrsf.Form();
718   
719   if ( Form != gp_Other) {
720     
721     // Alors, la GTrsf est en fait une Trsf. 
722     // La geometrie des courbes sera alors inchangee.
723
724     Handle(Geom2d_Curve) C = 
725       Handle(Geom2d_Curve)::DownCast(Curve->Transformed(GTrsf.Trsf2d()));
726     return C;
727   }
728   else { 
729     
730     // Alors, la GTrsf est une other Transformation.
731     // La geometrie des courbes est alors changee, et les conics devront
732     // etre converties en BSplines.
733     
734     Handle(Standard_Type) TheType = Curve->DynamicType();
735     
736     if ( TheType == STANDARD_TYPE(Geom2d_TrimmedCurve)) {
737       
738       // On va recurer sur la BasisCurve
739       
740       Handle(Geom2d_TrimmedCurve) C = 
741         Handle(Geom2d_TrimmedCurve)::DownCast(Curve->Copy());
742       
743       Handle(Standard_Type) TheBasisType = (C->BasisCurve())->DynamicType();
744       
745       if (TheBasisType == STANDARD_TYPE(Geom2d_BSplineCurve) ||
746           TheBasisType == STANDARD_TYPE(Geom2d_BezierCurve)    ) {
747         
748         // Dans ces cas le parametrage est conserve sur la courbe transformee
749         // on peut donc la trimmer avec les parametres de la courbe de base.
750         
751         Standard_Real U1 = C->FirstParameter();
752         Standard_Real U2 = C->LastParameter();
753         
754         Handle(Geom2d_TrimmedCurve) result = 
755           new Geom2d_TrimmedCurve(GTransform(C->BasisCurve(), GTrsf), U1,U2);
756         return result;
757       }
758       else if ( TheBasisType == STANDARD_TYPE(Geom2d_Line)) {
759         
760         // Dans ce cas, le parametrage n`est plus conserve.
761         // Il faut recalculer les parametres de Trimming sur la courbe 
762         // resultante. ( Calcul par projection ( ElCLib) des points debut 
763         // et fin transformes)
764         
765         Handle(Geom2d_Line) L = 
766           Handle(Geom2d_Line)::DownCast(GTransform(C->BasisCurve(), GTrsf));
767         gp_Lin2d Lin = L->Lin2d();
768         
769         gp_Pnt2d P1 = C->StartPoint();
770         gp_Pnt2d P2 = C->EndPoint();
771         P1.SetXY(GTrsf.Transformed(P1.XY()));
772         P2.SetXY(GTrsf.Transformed(P2.XY()));
773         Standard_Real U1 = ElCLib::Parameter(Lin,P1);
774         Standard_Real U2 = ElCLib::Parameter(Lin,P2);
775         
776         Handle(Geom2d_TrimmedCurve) result = 
777           new Geom2d_TrimmedCurve(L,U1,U2);
778         return result;
779       }
780       else if (TheBasisType == STANDARD_TYPE(Geom2d_Circle)   ||
781                TheBasisType == STANDARD_TYPE(Geom2d_Ellipse)  ||
782                TheBasisType == STANDARD_TYPE(Geom2d_Parabola) ||
783                TheBasisType == STANDARD_TYPE(Geom2d_Hyperbola)  ) {
784         
785         // Dans ces cas, la geometrie de la courbe n`est pas conservee
786         // on la convertir en BSpline avant de lui appliquer la Trsf.
787         
788         Handle(Geom2d_BSplineCurve) BS = 
789           Geom2dConvert::CurveToBSplineCurve(C);
790         return GTransform(BS,GTrsf);
791       }
792       else {
793         
794         // La transformee d`une OffsetCurve vaut ????? Sais pas faire !! 
795         
796         Handle(Geom2d_Curve) dummy;
797         return dummy;
798       }
799     }
800     else if ( TheType == STANDARD_TYPE(Geom2d_Line)) {
801       
802       Handle(Geom2d_Line) L = 
803         Handle(Geom2d_Line)::DownCast(Curve->Copy());
804       gp_Lin2d Lin = L->Lin2d();
805       gp_Pnt2d P  = Lin.Location();
806       gp_Pnt2d PP = L->Value(10.); // pourquoi pas !!
807       P.SetXY(GTrsf.Transformed(P.XY()));
808       PP.SetXY(GTrsf.Transformed(PP.XY()));
809       L->SetLocation(P);
810       gp_Vec2d V(P,PP);
811       L->SetDirection(gp_Dir2d(V));
812       return L;
813     }
814     else if ( TheType == STANDARD_TYPE(Geom2d_BezierCurve)) {
815       
816       // Les GTrsf etant des operation lineaires, la transformee d`une courbe
817       // a poles est la courbe dont les poles sont la transformee des poles
818       // de la courbe de base.
819       
820       Handle(Geom2d_BezierCurve) C = 
821         Handle(Geom2d_BezierCurve)::DownCast(Curve->Copy());
822       Standard_Integer NbPoles = C->NbPoles();
823       TColgp_Array1OfPnt2d Poles(1,NbPoles);
824       C->Poles(Poles);
825       for ( Standard_Integer i = 1; i <= NbPoles; i++) {
826         Poles(i).SetXY(GTrsf.Transformed(Poles(i).XY()));
827         C->SetPole(i,Poles(i));
828       }
829       return C;
830     }
831     else if ( TheType == STANDARD_TYPE(Geom2d_BSplineCurve)) {
832       
833       // Voir commentaire pour les Bezier.
834       
835       Handle(Geom2d_BSplineCurve) C = 
836         Handle(Geom2d_BSplineCurve)::DownCast(Curve->Copy());
837       Standard_Integer NbPoles = C->NbPoles();
838       TColgp_Array1OfPnt2d Poles(1,NbPoles);
839       C->Poles(Poles);
840       for ( Standard_Integer i = 1; i <= NbPoles; i++) {
841         Poles(i).SetXY(GTrsf.Transformed(Poles(i).XY()));
842         C->SetPole(i,Poles(i));
843       }
844       return C;
845     }
846     else if ( TheType == STANDARD_TYPE(Geom2d_Circle) ||
847               TheType == STANDARD_TYPE(Geom2d_Ellipse)  ) {
848       
849       // Dans ces cas, la geometrie de la courbe n`est pas conservee
850       // on la convertir en BSpline avant de lui appliquer la Trsf.
851       
852       Handle(Geom2d_BSplineCurve) C = 
853         Geom2dConvert::CurveToBSplineCurve(Curve);
854       return GTransform(C, GTrsf);
855     }
856     else if ( TheType == STANDARD_TYPE(Geom2d_Parabola)   ||
857               TheType == STANDARD_TYPE(Geom2d_Hyperbola)  ||
858               TheType == STANDARD_TYPE(Geom2d_OffsetCurve)  ) {
859       
860       // On ne sait pas faire : return a null Handle;
861       
862       Handle(Geom2d_Curve) dummy;
863       return dummy;
864     }
865   }
866
867   Handle(Geom2d_Curve) WNT__; // portage Windows.
868   return WNT__;
869 }
870
871
872 //=======================================================================
873 //function : SameRange
874 //purpose  : 
875 //=======================================================================
876 void GeomLib::SameRange(const Standard_Real         Tolerance,
877                         const Handle(Geom2d_Curve)& CurvePtr,
878                         const Standard_Real         FirstOnCurve,
879                         const Standard_Real         LastOnCurve,
880                         const Standard_Real         RequestedFirst,
881                         const Standard_Real         RequestedLast,
882                               Handle(Geom2d_Curve)& NewCurvePtr) 
883 {
884   if(CurvePtr.IsNull()) Standard_Failure::Raise();
885   if (Abs(LastOnCurve - RequestedLast) <= Tolerance &&
886       Abs(FirstOnCurve - RequestedFirst) <= Tolerance) { 
887     NewCurvePtr = CurvePtr;
888     return;
889   }
890
891   // the parametrisation lentgh  must at least be the same.
892   if (Abs(LastOnCurve - FirstOnCurve - RequestedLast + RequestedFirst) 
893       <= Tolerance) { 
894     if (CurvePtr->IsKind(STANDARD_TYPE(Geom2d_Line))) {
895       Handle(Geom2d_Line) Line =
896         Handle(Geom2d_Line)::DownCast(CurvePtr->Copy());
897       Standard_Real dU = FirstOnCurve - RequestedFirst;
898       gp_Dir2d D = Line->Direction() ;
899       Line->Translate(dU * gp_Vec2d(D));
900       NewCurvePtr = Line;
901     }
902     else if (CurvePtr->IsKind(STANDARD_TYPE(Geom2d_Circle))) {
903       gp_Trsf2d Trsf;
904       NewCurvePtr = Handle(Geom2d_Curve)::DownCast(CurvePtr->Copy()); 
905       Handle(Geom2d_Circle) Circ = 
906         Handle(Geom2d_Circle)::DownCast(NewCurvePtr);
907       gp_Pnt2d P = Circ->Location();
908       Standard_Real dU;
909       if (Circ->Circ2d().IsDirect()) {
910         dU = FirstOnCurve - RequestedFirst;
911       }
912       else {
913         dU = RequestedFirst - FirstOnCurve;
914       }
915       Trsf.SetRotation(P,dU);
916       NewCurvePtr->Transform(Trsf) ;
917     }
918     else if (CurvePtr->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve))) {
919       Handle(Geom2d_TrimmedCurve) TC = 
920         Handle(Geom2d_TrimmedCurve)::DownCast(CurvePtr);
921       GeomLib::SameRange(Tolerance,
922                          TC->BasisCurve(),
923                          FirstOnCurve  , LastOnCurve,
924                          RequestedFirst, RequestedLast,
925                          NewCurvePtr);
926       NewCurvePtr = new Geom2d_TrimmedCurve( NewCurvePtr, RequestedFirst, RequestedLast );
927     }
928 //
929 //  attention a des problemes de limitation : utiliser le MEME test que dans
930 //  Geom2d_TrimmedCurve::SetTrim car sinon comme on risque de relimite sur 
931 //  RequestedFirst et RequestedLast on aura un probleme
932 //
933 // 
934     else if (Abs(LastOnCurve - FirstOnCurve) > Precision::PConfusion() ||
935              Abs(RequestedLast + RequestedFirst) > Precision::PConfusion()) {
936       
937       Handle(Geom2d_TrimmedCurve) TC =
938         new Geom2d_TrimmedCurve(CurvePtr,FirstOnCurve,LastOnCurve);
939       
940       Handle(Geom2d_BSplineCurve) BS =
941         Geom2dConvert::CurveToBSplineCurve(TC);
942       TColStd_Array1OfReal Knots(1,BS->NbKnots());
943       BS->Knots(Knots);
944       
945       BSplCLib::Reparametrize(RequestedFirst,RequestedLast,Knots);
946       
947       BS->SetKnots(Knots);
948       NewCurvePtr = BS;
949     }
950   
951   }
952   else { // On segmente le resultat
953     Handle(Geom2d_TrimmedCurve) TC =
954       new Geom2d_TrimmedCurve( CurvePtr, FirstOnCurve, LastOnCurve );
955
956     Standard_Real newFirstOnCurve = TC->FirstParameter(), newLastOnCurve = TC->LastParameter();
957     
958     Handle(Geom2d_BSplineCurve) BS =
959       Geom2dConvert::CurveToBSplineCurve(TC);
960
961     if (BS->IsPeriodic()) 
962       BS->Segment( newFirstOnCurve, newLastOnCurve) ;
963     else 
964       BS->Segment( Max(newFirstOnCurve, BS->FirstParameter()),
965                    Min(newLastOnCurve,  BS->LastParameter()) );
966
967     TColStd_Array1OfReal Knots(1,BS->NbKnots());
968     BS->Knots(Knots);
969     
970     BSplCLib::Reparametrize(RequestedFirst,RequestedLast,Knots);
971     
972     BS->SetKnots(Knots);
973     NewCurvePtr = BS;
974   }
975 }
976
977 //=======================================================================
978 //class : GeomLib_CurveOnSurfaceEvaluator
979 //purpose: The evaluator for the Curve 3D building
980 //=======================================================================
981
982 class GeomLib_CurveOnSurfaceEvaluator : public AdvApprox_EvaluatorFunction
983 {
984  public:
985   GeomLib_CurveOnSurfaceEvaluator (Adaptor3d_CurveOnSurface& theCurveOnSurface,
986                                    Standard_Real theFirst, Standard_Real theLast)
987     : CurveOnSurface(theCurveOnSurface), FirstParam(theFirst), LastParam(theLast) {}
988   
989   virtual void Evaluate (Standard_Integer *Dimension,
990                          Standard_Real     StartEnd[2],
991                          Standard_Real    *Parameter,
992                          Standard_Integer *DerivativeRequest,
993                          Standard_Real    *Result, // [Dimension]
994                          Standard_Integer *ErrorCode);
995   
996  private:
997   Adaptor3d_CurveOnSurface& CurveOnSurface;
998   Standard_Real FirstParam;
999   Standard_Real LastParam; 
1000
1001   Handle(Adaptor3d_HCurve) TrimCurve;
1002 };
1003
1004 void GeomLib_CurveOnSurfaceEvaluator::Evaluate (Standard_Integer *,/*Dimension*/
1005                                                 Standard_Real     DebutFin[2],
1006                                                 Standard_Real    *Parameter,
1007                                                 Standard_Integer *DerivativeRequest,
1008                                                 Standard_Real    *Result,// [Dimension]
1009                                                 Standard_Integer *ReturnCode)
1010
1011   register Standard_Integer ii ;
1012   gp_Pnt Point ;
1013
1014   //Gestion des positionnements gauche / droite
1015   if ((DebutFin[0] != FirstParam) || (DebutFin[1] != LastParam)) 
1016     { 
1017       TrimCurve = CurveOnSurface.Trim(DebutFin[0], DebutFin[1], Precision::PConfusion());
1018       FirstParam = DebutFin[0];
1019       LastParam  = DebutFin[1];
1020     }
1021
1022   //Positionemment
1023   if (*DerivativeRequest == 0)
1024     {
1025      TrimCurve->D0((*Parameter), Point) ;
1026    
1027      for (ii = 0 ; ii < 3 ; ii++)
1028        Result[ii] = Point.Coord(ii + 1);
1029    }
1030   if (*DerivativeRequest == 1) 
1031     {
1032       gp_Vec Vector;
1033       TrimCurve->D1((*Parameter), Point, Vector);
1034       for (ii = 0 ; ii < 3 ; ii++)
1035         Result[ii] = Vector.Coord(ii + 1) ;
1036     }
1037   if (*DerivativeRequest == 2) 
1038     {
1039       gp_Vec Vector, VecBis;
1040       TrimCurve->D2((*Parameter), Point, VecBis, Vector);
1041       for (ii = 0 ; ii < 3 ; ii++)
1042         Result[ii] = Vector.Coord(ii + 1) ;
1043     }
1044   ReturnCode[0] = 0;
1045 }
1046
1047 //=======================================================================
1048 //function : BuildCurve3d
1049 //purpose  : 
1050 //=======================================================================
1051
1052 void GeomLib::BuildCurve3d(const Standard_Real           Tolerance,
1053                            Adaptor3d_CurveOnSurface&       Curve, 
1054                            const Standard_Real           FirstParameter,
1055                            const Standard_Real           LastParameter,
1056                            Handle(Geom_Curve)&            NewCurvePtr, 
1057                            Standard_Real&                MaxDeviation,
1058                            Standard_Real&                AverageDeviation,
1059                            const GeomAbs_Shape           Continuity,
1060                            const Standard_Integer        MaxDegree,
1061                            const Standard_Integer        MaxSegment) 
1062
1063 {
1064    
1065
1066   Standard_Integer curve_not_computed = 1 ;
1067   MaxDeviation     = 0.0e0 ;
1068   AverageDeviation = 0.0e0 ;
1069   const Handle(GeomAdaptor_HSurface) &     geom_adaptor_surface_ptr =
1070   Handle(GeomAdaptor_HSurface)::DownCast(Curve.GetSurface()) ;
1071   const Handle(Geom2dAdaptor_HCurve) &     geom_adaptor_curve_ptr =
1072   Handle(Geom2dAdaptor_HCurve)::DownCast(Curve.GetCurve()) ;
1073    
1074   if (! geom_adaptor_curve_ptr.IsNull() &&
1075       ! geom_adaptor_surface_ptr.IsNull()) {
1076      Handle(Geom_Plane) P ;
1077      const GeomAdaptor_Surface  &   geom_surface =
1078        * (GeomAdaptor_Surface *) &geom_adaptor_surface_ptr->Surface() ;
1079
1080     Handle(Geom_RectangularTrimmedSurface) RT = 
1081       Handle(Geom_RectangularTrimmedSurface)::
1082         DownCast(geom_surface.Surface());
1083     if ( RT.IsNull()) {
1084       P = Handle(Geom_Plane)::DownCast(geom_surface.Surface());
1085     }
1086     else {
1087       P = Handle(Geom_Plane)::DownCast(RT->BasisSurface());
1088     }
1089
1090    
1091     if (! P.IsNull()) {
1092       // compute the 3d curve
1093       gp_Ax2 axes = P->Position().Ax2();
1094       const Geom2dAdaptor_Curve & geom2d_curve =
1095         * (Geom2dAdaptor_Curve *) & geom_adaptor_curve_ptr->Curve2d() ;
1096       NewCurvePtr = 
1097         GeomLib::To3d(axes,
1098                       geom2d_curve.Curve());
1099      curve_not_computed = 0 ;
1100       
1101     }
1102   }
1103   if (curve_not_computed) {
1104
1105       //
1106       // Entree
1107       //
1108     Handle(TColStd_HArray1OfReal)   Tolerance1DPtr,Tolerance2DPtr; 
1109     Handle(TColStd_HArray1OfReal) Tolerance3DPtr =
1110       new TColStd_HArray1OfReal(1,1) ;
1111     Tolerance3DPtr->SetValue(1,Tolerance);
1112
1113      // Recherche des discontinuitees
1114      Standard_Integer NbIntervalC2 = Curve.NbIntervals(GeomAbs_C2);
1115      TColStd_Array1OfReal Param_de_decoupeC2 (1, NbIntervalC2+1);
1116      Curve.Intervals(Param_de_decoupeC2, GeomAbs_C2);
1117      
1118      Standard_Integer NbIntervalC3 = Curve.NbIntervals(GeomAbs_C3);
1119      TColStd_Array1OfReal Param_de_decoupeC3 (1, NbIntervalC3+1);
1120      Curve.Intervals(Param_de_decoupeC3, GeomAbs_C3);
1121
1122      // Note extension of the parameteric range  
1123      // Pour forcer le Trim au premier appel de l'evaluateur
1124      GeomLib_CurveOnSurfaceEvaluator ev (Curve, FirstParameter - 1., LastParameter  + 1.);
1125                                          
1126      // Approximation avec decoupe preferentiel
1127      AdvApprox_PrefAndRec Preferentiel(Param_de_decoupeC2,
1128                                        Param_de_decoupeC3);
1129      AdvApprox_ApproxAFunction  anApproximator(0,
1130                                               0,
1131                                               1,
1132                                               Tolerance1DPtr,
1133                                               Tolerance2DPtr,
1134                                               Tolerance3DPtr,
1135                                               FirstParameter,
1136                                               LastParameter,
1137                                               Continuity,
1138                                               MaxDegree,  
1139                                               MaxSegment,
1140                                               ev,
1141 //                                            CurveOnSurfaceEvaluator,
1142                                               Preferentiel) ;
1143     
1144     if (anApproximator.HasResult()) {
1145       GeomLib_MakeCurvefromApprox 
1146         aCurveBuilder(anApproximator) ;    
1147
1148       Handle(Geom_BSplineCurve) aCurvePtr = 
1149         aCurveBuilder.Curve(1) ;
1150       // On rend les resultats de l'approx
1151       MaxDeviation = anApproximator.MaxError(3,1) ;
1152       AverageDeviation = anApproximator.AverageError(3,1) ;
1153       NewCurvePtr = aCurvePtr ;
1154     }
1155   }  
1156  }
1157
1158 //=======================================================================
1159 //function :  AdjustExtremity
1160 //purpose  : 
1161 //=======================================================================
1162
1163 void GeomLib::AdjustExtremity(Handle(Geom_BoundedCurve)& Curve, 
1164                               const gp_Pnt& P1,
1165                               const gp_Pnt& P2,
1166                               const gp_Vec& T1,
1167                               const gp_Vec& T2)
1168 {
1169 // il faut Convertir l'entree (en preservant si possible le parametrage)
1170   Handle(Geom_BSplineCurve) aIn, aDef;  
1171   aIn = GeomConvert::CurveToBSplineCurve(Curve, Convert_QuasiAngular);
1172
1173   Standard_Integer ii, jj;
1174   gp_Pnt P;
1175   gp_Vec V, Vtan, DV;
1176   TColgp_Array1OfPnt PolesDef(1,4), Coeffs(1,4);
1177   TColStd_Array1OfReal FK(1, 8);
1178   TColStd_Array1OfReal Ti(1, 4);
1179   TColStd_Array1OfInteger Contact(1, 4);
1180
1181   Ti(1) = Ti(2) = aIn->FirstParameter();
1182   Ti(3) = Ti(4) = aIn->LastParameter();
1183   Contact(1) =  Contact(3) = 0;
1184   Contact(2) =  Contact(4) = 1;
1185   for (ii=1; ii<=4; ii++) {
1186     FK(ii) = aIn->FirstParameter();
1187     FK(ii) = aIn->LastParameter();
1188   }
1189
1190   // Calculs des contraintes de deformations
1191   aIn->D1(Ti(1), P, V);
1192   PolesDef(1).ChangeCoord() = P1.XYZ()-P.XYZ();
1193   Vtan = T1;
1194   Vtan.Normalize();
1195   DV = Vtan * (Vtan * V) - V;
1196   PolesDef(2).ChangeCoord() = (Ti(4)-Ti(1))*DV.XYZ();
1197
1198   aIn->D1(Ti(4), P, V);
1199   PolesDef(3).ChangeCoord() = P2.XYZ()-P.XYZ();
1200   Vtan = T2;
1201   Vtan.Normalize();
1202   DV = Vtan * (Vtan * V) - V;
1203   PolesDef(4).ChangeCoord() = (Ti(4)-Ti(1))* DV.XYZ();
1204  
1205   // Interpolation des contraintes
1206   math_Matrix Mat(1, 4, 1, 4);
1207   if (!PLib::HermiteCoefficients(0., 1., 1, 1, Mat)) 
1208     Standard_ConstructionError::Raise();
1209
1210   for (jj=1; jj<=4; jj++) {
1211     gp_XYZ aux(0.,0.,0.);
1212     for (ii=1; ii<=4; ii++) {
1213       aux.SetLinearForm(Mat(ii,jj), PolesDef(ii).XYZ(), aux);
1214     }
1215     Coeffs(jj).SetXYZ(aux);
1216   }
1217
1218   PLib::CoefficientsPoles(Coeffs, PLib::NoWeights(),
1219                           PolesDef,  PLib::NoWeights());
1220
1221   // Ajout de la deformation
1222   TColStd_Array1OfReal K(1, 2);
1223   TColStd_Array1OfInteger M(1, 2);
1224   K(1) = Ti(1);
1225   K(2) = Ti(4);
1226   M.Init(4);
1227
1228   aDef = new (Geom_BSplineCurve) (PolesDef, K, M, 3);
1229   if (aIn->Degree() < 3) aIn->IncreaseDegree(3);
1230   else aDef->IncreaseDegree(aIn->Degree());
1231
1232   for (ii=2; ii<aIn->NbKnots(); ii++) {
1233     aDef->InsertKnot(aIn->Knot(ii), aIn->Multiplicity(ii));
1234   }
1235
1236   if (aDef->NbPoles() != aIn->NbPoles()) 
1237     Standard_ConstructionError::Raise("Inconsistent poles's number");
1238
1239   for (ii=1; ii<=aDef->NbPoles(); ii++) {
1240     P = aIn->Pole(ii);
1241     P.ChangeCoord() += aDef->Pole(ii).XYZ();
1242     aIn->SetPole(ii, P);
1243   }
1244   Curve = aIn;
1245 }
1246 //=======================================================================
1247 //function : ExtendCurveToPoint
1248 //purpose  : 
1249 //=======================================================================
1250
1251 void GeomLib::ExtendCurveToPoint(Handle(Geom_BoundedCurve)& Curve, 
1252                                  const gp_Pnt& Point,
1253                                  const Standard_Integer Continuity,
1254                                  const Standard_Boolean After)
1255 {
1256   if(Continuity < 1 || Continuity > 3) return;
1257   Standard_Integer size = Continuity + 2;
1258   Standard_Real Ubord, Tol=1.e-6;
1259   math_Matrix  MatCoefs(1,size, 1,size);
1260   Standard_Real Lambda, L1;
1261   Standard_Integer ii, jj;
1262   gp_Vec d1, d2, d3;
1263   gp_Pnt p0;
1264 // il faut Convertir l'entree (en preservant si possible le parametrage)
1265   GeomConvert_CompCurveToBSplineCurve Concat(Curve, Convert_QuasiAngular);
1266
1267 // Les contraintes de constructions
1268   TColgp_Array1OfXYZ Cont(1,size);
1269   if (After) {
1270      Ubord = Curve->LastParameter();
1271     
1272    }
1273   else {
1274      Ubord = Curve->FirstParameter(); 
1275    }
1276   PLib::HermiteCoefficients(0, 1,           // Les Bornes
1277                             Continuity, 0,  // Les Ordres de contraintes
1278                             MatCoefs);
1279
1280   Curve->D3(Ubord, p0, d1, d2, d3);
1281   if (!After) { // Inversion du parametrage
1282     d1 *= -1;
1283     d3 *= -1;
1284   }
1285   
1286   L1 = p0.Distance(Point);
1287   if (L1 > Tol) {
1288     // Lambda est le ratio qu'il faut appliquer a la derive de la courbe
1289     // pour obtenir la derive du prolongement (fixe arbitrairement a la
1290     // longueur du segment bout de la courbe - point cible.
1291     // On essai d'avoir sur le prolongement la vitesse moyenne que l'on
1292     // a sur la courbe.
1293     gp_Vec daux;
1294     gp_Pnt pp;
1295     Standard_Real f= Curve->FirstParameter(), t, dt, norm; 
1296     dt = (Curve->LastParameter()-f)/9;
1297     norm = d1.Magnitude();
1298     for (ii=1, t=f+dt; ii<=8; ii++, t+=dt) {
1299       Curve->D1(t, pp, daux);
1300       norm += daux.Magnitude();
1301     }
1302     norm /= 9;
1303     dt = d1.Magnitude() / norm;
1304     if ((dt<1.5) && (dt>0.75)) { // Le bord est dans la moyenne on le garde
1305       Lambda = ((Standard_Real)1) / Max (d1.Magnitude() / L1, Tol);
1306     }
1307     else {
1308       Lambda = ((Standard_Real)1) / Max (norm / L1, Tol);
1309     }
1310   }
1311   else {
1312     return; // Pas d'extension
1313   }
1314
1315   // Optimisation du Lambda
1316   math_Matrix Cons(1, 3, 1, size);
1317   Cons(1,1) = p0.X();  Cons(2,1) = p0.Y(); Cons(3,1) = p0.Z();
1318   Cons(1,2) = d1.X();  Cons(2,2) = d1.Y(); Cons(3,2) = d1.Z();
1319   Cons(1,size) = Point.X();  Cons(2,size) = Point.Y(); Cons(3,size) = Point.Z();
1320   if (Continuity >= 2) {
1321      Cons(1,3) = d2.X();  Cons(2,3) = d2.Y(); Cons(3,3) = d2.Z(); 
1322   }
1323   if (Continuity >= 3) {
1324      Cons(1,4) = d3.X();  Cons(2,4) = d3.Y(); Cons(3,4) = d3.Z(); 
1325   }
1326   ComputeLambda(Cons, MatCoefs, L1, Lambda);
1327
1328   // Construction dans la Base Polynomiale
1329   Cont(1) = p0.XYZ();
1330   Cont(2) = d1.XYZ() * Lambda;
1331   if(Continuity >= 2) Cont(3) = d2.XYZ() * Pow(Lambda,2);
1332   if(Continuity >= 3) Cont(4) = d3.XYZ() * Pow(Lambda,3);
1333   Cont(size) = Point.XYZ();
1334     
1335
1336   TColgp_Array1OfPnt ExtrapPoles(1, size);
1337   TColgp_Array1OfPnt ExtraCoeffs(1, size);
1338
1339   gp_Pnt PNull(0.,0.,0.);
1340   ExtraCoeffs.Init(PNull);
1341   for (ii=1; ii<=size; ii++) {
1342     for (jj=1; jj<=size; jj++) {
1343       ExtraCoeffs(jj).ChangeCoord() += MatCoefs(ii,jj)*Cont(ii);
1344     }
1345   }
1346
1347   // Convertion Dans la Base de Bernstein
1348   PLib::CoefficientsPoles(ExtraCoeffs,  PLib::NoWeights(),
1349                           ExtrapPoles,  PLib::NoWeights());
1350   
1351   Handle(Geom_BezierCurve) Bezier = new (Geom_BezierCurve) (ExtrapPoles);
1352
1353   Standard_Real dist = ExtrapPoles(1).Distance(p0);
1354   Standard_Boolean Ok;
1355   Tol += dist;
1356
1357   // Concatenation
1358   Ok = Concat.Add(Bezier, Tol, After);
1359   if (!Ok) Standard_ConstructionError::Raise("ExtendCurveToPoint");
1360   
1361   Curve =  Concat.BSplineCurve();
1362 }
1363
1364
1365 //=======================================================================
1366 //function : ExtendKPart
1367 //purpose  : Extension par longueur des surfaces cannonique
1368 //=======================================================================
1369 static Standard_Boolean 
1370 ExtendKPart(Handle(Geom_RectangularTrimmedSurface)& Surface, 
1371             const Standard_Real Length,
1372             const Standard_Boolean InU,
1373             const Standard_Boolean After)
1374 {
1375
1376   if  (Surface.IsNull()) return Standard_False;
1377
1378   Standard_Boolean Ok=Standard_True;
1379   Standard_Real Uf, Ul, Vf, Vl;
1380   Handle(Geom_Surface) Support = Surface->BasisSurface();
1381   GeomAbs_SurfaceType Type;
1382
1383   Surface->Bounds(Uf, Ul, Vf, Vl);
1384   GeomAdaptor_Surface AS(Surface);
1385   Type = AS.GetType();
1386
1387   if (InU) {
1388     switch(Type) {
1389     case GeomAbs_Plane :
1390       {
1391         if (After) Ul+=Length;
1392         else       Uf-=Length;
1393         Surface = new (Geom_RectangularTrimmedSurface)
1394           (Support, Uf, Ul, Vf, Vl);
1395         break;
1396       }
1397
1398     default:
1399       Ok = Standard_False;
1400     }
1401   }
1402   else {
1403     switch(Type) {
1404     case GeomAbs_Plane :
1405     case GeomAbs_Cylinder :
1406     case GeomAbs_SurfaceOfExtrusion :
1407       {
1408         if (After) Vl+=Length;
1409         else       Vf-=Length;
1410         Surface = new (Geom_RectangularTrimmedSurface)
1411           (Support, Uf, Ul, Vf, Vl);
1412         break;
1413       }    
1414     default:
1415       Ok = Standard_False;
1416     }
1417   }
1418
1419   return Ok;
1420 }
1421
1422 //=======================================================================
1423 //function : ExtendSurfByLength
1424 //purpose  : 
1425 //=======================================================================
1426 void GeomLib::ExtendSurfByLength(Handle(Geom_BoundedSurface)& Surface, 
1427                                  const Standard_Real Length,
1428                                  const Standard_Integer Continuity,
1429                                  const Standard_Boolean InU,
1430                                  const Standard_Boolean After)
1431 {
1432   if(Continuity < 0 || Continuity > 3) return;
1433   Standard_Integer Cont = Continuity;
1434
1435   // Kpart ?
1436   Handle(Geom_RectangularTrimmedSurface) TS = 
1437     Handle(Geom_RectangularTrimmedSurface)::DownCast (Surface);
1438   if (ExtendKPart(TS,Length, InU, After) ) {
1439     Surface = TS;
1440     return;
1441   }
1442
1443 //  format BSplineSurface avec un degre suffisant pour la continuite voulue
1444   Handle(Geom_BSplineSurface) BS = 
1445     Handle(Geom_BSplineSurface)::DownCast (Surface);
1446   if (BS.IsNull()) {
1447     //BS = GeomConvert::SurfaceToBSplineSurface(Surface);
1448     Standard_Real Tol = Precision::Confusion(); //1.e-4;
1449     GeomAbs_Shape UCont = GeomAbs_C1, VCont = GeomAbs_C1;
1450     Standard_Integer degU = 14, degV = 14;
1451     Standard_Integer nmax = 16;
1452     Standard_Integer thePrec = 1;  
1453     GeomConvert_ApproxSurface theApprox(Surface,Tol,UCont,VCont,degU,degV,nmax,thePrec);
1454     if (theApprox.HasResult())
1455       BS = theApprox.Surface();
1456     else
1457       BS = GeomConvert::SurfaceToBSplineSurface(Surface);
1458   }
1459   if (InU&&(BS->UDegree()<Continuity+1)) 
1460     BS->IncreaseDegree(Continuity+1,BS->VDegree());      
1461   if (!InU&&(BS->VDegree()<Continuity+1))
1462     BS->IncreaseDegree(BS->UDegree(),Continuity+1);      
1463
1464   // si BS etait periodique dans le sens de l'extension, elle ne le sera plus
1465   if ( (InU&&(BS->IsUPeriodic())) || (!InU&&(BS->IsVPeriodic())) ) {
1466     Standard_Real U0,U1,V0,V1;
1467     BS->Bounds(U0,U1,V0,V1);
1468     BS->Segment(U0,U1,V0,V1);
1469   }     
1470
1471
1472 // IFV Fix OCC bug 0022694 - wrong result extrapolating rational surfaces
1473 //   Standard_Boolean rational = ( InU && BS->IsURational() ) 
1474 //                                   || ( !InU && BS->IsVRational() ) ;
1475   Standard_Boolean rational = (BS->IsURational() ||  BS->IsVRational());
1476   Standard_Boolean NullWeight;
1477    Standard_Real EpsW = 10*Precision::PConfusion();
1478   Standard_Integer gap = 3;
1479   if ( rational ) gap++;
1480
1481
1482         
1483   Standard_Integer Cdeg = 0, Cdim = 0, NbP = 0, Ksize = 0, Psize = 1;
1484   Standard_Integer ii, jj, ipole, Kount;  
1485   Standard_Real Tbord, lambmin=Length;
1486   Standard_Real * Padr = NULL;
1487   Standard_Boolean Ok;
1488   Handle(TColStd_HArray1OfReal)  FKnots, Point, lambda, Tgte, Poles;
1489
1490   
1491
1492
1493   for (Kount=0, Ok=Standard_False; Kount<=2 && !Ok; Kount++) {
1494     //  transformation de la surface en une BSpline non rationnelle a une variable
1495     //  de degre UDegree ou VDegree et de dimension 3 ou 4 x NbVpoles ou NbUpoles
1496     //  le nombre de poles egal a NbUpoles ou NbVpoles
1497     //  ATTENTION : dans le cas rationnel, un point de coordonnees (x,y,z)
1498     //              et de poids w devient un point de coordonnees (wx, wy, wz, w )
1499   
1500
1501     if (InU) {
1502       Cdeg = BS->UDegree();
1503       NbP = BS->NbUPoles();
1504       Cdim = BS->NbVPoles() * gap;
1505     }
1506     else {
1507       Cdeg = BS->VDegree();
1508       NbP = BS->NbVPoles();
1509       Cdim = BS->NbUPoles() * gap;
1510     }
1511
1512     //  les noeuds plats
1513     Ksize = NbP + Cdeg + 1;
1514     FKnots = new (TColStd_HArray1OfReal) (1,Ksize);
1515     if (InU) 
1516       BS->UKnotSequence(FKnots->ChangeArray1());
1517     else 
1518       BS->VKnotSequence(FKnots->ChangeArray1());
1519
1520     //  le parametre du noeud de raccord
1521     if (After)
1522       Tbord = FKnots->Value(FKnots->Upper()-Cdeg);
1523     else
1524       Tbord = FKnots->Value(FKnots->Lower()+Cdeg);
1525
1526     //  les poles
1527     Psize = Cdim * NbP;
1528     Poles = new (TColStd_HArray1OfReal) (1,Psize);
1529
1530     if (InU) {
1531       for (ii=1,ipole=1; ii<=NbP; ii++) {
1532         for (jj=1;jj<=BS->NbVPoles();jj++) {
1533           Poles->SetValue(ipole,   BS->Pole(ii,jj).X());
1534           Poles->SetValue(ipole+1, BS->Pole(ii,jj).Y());
1535           Poles->SetValue(ipole+2, BS->Pole(ii,jj).Z());
1536           if (rational) Poles->SetValue(ipole+3, BS->Weight(ii,jj));
1537           ipole+=gap;
1538         }
1539       }
1540     }
1541     else {
1542       for (jj=1,ipole=1; jj<=NbP; jj++) {
1543         for (ii=1;ii<=BS->NbUPoles();ii++) {
1544           Poles->SetValue(ipole,   BS->Pole(ii,jj).X());
1545           Poles->SetValue(ipole+1, BS->Pole(ii,jj).Y());
1546           Poles->SetValue(ipole+2, BS->Pole(ii,jj).Z());
1547           if (rational) Poles->SetValue(ipole+3, BS->Weight(ii,jj));
1548           ipole+=gap;
1549         }
1550       }
1551     }
1552     Padr = (Standard_Real *) &Poles->ChangeValue(1);
1553
1554     //  calcul du point de raccord et de la tangente
1555     Point = new (TColStd_HArray1OfReal)(1,Cdim);
1556     Tgte  = new (TColStd_HArray1OfReal)(1,Cdim);
1557     lambda = new (TColStd_HArray1OfReal)(1,Cdim);
1558
1559     Standard_Boolean  periodic_flag = Standard_False ;
1560     Standard_Integer extrap_mode[2], derivative_request = Max(Continuity,1);
1561     extrap_mode[0] = extrap_mode[1] = Cdeg;
1562     TColStd_Array1OfReal  Result(1, Cdim * (derivative_request+1)) ; 
1563     
1564     TColStd_Array1OfReal& tgte = Tgte->ChangeArray1();
1565     TColStd_Array1OfReal& point = Point->ChangeArray1();
1566     TColStd_Array1OfReal& lamb = lambda->ChangeArray1();
1567
1568     Standard_Real * Radr = (Standard_Real *) &Result(1) ;
1569
1570     BSplCLib::Eval(Tbord,periodic_flag,derivative_request,extrap_mode[0],
1571                    Cdeg,FKnots->Array1(),Cdim,*Padr,*Radr);
1572     Ok = Standard_True;
1573     for (ii=1;ii<=Cdim;ii++) {
1574       point(ii) = Result(ii);
1575       tgte(ii) = Result(ii+Cdim);
1576     }
1577   
1578     //  calcul de la contrainte a atteindre
1579
1580     gp_Vec CurT, OldT;
1581   
1582     Standard_Real NTgte, val, Tgtol = 1.e-12, OldN = 0.0;
1583     if (rational) {
1584       for (ii=gap;ii<=Cdim;ii+=gap) {
1585         tgte(ii) = 0.;
1586       }
1587       for (ii=gap;ii<=Cdim;ii+=gap) {
1588         CurT.SetCoord(tgte(ii-3),tgte(ii-2), tgte(ii-1)); 
1589         NTgte=CurT.Magnitude();
1590         if (NTgte>Tgtol) {
1591           val =  Length/NTgte;
1592           // Attentions aux Cas ou le segment donne par les poles 
1593           // est oppose au sens de la derive
1594           // Exemple: Certaine portions de tore.
1595           if ( (OldN > Tgtol) && (CurT.Angle(OldT) > 2)) {
1596             Ok = Standard_False;
1597           }
1598
1599           lamb(ii-1) = lamb(ii-2) = lamb(ii-3) = val;
1600           lamb(ii) = 0.;
1601           lambmin = Min(lambmin, val);
1602         }
1603         else {
1604           lamb(ii-1) = lamb(ii-2) = lamb(ii-3) = 0.;
1605           lamb(ii) = 0.;
1606         }
1607         OldT = CurT;
1608         OldN = NTgte;
1609       }
1610     }
1611     else {
1612       for (ii=gap;ii<=Cdim;ii+=gap) {
1613         CurT.SetCoord(tgte(ii-2),tgte(ii-1), tgte(ii)); 
1614         NTgte=CurT.Magnitude();
1615         if (NTgte>Tgtol) {
1616           val =  Length/NTgte;
1617           // Attentions aux Cas ou le segment donne par les poles 
1618           // est oppose au sens de la derive
1619           // Exemple: Certaine portion de tore.
1620           if ( (OldN > Tgtol) && (CurT.Angle(OldT) > 2)) {
1621              Ok = Standard_False;
1622           }
1623           lamb(ii) = lamb(ii-1) = lamb(ii-2) = val;
1624           lambmin = Min(lambmin, val);
1625         }
1626         else {
1627           lamb(ii) =lamb(ii-1) = lamb(ii-2) = 0.;
1628         }
1629         OldT = CurT;
1630         OldN = NTgte;
1631       }
1632     }
1633     if (!Ok && Kount<2) {
1634       // On augmente le degre de l'iso bord afin de rapprocher les poles de la surface
1635       // Et on ressaye
1636       if (InU) BS->IncreaseDegree(BS->UDegree(), BS->VDegree()+2);
1637       else     BS->IncreaseDegree(BS->UDegree()+2, BS->VDegree());
1638     }
1639   }
1640
1641
1642   TColStd_Array1OfReal ConstraintPoint(1,Cdim);
1643   if (After) {
1644     for (ii=1;ii<=Cdim;ii++) {
1645       ConstraintPoint(ii) = Point->Value(ii) + lambda->Value(ii)*Tgte->Value(ii);
1646     }
1647   }
1648   else {
1649     for (ii=1;ii<=Cdim;ii++) {
1650       ConstraintPoint(ii) = Point->Value(ii) - lambda->Value(ii)*Tgte->Value(ii);
1651     }
1652   }
1653
1654 //  cas particulier du rationnel
1655   if (rational) {
1656     for (ipole=1;ipole<=Psize;ipole+=gap) {
1657       Poles->ChangeValue(ipole) *= Poles->Value(ipole+3);
1658       Poles->ChangeValue(ipole+1) *= Poles->Value(ipole+3);
1659       Poles->ChangeValue(ipole+2) *= Poles->Value(ipole+3);
1660     }
1661     for (ii=1;ii<=Cdim;ii+=gap) {
1662       ConstraintPoint(ii) *= ConstraintPoint(ii+3);
1663       ConstraintPoint(ii+1) *= ConstraintPoint(ii+3);
1664       ConstraintPoint(ii+2) *= ConstraintPoint(ii+3);
1665     }
1666   }
1667   
1668 //  tableaux necessaires pour l'extension
1669   Standard_Integer Ksize2 = Ksize+Cdeg, NbPoles, NbKnots = 0;
1670   TColStd_Array1OfReal  FK(1, Ksize2) ; 
1671   Standard_Real * FKRadr = &FK(1);
1672
1673   Standard_Integer Psize2 = Psize+Cdeg*Cdim;
1674   TColStd_Array1OfReal  PRes(1, Psize2) ; 
1675   Standard_Real * PRadr = &PRes(1);
1676   Standard_Real ww;
1677   Standard_Boolean ExtOk = Standard_False;
1678   Handle(TColgp_HArray2OfPnt) NewPoles;
1679   Handle(TColStd_HArray2OfReal) NewWeights;
1680
1681
1682   for (Kount=1; Kount<=5 && !ExtOk; Kount++) {
1683     //  extension
1684     BSplCLib::TangExtendToConstraint(FKnots->Array1(),
1685                                      lambmin,NbP,*Padr,
1686                                      Cdim,Cdeg,
1687                                      ConstraintPoint, Cont, After,
1688                                      NbPoles, NbKnots,*FKRadr, *PRadr);
1689
1690     //  recopie des poles du resultat sous forme de points 3D et de poids
1691     Standard_Integer NU, NV, indice ;
1692     if (InU) {
1693       NU = NbPoles;
1694       NV = BS->NbVPoles();
1695     }
1696     else {
1697       NU = BS->NbUPoles();
1698       NV = NbPoles;
1699     }
1700
1701     NewPoles = new (TColgp_HArray2OfPnt)(1,NU,1,NV);
1702     TColgp_Array2OfPnt& NewP = NewPoles->ChangeArray2();
1703     NewWeights = new (TColStd_HArray2OfReal) (1,NU,1,NV);
1704     TColStd_Array2OfReal& NewW = NewWeights->ChangeArray2();
1705
1706     if (!rational) NewW.Init(1.);
1707     NullWeight= Standard_False;
1708
1709     if (InU) {
1710       for (ii=1; ii<=NU && !NullWeight; ii++) {
1711         for (jj=1; jj<=NV && !NullWeight; jj++) {
1712           indice = 1+(ii-1)*Cdim+(jj-1)*gap;
1713           NewP(ii,jj).SetCoord(1,PRes(indice));
1714           NewP(ii,jj).SetCoord(2,PRes(indice+1));
1715           NewP(ii,jj).SetCoord(3,PRes(indice+2));
1716           if (rational) {
1717             ww =  PRes(indice+3);
1718             if (ww < EpsW) {
1719               NullWeight = Standard_True;
1720             }
1721             else {
1722               NewW(ii,jj) = ww;
1723               NewP(ii,jj).ChangeCoord() /= ww;
1724             }
1725           }
1726         }
1727       }
1728     }
1729     else {
1730       for (jj=1; jj<=NV && !NullWeight; jj++) {
1731         for (ii=1; ii<=NU && !NullWeight; ii++) {
1732           indice = 1+(ii-1)*gap+(jj-1)*Cdim;
1733           NewP(ii,jj).SetCoord(1,PRes(indice));
1734           NewP(ii,jj).SetCoord(2,PRes(indice+1));
1735           NewP(ii,jj).SetCoord(3,PRes(indice+2));
1736           if (rational) {
1737             ww =  PRes(indice+3);
1738             if (ww < EpsW) {
1739               NullWeight = Standard_True;
1740             }
1741             else {
1742               NewW(ii,jj) = ww;
1743               NewP(ii,jj).ChangeCoord() /= ww;
1744             }
1745           }
1746         }
1747       }
1748     }
1749
1750     if (NullWeight) {
1751 #ifdef OCCT_DEBUG
1752       cout << "Echec de l'Extension rationnelle" << endl;    
1753 #endif
1754       lambmin /= 3.;
1755       NullWeight = Standard_False;
1756     }
1757     else {
1758       ExtOk = Standard_True;
1759     }
1760   }
1761     
1762
1763 // recopie des noeuds plats sous forme de noeuds avec leurs multiplicites
1764 // calcul des degres du resultat
1765   Standard_Integer Usize = BS->NbUKnots(), Vsize = BS->NbVKnots(), UDeg, VDeg;
1766   if (InU) 
1767     Usize++;
1768   else
1769     Vsize++;
1770   TColStd_Array1OfReal UKnots(1,Usize);
1771   TColStd_Array1OfReal VKnots(1,Vsize);
1772   TColStd_Array1OfInteger UMults(1,Usize);
1773   TColStd_Array1OfInteger VMults(1,Vsize);
1774   TColStd_Array1OfReal FKRes(1, NbKnots);
1775
1776   for (ii=1; ii<=NbKnots; ii++)
1777      FKRes(ii) = FK(ii);
1778
1779   if (InU) {
1780     BSplCLib::Knots(FKRes, UKnots, UMults);
1781     UDeg = Cdeg;
1782     UMults(Usize) = UDeg+1; // Petite verrue utile quand la continuite 
1783                              // n'est pas ok.
1784     BS->VKnots(VKnots);
1785     BS->VMultiplicities(VMults);
1786     VDeg = BS->VDegree();
1787   }
1788   else {
1789     BSplCLib::Knots(FKRes, VKnots, VMults);
1790     VDeg = Cdeg;
1791     VMults(Vsize) = VDeg+1;
1792     BS->UKnots(UKnots);
1793     BS->UMultiplicities(UMults);
1794     UDeg = BS->UDegree();
1795   }
1796
1797 //  construction de la surface BSpline resultat
1798   Handle(Geom_BSplineSurface) Res = 
1799     new (Geom_BSplineSurface) (NewPoles->Array2(),
1800                                NewWeights->Array2(),
1801                                UKnots,VKnots,
1802                                UMults,VMults,
1803                                UDeg,VDeg,
1804                                BS->IsUPeriodic(),
1805                                BS->IsVPeriodic());
1806   Surface = Res;
1807 }
1808
1809 //=======================================================================
1810 //function : Inertia
1811 //purpose  : 
1812 //=======================================================================
1813 void GeomLib::Inertia(const TColgp_Array1OfPnt& Points,
1814                       gp_Pnt& Bary,
1815                       gp_Dir& XDir,
1816                       gp_Dir& YDir,
1817                       Standard_Real& Xgap,
1818                       Standard_Real& Ygap,
1819                       Standard_Real& Zgap)
1820 {
1821   gp_XYZ GB(0., 0., 0.), Diff;
1822 //  gp_Vec A,B,C,D;
1823
1824   Standard_Integer i,nb=Points.Length();
1825   GB.SetCoord(0.,0.,0.);
1826   for (i=1; i<=nb; i++) 
1827       GB += Points(i).XYZ();
1828
1829   GB /= nb;
1830
1831   math_Matrix M (1, 3, 1, 3);
1832   M.Init(0.);
1833   for (i=1; i<=nb; i++) {
1834     Diff.SetLinearForm(-1, Points(i).XYZ(), GB);
1835     M(1,1) += Diff.X() *  Diff.X();
1836     M(2,2) += Diff.Y() *  Diff.Y();
1837     M(3,3) += Diff.Z() *  Diff.Z();
1838     M(1,2) += Diff.X() *  Diff.Y();
1839     M(1,3) += Diff.X() *  Diff.Z();
1840     M(2,3) += Diff.Y() *  Diff.Z();
1841   }
1842
1843   M(2,1)=M(1,2) ;
1844   M(3,1)=M(1,3) ;
1845   M(3,2)=M(2,3) ;
1846
1847   M /= nb;
1848
1849   math_Jacobi J(M);
1850   if (!J.IsDone()) {
1851 #ifdef OCCT_DEBUG
1852     cout << "Erreur dans Jacobbi" << endl;
1853     M.Dump(cout);
1854 #endif
1855   }
1856
1857   Standard_Real n1,n2,n3;
1858
1859   n1=J.Value(1);
1860   n2=J.Value(2);
1861   n3=J.Value(3);
1862
1863   Standard_Real r1 = Min(Min(n1,n2),n3), r2;
1864   Standard_Integer m1, m2, m3;
1865   if (r1==n1) {
1866     m1 = 1;
1867     r2 = Min(n2,n3);
1868     if (r2==n2) {
1869       m2 = 2;
1870       m3 = 3;
1871     }
1872     else {
1873       m2 = 3;
1874       m3 = 2;
1875     }
1876   }
1877   else {
1878     if (r1==n2) {
1879       m1 = 2 ;
1880       r2 = Min(n1,n3);
1881       if (r2==n1) {
1882         m2 = 1;
1883         m3 = 3;
1884       }
1885       else {
1886         m2 = 3;
1887         m3 = 1;
1888       }
1889     }
1890     else {
1891       m1 = 3 ;
1892       r2 = Min(n1,n2);
1893       if (r2==n1) {
1894         m2 = 1;
1895         m3 = 2;
1896       }
1897       else {
1898         m2 = 2;
1899         m3 = 1;
1900       }
1901     }
1902   }
1903
1904   math_Vector V2(1,3),V3(1,3);
1905   J.Vector(m2,V2);
1906   J.Vector(m3,V3);
1907   
1908   Bary.SetXYZ(GB);
1909   XDir.SetCoord(V3(1),V3(2),V3(3));
1910   YDir.SetCoord(V2(1),V2(2),V2(3));
1911
1912   Zgap = sqrt(Abs(J.Value(m1)));
1913   Ygap = sqrt(Abs(J.Value(m2)));
1914   Xgap = sqrt(Abs(J.Value(m3)));
1915 }
1916 //=======================================================================
1917 //function : AxeOfInertia
1918 //purpose  : 
1919 //=======================================================================
1920 void GeomLib::AxeOfInertia(const TColgp_Array1OfPnt& Points,
1921                            gp_Ax2& Axe,
1922                            Standard_Boolean& IsSingular,
1923                            const Standard_Real Tol)
1924 {
1925   gp_Pnt Bary;
1926   gp_Dir OX,OY,OZ;
1927   Standard_Real gx, gy, gz;
1928
1929   GeomLib::Inertia(Points, Bary, OX, OY, gx, gy, gz);
1930   
1931   if (gy*Points.Length()<=Tol) {
1932     gp_Ax2 axe (Bary, OX);
1933     OY = axe.XDirection();
1934     IsSingular = Standard_True;
1935   }
1936   else {
1937     IsSingular = Standard_False;
1938   }
1939
1940   OZ = OX^OY;
1941   gp_Ax2 TheAxe(Bary, OZ, OX);
1942   Axe = TheAxe;
1943 }
1944
1945 //=======================================================================
1946 //function : CanBeTreated
1947 //purpose  : indicates if the surface can be treated(if the conditions are
1948 //           filled) and need to be treated(if the surface hasn't been yet
1949 //           treated or if the surface is rationnal and non periodic)
1950 //=======================================================================
1951
1952 static Standard_Boolean CanBeTreated(Handle(Geom_BSplineSurface)& BSurf)
1953      
1954 {Standard_Integer i;
1955  Standard_Real    lambda;                                    //proportionnality coefficient
1956  Standard_Boolean AlreadyTreated=Standard_True;
1957  
1958  if (!BSurf->IsURational()||(BSurf->IsUPeriodic()))
1959    return Standard_False;
1960  else {
1961    lambda=(BSurf->Weight(1,1)/BSurf->Weight(BSurf->NbUPoles(),1));
1962    for (i=1;i<=BSurf->NbVPoles();i++)      //test of the proportionnality of the denominator on the boundaries
1963      if ((BSurf->Weight(1,i)/(lambda*BSurf->Weight(BSurf->NbUPoles(),i))<(1-Precision::Confusion()))||
1964          (BSurf->Weight(1,i)/(lambda*BSurf->Weight(BSurf->NbUPoles(),i))>(1+Precision::Confusion())))
1965        return Standard_False;
1966    i=1;
1967    while ((AlreadyTreated) && (i<=BSurf->NbVPoles())){        //tests if the surface has already been treated
1968      if (((BSurf->Weight(1,i)/(BSurf->Weight(2,i)))<(1-Precision::Confusion()))||
1969          ((BSurf->Weight(1,i)/(BSurf->Weight(2,i)))>(1+Precision::Confusion()))||
1970          ((BSurf->Weight(BSurf->NbUPoles()-1,i)/(BSurf->Weight(BSurf->NbUPoles(),i)))<(1-Precision::Confusion()))||
1971          ((BSurf->Weight(BSurf->NbUPoles()-1,i)/(BSurf->Weight(BSurf->NbUPoles(),i)))>(1+Precision::Confusion())))
1972        AlreadyTreated=Standard_False;
1973      i++;
1974    }
1975    if (AlreadyTreated)
1976      return Standard_False;
1977  }
1978  return Standard_True;  
1979 }
1980
1981 //=======================================================================
1982 //class   : law_evaluator
1983 //purpose : usefull to estimate the value of a function of 2 variables
1984 //=======================================================================
1985
1986 class law_evaluator : public BSplSLib_EvaluatorFunction
1987 {
1988
1989 public:
1990
1991   law_evaluator (const GeomLib_DenominatorMultiplierPtr theDenominatorPtr)
1992   : myDenominator (theDenominatorPtr) {}
1993
1994   virtual void Evaluate (const Standard_Integer theDerivativeRequest,
1995                          const Standard_Real    theUParameter,
1996                          const Standard_Real    theVParameter,
1997                          Standard_Real&         theResult,
1998                          Standard_Integer&      theErrorCode) const
1999   {
2000     if ((myDenominator != NULL) && (theDerivativeRequest == 0))
2001     {
2002       theResult = myDenominator->Value (theUParameter, theVParameter);
2003       theErrorCode = 0;
2004     }
2005     else
2006     {
2007       theErrorCode = 1;
2008     }
2009   }
2010
2011 private:
2012
2013   GeomLib_DenominatorMultiplierPtr myDenominator;
2014
2015 };
2016  
2017 //=======================================================================
2018 //function : CheckIfKnotExists
2019 //purpose  : true if the knot already exists in the knot sequence
2020 //=======================================================================
2021
2022 static Standard_Boolean CheckIfKnotExists(const TColStd_Array1OfReal&           surface_knots,
2023                                           const Standard_Real                   knot)
2024
2025 {Standard_Integer    i;
2026  for (i=1;i<=surface_knots.Length();i++)
2027    if ((surface_knots(i)-Precision::Confusion()<=knot)&&(surface_knots(i)+Precision::Confusion()>=knot))
2028      return Standard_True;
2029  return Standard_False;
2030 }
2031
2032 //=======================================================================
2033 //function : AddAKnot
2034 //purpose  : add a knot and its multiplicity to the knot sequence. This knot
2035 //           will be C2 and the degree is increased of deltasurface_degree 
2036 //=======================================================================
2037
2038 static void AddAKnot(const TColStd_Array1OfReal&           knots,
2039                      const TColStd_Array1OfInteger&        mults,
2040                      const Standard_Real                   knotinserted,
2041                      const Standard_Integer                deltasurface_degree,
2042                      const Standard_Integer                finalsurfacedegree,
2043                      Handle(TColStd_HArray1OfReal) &       newknots,
2044                      Handle(TColStd_HArray1OfInteger) &    newmults)
2045
2046 {Standard_Integer      i;
2047
2048  newknots=new TColStd_HArray1OfReal(1,knots.Length()+1);
2049  newmults=new TColStd_HArray1OfInteger(1,knots.Length()+1); 
2050  i=1;
2051  while (knots(i)<knotinserted){
2052    newknots->SetValue(i,knots(i));
2053    newmults->SetValue(i,mults(i)+deltasurface_degree);
2054    i++;
2055  }
2056  newknots->SetValue(i,knotinserted);                        //insertion of the new knot
2057  newmults->SetValue(i,finalsurfacedegree-2);
2058  i++;
2059  while (i<=newknots->Length()){
2060    newknots->SetValue(i,knots(i-1));
2061    newmults->SetValue(i,mults(i-1)+deltasurface_degree);
2062    i++;
2063  }
2064 }
2065
2066 //=======================================================================
2067 //function : Sort
2068 //purpose  : give the new flat knots(u or v) of the surface 
2069 //=======================================================================
2070
2071 static void BuildFlatKnot(const TColStd_Array1OfReal&           surface_knots,    
2072                  const TColStd_Array1OfInteger&        surface_mults,    
2073                  const Standard_Integer                deltasurface_degree, 
2074                  const Standard_Integer                finalsurface_degree, 
2075                  const Standard_Real                   knotmin,
2076                  const Standard_Real                   knotmax,
2077                  Handle(TColStd_HArray1OfReal)&        ResultKnots,
2078                  Handle(TColStd_HArray1OfInteger)&     ResultMults)
2079                  
2080 {
2081   Standard_Integer  i;
2082  
2083  if (CheckIfKnotExists(surface_knots,knotmin) &&
2084      CheckIfKnotExists(surface_knots,knotmax)){
2085    ResultKnots=new TColStd_HArray1OfReal(1,surface_knots.Length());
2086    ResultMults=new TColStd_HArray1OfInteger(1,surface_knots.Length());
2087    for (i=1;i<=surface_knots.Length();i++){
2088      ResultKnots->SetValue(i,surface_knots(i));
2089      ResultMults->SetValue(i,surface_mults(i)+deltasurface_degree);
2090    }
2091  }
2092  else{
2093    if ((CheckIfKnotExists(surface_knots,knotmin))&&(!CheckIfKnotExists(surface_knots,knotmax)))
2094      AddAKnot(surface_knots,surface_mults,knotmax,deltasurface_degree,finalsurface_degree,ResultKnots,ResultMults);
2095    else{
2096      if ((!CheckIfKnotExists(surface_knots,knotmin))&&(CheckIfKnotExists(surface_knots,knotmax)))
2097        AddAKnot(surface_knots,surface_mults,knotmin,deltasurface_degree,finalsurface_degree,ResultKnots,ResultMults);
2098      else{
2099        if ((!CheckIfKnotExists(surface_knots,knotmin))&&(!CheckIfKnotExists(surface_knots,knotmax))&&
2100            (knotmin==knotmax)){
2101          AddAKnot(surface_knots,surface_mults,knotmin,deltasurface_degree,finalsurface_degree,ResultKnots,ResultMults);
2102        }
2103        else{
2104          Handle(TColStd_HArray1OfReal)      IntermedKnots;
2105          Handle(TColStd_HArray1OfInteger)   IntermedMults;
2106          AddAKnot(surface_knots,surface_mults,knotmin,deltasurface_degree,finalsurface_degree,IntermedKnots,IntermedMults);
2107          AddAKnot(IntermedKnots->ChangeArray1(),IntermedMults->ChangeArray1(),knotmax,0,finalsurface_degree,ResultKnots,ResultMults);
2108        }
2109      }
2110    }
2111  }   
2112 }
2113
2114 //=======================================================================
2115 //function : FunctionMultiply 
2116 //purpose  : multiply the surface BSurf by a(u,v) (law_evaluator) on its
2117 //           numerator and denominator
2118 //=======================================================================
2119
2120 static void FunctionMultiply(Handle(Geom_BSplineSurface)&          BSurf,
2121                              const Standard_Real                   knotmin,
2122                              const Standard_Real                   knotmax)
2123      
2124 {TColStd_Array1OfReal      surface_u_knots(1,BSurf->NbUKnots()) ;
2125  TColStd_Array1OfInteger   surface_u_mults(1,BSurf->NbUKnots()) ;
2126  TColStd_Array1OfReal      surface_v_knots(1,BSurf->NbVKnots()) ;
2127  TColStd_Array1OfInteger   surface_v_mults(1,BSurf->NbVKnots()) ;
2128  TColgp_Array2OfPnt        surface_poles(1,BSurf->NbUPoles(),
2129                                          1,BSurf->NbVPoles()) ;
2130  TColStd_Array2OfReal      surface_weights(1,BSurf->NbUPoles(),
2131                                            1,BSurf->NbVPoles()) ;
2132  Standard_Integer          i,j,k,status,new_num_u_poles,new_num_v_poles,length=0;
2133  Handle(TColStd_HArray1OfReal)     newuknots,newvknots;
2134  Handle(TColStd_HArray1OfInteger)  newumults,newvmults;
2135
2136  BSurf->UKnots(surface_u_knots) ;
2137  BSurf->UMultiplicities(surface_u_mults) ;
2138  BSurf->VKnots(surface_v_knots) ;
2139  BSurf->VMultiplicities(surface_v_mults) ;
2140  BSurf->Poles(surface_poles) ;
2141  BSurf->Weights(surface_weights) ;
2142
2143  TColStd_Array1OfReal    Knots(1,2); 
2144  TColStd_Array1OfInteger Mults(1,2);
2145  Handle(TColStd_HArray1OfReal)      NewKnots;
2146  Handle(TColStd_HArray1OfInteger)   NewMults;
2147  
2148  Knots(1)=0;
2149  Knots(2)=1;
2150  Mults(1)=4;
2151  Mults(2)=4;
2152  BuildFlatKnot(Knots,Mults,0,3,knotmin,knotmax,NewKnots,NewMults);
2153
2154  for (i=1;i<=NewMults->Length();i++)
2155    length+=NewMults->Value(i);
2156  TColStd_Array1OfReal       FlatKnots(1,length);
2157  BSplCLib::KnotSequence(NewKnots->ChangeArray1(),NewMults->ChangeArray1(),FlatKnots);
2158
2159  GeomLib_DenominatorMultiplier aDenominator (BSurf, FlatKnots);
2160
2161  BuildFlatKnot(surface_u_knots,
2162                surface_u_mults,
2163                3,
2164                BSurf->UDegree()+3,
2165                knotmin,
2166                knotmax,
2167                newuknots,
2168                newumults);
2169  BuildFlatKnot(surface_v_knots,
2170                surface_v_mults,
2171                BSurf->VDegree(),
2172                2*(BSurf->VDegree()),
2173                1.0,
2174                0.0,
2175                newvknots,
2176                newvmults);
2177  length=0;
2178  for (i=1;i<=newumults->Length();i++)
2179    length+=newumults->Value(i);
2180  new_num_u_poles=(length-BSurf->UDegree()-3-1);
2181  TColStd_Array1OfReal       newuflatknots(1,length);
2182  length=0;
2183  for (i=1;i<=newvmults->Length();i++)
2184    length+=newvmults->Value(i);
2185  new_num_v_poles=(length-2*BSurf->VDegree()-1);
2186  TColStd_Array1OfReal       newvflatknots(1,length);
2187
2188  TColgp_Array2OfPnt        NewNumerator(1,new_num_u_poles,1,new_num_v_poles);
2189  TColStd_Array2OfReal      NewDenominator(1,new_num_u_poles,1,new_num_v_poles);
2190  
2191  BSplCLib::KnotSequence(newuknots->ChangeArray1(),newumults->ChangeArray1(),newuflatknots);
2192  BSplCLib::KnotSequence(newvknots->ChangeArray1(),newvmults->ChangeArray1(),newvflatknots);
2193 //POP pour WNT
2194  law_evaluator ev (&aDenominator);
2195 // BSplSLib::FunctionMultiply(law_evaluator,               //multiplication
2196  BSplSLib::FunctionMultiply(ev,               //multiplication
2197                             BSurf->UDegree(),
2198                             BSurf->VDegree(),
2199                             surface_u_knots,
2200                             surface_v_knots,
2201                             surface_u_mults,
2202                             surface_v_mults,
2203                             surface_poles,
2204                             surface_weights,
2205                             newuflatknots,
2206                             newvflatknots,
2207                             BSurf->UDegree()+3,
2208                             2*(BSurf->VDegree()),
2209                             NewNumerator,
2210                             NewDenominator,
2211                             status);
2212  if (status!=0)
2213    Standard_ConstructionError::Raise("GeomLib Multiplication Error") ;
2214  for (i = 1 ; i <= new_num_u_poles ; i++) {
2215       for (j = 1 ; j <= new_num_v_poles ; j++) {
2216         for (k = 1 ; k <= 3 ; k++) {
2217           NewNumerator(i,j).SetCoord(k,NewNumerator(i,j).Coord(k)/NewDenominator(i,j)) ;
2218         }
2219       }
2220     }
2221  BSurf= new Geom_BSplineSurface(NewNumerator,                  
2222                                 NewDenominator,
2223                                 newuknots->ChangeArray1(),
2224                                 newvknots->ChangeArray1(),
2225                                 newumults->ChangeArray1(),
2226                                 newvmults->ChangeArray1(),
2227                                 BSurf->UDegree()+3,
2228                                 2*(BSurf->VDegree()) );             
2229 }
2230
2231 //=======================================================================
2232 //function : CancelDenominatorDerivative1D
2233 //purpose  : cancel the denominator derivative in one direction
2234 //=======================================================================
2235
2236 static void CancelDenominatorDerivative1D(Handle(Geom_BSplineSurface) & BSurf)
2237      
2238 {Standard_Integer            i,j;
2239  Standard_Real               uknotmin=1.0,uknotmax=0.0,
2240                              x,y,
2241                              startu_value,
2242                              endu_value;
2243  TColStd_Array1OfReal        BSurf_u_knots(1,BSurf->NbUKnots()) ;
2244
2245  startu_value=BSurf->UKnot(1);
2246  endu_value=BSurf->UKnot(BSurf->NbUKnots());
2247  BSurf->UKnots(BSurf_u_knots) ;
2248  BSplCLib::Reparametrize(0.0,1.0,BSurf_u_knots);
2249  BSurf->SetUKnots(BSurf_u_knots);                             //reparametrisation of the surface
2250  Handle(Geom_BSplineCurve) BCurve;
2251  TColStd_Array1OfReal      BCurveWeights(1,BSurf->NbUPoles());
2252  TColgp_Array1OfPnt        BCurvePoles(1,BSurf->NbUPoles());
2253  TColStd_Array1OfReal      BCurveKnots(1,BSurf->NbUKnots());
2254  TColStd_Array1OfInteger   BCurveMults(1,BSurf->NbUKnots());
2255
2256  if (CanBeTreated(BSurf)){
2257    for (i=1;i<=BSurf->NbVPoles();i++){  //loop on each pole function
2258      x=1.0;y=0.0;
2259      for (j=1;j<=BSurf->NbUPoles();j++){
2260        BCurveWeights(j)=BSurf->Weight(j,i);
2261        BCurvePoles(j)=BSurf->Pole(j,i);
2262      }
2263      BSurf->UKnots(BCurveKnots);
2264      BSurf->UMultiplicities(BCurveMults);
2265      BCurve = new Geom_BSplineCurve(BCurvePoles, //building of a pole function 
2266                                     BCurveWeights,
2267                                     BCurveKnots,
2268                                     BCurveMults,
2269                                     BSurf->UDegree());
2270      Hermit::Solutionbis(BCurve,x,y,Precision::Confusion(),Precision::Confusion()); 
2271      if (x<uknotmin)
2272        uknotmin=x;    //uknotmin,uknotmax:extremal knots
2273      if ((x!=1.0)&&(x>uknotmax))
2274        uknotmax=x;
2275      if ((y!=0.0)&&(y<uknotmin))
2276        uknotmin=y;
2277      if (y>uknotmax)
2278        uknotmax=y;
2279    }
2280   
2281    FunctionMultiply(BSurf,uknotmin,uknotmax);                 //multiplication
2282
2283    BSurf->UKnots(BSurf_u_knots) ;
2284    BSplCLib::Reparametrize(startu_value,endu_value,BSurf_u_knots);
2285    BSurf->SetUKnots(BSurf_u_knots);
2286  }
2287 }
2288
2289 //=======================================================================
2290 //function : CancelDenominatorDerivative
2291 //purpose  : 
2292 //=======================================================================
2293
2294 void GeomLib::CancelDenominatorDerivative(Handle(Geom_BSplineSurface)         & BSurf,
2295                                           const Standard_Boolean              udirection,
2296                                           const Standard_Boolean              vdirection)
2297
2298 {if (udirection && !vdirection)
2299    CancelDenominatorDerivative1D(BSurf);
2300  else{
2301    if (!udirection && vdirection) {
2302      BSurf->ExchangeUV();
2303      CancelDenominatorDerivative1D(BSurf);
2304      BSurf->ExchangeUV();
2305    }
2306    else{
2307      if (udirection && vdirection){                            //optimize the treatment
2308        if (BSurf->UDegree()<=BSurf->VDegree()){
2309          CancelDenominatorDerivative1D(BSurf);
2310          BSurf->ExchangeUV();
2311          CancelDenominatorDerivative1D(BSurf);
2312          BSurf->ExchangeUV();
2313        }
2314        else{
2315          BSurf->ExchangeUV();
2316          CancelDenominatorDerivative1D(BSurf);
2317          BSurf->ExchangeUV();
2318          CancelDenominatorDerivative1D(BSurf);
2319        }
2320      }
2321    }
2322  }
2323 }
2324
2325 //=======================================================================
2326 //function : NormEstim
2327 //purpose  : 
2328 //=======================================================================
2329
2330 Standard_Integer GeomLib::NormEstim(const Handle(Geom_Surface)& S, 
2331                                     const gp_Pnt2d& UV, 
2332                                     const Standard_Real Tol, gp_Dir& N) 
2333 {
2334   gp_Vec DU, DV;
2335   gp_Pnt DummyPnt;
2336   Standard_Real aTol2 = Square(Tol);
2337
2338   S->D1(UV.X(), UV.Y(), DummyPnt, DU, DV);
2339
2340   Standard_Real MDU = DU.SquareMagnitude(), MDV = DV.SquareMagnitude();
2341
2342   if(MDU >= aTol2 && MDV >= aTol2) {
2343     gp_Vec Norm = DU^DV;
2344     Standard_Real Magn = Norm.SquareMagnitude();
2345     if(Magn < aTol2) return 3;
2346
2347     //Magn = sqrt(Magn);
2348     N.SetXYZ(Norm.XYZ());
2349
2350     return 0;
2351   }
2352   else {
2353     gp_Vec D2U, D2V, D2UV;
2354     Standard_Boolean isDone;
2355     CSLib_NormalStatus aStatus;
2356     gp_Dir aNormal;
2357
2358     S->D2(UV.X(), UV.Y(), DummyPnt, DU, DV, D2U, D2V, D2UV);
2359     CSLib::Normal(DU, DV, D2U, D2V, D2UV, Tol, isDone, aStatus, aNormal);
2360
2361     if (isDone) {
2362      Standard_Real Umin, Umax, Vmin, Vmax;
2363      Standard_Real step = 1.0e-5;
2364      Standard_Real eps = 1.0e-16;
2365      Standard_Real sign = -1.0;
2366      
2367      S->Bounds(Umin, Umax, Vmin, Vmax);
2368
2369      // check for cone apex singularity point
2370      if ((UV.Y() > Vmin + step) && (UV.Y() < Vmax - step))
2371      {
2372        gp_Dir aNormal1, aNormal2;
2373        Standard_Real aConeSingularityAngleEps = 1.0e-4;
2374        S->D1(UV.X(), UV.Y() - sign * step, DummyPnt, DU, DV);
2375        if ((DU.XYZ().SquareModulus() > eps) && (DV.XYZ().SquareModulus() > eps)) {
2376          aNormal1 = DU^DV;
2377          S->D1(UV.X(), UV.Y() + sign * step, DummyPnt, DU, DV);
2378          if ((DU.XYZ().SquareModulus() > eps) && (DV.XYZ().SquareModulus() > eps)) {
2379            aNormal2 = DU^DV;
2380            if (aNormal1.IsOpposite(aNormal2, aConeSingularityAngleEps))
2381              return 2;
2382          }
2383        }
2384      }
2385
2386      // Along V
2387      if(MDU < aTol2 && MDV >= aTol2) {
2388        if ((Vmax - UV.Y()) > (UV.Y() - Vmin))
2389          sign = 1.0;
2390        S->D1(UV.X(), UV.Y() + sign * step, DummyPnt, DU, DV);
2391        gp_Vec Norm = DU^DV;
2392        if (Norm.SquareMagnitude() < eps) {
2393          Standard_Real sign1 = -1.0;
2394          if ((Umax - UV.X()) > (UV.X() - Umin))
2395            sign1 = 1.0;
2396          S->D1(UV.X() + sign1 * step, UV.Y() + sign * step, DummyPnt, DU, DV);
2397          Norm = DU^DV;
2398        }
2399        if ((Norm.SquareMagnitude() >= eps) && (Norm.Dot(aNormal) < 0.0))
2400          aNormal.Reverse();
2401      }
2402
2403      // Along U
2404      if(MDV < aTol2 && MDU >= aTol2) {
2405        if ((Umax - UV.X()) > (UV.X() - Umin))
2406          sign = 1.0;
2407        S->D1(UV.X() + sign * step, UV.Y(), DummyPnt, DU, DV);
2408        gp_Vec Norm = DU^DV;
2409        if (Norm.SquareMagnitude() < eps) {
2410          Standard_Real sign1 = -1.0;
2411          if ((Vmax - UV.Y()) > (UV.Y() - Vmin))
2412            sign1 = 1.0;
2413          S->D1(UV.X() + sign * step, UV.Y() + sign1 * step, DummyPnt, DU, DV);
2414          Norm = DU^DV;
2415        }
2416        if ((Norm.SquareMagnitude() >= eps) && (Norm.Dot(aNormal) < 0.0))
2417          aNormal.Reverse();
2418      }
2419
2420       // quasysingular
2421       if ((aStatus == CSLib_D1NuIsNull) || (aStatus == CSLib_D1NvIsNull) || 
2422           (aStatus == CSLib_D1NuIsParallelD1Nv)) {
2423             N.SetXYZ(aNormal.XYZ());
2424             return 1;
2425       }
2426       // conical
2427       if (aStatus == CSLib_InfinityOfSolutions)
2428           return 2;
2429     }
2430     // computation is impossible
2431     else {
2432       // conical
2433       if (aStatus == CSLib_D1NIsNull) {
2434         return 2;
2435       }
2436       return 3;
2437     }
2438   }
2439   return 3;
2440 }
2441  
2442