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-2012 OPEN CASCADE SAS
   5 //
   6 // The content of this file is subject to the Open CASCADE Technology Public
   7 // License Version 6.5 (the "License"). You may not use the content of this file
   8 // except in compliance with the License. Please obtain a copy of the License
   9 // at http://www.opencascade.org and read it completely before using this file.
  10 //
  11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  13 //
  14 // The Original Code and all software distributed under the License is
  15 // distributed on an "AS IS" basis, without warranty of any kind, and the
  16 // Initial Developer hereby disclaims all such warranties, including without
  17 // limitation, any warranties of merchantability, fitness for a particular
  18 // purpose or non-infringement. Please see the License for the specific terms
  19 // and conditions governing the rights and limitations under the License.
  20
  21
  22
  23 #include <Extrema_ExtPRevS.ixx>
  24 #include <Adaptor3d_HCurve.hxx>
  25 #include <Extrema_ExtPElC.hxx>
  26 #include <Extrema_GenExtPS.hxx>
  27 #include <Extrema_POnCurv.hxx>
  28 #include <Extrema_POnSurf.hxx>
  29 #include <gp_Ax1.hxx>
  30 #include <gp_Ax2.hxx>
  31 #include <gp_Lin.hxx>
  32 #include <gp_Pln.hxx>
  33 #include <gp_Pnt.hxx>
  34 #include <gp_Trsf.hxx>
  35 #include <gp_Vec.hxx>
  36 #include <Precision.hxx>
  37 #include <ElCLib.hxx>
  38
  39 static gp_Ax2 GetPosition (const Adaptor3d_SurfaceOfRevolution& S)//const Handle(Adaptor_HCurve)& C)
  40 {
  41   Handle(Adaptor3d_HCurve) C = S.BasisCurve();
  42
  43   switch (C->GetType()) {
  44
  45   case GeomAbs_Line: {
  46     gp_Lin L = C->Line();
  47     gp_Dir N = S.AxeOfRevolution().Direction();
  48     if (N.IsParallel(L.Direction(), Precision::Angular())) {
  49       gp_Vec OO (L.Location(), S.AxeOfRevolution().Location());
  50       if (OO.Magnitude() <= gp::Resolution()) {
  51         OO = gp_Vec(L.Location(), ElCLib::Value(100,L));
  52         if (N.IsParallel(OO, Precision::Angular()))
  53           return gp_Ax2(); // Line and axe of revolution coinside
  54       }
  55       N ^= OO;
  56     }
  57     else {
  58       N ^= L.Direction();
  59     }
  60     return gp_Ax2 (L.Location(), N, L.Direction());
  61   }
  62   case GeomAbs_Circle:
  63     return C->Circle().Position();
  64   case GeomAbs_Ellipse:
  65     return C->Ellipse().Position();
  66   case GeomAbs_Hyperbola:
  67     return C->Hyperbola().Position();
  68   case GeomAbs_Parabola:
  69     return C->Parabola().Position();
  70   default:
  71     return gp_Ax2();
  72   }
  73 }
  74
  75 //=======================================================================
  76 //function : HasSingularity
  77 //purpose  :
  78 //=======================================================================
  79
  80 static Standard_Boolean HasSingularity(const Adaptor3d_SurfaceOfRevolution& S)
  81 {
  82
  83   const Handle(Adaptor3d_HCurve) C = S.BasisCurve();
  84   gp_Dir N = S.AxeOfRevolution().Direction();
  85   gp_Pnt P = S.AxeOfRevolution().Location();
  86
  87   gp_Lin L(P, N);
  88
  89   P = C->Value(C->FirstParameter());
  90
  91   if(L.SquareDistance(P) < Precision::Confusion() * Precision::Confusion()) return Standard_True;
  92
  93   P = C->Value(C->LastParameter());
  94
  95   if(L.SquareDistance(P) < Precision::Confusion() * Precision::Confusion()) return Standard_True;
  96
  97   return Standard_False;
  98 }
  99
 100 //=============================================================================
 101
 102 static void PerformExtPElC (Extrema_ExtPElC& E,
 103                             const gp_Pnt& P,
 104                             const Handle(Adaptor3d_HCurve)& C,
 105                             const Standard_Real Tol)
 106 {
 107   switch (C->GetType()) {
 108   case GeomAbs_Hyperbola:
 109     E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
 110     return;
 111   case GeomAbs_Line:
 112     E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
 113     return;
 114   case GeomAbs_Circle:
 115     E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
 116     return;
 117   case GeomAbs_Ellipse:
 118     E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
 119     return;
 120   case GeomAbs_Parabola:
 121     E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
 122     return;
 123 #ifndef DEB
 124   default:
 125     return ;
 126 #endif
 127   }
 128 }
 129
 130 //=======================================================================
 131 //function : IsCaseAnalyticallyComputable
 132 //purpose  :
 133 //=======================================================================
 134
 135 static Standard_Boolean IsCaseAnalyticallyComputable
 136   (const GeomAbs_CurveType& theType,
 137    const gp_Ax2& theCurvePos,
 138    const gp_Ax1& AxeOfRevolution)
 139 {
 140   // check type
 141   switch (theType) {
 142   case GeomAbs_Line:
 143   case GeomAbs_Circle:
 144   case GeomAbs_Ellipse:
 145   case GeomAbs_Hyperbola:
 146   case GeomAbs_Parabola:
 147     break;
 148   default:
 149     return  Standard_False;
 150   }
 151 //  the axe of revolution must be in the plane of the curve.
 152   gp_Pln pl(theCurvePos.Location(), theCurvePos.Direction());
 153   gp_Pnt p1 = AxeOfRevolution.Location();
 154   Standard_Real dist = 100., dist2 = dist * dist;
 155   Standard_Real aThreshold = Precision::Angular() * Precision::Angular() * dist2;
 156   gp_Pnt p2 = AxeOfRevolution.Location().XYZ() + dist * AxeOfRevolution.Direction().XYZ();
 157
 158   if((pl.SquareDistance(p1) < aThreshold) &&
 159      (pl.SquareDistance(p2) < aThreshold))
 160     return Standard_True;
 161   return Standard_False;
 162   //   gp_Vec V (AxeOfRevolution.Location(),theCurvePos.Location());
 163   //   if (Abs( V * theCurvePos.Direction()) <= gp::Resolution())
 164   //     return Standard_True;
 165   //   else
 166   //     return Standard_False;
 167 }
 168
 169 //=======================================================================
 170 //function : IsOriginalPnt
 171 //purpose  :
 172 //=======================================================================
 173
 174 static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
 175                                        const Extrema_POnSurf* Points,
 176                                        const Standard_Integer NbPoints)
 177 {
 178   for (Standard_Integer i=1; i<=NbPoints; i++) {
 179     if (Points[i].Value().IsEqual(P, Precision::Confusion())) {
 180       return Standard_False;
 181     }
 182   }
 183   return Standard_True;
 184 }
 185 //=======================================================================
 186 //function : IsExtremum
 187 //purpose  :
 188 //=======================================================================
 189
 190 static Standard_Boolean IsExtremum (const Standard_Real U, const Standard_Real V,
 191                                     const gp_Pnt& P, const Adaptor3d_SurfacePtr& S,
 192                                     gp_Pnt& E,        Standard_Real& Dist2,
 193                                     const Standard_Boolean IsVSup,
 194                                     const Standard_Boolean IsMin)
 195 {
 196   E = S->Value(U,V);
 197   Dist2 = P.SquareDistance(E);
 198   if (IsMin)
 199     return (Dist2 < P.SquareDistance(S->Value(U+1,V)) &&
 200             Dist2 < P.SquareDistance(S->Value(U-1,V)) &&
 201             Dist2 < P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
 202   else
 203     return (Dist2 > P.SquareDistance(S->Value(U+1,V)) &&
 204             Dist2 > P.SquareDistance(S->Value(U-1,V)) &&
 205             Dist2 > P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
 206 }
 207 //=======================================================================
 208 //function : Extrema_ExtPRevS
 209 //purpose  :
 210 //=======================================================================
 211
 212 Extrema_ExtPRevS::Extrema_ExtPRevS()
 213 {
 214   myS=NULL;
 215   myDone = Standard_False;
 216 }
 217 //=======================================================================
 218 //function : Extrema_ExtPRevS
 219 //purpose  :
 220 //=======================================================================
 221
 222 Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
 223                                    const Adaptor3d_SurfaceOfRevolution& S,
 224                                    const Standard_Real Umin,
 225                                    const Standard_Real Usup,
 226                                    const Standard_Real Vmin,
 227                                    const Standard_Real Vsup,
 228                                    const Standard_Real TolU,
 229                                    const Standard_Real TolV)
 230 {
 231   myS=NULL;
 232   Initialize (S,
 233               Umin, Usup, Vmin, Vsup,
 234               TolU, TolV);
 235   Perform(P);
 236 }
 237 //=======================================================================
 238 //function : Extrema_ExtPRevS
 239 //purpose  :
 240 //=======================================================================
 241
 242 Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
 243                                    const Adaptor3d_SurfaceOfRevolution& S,
 244                                    const Standard_Real TolU,
 245                                    const Standard_Real TolV)
 246 {
 247   myS=NULL;
 248   Initialize (S,
 249               S.FirstUParameter(), S.LastUParameter(),
 250               S.FirstVParameter(), S.LastVParameter(),
 251               TolU, TolV);
 252   Perform(P);
 253 }
 254 //=======================================================================
 255 //function : Initialize
 256 //purpose  :
 257 //=======================================================================
 258
 259 void Extrema_ExtPRevS::Initialize(const Adaptor3d_SurfaceOfRevolution& S,
 260                                   const Standard_Real Umin,
 261                                   const Standard_Real Usup,
 262                                   const Standard_Real Vmin,
 263                                   const Standard_Real Vsup,
 264                                   const Standard_Real TolU,
 265                                   const Standard_Real TolV)
 266 {
 267   myvinf=Vmin;
 268   myvsup=Vsup;
 269   mytolv=TolV;
 270
 271   Handle(Adaptor3d_HCurve) anACurve = S.BasisCurve();
 272   if (!myS || myS != (Adaptor3d_SurfacePtr)&S) {
 273     myS = Adaptor3d_SurfacePtr(&S);
 274     myPosition = GetPosition(S);
 275     myIsAnalyticallyComputable =
 276       IsCaseAnalyticallyComputable (anACurve->GetType(),myPosition,S.AxeOfRevolution());
 277   }
 278   if (!myIsAnalyticallyComputable) {
 279
 280     Standard_Integer nbu = 32, nbv = 32;
 281
 282     if(HasSingularity(S)) {nbv = 100;}
 283
 284     myExtPS.Initialize(S, nbu, nbv,
 285                        Umin, Usup, Vmin, Vsup,
 286                        TolU, TolV);
 287   }
 288 }
 289 //=======================================================================
 290 //function : Perform
 291 //purpose  :
 292 //=======================================================================
 293
 294 void Extrema_ExtPRevS::Perform(const gp_Pnt& P)
 295 {
 296   myDone = Standard_False;
 297   myNbExt = 0;
 298
 299   if (!myIsAnalyticallyComputable) {
 300
 301     myExtPS.Perform(P);
 302     myDone = myExtPS.IsDone();
 303     myNbExt = myExtPS.NbExt();
 304     return;
 305   }
 306
 307   Handle(Adaptor3d_HCurve) anACurve = myS->BasisCurve();
 308
 309   gp_Ax1 Ax = myS->AxeOfRevolution();
 310   gp_Vec Dir = Ax.Direction(), Z = myPosition.Direction();
 311   gp_Pnt O = Ax.Location();
 312
 313   Standard_Real OPdir = gp_Vec(O,P).Dot(Dir);
 314   gp_Pnt Pp = P.Translated(Dir.Multiplied(-OPdir));
 315   if (O.IsEqual(Pp,Precision::Confusion())) // P is on the AxeOfRevolution
 316     return;
 317
 318   Standard_Real U,V;
 319   gp_Pnt P1, Ppp;
 320   Standard_Real OPpz = gp_Vec(O,Pp).Dot(Z);
 321   if (Abs(OPpz) <= gp::Resolution()) {
 322     Ppp = Pp;
 323     U = 0;
 324   }
 325   else {
 326     Ppp = Pp.Translated(Z.Multiplied(-OPpz));
 327     if (O.IsEqual(Ppp,Precision::Confusion()))
 328       U = M_PI/2;
 329     else {
 330       U = gp_Vec(O,Ppp).AngleWithRef(gp_Vec(O,Pp),Dir);
 331     }
 332   }
 333
 334   gp_Vec OPpp (O,Ppp), OPq (O, myS->Value(M_PI/2,0));
 335   if (U != M_PI/2) {
 336     if (Abs(OPq.Magnitude()) <= gp::Resolution())
 337       OPq = gp_Vec(O, myS->Value(M_PI/2,anACurve->LastParameter()/10));
 338     if (OPpp.AngleWithRef(OPq,Dir) < 0)
 339       U += M_PI;
 340   }
 341
 342   gp_Trsf T;
 343   T.SetRotation(Ax, -U);
 344   P1 = P.Transformed(T);
 345
 346   gp_Pnt E;
 347   Standard_Real Dist2;
 348   Standard_Integer i;
 349
 350   Extrema_ExtPElC anExt;
 351   PerformExtPElC(anExt, P1, anACurve, mytolv);
 352
 353   if (anExt.IsDone()) {
 354     myDone = Standard_True;
 355     for (i=1; i<=anExt.NbExt(); i++) {
 356       Extrema_POnCurv POC=anExt.Point(i);
 357       V = POC.Parameter();
 358       if (V > myvsup) {
 359         //       if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
 360         //                         Standard_True, anExt.IsMin(i))) continue;
 361         Standard_Real newV = myvsup;
 362
 363         if((anACurve->GetType() == GeomAbs_Circle) ||
 364            (anACurve->GetType() == GeomAbs_Ellipse)) {
 365           newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
 366
 367           if (newV > myvsup) {
 368             newV = myvsup;
 369           }
 370         }
 371         V = newV;
 372
 373         if ( !IsExtremum (U, V, P, myS, E, Dist2,
 374                            Standard_True, anExt.IsMin(i))) {
 375           continue;
 376         }
 377       } else if (V < myvinf) {
 378         //      if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
 379         //                        Standard_False, anExt.IsMin(i))) continue;
 380
 381         Standard_Real newV = myvinf;
 382
 383         if((anACurve->GetType() == GeomAbs_Circle) ||
 384            (anACurve->GetType() == GeomAbs_Ellipse)) {
 385           newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
 386
 387           if(newV < myvinf)
 388             newV = myvinf;
 389         }
 390         V = newV;
 391
 392         if ( !IsExtremum (U, V, P, myS, E, Dist2,
 393                           Standard_False, anExt.IsMin(i))) continue;
 394       } else {
 395         E = myS->Value(U,V);
 396         Dist2 = P.SquareDistance(E);
 397       }
 398       if (IsOriginalPnt(E, myPoint, myNbExt)) {
 399         myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
 400         mySqDist[myNbExt] = Dist2;
 401       }
 402     }
 403   }
 404   T.SetRotation(Ax, M_PI);
 405   P1.Transform(T);
 406
 407   PerformExtPElC(anExt, P1, anACurve, mytolv);
 408   if (anExt.IsDone()) {
 409     myDone = Standard_True;
 410
 411     U += M_PI;
 412
 413     for (i=1; i<=anExt.NbExt(); i++) {
 414       Extrema_POnCurv POC=anExt.Point(i);
 415       V = POC.Parameter();
 416       if (V > myvsup) {
 417         //      if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
 418         //                         Standard_True, anExt.IsMin(i))) continue;
 419
 420         Standard_Real newV = myvsup;
 421
 422         if((anACurve->GetType() == GeomAbs_Circle) ||
 423            (anACurve->GetType() == GeomAbs_Ellipse)) {
 424           newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
 425
 426           if (newV > myvsup) {
 427             newV = myvsup;
 428           }
 429         }
 430         V = newV;
 431
 432         if ( !IsExtremum (U, V, P, myS, E, Dist2,
 433                            Standard_True, anExt.IsMin(i))) continue;
 434       } else if (V < myvinf) {
 435         //      if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
 436         //                        Standard_False, anExt.IsMin(i))) continue;
 437         Standard_Real newV = myvinf;
 438
 439         if((anACurve->GetType() == GeomAbs_Circle) ||
 440            (anACurve->GetType() == GeomAbs_Ellipse)) {
 441           newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
 442
 443           if(newV < myvinf)
 444             newV = myvinf;
 445         }
 446         V = newV;
 447
 448         if ( !IsExtremum (U, V, P, myS, E, Dist2,
 449                           Standard_False, anExt.IsMin(i))) continue;
 450       } else {
 451         E = myS->Value(U,V);
 452         Dist2 = P.SquareDistance(E);
 453       }
 454       if (IsOriginalPnt(E, myPoint, myNbExt)) {
 455         myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
 456         mySqDist[myNbExt] = Dist2;
 457
 458       }
 459     }
 460   }
 461 }
 462
 463
 464 //=======================================================================
 465 //function : IsDone
 466 //purpose  :
 467 //=======================================================================
 468
 469 Standard_Boolean Extrema_ExtPRevS::IsDone() const
 470 {
 471   return myDone;
 472 }
 473
 474
 475 //=======================================================================
 476 //function : NbExt
 477 //purpose  :
 478 //=======================================================================
 479
 480 Standard_Integer Extrema_ExtPRevS::NbExt() const
 481 {
 482   if (!IsDone()) { StdFail_NotDone::Raise(); }
 483   return myNbExt;
 484 }
 485
 486
 487 //=======================================================================
 488 //function : Value
 489 //purpose  :
 490 //=======================================================================
 491
 492 Standard_Real Extrema_ExtPRevS::SquareDistance(const Standard_Integer N) const
 493 {
 494   if (!IsDone()) { StdFail_NotDone::Raise(); }
 495   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
 496   if (myIsAnalyticallyComputable)
 497     return mySqDist[N];
 498   else
 499     return myExtPS.SquareDistance(N);
 500 }
 501 //=======================================================================
 502 //function : Point
 503 //purpose  :
 504 //=======================================================================
 505
 506 Extrema_POnSurf Extrema_ExtPRevS::Point(const Standard_Integer N) const
 507 {
 508   if (!IsDone()) { StdFail_NotDone::Raise(); }
 509   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
 510   if (myIsAnalyticallyComputable)
 511     return myPoint[N];
 512   else
 513     return myExtPS.Point(N);
 514 }
 515
 516