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