src/GeomLib/GeomLib.cxx

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