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