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           Standard_Real t, delta;
 287           if (myAdpSection.GetType() == GeomAbs_BSplineCurve)
 288             {
 289               Handle(Geom_BSplineCurve) BC =
 290                 Handle(Geom_BSplineCurve)::DownCast(myAdpSection.Curve());
 291               Standard_Integer I1, I2, I3, I4;
 292               BC->LocateU( first, Precision::Confusion(), I1, I2 );
 293               BC->LocateU( last,  Precision::Confusion(), I3, I4 );
 294               Standard_Integer NbKnots = I3 - I2 + 1;
 295
 296               Standard_Integer NbLocalPnts = 10;
 297               Standard_Integer NbPnts = (NbKnots-1) * NbLocalPnts;
 298               if (I1 != I2)
 299                 NbPnts += NbLocalPnts;
 300               if (I3 != I4 && first < BC->Knot(I3))
 301                 NbPnts += NbLocalPnts;
 302               if (!myAdpSection.IsClosed())
 303                 NbPnts++;
 304               Pnts = new TColgp_HArray1OfPnt(1, NbPnts);
 305               Standard_Integer nb = 1;
 306               if (I1 != I2)
 307                 {
 308                   Standard_Real locallast = (BC->Knot(I2) < last)? BC->Knot(I2) : last;
 309                   delta = (locallast - first) / NbLocalPnts;
 310                   for (j = 0; j < NbLocalPnts; j++)
 311                     {
 312                       t = first + j*delta;
 313                       Pnts->SetValue( nb++, myAdpSection.Value(t) );
 314                     }
 315                 }
 316               for (i = I2; i < I3; i++)
 317                 {
 318                   t = BC->Knot(i);
 319                   delta = (BC->Knot(i+1) - t) / NbLocalPnts;
 320                   for (j = 0; j < NbLocalPnts; j++)
 321                     {
 322                       Pnts->SetValue( nb++, myAdpSection.Value(t) );
 323                       t += delta;
 324                     }
 325                 }
 326               if (I3 != I4 && first < BC->Knot(I3))
 327                 {
 328                   t = BC->Knot(I3);
 329                   delta = (last - t) / NbLocalPnts;
 330                   for (j = 0; j < NbLocalPnts; j++)
 331                     {
 332                       Pnts->SetValue( nb++, myAdpSection.Value(t) );
 333                       t += delta;
 334                     }
 335                 }
 336               if (!myAdpSection.IsClosed())
 337                 Pnts->SetValue( nb, myAdpSection.Value(last) );
 338             }
 339           else // other type
 340             {
 341               Standard_Integer NbPnts = NbPoles-1;
 342               if (!myAdpSection.IsClosed())
 343                 NbPnts++;
 344               Pnts = new TColgp_HArray1OfPnt(1, NbPnts);
 345               delta = (last - first) / (NbPoles-1);
 346               for (i = 0; i < NbPoles-1; i++)
 347                 {
 348                   t = first + i*delta;
 349                   Pnts->SetValue( i+1, myAdpSection.Value(t) );
 350                 }
 351               if (!myAdpSection.IsClosed())
 352                 Pnts->SetValue( NbPnts, myAdpSection.Value(last) );
 353             }
 354
 355           Standard_Boolean issing;
 356           gp_Ax2 axe;
 357           GeomLib::AxeOfInertia(Pnts->Array1(), axe, issing, Precision::Confusion());
 358           if (!issing) {
 359             isplan = Standard_True;
 360             TheAxe.SetLocation(axe.Location());
 361             TheAxe.SetDirection(axe.Direction());
 362           }
 363         }
 364
 365
 366       myExt.Initialize(myAdpSection,
 367                        myAdpSection.FirstParameter(),
 368                        myAdpSection.LastParameter(),
 369                        Precision::Confusion());
 370     }
 371 }
 372
 373 //===============================================================
 374 // Function :SetLocation
 375 // Purpose :
 376 //===============================================================
 377 void GeomFill_SectionPlacement::
 378 SetLocation(const Handle(GeomFill_LocationLaw)& L)
 379 {
 380   myLaw = L;
 381 }
 382
 383 //===============================================================
 384 // Function : Perform
 385 // Purpose : Le plus simple
 386 //===============================================================
 387 void GeomFill_SectionPlacement::Perform(const Standard_Real Tol)
 388 {
 389   Handle(Adaptor3d_HCurve) Path;
 390   Path = myLaw->GetCurve();
 391   Perform(Path, Tol);
 392 }
 393
 394 //===============================================================
 395 // Function :Perform
 396 // Purpose : Recherche automatique
 397 //===============================================================
 398 void GeomFill_SectionPlacement::Perform(const Handle(Adaptor3d_HCurve)& Path,
 399                                         const Standard_Real Tol)
 400 {
 401   Standard_Real IntTol = 1.e-5;
 402   Standard_Real DistCenter = Precision::Infinite();
 403
 404   if (myIsPoint)
 405     {
 406       Extrema_ExtPC Projector(myPoint, Path->Curve(), Precision::Confusion());
 407       DistMini( Projector, Path->Curve(), Dist, PathParam );
 408       AngleMax = M_PI/2;
 409     }
 410   else
 411     {
 412       PathParam = Path->FirstParameter();
 413       SecParam =  myAdpSection.FirstParameter();
 414
 415       Standard_Real distaux, taux, alpha;
 416       gp_Pnt PonPath, PonSec, P;
 417       gp_Vec VRef, dp1;
 418       VRef.SetXYZ(TheAxe.Direction().XYZ());
 419
 420       Tangente( Path->Curve(), PathParam, PonPath, dp1);
 421       PonSec = myAdpSection.Value(SecParam);
 422       Dist =  PonPath.Distance(PonSec);
 423       if (Dist > Tol) { // On Cherche un meilleur point sur la section
 424         myExt.Perform(PonPath);
 425         if ( myExt.IsDone() ) {
 426           DistMini(myExt, myAdpSection, Dist, SecParam);
 427           PonSec = myAdpSection.Value(SecParam);
 428         }
 429       }
 430       AngleMax = EvalAngle(VRef, dp1);
 431       if (isplan) AngleMax = M_PI/2 - AngleMax;
 432
 433       Standard_Boolean Trouve = Standard_False;
 434       Standard_Integer ii;
 435
 436       if (isplan) {
 437         // (1.1) Distances Point-Plan
 438         Standard_Real DistPlan;
 439         gp_Vec V1 (PonPath, TheAxe.Location());
 440         DistPlan = Abs(V1.Dot(VRef));
 441         if (DistPlan <= IntTol)
 442           DistCenter = V1.Magnitude();
 443
 444         gp_Pnt Plast = Path->Value( Path->LastParameter() );
 445         V1.SetXYZ( TheAxe.Location().XYZ() - Plast.XYZ() );
 446         DistPlan = Abs(V1.Dot(VRef));
 447         if (DistPlan <= IntTol)
 448           {
 449             Standard_Real aDist = V1.Magnitude();
 450             if (aDist < DistCenter)
 451               {
 452                 DistCenter = aDist;
 453                 PonPath = Plast;
 454                 PathParam = Path->LastParameter();
 455               }
 456           }
 457
 458         // (1.2) Intersection Plan-courbe
 459         gp_Ax3 axe (TheAxe.Location(), TheAxe.Direction());
 460         Handle(Geom_Plane) plan = new (Geom_Plane)(axe);
 461         Handle(GeomAdaptor_HSurface) adplan =
 462           new (GeomAdaptor_HSurface)(plan);
 463         IntCurveSurface_HInter Intersector;
 464         Intersector.Perform(Path, adplan);
 465         if (Intersector.IsDone())
 466           {
 467             Standard_Real w;
 468             Standard_Real aDist;
 469             for (ii=1; ii<=Intersector.NbPoints(); ii++)
 470               {
 471                 w = Intersector.Point(ii).W();
 472                 P = Path->Value( w );
 473                 aDist = P.Distance( TheAxe.Location() );
 474                 if (aDist < DistCenter)
 475                   {
 476                     DistCenter = aDist;
 477                     PonPath = P;
 478                     PathParam = w;
 479                   }
 480               }
 481           }
 482         if (!Intersector.IsDone() || Intersector.NbPoints() == 0)
 483           {
 484             Standard_Integer NbPnts = NbSamples( mySection );
 485             TColgp_Array1OfPnt Pnts( 1, NbPnts+1 );
 486             Standard_Real delta = (mySection->LastParameter()-mySection->FirstParameter())/NbPnts;
 487             for (ii = 0; ii <= NbPnts; ii++)
 488               Pnts(ii+1) = mySection->Value( mySection->FirstParameter() + ii*delta );
 489
 490             gp_Pnt BaryCenter;
 491             gp_Dir Xdir, Ydir;
 492             Standard_Real Xgap, Ygap, Zgap;
 493             GeomLib::Inertia( Pnts, BaryCenter, Xdir, Ydir, Xgap, Ygap, Zgap );
 494
 495             gp_Pnt Pfirst = Path->Value( Path->FirstParameter() );
 496             if (Pfirst.Distance(BaryCenter) < Plast.Distance(BaryCenter))
 497               PathParam = Path->FirstParameter();
 498             else
 499               {
 500                 PathParam = Path->LastParameter();
 501                 Tangente( Path->Curve(), PathParam, PonPath, dp1);
 502                 PonSec = myAdpSection.Value(SecParam);
 503                 Dist =  PonPath.Distance(PonSec);
 504                 if (Dist > Tol) { // On Cherche un meilleur point sur la section
 505                   myExt.Perform(PonPath);
 506                   if ( myExt.IsDone() ) {
 507                     DistMini(myExt, myAdpSection, Dist, SecParam);
 508                     PonSec = myAdpSection.Value(SecParam);
 509                   }
 510                 }
 511                 AngleMax = EvalAngle(VRef, dp1);
 512                 AngleMax = M_PI/2 - AngleMax;
 513               }
 514           }
 515
 516         /*
 517            // (1.1) Distances Point-Plan
 518            Standard_Real DistPlan;
 519            gp_Vec V1 (PonPath, TheAxe.Location());
 520            DistPlan = Abs(V1.Dot(VRef));
 521
 522            // On examine l'autre extremite
 523            gp_Pnt P;
 524            Tangente(Path->Curve(), Path->LastParameter(), P, dp1);
 525            V1.SetXYZ(TheAxe.Location().XYZ()-P.XYZ());
 526            if  (Abs(V1.Dot(VRef)) <= DistPlan ) { // On prend l'autre extremite
 527            alpha =  M_PI/2 - EvalAngle(VRef, dp1);
 528            distaux =  PonPath.Distance(PonSec);
 529            if (distaux > Tol) {
 530            myExt.Perform(P);
 531            if ( myExt.IsDone() )
 532            DistMini(myExt, myAdpSection, distaux, taux);
 533            }
 534            else
 535            taux = SecParam;
 536
 537            if (Choix(distaux, alpha)) {
 538            Dist = distaux;
 539            AngleMax = alpha;
 540            PonPath = P;
 541            PathParam = Path->LastParameter();
 542            }
 543            }
 544
 545            // (1.2) Intersection Plan-courbe
 546            gp_Ax3 axe (TheAxe.Location(), TheAxe.Direction());
 547            Handle(Geom_Plane) plan = new (Geom_Plane)(axe);
 548            Handle(GeomAdaptor_HSurface) adplan =
 549            new (GeomAdaptor_HSurface)(plan);
 550            IntCurveSurface_HInter Intersector;
 551            Intersector.Perform(Path, adplan);
 552            if (Intersector.IsDone()) {
 553            Standard_Real w;
 554            gp_Vec V;
 555            for (ii=1; ii<=Intersector.NbPoints(); ii++) {
 556            w = Intersector.Point(ii).W();
 557            //(1.3) test d'angle
 558            Tangente( Path->Curve(), w, P, V);
 559            alpha = M_PI/2 - EvalAngle(V, VRef);
 560            //(1.4) Test de distance Point-Courbe
 561            myExt.Perform(P);
 562            if ( myExt.IsDone() ) {
 563            DistMini(myExt, myAdpSection, distaux, taux);
 564            if (Choix(distaux, alpha)) {
 565            Dist = distaux;
 566            SecParam  = taux;
 567            AngleMax = alpha;
 568            PonPath = P;
 569            PathParam = w;
 570            PonSec = myAdpSection.Value(SecParam);
 571            }
 572            }
 573            else {
 574            distaux = P.Distance(PonSec);
 575            if (Choix(distaux, alpha)) {
 576            Dist = distaux;
 577            AngleMax = alpha;
 578            PonPath = P;
 579            PathParam = w;
 580            }
 581            }
 582            }
 583            }
 584     */
 585 #ifdef OCCT_DEBUG
 586         if (Intersector.NbPoints() == 0) {
 587           Intersector.Dump();
 588         }
 589 #endif
 590       }
 591
 592       // Cas General
 593       if (! isplan) {
 594         // (2.1) Distance avec les extremites ...
 595         myExt.Perform(PonPath);
 596         if ( myExt.IsDone() ) {
 597           DistMini(myExt, myAdpSection, distaux, taux);
 598           if (distaux < Dist) {
 599             Dist = distaux;
 600             SecParam  = taux;
 601           }
 602         }
 603         Trouve = (Dist <= Tol);
 604         if (!Trouve) {
 605           Tangente( Path->Curve(), Path->LastParameter(), P, dp1);
 606           alpha = EvalAngle(VRef, dp1);
 607           myExt.Perform(P);
 608           if ( myExt.IsDone() ) {
 609             if ( myExt.IsDone() ) {
 610               DistMini(myExt, myAdpSection, distaux, taux);
 611               if (Choix(distaux, alpha)) {
 612                 Dist = distaux;
 613                 SecParam  = taux;
 614                 AngleMax = alpha;
 615                 PonPath = P;
 616                 PathParam = Path->LastParameter();
 617               }
 618             }
 619           }
 620           Trouve = (Dist <= Tol);
 621         }
 622
 623         // (2.2) Distance courbe-courbe
 624         if (!Trouve) {
 625           Extrema_ExtCC Ext(Path->Curve(), myAdpSection,
 626                             Path->FirstParameter(), Path->LastParameter(),
 627                             myAdpSection.FirstParameter(),
 628                             myAdpSection.LastParameter(),
 629                             Path->Resolution(Tol/100),
 630                             myAdpSection.Resolution(Tol/100));
 631           if (Ext.IsDone()) {
 632             Extrema_POnCurv P1, P2;
 633             for (ii=1; ii<=Ext.NbExt(); ii++) {
 634               distaux = sqrt (Ext.SquareDistance(ii));
 635               Ext.Points(ii, P1, P2);
 636               Tangente(Path->Curve(), P1.Parameter(), P, dp1);
 637               alpha =  EvalAngle(VRef, dp1);
 638               if (Choix(distaux, alpha)) {
 639                 Trouve = Standard_True;
 640                 Dist = distaux;
 641                 PathParam = P1.Parameter();
 642                 SecParam = P2.Parameter();
 643                 PonSec =  P2.Value();
 644                 PonPath = P;
 645                 AngleMax = alpha;
 646               }
 647             }
 648           }
 649           if (!Trouve){
 650             // Si l'on a toujours rien, on essai une distance point/path
 651             // c'est la derniere chance.
 652             Extrema_ExtPC PExt;
 653             PExt.Initialize(Path->Curve(),
 654                             Path->FirstParameter(),
 655                             Path->LastParameter(),
 656                             Precision::Confusion());
 657             PExt.Perform(PonSec);
 658             if ( PExt.IsDone() ) {
 659               // modified for OCC13595  Tue Oct 17 15:00:08 2006.BEGIN
 660               //              DistMini(PExt, myAdpSection, distaux, taux);
 661               DistMini(PExt, Path->Curve(), distaux, taux);
 662               // modified for OCC13595  Tue Oct 17 15:00:11 2006.END
 663               Tangente(Path->Curve(), taux, P, dp1);
 664               alpha = EvalAngle(VRef, dp1);
 665               if (Choix(distaux, alpha)) {
 666                 Dist = distaux;
 667                 PonPath = P;
 668                 AngleMax = alpha;
 669                 PathParam = taux;
 670               }
 671             }
 672           }
 673         }
 674       }
 675     }
 676
 677   done = Standard_True;
 678 }
 679
 680
 681 //===============================================================
 682 // Function : Perform
 683 // Purpose : Calcul le placement pour un parametre donne.
 684 //===============================================================
 685 void GeomFill_SectionPlacement::Perform(const Standard_Real Param,
 686                                         const Standard_Real Tol)
 687 {
 688   done = Standard_True;
 689   Handle(Adaptor3d_HCurve) Path;
 690   Path = myLaw->GetCurve();
 691
 692   PathParam = Param;
 693   if (myIsPoint)
 694     {
 695       gp_Pnt PonPath = Path->Value( PathParam );
 696       Dist = PonPath.Distance(myPoint);
 697       AngleMax = M_PI/2;
 698     }
 699   else
 700     {
 701       SecParam =  myAdpSection.FirstParameter();
 702
 703       //  Standard_Real distaux, taux, alpha;
 704       //  gp_Pnt PonPath, PonSec, P;
 705       gp_Pnt PonPath, PonSec;
 706       gp_Vec VRef, dp1;
 707       VRef.SetXYZ(TheAxe.Direction().XYZ());
 708
 709       Tangente( Path->Curve(), PathParam, PonPath, dp1);
 710       PonSec = myAdpSection.Value(SecParam);
 711       Dist =  PonPath.Distance(PonSec);
 712       if (Dist > Tol) { // On Cherche un meilleur point sur la section
 713         myExt.Perform(PonPath);
 714         if ( myExt.IsDone() ) {
 715           DistMini(myExt, myAdpSection, Dist, SecParam);
 716           PonSec = myAdpSection.Value(SecParam);
 717         }
 718       }
 719       AngleMax = EvalAngle(VRef, dp1);
 720       if (isplan) AngleMax = M_PI/2 - AngleMax;
 721     }
 722
 723   done = Standard_True;
 724 }
 725
 726 //===============================================================
 727 // Function :IsDone
 728 // Purpose :
 729 //===============================================================
 730  Standard_Boolean GeomFill_SectionPlacement::IsDone() const
 731 {
 732   return done;
 733 }
 734
 735 //===============================================================
 736 // Function : ParameterOnPath
 737 // Purpose :
 738 //===============================================================
 739  Standard_Real GeomFill_SectionPlacement::ParameterOnPath() const
 740 {
 741   return PathParam;
 742 }
 743
 744 //===============================================================
 745 // Function : ParameterOnSection
 746 // Purpose :
 747 //===============================================================
 748  Standard_Real GeomFill_SectionPlacement::ParameterOnSection() const
 749 {
 750  return SecParam;
 751 }
 752
 753 //===============================================================
 754 // Function : Distance
 755 // Purpose :
 756 //===============================================================
 757  Standard_Real GeomFill_SectionPlacement::Distance() const
 758 {
 759   return Dist;
 760 }
 761
 762 //===============================================================
 763 // Function : Angle
 764 // Purpose :
 765 //===============================================================
 766  Standard_Real GeomFill_SectionPlacement::Angle() const
 767 {
 768   return AngleMax;
 769 }
 770
 771 //===============================================================
 772 // Function : Transformation
 773 // Purpose :
 774 //===============================================================
 775  gp_Trsf GeomFill_SectionPlacement::Transformation(
 776          const Standard_Boolean WithTranslation,
 777          const Standard_Boolean WithCorrection) const
 778 {
 779   gp_Vec V;
 780   gp_Mat M;
 781   gp_Dir DT, DN, D;
 782 // modified by NIZHNY-MKK  Fri Oct 17 15:27:07 2003
 783   gp_Pnt P(0., 0., 0.), PSection(0., 0., 0.);
 784
 785   // Calcul des reperes
 786   myLaw->D0(PathParam, M, V);
 787
 788   P.SetXYZ(V.XYZ());
 789   D.SetXYZ (M.Column(3));
 790   DN.SetXYZ(M.Column(1));
 791   gp_Ax3 Paxe(P, D, DN);
 792
 793
 794   if (WithTranslation || WithCorrection) {
 795     if (myIsPoint)
 796       PSection = myPoint;
 797     else
 798       PSection = mySection->Value(SecParam);
 799   }
 800   // modified by NIZHNY-MKK  Wed Oct  8 15:03:19 2003.BEGIN
 801   gp_Trsf Rot;
 802
 803   if (WithCorrection && !myIsPoint) {
 804     Standard_Real angle;
 805
 806     if (!isplan)
 807       Standard_Failure::Raise("Illegal usage: can't rotate non-planar profile");
 808
 809     gp_Dir ProfileNormal = TheAxe.Direction();
 810     gp_Dir SpineStartDir = Paxe.Direction();
 811
 812     if (! ProfileNormal.IsParallel( SpineStartDir, Precision::Angular() )) {
 813       gp_Dir DirAxeOfRotation = ProfileNormal ^ SpineStartDir;
 814       angle = ProfileNormal.AngleWithRef( SpineStartDir, DirAxeOfRotation );
 815       gp_Ax1 AxeOfRotation( TheAxe.Location(), DirAxeOfRotation );
 816       Rot.SetRotation( AxeOfRotation, angle );
 817     }
 818     PSection.Transform(Rot);
 819   }
 820   // modified by NIZHNY-MKK  Wed Oct  8 15:03:21 2003.END
 821
 822   if (WithTranslation) {
 823     P.ChangeCoord().SetLinearForm(-1,PSection.XYZ(),
 824                                   V.XYZ());
 825   }
 826   else {
 827     P.SetCoord(0., 0., 0.);
 828   }
 829
 830   gp_Ax3 Saxe(P, gp::DZ(), gp::DX());
 831
 832   // Transfo
 833   gp_Trsf Tf;
 834   Tf.SetTransformation(Saxe, Paxe);
 835
 836   if (WithCorrection) {
 837     // modified by NIZHNY-MKK  Fri Oct 17 15:26:36 2003.BEGIN
 838     //     Standard_Real angle;
 839     //     gp_Trsf Rot;
 840
 841     //     if (!isplan)
 842     //       Standard_Failure::Raise("Illegal usage: can't rotate non-planar profile");
 843
 844     //     gp_Dir ProfileNormal = TheAxe.Direction();
 845     //     gp_Dir SpineStartDir = Paxe.Direction();
 846     //     if (! ProfileNormal.IsParallel( SpineStartDir, Precision::Angular() ))
 847     //       {
 848     //  gp_Dir DirAxeOfRotation = ProfileNormal ^ SpineStartDir;
 849     //  angle = ProfileNormal.AngleWithRef( SpineStartDir, DirAxeOfRotation );
 850     //  gp_Ax1 AxeOfRotation( TheAxe.Location(), DirAxeOfRotation );
 851     //  Rot.SetRotation( AxeOfRotation, angle );
 852     //       }
 853     // modified by NIZHNY-MKK  Fri Oct 17 15:26:42 2003.END
 854
 855     Tf *= Rot;
 856   }
 857
 858   return Tf;
 859 }
 860
 861 //===============================================================
 862 // Function : Section
 863 // Purpose :
 864 //===============================================================
 865  Handle(Geom_Curve) GeomFill_SectionPlacement::
 866  Section(const Standard_Boolean WithTranslation) const
 867 {
 868   Handle(Geom_Curve) TheSection =
 869     Handle(Geom_Curve)::DownCast(mySection->Copy());
 870   TheSection->Transform(Transformation(WithTranslation, Standard_False));
 871   return TheSection;
 872 }
 873
 874
 875 //===============================================================
 876 // Function :
 877 // Purpose :
 878 //===============================================================
 879  Handle(Geom_Curve) GeomFill_SectionPlacement::
 880  ModifiedSection(const Standard_Boolean WithTranslation) const
 881 {
 882   Handle(Geom_Curve) TheSection =
 883     Handle(Geom_Curve)::DownCast(mySection->Copy());
 884   TheSection->Transform(Transformation(WithTranslation, Standard_True));
 885   return TheSection;
 886 }
 887
 888 //===============================================================
 889 // Function :SectionAxis
 890 // Purpose :
 891 //===============================================================
 892  void GeomFill_SectionPlacement::SectionAxis(const gp_Mat& M,
 893                                              gp_Vec& T,
 894                                              gp_Vec& N,
 895                                              gp_Vec& BN) const
 896 {
 897   Standard_Real Eps = 1.e-10;
 898   gp_Dir D;
 899   gp_Vec PathNormal;
 900   GeomLProp_CLProps CP(mySection, SecParam, 2, Eps);
 901   if (CP.IsTangentDefined()) {
 902     CP.Tangent(D);
 903     T.SetXYZ(D.XYZ());
 904     T.Normalize();
 905     if (CP.Curvature() > Eps) {
 906       CP.Normal(D);
 907       N.SetXYZ(D.XYZ());
 908     }
 909     else {
 910       // Cas ambigu, on essai de recuperer la normal a la trajectoire
 911       // Reste le probleme des points d'inflexions, qui n'est pas
 912       // bien traiter par LProp (pas de recuperation de la normal) :
 913       // A voir...
 914       PathNormal.SetXYZ(M.Column(1));
 915       PathNormal.Normalize();
 916       BN = T ^ PathNormal;
 917       if (BN.Magnitude() > Eps) {
 918         BN.Normalize();
 919         //N = BN ^ T;
 920       }
 921       N = BN ^ T;
 922     }
 923   }
 924   else { // Cas indefinie, on prend le triedre
 925         // complet sur la trajectoire
 926     T.SetXYZ(M.Column(3));
 927     N.SetXYZ(M.Column(2));
 928   }
 929   BN = T ^ N;
 930 }
 931
 932
 933 //===============================================================
 934 // Function :Choix
 935 // Purpose : Decide si le couple (dist, angle) est "meilleur" que
 936 // couple courrant.
 937 //===============================================================
 938  Standard_Boolean GeomFill_SectionPlacement::Choix(const Standard_Real dist,
 939                                                    const Standard_Real  angle) const
 940 {
 941   Standard_Real evoldist, evolangle;
 942   evoldist = dist - Dist;
 943   evolangle = angle - AngleMax;
 944   // (1) Si la gain en distance est > que le gabarit, on prend
 945   if (evoldist < - Gabarit)
 946     return Standard_True;
 947
 948  //  (2) si l'ecart en distance est de l'ordre du gabarit
 949   if (Abs(evoldist) < Gabarit) {
 950 //  (2.1) si le gain en angle est important on garde
 951     if  (evolangle > 0.5)
 952        return Standard_True;
 953  //  (2.2) si la variation d'angle est moderee on evalue
 954  //  une fonction de penalite
 955     if (Penalite(angle, dist/Gabarit) < Penalite(AngleMax, Dist/Gabarit) )
 956       return Standard_True;
 957   }
 958
 959   return Standard_False;
 960 }
 961
 962
 963
 964
 965
 966
 967
 968