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