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