src/Extrema/Extrema_ExtCS.cxx

   1 // Created on: 1995-07-18
   2 // Created by: Modelistation
   3 // Copyright (c) 1995-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 //  Modified by skv - Thu Jul  7 12:29:34 2005 OCC9134
  18
  19 #include <Adaptor3d_Curve.hxx>
  20 #include <Adaptor3d_Surface.hxx>
  21 #include <Bnd_Box.hxx>
  22 #include <BndLib_AddSurface.hxx>
  23 #include <ElCLib.hxx>
  24 #include <ElSLib.hxx>
  25 #include <Extrema_ExtCS.hxx>
  26 #include <Extrema_ExtPElC.hxx>
  27 #include <Extrema_ExtPElS.hxx>
  28 #include <Extrema_ExtPS.hxx>
  29 #include <Extrema_GenExtCS.hxx>
  30 #include <Extrema_POnCurv.hxx>
  31 #include <Extrema_POnSurf.hxx>
  32 #include <GeomAbs_CurveType.hxx>
  33 #include <gp_Cone.hxx>
  34 #include <gp_Cylinder.hxx>
  35 #include <gp_Lin.hxx>
  36 #include <gp_Pln.hxx>
  37 #include <gp_Pnt.hxx>
  38 #include <gp_Sphere.hxx>
  39 #include <gp_Torus.hxx>
  40 #include <Precision.hxx>
  41 #include <Standard_NotImplemented.hxx>
  42 #include <Standard_OutOfRange.hxx>
  43 #include <Standard_TypeMismatch.hxx>
  44 #include <StdFail_InfiniteSolutions.hxx>
  45 #include <StdFail_NotDone.hxx>
  46 #include <TColStd_Array1OfReal.hxx>
  47
  48 Extrema_ExtCS::Extrema_ExtCS()
  49 {
  50   myDone = Standard_False;
  51 }
  52
  53 Extrema_ExtCS::Extrema_ExtCS(const Adaptor3d_Curve&   C,
  54   const Adaptor3d_Surface& S,
  55   const Standard_Real    TolC,
  56   const Standard_Real    TolS)
  57
  58 {
  59   Initialize(S, S.FirstUParameter(), S.LastUParameter(),
  60     S.FirstVParameter(), S.LastVParameter(),
  61     TolC, TolS);
  62   Perform(C, C.FirstParameter(), C.LastParameter());
  63 }
  64
  65 Extrema_ExtCS::Extrema_ExtCS(const Adaptor3d_Curve&   C,
  66   const Adaptor3d_Surface& S,
  67   const Standard_Real    UCinf,
  68   const Standard_Real    UCsup,
  69   const Standard_Real    Uinf,
  70   const Standard_Real    Usup,
  71   const Standard_Real    Vinf,
  72   const Standard_Real    Vsup,
  73   const Standard_Real    TolC,
  74   const Standard_Real    TolS)
  75
  76 {
  77   Initialize(S, Uinf, Usup, Vinf, Vsup, TolC, TolS);
  78   Perform(C, UCinf, UCsup);
  79 }
  80
  81
  82 void Extrema_ExtCS::Initialize(const Adaptor3d_Surface& S,
  83   const Standard_Real    Uinf,
  84   const Standard_Real    Usup,
  85   const Standard_Real    Vinf,
  86   const Standard_Real    Vsup,
  87   const Standard_Real    TolC,
  88   const Standard_Real    TolS)
  89 {
  90   myS = (Adaptor3d_SurfacePtr)&S;
  91   myIsPar = Standard_False;
  92   myuinf  = Uinf;
  93   myusup  = Usup;
  94   myvinf  = Vinf;
  95   myvsup  = Vsup;
  96   mytolC  = TolC;
  97   mytolS  = TolS;
  98   myStype  = myS->GetType();
  99 }
 100
 101
 102 void Extrema_ExtCS::Perform(const Adaptor3d_Curve& C,
 103   const Standard_Real  Uinf,
 104   const Standard_Real  Usup)
 105 {
 106   myucinf = Uinf;
 107   myucsup = Usup;
 108   myPOnS.Clear();
 109   myPOnC.Clear();
 110   mySqDist.Clear();
 111   Standard_Integer i, j;
 112   Standard_Integer NbT, NbU, NbV;
 113   NbT = NbU = NbV = 10;
 114   GeomAbs_CurveType myCtype  = C.GetType();
 115
 116   myDone = Standard_False;
 117   // Try analytic computation of extrema
 118   Standard_Boolean isComputeAnalytic = Standard_True;
 119
 120   switch(myCtype) {
 121
 122   case GeomAbs_Line:
 123     {
 124
 125       switch(myStype) {
 126       case GeomAbs_Sphere:
 127         myExtElCS.Perform(C.Line(), myS->Sphere());
 128         break;
 129       case GeomAbs_Cylinder:
 130         myExtElCS.Perform(C.Line(), myS->Cylinder());
 131         break;
 132       case GeomAbs_Plane:
 133         myExtElCS.Perform(C.Line(), myS->Plane());
 134         if (myExtElCS.IsParallel())   break;
 135         Standard_FALLTHROUGH
 136
 137       case GeomAbs_Torus:
 138       case GeomAbs_Cone:
 139       case GeomAbs_BezierSurface:
 140       case GeomAbs_BSplineSurface:
 141       case GeomAbs_SurfaceOfRevolution:
 142       case GeomAbs_SurfaceOfExtrusion:
 143       case GeomAbs_OffsetSurface:
 144       case GeomAbs_OtherSurface:
 145         {
 146           Standard_Real cfirst = myucinf, clast = myucsup;
 147           Standard_Real ufirst = myS->FirstUParameter(), ulast = myS->LastUParameter(),
 148             vfirst = myS->FirstVParameter(), vlast = myS->LastVParameter();
 149
 150           if(Precision::IsInfinite(Abs(cfirst)) || Precision::IsInfinite(Abs(clast))) {
 151
 152             Bnd_Box aSurfBox;
 153             BndLib_AddSurface::Add(*myS, ufirst, ulast, vfirst, vlast, Precision::Confusion(), aSurfBox);
 154             Standard_Real xmin, ymin, zmin, xmax, ymax, zmax;
 155             aSurfBox.Get(xmin, ymin, zmin, xmax, ymax, zmax);
 156             Standard_Real tmin = Precision::Infinite(), tmax = -tmin;
 157             gp_Lin aLin = C.Line();
 158
 159
 160             if(!( Precision::IsInfinite(Abs(xmin)) || Precision::IsInfinite(Abs(xmax)) ||
 161               Precision::IsInfinite(Abs(ymin)) || Precision::IsInfinite(Abs(ymax)) ||
 162               Precision::IsInfinite(Abs(zmin)) || Precision::IsInfinite(Abs(zmax)))  ) {
 163
 164                 Extrema_ExtPElC anExt;
 165                 Extrema_POnCurv aPntOnLin;
 166                 Standard_Real aParOnLin;
 167                 Standard_Real lim = Precision::Infinite();
 168                 gp_Pnt aLimPntArray[8];
 169
 170                 aLimPntArray[0].SetCoord(xmin, ymin, zmin);
 171                 aLimPntArray[1].SetCoord(xmax, ymin, zmin);
 172                 aLimPntArray[2].SetCoord(xmin, ymax, zmin);
 173                 aLimPntArray[3].SetCoord(xmax, ymax, zmin);
 174                 aLimPntArray[4].SetCoord(xmin, ymin, zmax);
 175                 aLimPntArray[5].SetCoord(xmax, ymin, zmax);
 176                 aLimPntArray[6].SetCoord(xmin, ymax, zmax);
 177                 aLimPntArray[7].SetCoord(xmax, ymax, zmax);
 178
 179                 for(i = 0; i <= 7; i++) {
 180                   anExt.Perform(aLimPntArray[i], aLin, Precision::Confusion(), -lim, lim);
 181                   aPntOnLin = anExt.Point(1);
 182                   aParOnLin = aPntOnLin.Parameter();
 183                   tmin = Min(aParOnLin, tmin);
 184                   tmax = Max(aParOnLin, tmax);
 185                 }
 186
 187             }
 188             else {
 189               tmin = -1.e+50;
 190               tmax =  1.e+50;
 191             }
 192
 193
 194             cfirst = Max(cfirst, tmin);
 195             clast  = Min(clast,  tmax);
 196
 197           }
 198
 199           if (myS->IsUPeriodic())
 200             NbU = 13;
 201           if (myS->IsVPeriodic())
 202             NbV = 13;
 203
 204           Extrema_GenExtCS Ext(C, *myS, NbT, NbU, NbV, cfirst, clast, ufirst, ulast,
 205             vfirst, vlast, mytolC, mytolS);
 206
 207           myDone = Ext.IsDone();
 208           if (myDone) {
 209             Standard_Integer NbExt = Ext.NbExt();
 210             Standard_Real T,U,V;
 211             Extrema_POnCurv PC;
 212             Extrema_POnSurf PS;
 213             for (i = 1; i <= NbExt; i++) {
 214               PC = Ext.PointOnCurve(i);
 215               PS = Ext.PointOnSurface(i);
 216               T = PC.Parameter();
 217               PS.Parameter(U, V);
 218               AddSolution(C, T, U, V, PC.Value(), PS.Value(), Ext.SquareDistance(i));
 219             }
 220           }
 221           return;
 222
 223         }
 224       }
 225       break;
 226     }
 227     //  Modified by skv - Thu Jul  7 12:29:34 2005 OCC9134 Begin
 228   case GeomAbs_Circle:
 229     {
 230       if(myStype == GeomAbs_Cylinder) {
 231         myExtElCS.Perform(C.Circle(), myS->Cylinder());
 232         break;
 233       }
 234       else if(myStype == GeomAbs_Plane)
 235       {
 236         myExtElCS.Perform(C.Circle(), myS->Plane());
 237         break;
 238       }
 239       else if (myStype == GeomAbs_Sphere)
 240       {
 241         myExtElCS.Perform(C.Circle(), myS->Sphere());
 242         break;
 243       }
 244     }
 245     Standard_FALLTHROUGH
 246   case GeomAbs_Hyperbola:
 247     {
 248       if(myCtype == GeomAbs_Hyperbola && myStype == GeomAbs_Plane) {
 249         //  Modified by skv - Thu Jul  7 12:29:34 2005 OCC9134 End
 250         myExtElCS.Perform(C.Hyperbola(), myS->Plane());
 251         break;
 252       }
 253     }
 254     Standard_FALLTHROUGH
 255   default:
 256     {
 257       isComputeAnalytic = Standard_False;
 258       break;
 259     }
 260   }
 261
 262   if (isComputeAnalytic)
 263   {
 264     if (myExtElCS.IsDone())
 265     {
 266       myDone = Standard_True;
 267       myIsPar = myExtElCS.IsParallel();
 268       if (myIsPar)
 269       {
 270         mySqDist.Append(myExtElCS.SquareDistance(1));
 271       }
 272       else
 273       {
 274         Standard_Integer NbExt = myExtElCS.NbExt();
 275         for (i = 1; i <= NbExt; i++)
 276         {
 277           Extrema_POnCurv PC;
 278           Extrema_POnSurf PS;
 279           myExtElCS.Points(i, PC, PS);
 280           Standard_Real Ucurve = PC.Parameter();
 281           Standard_Real U, V;
 282           PS.Parameter(U, V);
 283           AddSolution(C, Ucurve, U, V, PC.Value(), PS.Value(), myExtElCS.SquareDistance(i));
 284         }
 285
 286         if (mySqDist.Length() == 0 && NbExt > 0)
 287         {
 288           // Analytical extrema seem to be out of curve/surface boundaries.
 289           // Try extremity points of curve.
 290           gp_Pnt aPOnC[2], aPOnS[2];
 291           Standard_Real aT[2] = { myucinf, myucsup }, U[2], V[2];
 292           Standard_Real aDist[2] = { -1, -1 };
 293           for (i = 0; i < 2; ++i)
 294           {
 295             if (Precision::IsInfinite(aT[i]))
 296               continue;
 297
 298             aPOnC[i] = C.Value(aT[i]);
 299             switch (myStype)
 300             {
 301               case GeomAbs_Plane:
 302               {
 303                 ElSLib::Parameters(myS->Plane(), aPOnC[i], U[i], V[i]);
 304                 aPOnS[i] = ElSLib::Value(U[i], V[i], myS->Plane());
 305                 break;
 306               }
 307               case GeomAbs_Sphere:
 308               {
 309                 ElSLib::Parameters(myS->Sphere(), aPOnC[i], U[i], V[i]);
 310                 aPOnS[i] = ElSLib::Value(U[i], V[i], myS->Sphere());
 311                 break;
 312               }
 313               case GeomAbs_Cylinder:
 314               {
 315                 ElSLib::Parameters(myS->Cylinder(), aPOnC[i], U[i], V[i]);
 316                 aPOnS[i] = ElSLib::Value(U[i], V[i], myS->Cylinder());
 317                 break;
 318               }
 319               case GeomAbs_Torus:
 320               {
 321                 ElSLib::Parameters(myS->Torus(), aPOnC[i], U[i], V[i]);
 322                 aPOnS[i] = ElSLib::Value(U[i], V[i], myS->Torus());
 323                 break;
 324               }
 325               case GeomAbs_Cone:
 326               {
 327                 ElSLib::Parameters(myS->Cone(), aPOnC[i], U[i], V[i]);
 328                 aPOnS[i] = ElSLib::Value(U[i], V[i], myS->Cone());
 329                 break;
 330               }
 331               default:
 332                 continue;
 333             }
 334
 335             aDist[i] = aPOnC[i].SquareDistance(aPOnS[i]);
 336           }
 337
 338           Standard_Boolean bAdd[2] = {Standard_False, Standard_False};
 339
 340           // Choose solution to add
 341           if (aDist[0] >= 0. && aDist[1] >= 0.)
 342           {
 343             Standard_Real aDiff = aDist[0] - aDist[1];
 344             // Both computed -> take only minimal
 345             if (Abs(aDiff) < Precision::Confusion())
 346               // Add both
 347               bAdd[0] = bAdd[1] = Standard_True;
 348             else if (aDiff < 0)
 349               // Add first
 350               bAdd[0] = Standard_True;
 351             else
 352               // Add second
 353               bAdd[1] = Standard_True;
 354           }
 355           else if (aDist[0] >= 0.)
 356             // Add first
 357             bAdd[0] = Standard_True;
 358           else if (aDist[1] >= 0.)
 359             // Add second
 360             bAdd[1] = Standard_True;
 361
 362           for (i = 0; i < 2; ++i)
 363           {
 364             if (bAdd[i])
 365               AddSolution(C, aT[i], U[i], V[i], aPOnC[i], aPOnS[i], aDist[i]);
 366           }
 367         }
 368       }
 369       return;
 370     }
 371   }
 372
 373   // Elementary extrema is not done, try generic solution
 374   Extrema_GenExtCS Ext;
 375   Ext.Initialize(*myS, NbU, NbV, mytolS);
 376   if (myCtype == GeomAbs_Hyperbola) {
 377     Standard_Real tmin = Max(-20., C.FirstParameter());
 378     Standard_Real tmax = Min(20., C.LastParameter());
 379     Ext.Perform(C, NbT, tmin, tmax, mytolC); // to avoid overflow
 380   }
 381   else {
 382     if ((myCtype == GeomAbs_Circle       && NbT < 13) ||
 383       (myCtype == GeomAbs_BSplineCurve && NbT < 13))
 384     {
 385       NbT = 13;
 386     }
 387     Ext.Perform(C, NbT, mytolC);
 388   }
 389
 390   myDone = Ext.IsDone();
 391   if (myDone) {
 392     Standard_Integer NbExt = Ext.NbExt();
 393     Standard_Real T, U, V;
 394     Extrema_POnCurv PC;
 395     Extrema_POnSurf PS;
 396     for (i = 1; i <= NbExt; i++) {
 397       PC = Ext.PointOnCurve(i);
 398       PS = Ext.PointOnSurface(i);
 399       T = PC.Parameter();
 400       PS.Parameter(U, V);
 401       AddSolution(C, T, U, V, PC.Value(), PS.Value(), Ext.SquareDistance(i));
 402     }
 403
 404     //Add sharp points
 405     Standard_Integer SolNumber = mySqDist.Length();
 406     Standard_Address CopyC = (Standard_Address)&C;
 407     Adaptor3d_Curve& aC = *(Adaptor3d_Curve*)CopyC;
 408     Standard_Integer NbIntervals = aC.NbIntervals(GeomAbs_C1);
 409     TColStd_Array1OfReal SharpPoints(1, NbIntervals + 1);
 410     aC.Intervals(SharpPoints, GeomAbs_C1);
 411
 412     Extrema_ExtPS aProjPS;
 413     aProjPS.Initialize(*myS,
 414       myS->FirstUParameter(),
 415       myS->LastUParameter(),
 416       myS->FirstVParameter(),
 417       myS->LastVParameter(),
 418       mytolS,
 419       mytolS);
 420
 421     for (i = 2; i < SharpPoints.Upper(); ++i)
 422     {
 423       T = SharpPoints(i);
 424       gp_Pnt aPnt = C.Value(T);
 425       aProjPS.Perform(aPnt);
 426       if (!aProjPS.IsDone())
 427         continue;
 428       Standard_Integer NbProj = aProjPS.NbExt(), jmin = 0;
 429       Standard_Real MinSqDist = RealLast();
 430       for (j = 1; j <= NbProj; j++)
 431       {
 432         Standard_Real aSqDist = aProjPS.SquareDistance(j);
 433         if (aSqDist < MinSqDist)
 434         {
 435           MinSqDist = aSqDist;
 436           jmin = j;
 437         }
 438       }
 439       if (jmin != 0)
 440       {
 441         aProjPS.Point(jmin).Parameter(U, V);
 442         AddSolution(C, T, U, V,
 443           aPnt, aProjPS.Point(jmin).Value(), MinSqDist);
 444       }
 445     }
 446     //Cut sharp solutions to keep only minimum and maximum
 447     Standard_Integer imin = SolNumber + 1, imax = mySqDist.Length();
 448     for (i = SolNumber + 1; i <= mySqDist.Length(); i++)
 449     {
 450       if (mySqDist(i) < mySqDist(imin))
 451         imin = i;
 452       if (mySqDist(i) > mySqDist(imax))
 453         imax = i;
 454     }
 455     if (mySqDist.Length() > SolNumber + 2)
 456     {
 457       Standard_Real MinSqDist = mySqDist(imin);
 458       Standard_Real MaxSqDist = mySqDist(imax);
 459       Extrema_POnCurv MinPC = myPOnC(imin);
 460       Extrema_POnCurv MaxPC = myPOnC(imax);
 461       Extrema_POnSurf MinPS = myPOnS(imin);
 462       Extrema_POnSurf MaxPS = myPOnS(imax);
 463
 464       mySqDist.Remove(SolNumber + 1, mySqDist.Length());
 465       myPOnC.Remove(SolNumber + 1, myPOnC.Length());
 466       myPOnS.Remove(SolNumber + 1, myPOnS.Length());
 467
 468       mySqDist.Append(MinSqDist);
 469       myPOnC.Append(MinPC);
 470       myPOnS.Append(MinPS);
 471       mySqDist.Append(MaxSqDist);
 472       myPOnC.Append(MaxPC);
 473       myPOnS.Append(MaxPS);
 474     }
 475   }
 476 }
 477
 478
 479 Standard_Boolean Extrema_ExtCS::IsDone() const
 480 {
 481   return myDone;
 482 }
 483
 484 Standard_Boolean Extrema_ExtCS::IsParallel() const
 485 {
 486   return myIsPar;
 487 }
 488
 489
 490 Standard_Real Extrema_ExtCS::SquareDistance(const Standard_Integer N) const
 491 {
 492   if(!myDone) throw StdFail_NotDone();
 493   if (myIsPar && N != 1) throw StdFail_InfiniteSolutions();
 494   if ((N < 1) || (N > mySqDist.Length())) throw Standard_OutOfRange();
 495   return mySqDist.Value(N);
 496 }
 497
 498
 499 Standard_Integer Extrema_ExtCS::NbExt() const
 500 {
 501   if(!myDone) throw StdFail_NotDone();
 502   return myPOnC.Length();
 503 }
 504
 505
 506
 507 void Extrema_ExtCS::Points(const Standard_Integer N,
 508   Extrema_POnCurv&       P1,
 509   Extrema_POnSurf&       P2) const
 510 {
 511   if(!myDone) throw StdFail_NotDone();
 512   P1 = myPOnC.Value(N);
 513   P2 = myPOnS.Value(N);
 514 }
 515
 516 Standard_Boolean Extrema_ExtCS::AddSolution(const Adaptor3d_Curve& theCurve,
 517   const Standard_Real aT,
 518   const Standard_Real aU,
 519   const Standard_Real aV,
 520   const gp_Pnt& PointOnCurve,
 521   const gp_Pnt& PointOnSurf,
 522   const Standard_Real SquareDist)
 523 {
 524   Standard_Boolean Added = Standard_False;
 525
 526   Standard_Real T = aT, U = aU, V = aV;
 527
 528   if (theCurve.IsPeriodic())
 529     T = ElCLib::InPeriod(T, myucinf, myucinf + theCurve.Period());
 530   if (myS->IsUPeriodic())
 531     U = ElCLib::InPeriod(U, myuinf, myuinf + myS->UPeriod());
 532   if (myS->IsVPeriodic())
 533     V = ElCLib::InPeriod(V, myvinf, myvinf + myS->VPeriod());
 534
 535   Extrema_POnCurv aPC;
 536   Extrema_POnSurf aPS;
 537   if ((myucinf-T) <= mytolC && (T-myucsup) <= mytolC &&
 538     (myuinf-U) <= mytolS && (U-myusup) <= mytolS &&
 539     (myvinf-V) <= mytolS && (V-myvsup) <= mytolS)
 540   {
 541     Standard_Boolean IsNewSolution = Standard_True;
 542     for (Standard_Integer j = 1; j <= mySqDist.Length(); j++)
 543     {
 544       aPC = myPOnC(j);
 545       aPS = myPOnS(j);
 546       Standard_Real Tj = aPC.Parameter();
 547       Standard_Real Uj, Vj;
 548       aPS.Parameter(Uj, Vj);
 549       if (Abs(T - Tj) <= mytolC &&
 550         Abs(U - Uj) <= mytolS &&
 551         Abs(V - Vj) <= mytolS)
 552       {
 553         IsNewSolution = Standard_False;
 554         break;
 555       }
 556     }
 557     if (IsNewSolution)
 558     {
 559       mySqDist.Append(SquareDist);
 560       aPC.SetValues(T, PointOnCurve);
 561       myPOnC.Append(aPC);
 562       myPOnS.Append(Extrema_POnSurf(U, V, PointOnSurf));
 563       Added = Standard_True;
 564     }
 565   }
 566   return Added;
 567 }