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