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::SquareConfusion()) return Standard_True;
  92
  93   P = C->Value(C->LastParameter());
  94
  95   if(L.SquareDistance(P) < Precision::SquareConfusion()) 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   default:
 124     return ;
 125   }
 126 }
 127
 128 //=======================================================================
 129 //function : IsCaseAnalyticallyComputable
 130 //purpose  :
 131 //=======================================================================
 132
 133 static Standard_Boolean IsCaseAnalyticallyComputable
 134   (const GeomAbs_CurveType& theType,
 135    const gp_Ax2& theCurvePos,
 136    const gp_Ax1& AxeOfRevolution)
 137 {
 138   // check type
 139   switch (theType) {
 140   case GeomAbs_Line:
 141   case GeomAbs_Circle:
 142   case GeomAbs_Ellipse:
 143   case GeomAbs_Hyperbola:
 144   case GeomAbs_Parabola:
 145     break;
 146   default:
 147     return  Standard_False;
 148   }
 149 //  the axe of revolution must be in the plane of the curve.
 150   gp_Pln pl(theCurvePos.Location(), theCurvePos.Direction());
 151   gp_Pnt p1 = AxeOfRevolution.Location();
 152   Standard_Real dist = 100., dist2 = dist * dist;
 153   Standard_Real aThreshold = Precision::Angular() * Precision::Angular() * dist2;
 154   gp_Pnt p2 = AxeOfRevolution.Location().XYZ() + dist * AxeOfRevolution.Direction().XYZ();
 155
 156   if((pl.SquareDistance(p1) < aThreshold) &&
 157      (pl.SquareDistance(p2) < aThreshold))
 158     return Standard_True;
 159   return Standard_False;
 160   //   gp_Vec V (AxeOfRevolution.Location(),theCurvePos.Location());
 161   //   if (Abs( V * theCurvePos.Direction()) <= gp::Resolution())
 162   //     return Standard_True;
 163   //   else
 164   //     return Standard_False;
 165 }
 166
 167 //=======================================================================
 168 //function : IsOriginalPnt
 169 //purpose  :
 170 //=======================================================================
 171
 172 static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
 173                                        const Extrema_POnSurf* Points,
 174                                        const Standard_Integer NbPoints)
 175 {
 176   for (Standard_Integer i=1; i<=NbPoints; i++) {
 177     if (Points[i].Value().IsEqual(P, Precision::Confusion())) {
 178       return Standard_False;
 179     }
 180   }
 181   return Standard_True;
 182 }
 183 //=======================================================================
 184 //function : IsExtremum
 185 //purpose  :
 186 //=======================================================================
 187
 188 static Standard_Boolean IsExtremum (const Standard_Real U, const Standard_Real V,
 189                                     const gp_Pnt& P, const Adaptor3d_SurfacePtr& S,
 190                                     gp_Pnt& E,        Standard_Real& Dist2,
 191                                     const Standard_Boolean IsVSup,
 192                                     const Standard_Boolean IsMin)
 193 {
 194   E = S->Value(U,V);
 195   Dist2 = P.SquareDistance(E);
 196   if (IsMin)
 197     return (Dist2 < P.SquareDistance(S->Value(U+1,V)) &&
 198             Dist2 < P.SquareDistance(S->Value(U-1,V)) &&
 199             Dist2 < P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
 200   else
 201     return (Dist2 > P.SquareDistance(S->Value(U+1,V)) &&
 202             Dist2 > P.SquareDistance(S->Value(U-1,V)) &&
 203             Dist2 > P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
 204 }
 205 //=======================================================================
 206 //function : Extrema_ExtPRevS
 207 //purpose  :
 208 //=======================================================================
 209
 210 Extrema_ExtPRevS::Extrema_ExtPRevS()
 211 {
 212   myS=NULL;
 213   myDone = Standard_False;
 214 }
 215 //=======================================================================
 216 //function : Extrema_ExtPRevS
 217 //purpose  :
 218 //=======================================================================
 219
 220 Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
 221                                    const Adaptor3d_SurfaceOfRevolution& S,
 222                                    const Standard_Real Umin,
 223                                    const Standard_Real Usup,
 224                                    const Standard_Real Vmin,
 225                                    const Standard_Real Vsup,
 226                                    const Standard_Real TolU,
 227                                    const Standard_Real TolV)
 228 {
 229   myS=NULL;
 230   Initialize (S,
 231               Umin, Usup, Vmin, Vsup,
 232               TolU, TolV);
 233   Perform(P);
 234 }
 235 //=======================================================================
 236 //function : Extrema_ExtPRevS
 237 //purpose  :
 238 //=======================================================================
 239
 240 Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
 241                                    const Adaptor3d_SurfaceOfRevolution& S,
 242                                    const Standard_Real TolU,
 243                                    const Standard_Real TolV)
 244 {
 245   myS=NULL;
 246   Initialize (S,
 247               S.FirstUParameter(), S.LastUParameter(),
 248               S.FirstVParameter(), S.LastVParameter(),
 249               TolU, TolV);
 250   Perform(P);
 251 }
 252 //=======================================================================
 253 //function : Initialize
 254 //purpose  :
 255 //=======================================================================
 256
 257 void Extrema_ExtPRevS::Initialize(const Adaptor3d_SurfaceOfRevolution& S,
 258                                   const Standard_Real Umin,
 259                                   const Standard_Real Usup,
 260                                   const Standard_Real Vmin,
 261                                   const Standard_Real Vsup,
 262                                   const Standard_Real TolU,
 263                                   const Standard_Real TolV)
 264 {
 265   myvinf=Vmin;
 266   myvsup=Vsup;
 267   mytolv=TolV;
 268
 269   Handle(Adaptor3d_HCurve) anACurve = S.BasisCurve();
 270   if (!myS || myS != (Adaptor3d_SurfacePtr)&S) {
 271     myS = Adaptor3d_SurfacePtr(&S);
 272     myPosition = GetPosition(S);
 273     myIsAnalyticallyComputable =
 274       IsCaseAnalyticallyComputable (anACurve->GetType(),myPosition,S.AxeOfRevolution());
 275   }
 276   if (!myIsAnalyticallyComputable) {
 277
 278     Standard_Integer nbu = 32, nbv = 32;
 279
 280     if(HasSingularity(S)) {nbv = 100;}
 281
 282     myExtPS.Initialize(S, nbu, nbv,
 283                        Umin, Usup, Vmin, Vsup,
 284                        TolU, TolV);
 285   }
 286 }
 287 //=======================================================================
 288 //function : Perform
 289 //purpose  :
 290 //=======================================================================
 291
 292 void Extrema_ExtPRevS::Perform(const gp_Pnt& P)
 293 {
 294   myDone = Standard_False;
 295   myNbExt = 0;
 296
 297   if (!myIsAnalyticallyComputable) {
 298
 299     myExtPS.Perform(P);
 300     myDone = myExtPS.IsDone();
 301     myNbExt = myExtPS.NbExt();
 302     return;
 303   }
 304
 305   Handle(Adaptor3d_HCurve) anACurve = myS->BasisCurve();
 306
 307   gp_Ax1 Ax = myS->AxeOfRevolution();
 308   gp_Vec Dir = Ax.Direction(), Z = myPosition.Direction();
 309   gp_Pnt O = Ax.Location();
 310
 311   Standard_Real OPdir = gp_Vec(O,P).Dot(Dir);
 312   gp_Pnt Pp = P.Translated(Dir.Multiplied(-OPdir));
 313   if (O.IsEqual(Pp,Precision::Confusion())) // P is on the AxeOfRevolution
 314     return;
 315
 316   Standard_Real U,V;
 317   gp_Pnt P1, Ppp;
 318   Standard_Real OPpz = gp_Vec(O,Pp).Dot(Z);
 319   if (Abs(OPpz) <= gp::Resolution()) {
 320     Ppp = Pp;
 321     U = 0;
 322   }
 323   else {
 324     Ppp = Pp.Translated(Z.Multiplied(-OPpz));
 325     if (O.IsEqual(Ppp,Precision::Confusion()))
 326       U = M_PI/2;
 327     else {
 328       U = gp_Vec(O,Ppp).AngleWithRef(gp_Vec(O,Pp),Dir);
 329     }
 330   }
 331
 332   gp_Vec OPpp (O,Ppp), OPq (O, myS->Value(M_PI/2,0));
 333   if (U != M_PI/2) {
 334     if (Abs(OPq.Magnitude()) <= gp::Resolution())
 335       OPq = gp_Vec(O, myS->Value(M_PI/2,anACurve->LastParameter()/10));
 336     if (OPpp.AngleWithRef(OPq,Dir) < 0)
 337       U += M_PI;
 338   }
 339
 340   gp_Trsf T;
 341   T.SetRotation(Ax, -U);
 342   P1 = P.Transformed(T);
 343
 344   gp_Pnt E;
 345   Standard_Real Dist2;
 346   Standard_Integer i;
 347
 348   Extrema_ExtPElC anExt;
 349   PerformExtPElC(anExt, P1, anACurve, mytolv);
 350
 351   if (anExt.IsDone()) {
 352     myDone = Standard_True;
 353     for (i=1; i<=anExt.NbExt(); i++) {
 354       Extrema_POnCurv POC=anExt.Point(i);
 355       V = POC.Parameter();
 356       if (V > myvsup) {
 357         //       if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
 358         //                         Standard_True, anExt.IsMin(i))) continue;
 359         Standard_Real newV = myvsup;
 360
 361         if((anACurve->GetType() == GeomAbs_Circle) ||
 362            (anACurve->GetType() == GeomAbs_Ellipse)) {
 363           newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
 364
 365           if (newV > myvsup) {
 366             newV -= 2. * M_PI;
 367
 368             if (newV + mytolv < myvinf) {
 369               newV = myvsup;
 370             } else if (newV < myvinf) {
 371               newV = myvinf;
 372             }
 373           }
 374         }
 375         V = newV;
 376
 377         if ( !IsExtremum (U, V, P, myS, E, Dist2,
 378                            Standard_True, anExt.IsMin(i))) {
 379           continue;
 380         }
 381       } else if (V < myvinf) {
 382         //      if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
 383         //                        Standard_False, anExt.IsMin(i))) continue;
 384
 385         Standard_Real newV = myvinf;
 386
 387         if((anACurve->GetType() == GeomAbs_Circle) ||
 388            (anACurve->GetType() == GeomAbs_Ellipse)) {
 389           newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
 390
 391           if(newV < myvinf) {
 392             newV += 2. * M_PI;
 393
 394             if (newV - mytolv > myvsup) {
 395               newV = myvinf;
 396             } else if (newV > myvsup) {
 397               newV = myvsup;
 398             }
 399           }
 400         }
 401         V = newV;
 402
 403         if ( !IsExtremum (U, V, P, myS, E, Dist2,
 404                           Standard_False, anExt.IsMin(i))) continue;
 405       } else {
 406         E = myS->Value(U,V);
 407         Dist2 = P.SquareDistance(E);
 408       }
 409       if (IsOriginalPnt(E, myPoint, myNbExt)) {
 410         myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
 411         mySqDist[myNbExt] = Dist2;
 412       }
 413     }
 414   }
 415   T.SetRotation(Ax, M_PI);
 416   P1.Transform(T);
 417
 418   PerformExtPElC(anExt, P1, anACurve, mytolv);
 419   if (anExt.IsDone()) {
 420     myDone = Standard_True;
 421
 422     U += M_PI;
 423
 424     for (i=1; i<=anExt.NbExt(); i++) {
 425       Extrema_POnCurv POC=anExt.Point(i);
 426       V = POC.Parameter();
 427       if (V > myvsup) {
 428         //      if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
 429         //                         Standard_True, anExt.IsMin(i))) continue;
 430
 431         Standard_Real newV = myvsup;
 432
 433         if((anACurve->GetType() == GeomAbs_Circle) ||
 434            (anACurve->GetType() == GeomAbs_Ellipse)) {
 435           newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
 436
 437           if (newV > myvsup) {
 438             newV -= 2. * M_PI;
 439
 440             if (newV + mytolv < myvinf) {
 441               newV = myvsup;
 442             } else if (newV < myvinf) {
 443               newV = myvinf;
 444             }
 445           }
 446         }
 447         V = newV;
 448
 449         if ( !IsExtremum (U, V, P, myS, E, Dist2,
 450                            Standard_True, anExt.IsMin(i))) continue;
 451       } else if (V < myvinf) {
 452         //      if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
 453         //                        Standard_False, anExt.IsMin(i))) continue;
 454         Standard_Real newV = myvinf;
 455
 456         if((anACurve->GetType() == GeomAbs_Circle) ||
 457            (anACurve->GetType() == GeomAbs_Ellipse)) {
 458           newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
 459
 460           if(newV < myvinf) {
 461             newV += 2. * M_PI;
 462
 463             if (newV - mytolv > myvsup) {
 464               newV = myvinf;
 465             } else if (newV > myvsup) {
 466               newV = myvsup;
 467             }
 468           }
 469         }
 470         V = newV;
 471
 472         if ( !IsExtremum (U, V, P, myS, E, Dist2,
 473                           Standard_False, anExt.IsMin(i))) continue;
 474       } else {
 475         E = myS->Value(U,V);
 476         Dist2 = P.SquareDistance(E);
 477       }
 478       if (IsOriginalPnt(E, myPoint, myNbExt)) {
 479         myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
 480         mySqDist[myNbExt] = Dist2;
 481
 482       }
 483     }
 484   }
 485 }
 486
 487
 488 //=======================================================================
 489 //function : IsDone
 490 //purpose  :
 491 //=======================================================================
 492
 493 Standard_Boolean Extrema_ExtPRevS::IsDone() const
 494 {
 495   return myDone;
 496 }
 497
 498
 499 //=======================================================================
 500 //function : NbExt
 501 //purpose  :
 502 //=======================================================================
 503
 504 Standard_Integer Extrema_ExtPRevS::NbExt() const
 505 {
 506   if (!IsDone()) { StdFail_NotDone::Raise(); }
 507   return myNbExt;
 508 }
 509
 510
 511 //=======================================================================
 512 //function : Value
 513 //purpose  :
 514 //=======================================================================
 515
 516 Standard_Real Extrema_ExtPRevS::SquareDistance(const Standard_Integer N) const
 517 {
 518   if (!IsDone()) { StdFail_NotDone::Raise(); }
 519   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
 520   if (myIsAnalyticallyComputable)
 521     return mySqDist[N];
 522   else
 523     return myExtPS.SquareDistance(N);
 524 }
 525 //=======================================================================
 526 //function : Point
 527 //purpose  :
 528 //=======================================================================
 529
 530 Extrema_POnSurf Extrema_ExtPRevS::Point(const Standard_Integer N) const
 531 {
 532   if (!IsDone()) { StdFail_NotDone::Raise(); }
 533   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
 534   if (myIsAnalyticallyComputable)
 535     return myPoint[N];
 536   else
 537     return myExtPS.Point(N);
 538 }
 539
 540