src/GeomFill/GeomFill_SectionPlacement.cxx

   1 // Created on: 1997-12-15
   2 // Created by: Philippe MANGIN
   3 // Copyright (c) 1997-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
  18 #include <Adaptor3d_HCurve.hxx>
  19 #include <Bnd_Box.hxx>
  20 #include <BndLib_Add3dCurve.hxx>
  21 #include <Extrema_ExtCC.hxx>
  22 #include <Extrema_POnCurv.hxx>
  23 #include <Geom_BezierCurve.hxx>
  24 #include <Geom_BSplineCurve.hxx>
  25 #include <Geom_CartesianPoint.hxx>
  26 #include <Geom_Conic.hxx>
  27 #include <Geom_Curve.hxx>
  28 #include <Geom_Geometry.hxx>
  29 #include <Geom_Line.hxx>
  30 #include <Geom_Plane.hxx>
  31 #include <Geom_TrimmedCurve.hxx>
  32 #include <GeomAbs_CurveType.hxx>
  33 #include <GeomAdaptor_HSurface.hxx>
  34 #include <GeomFill_LocationLaw.hxx>
  35 #include <GeomFill_SectionPlacement.hxx>
  36 #include <GeomLib.hxx>
  37 #include <GeomLProp_CLProps.hxx>
  38 #include <gp.hxx>
  39 #include <gp_Ax2.hxx>
  40 #include <gp_Ax3.hxx>
  41 #include <gp_Dir.hxx>
  42 #include <gp_Mat.hxx>
  43 #include <gp_Pnt.hxx>
  44 #include <gp_Trsf.hxx>
  45 #include <gp_Vec.hxx>
  46 #include <IntCurveSurface_HInter.hxx>
  47 #include <IntCurveSurface_IntersectionPoint.hxx>
  48 #include <Precision.hxx>
  49 #include <StdFail_NotDone.hxx>
  50 #include <TColgp_HArray1OfPnt.hxx>
  51
  52 //===============================================================
  53 // Function :
  54 // Purpose :
  55 //===============================================================
  56 static void Tangente(const Adaptor3d_Curve& Path,
  57                      const Standard_Real Param,
  58                      gp_Pnt& P,
  59                      gp_Vec& Tang)
  60 {
  61   Path.D1(Param, P, Tang);
  62   Standard_Real Norm = Tang.Magnitude();
  63
  64   for (Standard_Integer ii=2; (ii<12) && (Norm < Precision::Confusion());
  65        ii++) {
  66     Tang =  Path.DN(Param, ii);
  67     Norm = Tang.Magnitude();
  68   }
  69
  70   if (Norm > 100*gp::Resolution()) Tang /= Norm;
  71 }
  72
  73
  74 static Standard_Real Penalite(const Standard_Real angle,
  75                               const Standard_Real dist)
  76 {
  77   Standard_Real penal;
  78
  79   if (dist < 1)
  80     penal = Sqrt(dist);
  81   else if (dist <2)
  82     penal = Pow(dist, 2);
  83   else
  84     penal = dist + 2;
  85
  86   if (angle > 1.e-3) {
  87     penal += 1./angle -2./M_PI;
  88   }
  89   else {
  90     penal += 1.e3;
  91   }
  92
  93   return penal;
  94 }
  95
  96 static Standard_Real EvalAngle(const gp_Vec& V1,
  97                                const gp_Vec& V2)
  98 {
  99  Standard_Real angle;
 100  angle = V1.Angle(V2);
 101  if (angle > M_PI/2) angle = M_PI - angle;
 102  return angle;
 103 }
 104
 105 static Standard_Integer NbSamples(const Handle(Geom_Curve)& aCurve)
 106 {
 107   Standard_Real nbs = 100.; //on default
 108
 109   Handle(Geom_Curve) theCurve = aCurve;
 110   if (aCurve->IsInstance(STANDARD_TYPE(Geom_TrimmedCurve)))
 111     theCurve = (Handle(Geom_TrimmedCurve)::DownCast(aCurve))->BasisCurve();
 112
 113   if (theCurve->IsInstance(STANDARD_TYPE(Geom_Line)))
 114     nbs = 1;
 115   else if (theCurve->IsKind(STANDARD_TYPE(Geom_Conic)))
 116     nbs = 4;
 117   else if (theCurve->IsInstance(STANDARD_TYPE(Geom_BezierCurve)))
 118     {
 119       Handle(Geom_BezierCurve) BC = Handle(Geom_BezierCurve)::DownCast(theCurve);
 120       nbs = 3 + BC->NbPoles();
 121     }
 122   else if (theCurve->IsInstance(STANDARD_TYPE(Geom_BSplineCurve)))
 123     {
 124       Handle(Geom_BSplineCurve) BC = Handle(Geom_BSplineCurve)::DownCast(theCurve);
 125       nbs  = BC->NbKnots();
 126       nbs *= BC->Degree();
 127       Standard_Real ratio =
 128         (aCurve->LastParameter() - aCurve->FirstParameter())/(BC->LastParameter() - BC->FirstParameter());
 129       nbs *= ratio;
 130       if(nbs < 4.0)
 131         nbs = 4;
 132     }
 133
 134   if (nbs > 300.)
 135     nbs = 300;
 136
 137   return ((Standard_Integer)nbs);
 138 }
 139
 140 //===============================================================
 141 // Function :DistMini
 142 // Purpose : Examine un extrema pour updater <Dist> & <Param>
 143 //===============================================================
 144 static void DistMini(const Extrema_ExtPC& Ext,
 145                      const Adaptor3d_Curve& C,
 146                      Standard_Real& Dist,
 147                      Standard_Real& Param)
 148 {
 149   Standard_Real dist1, dist2;
 150   Standard_Integer ii;
 151   gp_Pnt P1, P2;
 152   Standard_Real Dist2 = RealLast();
 153
 154   Ext.TrimmedSquareDistances(dist1, dist2, P1, P2);
 155   if ( (dist1<Dist2) || (dist2<Dist2) ) {
 156     if (dist1 < dist2) {
 157       Dist2 = dist1;
 158       Param  = C.FirstParameter();
 159     }
 160     else {
 161       Dist2 = dist2;
 162       Param  = C.LastParameter();
 163     }
 164   }
 165
 166   if (Ext.IsDone())
 167     {
 168       for (ii=1; ii<= Ext.NbExt(); ii++) {
 169         if (Ext.SquareDistance(ii) < Dist2) {
 170           Dist2 = Ext.SquareDistance(ii);
 171           Param  =  Ext.Point(ii).Parameter();
 172         }
 173       }
 174     }
 175   Dist = sqrt (Dist2);
 176 }
 177
 178
 179 //===============================================================
 180 // Function : Constructeur
 181 // Purpose :
 182 //===============================================================
 183 GeomFill_SectionPlacement::
 184 GeomFill_SectionPlacement(const Handle(GeomFill_LocationLaw)& L,
 185                           const Handle(Geom_Geometry)& Section) :
 186                           myLaw(L), /* myAdpSection(Section),  mySection(Section), */
 187                           Dist(RealLast()), AngleMax(0.)
 188 {
 189
 190   done   = Standard_False;
 191   isplan = Standard_False;
 192   myIsPoint = Standard_False;
 193
 194   if (Section->IsInstance(STANDARD_TYPE(Geom_CartesianPoint)))
 195     {
 196       myIsPoint = Standard_True;
 197       myPoint   = Handle(Geom_CartesianPoint)::DownCast(Section)->Pnt();
 198       isplan    = Standard_True;
 199     }
 200   else
 201     {
 202       Handle(Geom_Curve) CurveSection = Handle(Geom_Curve)::DownCast(Section);
 203       myAdpSection.Load(CurveSection);
 204       mySection = CurveSection;
 205     }
 206
 207   Standard_Integer i, j, NbPoles=0;
 208   Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
 209
 210   // Boite d'encombrement de la section pour en deduire le gabarit
 211   Bnd_Box box;
 212   if (myIsPoint)
 213     box.Add(myPoint);
 214   else
 215     BndLib_Add3dCurve::Add(myAdpSection, 1.e-4, box);
 216   box.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 217
 218   Standard_Real DX = aXmax - aXmin ;
 219   Standard_Real DY = aYmax - aYmin ;
 220   Standard_Real DZ = aZmax - aZmin ;
 221   Gabarit = Sqrt( DX*DX+ DY*DY + DZ*DZ )/2. ;
 222
 223   Gabarit += Precision::Confusion(); // Cas des toute petite
 224
 225
 226   // Initialisation de TheAxe pour les cas singulier
 227   if (!myIsPoint)
 228     {
 229       gp_Pnt P;
 230       gp_Vec V;
 231       Tangente(myAdpSection,
 232                (myAdpSection.FirstParameter()+myAdpSection.LastParameter())/2,
 233                P, V);
 234       TheAxe.SetLocation(P);
 235       TheAxe.SetDirection(V);
 236
 237       // y a t'il un Plan moyen ?
 238       GeomAbs_CurveType TheType = myAdpSection.GetType();
 239       switch  (TheType) {
 240       case GeomAbs_Circle:
 241         {
 242           isplan = Standard_True;
 243           TheAxe =  myAdpSection.Circle().Axis();
 244           break;
 245         }
 246       case GeomAbs_Ellipse:
 247         {
 248           isplan = Standard_True;
 249           TheAxe =  myAdpSection.Ellipse().Axis();
 250           break;
 251         }
 252       case GeomAbs_Hyperbola:
 253         {
 254           isplan = Standard_True;
 255           TheAxe =  myAdpSection.Hyperbola().Axis();
 256           break;
 257         }
 258       case GeomAbs_Parabola:
 259         {
 260           isplan = Standard_True;
 261           TheAxe =  myAdpSection.Parabola().Axis();
 262           break;
 263         }
 264       case GeomAbs_Line:
 265         {
 266           NbPoles = 0; // Pas de Plan !!
 267           break;
 268         }
 269       case GeomAbs_BezierCurve:
 270       case GeomAbs_BSplineCurve:
 271         {
 272           NbPoles = myAdpSection.NbPoles();
 273           break;
 274         }
 275       default:
 276         NbPoles = 21;
 277       }
 278
 279
 280       if (!isplan && NbPoles>2)
 281         {
 282           // Calcul d'un plan moyen.
 283           Handle(TColgp_HArray1OfPnt) Pnts;
 284           Standard_Real first = myAdpSection.FirstParameter();
 285           Standard_Real last =  myAdpSection.LastParameter();
 286           if (myAdpSection.IsPeriodic())
 287           {
 288             //Correct boundaries to avoid mistake of LocateU
 289             Handle(Geom_Curve) aCurve = myAdpSection.Curve();
 290             if (aCurve->IsInstance(STANDARD_TYPE(Geom_TrimmedCurve)))
 291               aCurve = (Handle(Geom_TrimmedCurve)::DownCast(aCurve))->BasisCurve();
 292             Standard_Real Ufirst = aCurve->FirstParameter();
 293             Standard_Real aPeriod = aCurve->Period();
 294             Standard_Real U1 = Ufirst + Floor((first - Ufirst)/aPeriod) * aPeriod;
 295             Standard_Real U2 = U1 + aPeriod;
 296             if (Abs(first - U1) <= Precision::PConfusion())
 297               first = U1;
 298             if (Abs(last - U2) <= Precision::PConfusion())
 299               last = U2;
 300           }
 301           Standard_Real t, delta;
 302           if (myAdpSection.GetType() == GeomAbs_BSplineCurve)
 303             {
 304               Handle(Geom_BSplineCurve) BC =
 305                 Handle(Geom_BSplineCurve)::DownCast(myAdpSection.Curve());
 306               Standard_Integer I1, I2, I3, I4;
 307               BC->LocateU( first, Precision::Confusion(), I1, I2 );
 308               BC->LocateU( last,  Precision::Confusion(), I3, I4 );
 309               Standard_Integer NbKnots = I3 - I2 + 1;
 310
 311               Standard_Integer NbLocalPnts = 10;
 312               Standard_Integer NbPnts = (NbKnots-1) * NbLocalPnts;
 313               if (I1 != I2)
 314                 NbPnts += NbLocalPnts;
 315               if (I3 != I4 && first < BC->Knot(I3))
 316                 NbPnts += NbLocalPnts;
 317               if (!myAdpSection.IsClosed())
 318                 NbPnts++;
 319               Pnts = new TColgp_HArray1OfPnt(1, NbPnts);
 320               Standard_Integer nb = 1;
 321               if (I1 != I2)
 322                 {
 323                   Standard_Real locallast = (BC->Knot(I2) < last)? BC->Knot(I2) : last;
 324                   delta = (locallast - first) / NbLocalPnts;
 325                   for (j = 0; j < NbLocalPnts; j++)
 326                     {
 327                       t = first + j*delta;
 328                       Pnts->SetValue( nb++, myAdpSection.Value(t) );
 329                     }
 330                 }
 331               for (i = I2; i < I3; i++)
 332                 {
 333                   t = BC->Knot(i);
 334                   delta = (BC->Knot(i+1) - t) / NbLocalPnts;
 335                   for (j = 0; j < NbLocalPnts; j++)
 336                     {
 337                       Pnts->SetValue( nb++, myAdpSection.Value(t) );
 338                       t += delta;
 339                     }
 340                 }
 341               if (I3 != I4 && first < BC->Knot(I3))
 342                 {
 343                   t = BC->Knot(I3);
 344                   delta = (last - t) / NbLocalPnts;
 345                   for (j = 0; j < NbLocalPnts; j++)
 346                     {
 347                       Pnts->SetValue( nb++, myAdpSection.Value(t) );
 348                       t += delta;
 349                     }
 350                 }
 351               if (!myAdpSection.IsClosed())
 352                 Pnts->SetValue( nb, myAdpSection.Value(last) );
 353             }
 354           else // other type
 355             {
 356               Standard_Integer NbPnts = NbPoles-1;
 357               if (!myAdpSection.IsClosed())
 358                 NbPnts++;
 359               Pnts = new TColgp_HArray1OfPnt(1, NbPnts);
 360               delta = (last - first) / (NbPoles-1);
 361               for (i = 0; i < NbPoles-1; i++)
 362                 {
 363                   t = first + i*delta;
 364                   Pnts->SetValue( i+1, myAdpSection.Value(t) );
 365                 }
 366               if (!myAdpSection.IsClosed())
 367                 Pnts->SetValue( NbPnts, myAdpSection.Value(last) );
 368             }
 369
 370           Standard_Boolean issing;
 371           gp_Ax2 axe;
 372           GeomLib::AxeOfInertia(Pnts->Array1(), axe, issing, Precision::Confusion());
 373           if (!issing) {
 374             isplan = Standard_True;
 375             TheAxe.SetLocation(axe.Location());
 376             TheAxe.SetDirection(axe.Direction());
 377           }
 378         }
 379
 380
 381       myExt.Initialize(myAdpSection,
 382                        myAdpSection.FirstParameter(),
 383                        myAdpSection.LastParameter(),
 384                        Precision::Confusion());
 385     }
 386 }
 387
 388 //===============================================================
 389 // Function :SetLocation
 390 // Purpose :
 391 //===============================================================
 392 void GeomFill_SectionPlacement::
 393 SetLocation(const Handle(GeomFill_LocationLaw)& L)
 394 {
 395   myLaw = L;
 396 }
 397
 398 //===============================================================
 399 // Function : Perform
 400 // Purpose : Le plus simple
 401 //===============================================================
 402 void GeomFill_SectionPlacement::Perform(const Standard_Real Tol)
 403 {
 404   Handle(Adaptor3d_HCurve) Path;
 405   Path = myLaw->GetCurve();
 406   Perform(Path, Tol);
 407 }
 408
 409 //===============================================================
 410 // Function :Perform
 411 // Purpose : Recherche automatique
 412 //===============================================================
 413 void GeomFill_SectionPlacement::Perform(const Handle(Adaptor3d_HCurve)& Path,
 414                                         const Standard_Real Tol)
 415 {
 416   Standard_Real IntTol = 1.e-5;
 417   Standard_Real DistCenter = Precision::Infinite();
 418
 419   if (myIsPoint)
 420     {
 421       Extrema_ExtPC Projector(myPoint, Path->Curve(), Precision::Confusion());
 422       DistMini( Projector, Path->Curve(), Dist, PathParam );
 423       AngleMax = M_PI/2;
 424     }
 425   else
 426     {
 427       PathParam = Path->FirstParameter();
 428       SecParam =  myAdpSection.FirstParameter();
 429
 430       Standard_Real distaux, taux = 0.0, alpha;
 431       gp_Pnt PonPath, PonSec, P;
 432       gp_Vec VRef, dp1;
 433       VRef.SetXYZ(TheAxe.Direction().XYZ());
 434
 435       Tangente( Path->Curve(), PathParam, PonPath, dp1);
 436       PonSec = myAdpSection.Value(SecParam);
 437       Dist =  PonPath.Distance(PonSec);
 438       if (Dist > Tol) { // On Cherche un meilleur point sur la section
 439         myExt.Perform(PonPath);
 440         if ( myExt.IsDone() ) {
 441           DistMini(myExt, myAdpSection, Dist, SecParam);
 442           PonSec = myAdpSection.Value(SecParam);
 443         }
 444       }
 445       AngleMax = EvalAngle(VRef, dp1);
 446       if (isplan) AngleMax = M_PI/2 - AngleMax;
 447
 448       Standard_Boolean Trouve = Standard_False;
 449       Standard_Integer ii;
 450
 451       if (isplan) {
 452         // (1.1) Distances Point-Plan
 453         Standard_Real DistPlan;
 454         gp_Vec V1 (PonPath, TheAxe.Location());
 455         DistPlan = Abs(V1.Dot(VRef));
 456         if (DistPlan <= IntTol)
 457           DistCenter = V1.Magnitude();
 458
 459         gp_Pnt Plast = Path->Value( Path->LastParameter() );
 460         V1.SetXYZ( TheAxe.Location().XYZ() - Plast.XYZ() );
 461         DistPlan = Abs(V1.Dot(VRef));
 462         if (DistPlan <= IntTol)
 463           {
 464             Standard_Real aDist = V1.Magnitude();
 465             if (aDist < DistCenter)
 466               {
 467                 DistCenter = aDist;
 468                 PonPath = Plast;
 469                 PathParam = Path->LastParameter();
 470               }
 471           }
 472
 473         // (1.2) Intersection Plan-courbe
 474         gp_Ax3 axe (TheAxe.Location(), TheAxe.Direction());
 475         Handle(Geom_Plane) plan = new (Geom_Plane)(axe);
 476         Handle(GeomAdaptor_HSurface) adplan =
 477           new (GeomAdaptor_HSurface)(plan);
 478         IntCurveSurface_HInter Intersector;
 479         Intersector.Perform(Path, adplan);
 480         if (Intersector.IsDone())
 481           {
 482             Standard_Real w;
 483             Standard_Real aDist;
 484             for (ii=1; ii<=Intersector.NbPoints(); ii++)
 485               {
 486                 w = Intersector.Point(ii).W();
 487                 P = Path->Value( w );
 488                 aDist = P.Distance( TheAxe.Location() );
 489                 if (aDist < DistCenter)
 490                   {
 491                     DistCenter = aDist;
 492                     PonPath = P;
 493                     PathParam = w;
 494                   }
 495               }
 496           }
 497         if (!Intersector.IsDone() || Intersector.NbPoints() == 0)
 498           {
 499             Standard_Integer NbPnts = NbSamples( mySection );
 500             TColgp_Array1OfPnt Pnts( 1, NbPnts+1 );
 501             Standard_Real delta = (mySection->LastParameter()-mySection->FirstParameter())/NbPnts;
 502             for (ii = 0; ii <= NbPnts; ii++)
 503               Pnts(ii+1) = mySection->Value( mySection->FirstParameter() + ii*delta );
 504
 505             gp_Pnt BaryCenter;
 506             gp_Dir Xdir, Ydir;
 507             Standard_Real Xgap, Ygap, Zgap;
 508             GeomLib::Inertia( Pnts, BaryCenter, Xdir, Ydir, Xgap, Ygap, Zgap );
 509
 510             gp_Pnt Pfirst = Path->Value( Path->FirstParameter() );
 511             if (Pfirst.Distance(BaryCenter) < Plast.Distance(BaryCenter))
 512               PathParam = Path->FirstParameter();
 513             else
 514               {
 515                 PathParam = Path->LastParameter();
 516                 Tangente( Path->Curve(), PathParam, PonPath, dp1);
 517                 PonSec = myAdpSection.Value(SecParam);
 518                 Dist =  PonPath.Distance(PonSec);
 519                 if (Dist > Tol) { // On Cherche un meilleur point sur la section
 520                   myExt.Perform(PonPath);
 521                   if ( myExt.IsDone() ) {
 522                     DistMini(myExt, myAdpSection, Dist, SecParam);
 523                     PonSec = myAdpSection.Value(SecParam);
 524                   }
 525                 }
 526                 AngleMax = EvalAngle(VRef, dp1);
 527                 AngleMax = M_PI/2 - AngleMax;
 528               }
 529           }
 530
 531         /*
 532            // (1.1) Distances Point-Plan
 533            Standard_Real DistPlan;
 534            gp_Vec V1 (PonPath, TheAxe.Location());
 535            DistPlan = Abs(V1.Dot(VRef));
 536
 537            // On examine l'autre extremite
 538            gp_Pnt P;
 539            Tangente(Path->Curve(), Path->LastParameter(), P, dp1);
 540            V1.SetXYZ(TheAxe.Location().XYZ()-P.XYZ());
 541            if  (Abs(V1.Dot(VRef)) <= DistPlan ) { // On prend l'autre extremite
 542            alpha =  M_PI/2 - EvalAngle(VRef, dp1);
 543            distaux =  PonPath.Distance(PonSec);
 544            if (distaux > Tol) {
 545            myExt.Perform(P);
 546            if ( myExt.IsDone() )
 547            DistMini(myExt, myAdpSection, distaux, taux);
 548            }
 549            else
 550            taux = SecParam;
 551
 552            if (Choix(distaux, alpha)) {
 553            Dist = distaux;
 554            AngleMax = alpha;
 555            PonPath = P;
 556            PathParam = Path->LastParameter();
 557            }
 558            }
 559
 560            // (1.2) Intersection Plan-courbe
 561            gp_Ax3 axe (TheAxe.Location(), TheAxe.Direction());
 562            Handle(Geom_Plane) plan = new (Geom_Plane)(axe);
 563            Handle(GeomAdaptor_HSurface) adplan =
 564            new (GeomAdaptor_HSurface)(plan);
 565            IntCurveSurface_HInter Intersector;
 566            Intersector.Perform(Path, adplan);
 567            if (Intersector.IsDone()) {
 568            Standard_Real w;
 569            gp_Vec V;
 570            for (ii=1; ii<=Intersector.NbPoints(); ii++) {
 571            w = Intersector.Point(ii).W();
 572            //(1.3) test d'angle
 573            Tangente( Path->Curve(), w, P, V);
 574            alpha = M_PI/2 - EvalAngle(V, VRef);
 575            //(1.4) Test de distance Point-Courbe
 576            myExt.Perform(P);
 577            if ( myExt.IsDone() ) {
 578            DistMini(myExt, myAdpSection, distaux, taux);
 579            if (Choix(distaux, alpha)) {
 580            Dist = distaux;
 581            SecParam  = taux;
 582            AngleMax = alpha;
 583            PonPath = P;
 584            PathParam = w;
 585            PonSec = myAdpSection.Value(SecParam);
 586            }
 587            }
 588            else {
 589            distaux = P.Distance(PonSec);
 590            if (Choix(distaux, alpha)) {
 591            Dist = distaux;
 592            AngleMax = alpha;
 593            PonPath = P;
 594            PathParam = w;
 595            }
 596            }
 597            }
 598            }
 599     */
 600 #ifdef OCCT_DEBUG
 601         if (Intersector.NbPoints() == 0) {
 602           Intersector.Dump();
 603         }
 604 #endif
 605       }
 606
 607       // Cas General
 608       if (! isplan) {
 609         // (2.1) Distance avec les extremites ...
 610         myExt.Perform(PonPath);
 611         if ( myExt.IsDone() ) {
 612           DistMini(myExt, myAdpSection, distaux, taux);
 613           if (distaux < Dist) {
 614             Dist = distaux;
 615             SecParam  = taux;
 616           }
 617         }
 618         Trouve = (Dist <= Tol);
 619         if (!Trouve) {
 620           Tangente( Path->Curve(), Path->LastParameter(), P, dp1);
 621           alpha = EvalAngle(VRef, dp1);
 622           myExt.Perform(P);
 623           if ( myExt.IsDone() ) {
 624             if ( myExt.IsDone() ) {
 625               DistMini(myExt, myAdpSection, distaux, taux);
 626               if (Choix(distaux, alpha)) {
 627                 Dist = distaux;
 628                 SecParam  = taux;
 629                 AngleMax = alpha;
 630                 PonPath = P;
 631                 PathParam = Path->LastParameter();
 632               }
 633             }
 634           }
 635           Trouve = (Dist <= Tol);
 636         }
 637
 638         // (2.2) Distance courbe-courbe
 639         if (!Trouve) {
 640           Extrema_ExtCC Ext(Path->Curve(), myAdpSection,
 641                             Path->FirstParameter(), Path->LastParameter(),
 642                             myAdpSection.FirstParameter(),
 643                             myAdpSection.LastParameter(),
 644                             Path->Resolution(Tol/100),
 645                             myAdpSection.Resolution(Tol/100));
 646           if (Ext.IsDone()) {
 647             Extrema_POnCurv P1, P2;
 648             for (ii=1; ii<=Ext.NbExt(); ii++) {
 649               distaux = sqrt (Ext.SquareDistance(ii));
 650               Ext.Points(ii, P1, P2);
 651               Tangente(Path->Curve(), P1.Parameter(), P, dp1);
 652               alpha =  EvalAngle(VRef, dp1);
 653               if (Choix(distaux, alpha)) {
 654                 Trouve = Standard_True;
 655                 Dist = distaux;
 656                 PathParam = P1.Parameter();
 657                 SecParam = P2.Parameter();
 658                 PonSec =  P2.Value();
 659                 PonPath = P;
 660                 AngleMax = alpha;
 661               }
 662             }
 663           }
 664           if (!Trouve){
 665             // Si l'on a toujours rien, on essai une distance point/path
 666             // c'est la derniere chance.
 667             Extrema_ExtPC PExt;
 668             PExt.Initialize(Path->Curve(),
 669                             Path->FirstParameter(),
 670                             Path->LastParameter(),
 671                             Precision::Confusion());
 672             PExt.Perform(PonSec);
 673             if ( PExt.IsDone() ) {
 674               // modified for OCC13595  Tue Oct 17 15:00:08 2006.BEGIN
 675               //              DistMini(PExt, myAdpSection, distaux, taux);
 676               DistMini(PExt, Path->Curve(), distaux, taux);
 677               // modified for OCC13595  Tue Oct 17 15:00:11 2006.END
 678               Tangente(Path->Curve(), taux, P, dp1);
 679               alpha = EvalAngle(VRef, dp1);
 680               if (Choix(distaux, alpha)) {
 681                 Dist = distaux;
 682                 PonPath = P;
 683                 AngleMax = alpha;
 684                 PathParam = taux;
 685               }
 686             }
 687           }
 688         }
 689       }
 690     }
 691
 692   done = Standard_True;
 693 }
 694
 695
 696 //===============================================================
 697 // Function : Perform
 698 // Purpose : Calcul le placement pour un parametre donne.
 699 //===============================================================
 700 void GeomFill_SectionPlacement::Perform(const Standard_Real Param,
 701                                         const Standard_Real Tol)
 702 {
 703   done = Standard_True;
 704   Handle(Adaptor3d_HCurve) Path;
 705   Path = myLaw->GetCurve();
 706
 707   PathParam = Param;
 708   if (myIsPoint)
 709     {
 710       gp_Pnt PonPath = Path->Value( PathParam );
 711       Dist = PonPath.Distance(myPoint);
 712       AngleMax = M_PI/2;
 713     }
 714   else
 715     {
 716       SecParam =  myAdpSection.FirstParameter();
 717
 718       //  Standard_Real distaux, taux, alpha;
 719       //  gp_Pnt PonPath, PonSec, P;
 720       gp_Pnt PonPath, PonSec;
 721       gp_Vec VRef, dp1;
 722       VRef.SetXYZ(TheAxe.Direction().XYZ());
 723
 724       Tangente( Path->Curve(), PathParam, PonPath, dp1);
 725       PonSec = myAdpSection.Value(SecParam);
 726       Dist =  PonPath.Distance(PonSec);
 727       if (Dist > Tol) { // On Cherche un meilleur point sur la section
 728         myExt.Perform(PonPath);
 729         if ( myExt.IsDone() ) {
 730           DistMini(myExt, myAdpSection, Dist, SecParam);
 731           PonSec = myAdpSection.Value(SecParam);
 732         }
 733       }
 734       AngleMax = EvalAngle(VRef, dp1);
 735       if (isplan) AngleMax = M_PI/2 - AngleMax;
 736     }
 737
 738   done = Standard_True;
 739 }
 740
 741 //===============================================================
 742 // Function :IsDone
 743 // Purpose :
 744 //===============================================================
 745  Standard_Boolean GeomFill_SectionPlacement::IsDone() const
 746 {
 747   return done;
 748 }
 749
 750 //===============================================================
 751 // Function : ParameterOnPath
 752 // Purpose :
 753 //===============================================================
 754  Standard_Real GeomFill_SectionPlacement::ParameterOnPath() const
 755 {
 756   return PathParam;
 757 }
 758
 759 //===============================================================
 760 // Function : ParameterOnSection
 761 // Purpose :
 762 //===============================================================
 763  Standard_Real GeomFill_SectionPlacement::ParameterOnSection() const
 764 {
 765  return SecParam;
 766 }
 767
 768 //===============================================================
 769 // Function : Distance
 770 // Purpose :
 771 //===============================================================
 772  Standard_Real GeomFill_SectionPlacement::Distance() const
 773 {
 774   return Dist;
 775 }
 776
 777 //===============================================================
 778 // Function : Angle
 779 // Purpose :
 780 //===============================================================
 781  Standard_Real GeomFill_SectionPlacement::Angle() const
 782 {
 783   return AngleMax;
 784 }
 785
 786 //===============================================================
 787 // Function : Transformation
 788 // Purpose :
 789 //===============================================================
 790  gp_Trsf GeomFill_SectionPlacement::Transformation(
 791          const Standard_Boolean WithTranslation,
 792          const Standard_Boolean WithCorrection) const
 793 {
 794   gp_Vec V;
 795   gp_Mat M;
 796   gp_Dir DT, DN, D;
 797 // modified by NIZHNY-MKK  Fri Oct 17 15:27:07 2003
 798   gp_Pnt P(0., 0., 0.), PSection(0., 0., 0.);
 799
 800   // Calcul des reperes
 801   myLaw->D0(PathParam, M, V);
 802
 803   P.SetXYZ(V.XYZ());
 804   D.SetXYZ (M.Column(3));
 805   DN.SetXYZ(M.Column(1));
 806   gp_Ax3 Paxe(P, D, DN);
 807
 808
 809   if (WithTranslation || WithCorrection) {
 810     if (myIsPoint)
 811       PSection = myPoint;
 812     else
 813       PSection = mySection->Value(SecParam);
 814   }
 815   // modified by NIZHNY-MKK  Wed Oct  8 15:03:19 2003.BEGIN
 816   gp_Trsf Rot;
 817
 818   if (WithCorrection && !myIsPoint) {
 819     Standard_Real angle;
 820
 821     if (!isplan)
 822       throw Standard_Failure("Illegal usage: can't rotate non-planar profile");
 823
 824     gp_Dir ProfileNormal = TheAxe.Direction();
 825     gp_Dir SpineStartDir = Paxe.Direction();
 826
 827     if (! ProfileNormal.IsParallel( SpineStartDir, Precision::Angular() )) {
 828       gp_Dir DirAxeOfRotation = ProfileNormal ^ SpineStartDir;
 829       angle = ProfileNormal.AngleWithRef( SpineStartDir, DirAxeOfRotation );
 830       gp_Ax1 AxeOfRotation( TheAxe.Location(), DirAxeOfRotation );
 831       Rot.SetRotation( AxeOfRotation, angle );
 832     }
 833     PSection.Transform(Rot);
 834   }
 835   // modified by NIZHNY-MKK  Wed Oct  8 15:03:21 2003.END
 836
 837   if (WithTranslation) {
 838     P.ChangeCoord().SetLinearForm(-1,PSection.XYZ(),
 839                                   V.XYZ());
 840   }
 841   else {
 842     P.SetCoord(0., 0., 0.);
 843   }
 844
 845   gp_Ax3 Saxe(P, gp::DZ(), gp::DX());
 846
 847   // Transfo
 848   gp_Trsf Tf;
 849   Tf.SetTransformation(Saxe, Paxe);
 850
 851   if (WithCorrection) {
 852     // modified by NIZHNY-MKK  Fri Oct 17 15:26:36 2003.BEGIN
 853     //     Standard_Real angle;
 854     //     gp_Trsf Rot;
 855
 856     //     if (!isplan)
 857     //       throw Standard_Failure("Illegal usage: can't rotate non-planar profile");
 858
 859     //     gp_Dir ProfileNormal = TheAxe.Direction();
 860     //     gp_Dir SpineStartDir = Paxe.Direction();
 861     //     if (! ProfileNormal.IsParallel( SpineStartDir, Precision::Angular() ))
 862     //       {
 863     //  gp_Dir DirAxeOfRotation = ProfileNormal ^ SpineStartDir;
 864     //  angle = ProfileNormal.AngleWithRef( SpineStartDir, DirAxeOfRotation );
 865     //  gp_Ax1 AxeOfRotation( TheAxe.Location(), DirAxeOfRotation );
 866     //  Rot.SetRotation( AxeOfRotation, angle );
 867     //       }
 868     // modified by NIZHNY-MKK  Fri Oct 17 15:26:42 2003.END
 869
 870     Tf *= Rot;
 871   }
 872
 873   return Tf;
 874 }
 875
 876 //===============================================================
 877 // Function : Section
 878 // Purpose :
 879 //===============================================================
 880  Handle(Geom_Curve) GeomFill_SectionPlacement::
 881  Section(const Standard_Boolean WithTranslation) const
 882 {
 883   Handle(Geom_Curve) TheSection =
 884     Handle(Geom_Curve)::DownCast(mySection->Copy());
 885   TheSection->Transform(Transformation(WithTranslation, Standard_False));
 886   return TheSection;
 887 }
 888
 889
 890 //===============================================================
 891 // Function :
 892 // Purpose :
 893 //===============================================================
 894  Handle(Geom_Curve) GeomFill_SectionPlacement::
 895  ModifiedSection(const Standard_Boolean WithTranslation) const
 896 {
 897   Handle(Geom_Curve) TheSection =
 898     Handle(Geom_Curve)::DownCast(mySection->Copy());
 899   TheSection->Transform(Transformation(WithTranslation, Standard_True));
 900   return TheSection;
 901 }
 902
 903 //===============================================================
 904 // Function :SectionAxis
 905 // Purpose :
 906 //===============================================================
 907  void GeomFill_SectionPlacement::SectionAxis(const gp_Mat& M,
 908                                              gp_Vec& T,
 909                                              gp_Vec& N,
 910                                              gp_Vec& BN) const
 911 {
 912   Standard_Real Eps = 1.e-10;
 913   gp_Dir D;
 914   gp_Vec PathNormal;
 915   GeomLProp_CLProps CP(mySection, SecParam, 2, Eps);
 916   if (CP.IsTangentDefined()) {
 917     CP.Tangent(D);
 918     T.SetXYZ(D.XYZ());
 919     T.Normalize();
 920     if (CP.Curvature() > Eps) {
 921       CP.Normal(D);
 922       N.SetXYZ(D.XYZ());
 923     }
 924     else {
 925       // Cas ambigu, on essai de recuperer la normal a la trajectoire
 926       // Reste le probleme des points d'inflexions, qui n'est pas
 927       // bien traiter par LProp (pas de recuperation de la normal) :
 928       // A voir...
 929       PathNormal.SetXYZ(M.Column(1));
 930       PathNormal.Normalize();
 931       BN = T ^ PathNormal;
 932       if (BN.Magnitude() > Eps) {
 933         BN.Normalize();
 934         //N = BN ^ T;
 935       }
 936       N = BN ^ T;
 937     }
 938   }
 939   else { // Cas indefinie, on prend le triedre
 940         // complet sur la trajectoire
 941     T.SetXYZ(M.Column(3));
 942     N.SetXYZ(M.Column(2));
 943   }
 944   BN = T ^ N;
 945 }
 946
 947
 948 //===============================================================
 949 // Function :Choix
 950 // Purpose : Decide si le couple (dist, angle) est "meilleur" que
 951 // couple courrant.
 952 //===============================================================
 953  Standard_Boolean GeomFill_SectionPlacement::Choix(const Standard_Real dist,
 954                                                    const Standard_Real  angle) const
 955 {
 956   Standard_Real evoldist, evolangle;
 957   evoldist = dist - Dist;
 958   evolangle = angle - AngleMax;
 959   // (1) Si la gain en distance est > que le gabarit, on prend
 960   if (evoldist < - Gabarit)
 961     return Standard_True;
 962
 963  //  (2) si l'ecart en distance est de l'ordre du gabarit
 964   if (Abs(evoldist) < Gabarit) {
 965 //  (2.1) si le gain en angle est important on garde
 966     if  (evolangle > 0.5)
 967        return Standard_True;
 968  //  (2.2) si la variation d'angle est moderee on evalue
 969  //  une fonction de penalite
 970     if (Penalite(angle, dist/Gabarit) < Penalite(AngleMax, Dist/Gabarit) )
 971       return Standard_True;
 972   }
 973
 974   return Standard_False;
 975 }
 976
 977
 978
 979
 980
 981
 982
 983