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