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