src/GeomFill/GeomFill_SectionPlacement.cxx

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