src/Extrema/Extrema_ExtPRevS.cxx

   1 // Created on: 1999-09-21
   2 // Created by: Edward AGAPOV
   3 // Copyright (c) 1999-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 <GeomAdaptor_HSurfaceOfRevolution.hxx>
  20 #include <ElCLib.hxx>
  21 #include <Extrema_ExtPElC.hxx>
  22 #include <Extrema_ExtPRevS.hxx>
  23 #include <Extrema_GenExtPS.hxx>
  24 #include <Extrema_POnCurv.hxx>
  25 #include <Extrema_POnSurf.hxx>
  26 #include <gp_Ax1.hxx>
  27 #include <gp_Ax2.hxx>
  28 #include <gp_Lin.hxx>
  29 #include <gp_Pln.hxx>
  30 #include <gp_Pnt.hxx>
  31 #include <gp_Trsf.hxx>
  32 #include <gp_Vec.hxx>
  33 #include <Precision.hxx>
  34 #include <Standard_OutOfRange.hxx>
  35 #include <Standard_Type.hxx>
  36 #include <StdFail_NotDone.hxx>
  37
  38 IMPLEMENT_STANDARD_RTTIEXT(Extrema_ExtPRevS,Standard_Transient)
  39
  40 static gp_Ax2 GetPosition (const GeomAdaptor_SurfaceOfRevolution& S)//const Handle(Adaptor_HCurve)& C)
  41 {
  42   Handle(Adaptor3d_HCurve) C = S.BasisCurve();
  43
  44   switch (C->GetType()) {
  45
  46   case GeomAbs_Line: {
  47     gp_Lin L = C->Line();
  48     gp_Dir N = S.AxeOfRevolution().Direction();
  49     if (N.IsParallel(L.Direction(), Precision::Angular())) {
  50       gp_Vec OO (L.Location(), S.AxeOfRevolution().Location());
  51       if (OO.Magnitude() <= gp::Resolution()) {
  52         OO = gp_Vec(L.Location(), ElCLib::Value(100,L));
  53         if (N.IsParallel(OO, Precision::Angular()))
  54           return gp_Ax2(); // Line and axe of revolution coinside
  55       }
  56       N ^= OO;
  57     }
  58     else {
  59       N ^= L.Direction();
  60     }
  61     return gp_Ax2 (L.Location(), N, L.Direction());
  62   }
  63   case GeomAbs_Circle:
  64     return C->Circle().Position();
  65   case GeomAbs_Ellipse:
  66     return C->Ellipse().Position();
  67   case GeomAbs_Hyperbola:
  68     return C->Hyperbola().Position();
  69   case GeomAbs_Parabola:
  70     return C->Parabola().Position();
  71   default:
  72     return gp_Ax2();
  73   }
  74 }
  75
  76 //=======================================================================
  77 //function : HasSingularity
  78 //purpose  :
  79 //=======================================================================
  80
  81 static Standard_Boolean HasSingularity(const GeomAdaptor_SurfaceOfRevolution& S)
  82 {
  83
  84   const Handle(Adaptor3d_HCurve) C = S.BasisCurve();
  85   gp_Dir N = S.AxeOfRevolution().Direction();
  86   gp_Pnt P = S.AxeOfRevolution().Location();
  87
  88   gp_Lin L(P, N);
  89
  90   P = C->Value(C->FirstParameter());
  91
  92   if(L.SquareDistance(P) < Precision::SquareConfusion()) return Standard_True;
  93
  94   P = C->Value(C->LastParameter());
  95
  96   if(L.SquareDistance(P) < Precision::SquareConfusion()) return Standard_True;
  97
  98   return Standard_False;
  99 }
 100
 101 //=============================================================================
 102
 103 static void PerformExtPElC (Extrema_ExtPElC& E,
 104                             const gp_Pnt& P,
 105                             const Handle(Adaptor3d_HCurve)& C,
 106                             const Standard_Real Tol)
 107 {
 108   switch (C->GetType()) {
 109   case GeomAbs_Hyperbola:
 110     E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
 111     return;
 112   case GeomAbs_Line:
 113     E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
 114     return;
 115   case GeomAbs_Circle:
 116     E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
 117     return;
 118   case GeomAbs_Ellipse:
 119     E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
 120     return;
 121   case GeomAbs_Parabola:
 122     E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
 123     return;
 124   default:
 125     return ;
 126   }
 127 }
 128
 129 //=======================================================================
 130 //function : IsCaseAnalyticallyComputable
 131 //purpose  :
 132 //=======================================================================
 133
 134 static Standard_Boolean IsCaseAnalyticallyComputable
 135   (const GeomAbs_CurveType& theType,
 136    const gp_Ax2& theCurvePos,
 137    const gp_Ax1& AxeOfRevolution)
 138 {
 139   // check type
 140   switch (theType) {
 141   case GeomAbs_Line:
 142   case GeomAbs_Circle:
 143   case GeomAbs_Ellipse:
 144   case GeomAbs_Hyperbola:
 145   case GeomAbs_Parabola:
 146     break;
 147   default:
 148     return  Standard_False;
 149   }
 150 //  the axe of revolution must be in the plane of the curve.
 151   gp_Pln pl(theCurvePos.Location(), theCurvePos.Direction());
 152   gp_Pnt p1 = AxeOfRevolution.Location();
 153   Standard_Real dist = 100., dist2 = dist * dist;
 154   Standard_Real aThreshold = Precision::Angular() * Precision::Angular() * dist2;
 155   gp_Pnt p2 = AxeOfRevolution.Location().XYZ() + dist * AxeOfRevolution.Direction().XYZ();
 156
 157   if((pl.SquareDistance(p1) < aThreshold) &&
 158      (pl.SquareDistance(p2) < aThreshold))
 159     return Standard_True;
 160   return Standard_False;
 161   //   gp_Vec V (AxeOfRevolution.Location(),theCurvePos.Location());
 162   //   if (Abs( V * theCurvePos.Direction()) <= gp::Resolution())
 163   //     return Standard_True;
 164   //   else
 165   //     return Standard_False;
 166 }
 167
 168 //=======================================================================
 169 //function : IsOriginalPnt
 170 //purpose  :
 171 //=======================================================================
 172
 173 static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
 174                                        const Extrema_POnSurf* Points,
 175                                        const Standard_Integer NbPoints)
 176 {
 177   for (Standard_Integer i=1; i<=NbPoints; i++) {
 178     if (Points[i-1].Value().IsEqual(P, Precision::Confusion())) {
 179       return Standard_False;
 180     }
 181   }
 182   return Standard_True;
 183 }
 184 //=======================================================================
 185 //function : IsExtremum
 186 //purpose  :
 187 //=======================================================================
 188
 189 static Standard_Boolean IsExtremum (const Standard_Real U, const Standard_Real V,
 190                                     const gp_Pnt& P, const Adaptor3d_SurfacePtr& S,
 191                                     gp_Pnt& E,        Standard_Real& Dist2,
 192                                     const Standard_Boolean IsVSup,
 193                                     const Standard_Boolean IsMin)
 194 {
 195   E = S->Value(U,V);
 196   Dist2 = P.SquareDistance(E);
 197   if (IsMin)
 198     return (Dist2 < P.SquareDistance(S->Value(U+1,V)) &&
 199             Dist2 < P.SquareDistance(S->Value(U-1,V)) &&
 200             Dist2 < P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
 201   else
 202     return (Dist2 > P.SquareDistance(S->Value(U+1,V)) &&
 203             Dist2 > P.SquareDistance(S->Value(U-1,V)) &&
 204             Dist2 > P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
 205 }
 206 //=======================================================================
 207 //function : Extrema_ExtPRevS
 208 //purpose  :
 209 //=======================================================================
 210
 211 Extrema_ExtPRevS::Extrema_ExtPRevS()
 212 {
 213   myDone = Standard_False;
 214 }
 215 //=======================================================================
 216 //function : Extrema_ExtPRevS
 217 //purpose  :
 218 //=======================================================================
 219
 220 Extrema_ExtPRevS::Extrema_ExtPRevS (const gp_Pnt&                                 theP,
 221                                     const Handle(GeomAdaptor_HSurfaceOfRevolution)& theS,
 222                                     const Standard_Real                           theUmin,
 223                                     const Standard_Real                           theUsup,
 224                                     const Standard_Real                           theVmin,
 225                                     const Standard_Real                           theVsup,
 226                                     const Standard_Real                           theTolU,
 227                                     const Standard_Real                           theTolV)
 228 {
 229   Initialize (theS,
 230               theUmin,
 231               theUsup,
 232               theVmin,
 233               theVsup,
 234               theTolU,
 235               theTolV);
 236
 237   Perform (theP);
 238 }
 239 //=======================================================================
 240 //function : Extrema_ExtPRevS
 241 //purpose  :
 242 //=======================================================================
 243
 244 Extrema_ExtPRevS::Extrema_ExtPRevS (const gp_Pnt&                                 theP,
 245                                     const Handle(GeomAdaptor_HSurfaceOfRevolution)& theS,
 246                                     const Standard_Real                           theTolU,
 247                                     const Standard_Real                           theTolV)
 248 {
 249   Initialize (theS,
 250               theS->FirstUParameter(),
 251               theS->LastUParameter(),
 252               theS->FirstVParameter(),
 253               theS->LastVParameter(),
 254               theTolU,
 255               theTolV);
 256
 257   Perform (theP);
 258 }
 259 //=======================================================================
 260 //function : Initialize
 261 //purpose  :
 262 //=======================================================================
 263
 264 void Extrema_ExtPRevS::Initialize (const Handle(GeomAdaptor_HSurfaceOfRevolution)& theS,
 265                                    const Standard_Real                           theUmin,
 266                                    const Standard_Real                           theUsup,
 267                                    const Standard_Real                           theVmin,
 268                                    const Standard_Real                           theVsup,
 269                                    const Standard_Real                           theTolU,
 270                                    const Standard_Real                           theTolV)
 271 {
 272   myvinf = theVmin;
 273   myvsup = theVsup;
 274   mytolv = theTolV;
 275
 276   Handle(Adaptor3d_HCurve) anACurve = theS->BasisCurve();
 277
 278   if (myS != theS)
 279   {
 280     myS = theS;
 281     myPosition = GetPosition (theS->ChangeSurface());
 282     myIsAnalyticallyComputable =
 283       IsCaseAnalyticallyComputable (anACurve->GetType(), myPosition, theS->AxeOfRevolution());
 284   }
 285
 286   if (!myIsAnalyticallyComputable)
 287   {
 288     Standard_Integer aNbu = 32, aNbv = 32;
 289
 290     if (HasSingularity (theS->ChangeSurface()))
 291     {
 292       aNbv = 100;
 293     }
 294
 295     myExtPS.Initialize (theS->ChangeSurface(),
 296                         aNbu,
 297                         aNbv,
 298                         theUmin,
 299                         theUsup,
 300                         theVmin,
 301                         theVsup,
 302                         theTolU,
 303                         theTolV);
 304   }
 305 }
 306 //=======================================================================
 307 //function : Perform
 308 //purpose  :
 309 //=======================================================================
 310
 311 void Extrema_ExtPRevS::Perform(const gp_Pnt& P)
 312 {
 313   myDone = Standard_False;
 314   myNbExt = 0;
 315
 316   if (!myIsAnalyticallyComputable) {
 317
 318     myExtPS.Perform(P);
 319     myDone = myExtPS.IsDone();
 320     myNbExt = myExtPS.NbExt();
 321     return;
 322   }
 323
 324   Handle(Adaptor3d_HCurve) anACurve = myS->BasisCurve();
 325
 326   gp_Ax1 Ax = myS->AxeOfRevolution();
 327   gp_Vec Dir = Ax.Direction(), Z = myPosition.Direction();
 328   gp_Pnt O = Ax.Location();
 329
 330   Standard_Real OPdir = gp_Vec(O,P).Dot(Dir);
 331   gp_Pnt Pp = P.Translated(Dir.Multiplied(-OPdir));
 332   if (O.IsEqual(Pp,Precision::Confusion())) // P is on the AxeOfRevolution
 333     return;
 334
 335   Standard_Real U,V;
 336   gp_Pnt P1, Ppp;
 337   Standard_Real OPpz = gp_Vec(O,Pp).Dot(Z);
 338   if (Abs(OPpz) <= gp::Resolution()) {
 339     Ppp = Pp;
 340     U = 0;
 341   }
 342   else {
 343     Ppp = Pp.Translated(Z.Multiplied(-OPpz));
 344     if (O.IsEqual(Ppp,Precision::Confusion()))
 345       U = M_PI/2;
 346     else {
 347       U = gp_Vec(O,Ppp).AngleWithRef(gp_Vec(O,Pp),Dir);
 348     }
 349   }
 350
 351   gp_Vec OPpp (O,Ppp), OPq (O, myS->Value(M_PI/2,0));
 352   if (U != M_PI/2) {
 353     if (Abs(OPq.Magnitude()) <= gp::Resolution())
 354       OPq = gp_Vec(O, myS->Value(M_PI/2,anACurve->LastParameter()/10));
 355     if (OPpp.AngleWithRef(OPq,Dir) < 0)
 356       U += M_PI;
 357   }
 358
 359   gp_Trsf T;
 360   T.SetRotation(Ax, -U);
 361   P1 = P.Transformed(T);
 362
 363   gp_Pnt E;
 364   Standard_Real Dist2;
 365   Standard_Integer i;
 366
 367   Extrema_ExtPElC anExt;
 368   PerformExtPElC(anExt, P1, anACurve, mytolv);
 369
 370   if (anExt.IsDone()) {
 371     myDone = Standard_True;
 372     for (i=1; i<=anExt.NbExt(); i++) {
 373       Extrema_POnCurv POC=anExt.Point(i);
 374       V = POC.Parameter();
 375       if (V > myvsup) {
 376         //       if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
 377         //                         Standard_True, anExt.IsMin(i))) continue;
 378         Standard_Real newV = myvsup;
 379
 380         if((anACurve->GetType() == GeomAbs_Circle) ||
 381            (anACurve->GetType() == GeomAbs_Ellipse)) {
 382           newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
 383
 384           if (newV > myvsup) {
 385             newV -= 2. * M_PI;
 386
 387             if (newV + mytolv < myvinf) {
 388               newV = myvsup;
 389             } else if (newV < myvinf) {
 390               newV = myvinf;
 391             }
 392           }
 393         }
 394         V = newV;
 395
 396         if ( !IsExtremum (U, V, P, &(myS->ChangeSurface()), E, Dist2,
 397                            Standard_True, anExt.IsMin(i))) {
 398           continue;
 399         }
 400       } else if (V < myvinf) {
 401         //      if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
 402         //                        Standard_False, anExt.IsMin(i))) continue;
 403
 404         Standard_Real newV = myvinf;
 405
 406         if((anACurve->GetType() == GeomAbs_Circle) ||
 407            (anACurve->GetType() == GeomAbs_Ellipse)) {
 408           newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
 409
 410           if(newV < myvinf) {
 411             newV += 2. * M_PI;
 412
 413             if (newV - mytolv > myvsup) {
 414               newV = myvinf;
 415             } else if (newV > myvsup) {
 416               newV = myvsup;
 417             }
 418           }
 419         }
 420         V = newV;
 421
 422         if ( !IsExtremum (U, V, P, &(myS->ChangeSurface()), E, Dist2,
 423                           Standard_False, anExt.IsMin(i))) continue;
 424       } else {
 425         E = myS->Value(U,V);
 426         Dist2 = P.SquareDistance(E);
 427       }
 428       if (IsOriginalPnt(E, myPoint, myNbExt)) {
 429         myPoint[myNbExt] = Extrema_POnSurf(U,V,E);
 430         mySqDist[myNbExt] = Dist2;
 431         myNbExt++;
 432       }
 433     }
 434   }
 435   T.SetRotation(Ax, M_PI);
 436   P1.Transform(T);
 437
 438   PerformExtPElC(anExt, P1, anACurve, mytolv);
 439   if (anExt.IsDone()) {
 440     myDone = Standard_True;
 441
 442     U += M_PI;
 443
 444     for (i=1; i<=anExt.NbExt(); i++) {
 445       Extrema_POnCurv POC=anExt.Point(i);
 446       V = POC.Parameter();
 447       if (V > myvsup) {
 448         //      if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
 449         //                         Standard_True, anExt.IsMin(i))) continue;
 450
 451         Standard_Real newV = myvsup;
 452
 453         if((anACurve->GetType() == GeomAbs_Circle) ||
 454            (anACurve->GetType() == GeomAbs_Ellipse)) {
 455           newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
 456
 457           if (newV > myvsup) {
 458             newV -= 2. * M_PI;
 459
 460             if (newV + mytolv < myvinf) {
 461               newV = myvsup;
 462             } else if (newV < myvinf) {
 463               newV = myvinf;
 464             }
 465           }
 466         }
 467         V = newV;
 468
 469         if ( !IsExtremum (U, V, P, &(myS->ChangeSurface()), E, Dist2,
 470                            Standard_True, anExt.IsMin(i))) continue;
 471       } else if (V < myvinf) {
 472         //      if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
 473         //                        Standard_False, anExt.IsMin(i))) continue;
 474         Standard_Real newV = myvinf;
 475
 476         if((anACurve->GetType() == GeomAbs_Circle) ||
 477            (anACurve->GetType() == GeomAbs_Ellipse)) {
 478           newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
 479
 480           if(newV < myvinf) {
 481             newV += 2. * M_PI;
 482
 483             if (newV - mytolv > myvsup) {
 484               newV = myvinf;
 485             } else if (newV > myvsup) {
 486               newV = myvsup;
 487             }
 488           }
 489         }
 490         V = newV;
 491
 492         if ( !IsExtremum (U, V, P, &(myS->ChangeSurface()), E, Dist2,
 493                           Standard_False, anExt.IsMin(i))) continue;
 494       } else {
 495         E = myS->Value(U,V);
 496         Dist2 = P.SquareDistance(E);
 497       }
 498       if (IsOriginalPnt(E, myPoint, myNbExt)) {
 499         myPoint[myNbExt] = Extrema_POnSurf(U,V,E);
 500         mySqDist[myNbExt] = Dist2;
 501         myNbExt++;
 502       }
 503     }
 504   }
 505 }
 506
 507
 508 //=======================================================================
 509 //function : IsDone
 510 //purpose  :
 511 //=======================================================================
 512
 513 Standard_Boolean Extrema_ExtPRevS::IsDone() const
 514 {
 515   return myDone;
 516 }
 517
 518
 519 //=======================================================================
 520 //function : NbExt
 521 //purpose  :
 522 //=======================================================================
 523
 524 Standard_Integer Extrema_ExtPRevS::NbExt() const
 525 {
 526   if (!IsDone()) { throw StdFail_NotDone(); }
 527   return myNbExt;
 528 }
 529
 530
 531 //=======================================================================
 532 //function : Value
 533 //purpose  :
 534 //=======================================================================
 535
 536 Standard_Real Extrema_ExtPRevS::SquareDistance(const Standard_Integer N) const
 537 {
 538   if (!IsDone()) { throw StdFail_NotDone(); }
 539   if ((N < 1) || (N > myNbExt)) { throw Standard_OutOfRange(); }
 540   if (myIsAnalyticallyComputable)
 541     return mySqDist[N-1];
 542   else
 543     return myExtPS.SquareDistance(N);
 544 }
 545 //=======================================================================
 546 //function : Point
 547 //purpose  :
 548 //=======================================================================
 549
 550 const Extrema_POnSurf& Extrema_ExtPRevS::Point(const Standard_Integer N) const
 551 {
 552   if (!IsDone()) { throw StdFail_NotDone(); }
 553   if ((N < 1) || (N > myNbExt)) { throw Standard_OutOfRange(); }
 554   if (myIsAnalyticallyComputable)
 555     return myPoint[N-1];
 556   else
 557     return myExtPS.Point(N);
 558 }
 559
 560