src/Extrema/Extrema_ExtPExtS.cxx

   1 // Created on: 1999-09-16
   2 // Created by: Edward AGAPOV
   3 // Copyright (c) 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 #include <Standard_NotImplemented.hxx>
  23 #include <Standard_OutOfRange.hxx>
  24 #include <StdFail_NotDone.hxx>
  25 #include <Adaptor3d_HCurve.hxx>
  26 #include <ElCLib.hxx>
  27 #include <Extrema_ExtPElC.hxx>
  28 #include <Extrema_ExtPExtS.hxx>
  29 #include <Extrema_POnCurv.hxx>
  30 #include <Extrema_POnSurf.hxx>
  31 #include <Precision.hxx>
  32 #include <gp.hxx>
  33 #include <gp_Ax2.hxx>
  34 #include <gp_Dir.hxx>
  35 #include <gp_Lin.hxx>
  36 #include <gp_Pln.hxx>
  37 #include <gp_Pnt.hxx>
  38 #include <gp_Vec.hxx>
  39 #include <math_FunctionSetRoot.hxx>
  40 #include <math_Vector.hxx>
  41 #include <Adaptor3d_SurfaceOfLinearExtrusion.hxx>
  42
  43 static gp_Ax2 GetPosition (const Handle(Adaptor3d_HCurve)& C);
  44
  45 static void PerformExtPElC (Extrema_ExtPElC& E,
  46                             const gp_Pnt& P,
  47                             const Handle(Adaptor3d_HCurve)& C,
  48                             const Standard_Real Tol);
  49
  50 static Standard_Boolean
  51   IsCaseAnalyticallyComputable (const GeomAbs_CurveType& theType,
  52                                 const gp_Ax2& theCurvePos,
  53                                 const gp_Dir& theSurfaceDirection) ;
  54
  55 static gp_Pnt GetValue(const Standard_Real U,
  56                        const Handle(Adaptor3d_HCurve)& C);
  57 //=======================================================================
  58 //function : Project
  59 //purpose  : Returns the projection of a point <Point> on a plane
  60 //           <ThePlane>  along a direction <TheDir>.
  61 //=======================================================================
  62 static gp_Pnt ProjectPnt(const gp_Ax2& ThePlane,
  63                          const gp_Dir& TheDir,
  64                          const gp_Pnt& Point)
  65 {
  66   gp_Vec PO(Point,ThePlane.Location());
  67   Standard_Real Alpha = PO * gp_Vec(ThePlane.Direction());
  68   Alpha /= TheDir * ThePlane.Direction();
  69   gp_Pnt P;
  70   P.SetXYZ(Point.XYZ() + Alpha * TheDir.XYZ());
  71   return P;
  72 }
  73
  74 //=======================================================================
  75 //function : IsOriginalPnt
  76 //purpose  :
  77 //=======================================================================
  78
  79 static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
  80                                        const Extrema_POnSurf* Points,
  81                                        const Standard_Integer NbPoints)
  82 {
  83   for (Standard_Integer i=1; i<=NbPoints; i++) {
  84     if (Points[i-1].Value().IsEqual(P, Precision::Confusion())) {
  85       return Standard_False;
  86     }
  87   }
  88   return Standard_True;
  89 }
  90
  91 //=======================================================================
  92 //function : MakePreciser
  93 //purpose  :
  94 //=======================================================================
  95
  96 void Extrema_ExtPExtS::MakePreciser (Standard_Real& U,
  97                                      const gp_Pnt& P,
  98                                      const Standard_Boolean isMin,
  99                                      const gp_Ax2& OrtogSection) const
 100 {
 101   if (U>myusup) {
 102     U = myusup;
 103   } else if (U<myuinf) {
 104     U = myuinf;
 105   } else {
 106
 107     Standard_Real step = (myusup - myuinf) / 30, D2e, D2next ,D2prev;
 108     gp_Pnt
 109       Pe = ProjectPnt (OrtogSection, myDirection, GetValue(U,myC)),
 110       Pprev = ProjectPnt (OrtogSection, myDirection, GetValue(U-step, myC)),
 111       Pnext = ProjectPnt (OrtogSection, myDirection, GetValue(U+step, myC));
 112     D2e = P.SquareDistance(Pe),
 113     D2next = P.SquareDistance(Pnext),
 114     D2prev = P.SquareDistance(Pprev);
 115     Standard_Boolean notFound;
 116     if (isMin)
 117       notFound =  (D2e > D2prev || D2e > D2next);
 118     else
 119       notFound = (D2e < D2prev || D2e < D2next);
 120
 121     if (notFound && (D2e < D2next && isMin)) {
 122       step = -step;
 123       D2next = D2prev;
 124       Pnext = Pprev;
 125     }
 126     while (notFound) {
 127       U = U + step;
 128       if (U > myusup) {
 129         U = myusup;
 130         break;
 131       }
 132       if (U < myuinf) {
 133         U = myuinf;
 134         break;
 135       }
 136       D2e = D2next;
 137       Pe = Pnext;
 138       Pnext = ProjectPnt (OrtogSection, myDirection, GetValue(U+step, myC));
 139       D2next = P.SquareDistance(Pnext);
 140       if (isMin)
 141         notFound = D2e > D2next;
 142       else
 143         notFound = D2e < D2next;
 144     }
 145   }
 146 }
 147 //=============================================================================
 148
 149 Extrema_ExtPExtS::Extrema_ExtPExtS ()
 150 {
 151   myDone = Standard_False;
 152 }
 153
 154 //=============================================================================
 155
 156 Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt&       P,
 157                                     const Adaptor3d_SurfaceOfLinearExtrusion&  S,
 158                                     const Standard_Real    Umin,
 159                                     const Standard_Real    Usup,
 160                                     const Standard_Real    Vmin,
 161                                     const Standard_Real    Vsup,
 162                                     const Standard_Real    TolU,
 163                                     const Standard_Real    TolV)
 164 {
 165   Initialize (S,
 166               Umin, Usup, Vmin, Vsup,
 167               TolU, TolV);
 168   Perform(P);
 169 }
 170 //=============================================================================
 171
 172 Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt&       P,
 173                                     const Adaptor3d_SurfaceOfLinearExtrusion&  S,
 174                                     const Standard_Real    TolU,
 175                                     const Standard_Real    TolV)
 176 {
 177   Initialize (S,
 178               S.FirstUParameter(), S.LastUParameter(),
 179               S.FirstVParameter(), S.LastVParameter(),
 180               TolU, TolV);
 181   Perform(P);
 182 }
 183 //=======================================================================
 184 //function : Initialize
 185 //purpose  :
 186 //=======================================================================
 187
 188 void Extrema_ExtPExtS::Initialize(const Adaptor3d_SurfaceOfLinearExtrusion& S,
 189                                   const Standard_Real Uinf,
 190                                   const Standard_Real Usup,
 191                                   const Standard_Real Vinf,
 192                                   const Standard_Real Vsup,
 193                                   const Standard_Real TolU,
 194                                   const Standard_Real TolV)
 195 {
 196   myuinf=Uinf;
 197   myusup=Usup;
 198   mytolu=TolU;
 199
 200   myvinf=Vinf;
 201   myvsup=Vsup;
 202   mytolv=TolV;
 203
 204   Handle(Adaptor3d_HCurve) anACurve = S.BasisCurve();
 205
 206   myF.Initialize(S);
 207   myC = anACurve;
 208   myPosition = GetPosition(myC);
 209   myDirection = S.Direction();
 210   myIsAnalyticallyComputable = //Standard_False;
 211     IsCaseAnalyticallyComputable (myC->GetType(),myPosition,myDirection);
 212
 213   if (!myIsAnalyticallyComputable)
 214
 215     myExtPS.Initialize(S, 32, 32,
 216                        Uinf, Usup, Vinf, Vsup,
 217                        TolU, TolV);
 218 }
 219
 220
 221 //=======================================================================
 222 //function : Perform
 223 //purpose  :
 224 //=======================================================================
 225
 226 void Extrema_ExtPExtS::Perform (const gp_Pnt& P)
 227 {
 228   myDone = Standard_False;
 229   myNbExt = 0;
 230
 231   if (!myIsAnalyticallyComputable) {
 232     myExtPS.Perform(P);
 233     myDone = myExtPS.IsDone();
 234 //  modified by NIZHNY-EAP Wed Nov 17 12:59:08 1999 ___BEGIN___
 235     myNbExt = myExtPS.NbExt();
 236 //  modified by NIZHNY-EAP Wed Nov 17 12:59:09 1999 ___END___
 237     return;
 238   }
 239
 240   gp_Pnt Pe, Pp = ProjectPnt(myPosition,myDirection,P);
 241   Extrema_ExtPElC anExt;
 242   PerformExtPElC(anExt, Pp, myC, mytolu);
 243   if (!anExt.IsDone()) return;
 244
 245   gp_Ax2 anOrtogSection (P, myDirection);
 246   Standard_Real U,V;
 247   Standard_Boolean
 248     isMin,
 249     isSimpleCase =
 250       myDirection.IsParallel(myPosition.Direction(),Precision::Angular());
 251   Standard_Integer i, aNbExt = anExt.NbExt();
 252   math_Vector UV(1,2), Tol(1,2), UVinf(1,2), UVsup(1,2);
 253   Tol(1) = mytolu;   Tol(2) = mytolv;
 254   UVinf(1) = myuinf; UVinf(2) = myvinf;
 255   UVsup(1) = myusup; UVsup(2) = myvsup;
 256
 257
 258   for (i=1; i<=aNbExt; i++) {
 259     Extrema_POnCurv POC=anExt.Point(i);
 260     U = POC.Parameter();
 261     //// modified by jgv, 23.12.2008 for OCC17194 ////
 262     if (myC->IsPeriodic())
 263       {
 264         Standard_Real U2 = U;
 265         ElCLib::AdjustPeriodic(myuinf, myuinf + 2.*M_PI, Precision::PConfusion(), U, U2);
 266       }
 267     //////////////////////////////////////////////////
 268     gp_Pnt E = POC.Value();
 269     Pe = ProjectPnt(anOrtogSection, myDirection, E);
 270
 271     if (isSimpleCase) {
 272       V = gp_Vec(E,Pe) * gp_Vec(myDirection);
 273       // modified by NIZHNY-MKK  Thu Sep 18 14:46:14 2003.BEGIN
 274       //       myPoint[++myNbExt] = Extrema_POnSurf(U, V, Pe);
 275       //       myValue[myNbExt] = anExt.Value(i);
 276       myPoint[myNbExt] = Extrema_POnSurf(U, V, Pe);
 277       mySqDist[myNbExt] = anExt.SquareDistance(i);
 278       myNbExt++;
 279       // modified by NIZHNY-MKK  Thu Sep 18 14:46:18 2003.END
 280     }
 281     else {
 282       myF.SetPoint(P);
 283       isMin = anExt.IsMin(i);//( Pp.Distance(GetValue(U+10,myC)) > anExt.Value(i) );
 284
 285       MakePreciser(U, P, isMin, anOrtogSection);
 286       E = GetValue(U, myC);
 287       Pe = ProjectPnt (anOrtogSection, myDirection, E),
 288       V = gp_Vec(E,Pe) * gp_Vec(myDirection);
 289       UV(1) = U;         UV(2) = V;
 290       math_FunctionSetRoot aFSR (myF,UV,Tol,UVinf,UVsup);
 291 //      for (Standard_Integer k=1 ; k <= myF.NbExt();
 292       Standard_Integer k;
 293       for ( k=1 ; k <= myF.NbExt(); k++) {
 294         if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
 295           // modified by NIZHNY-MKK  Thu Sep 18 14:46:41 2003.BEGIN
 296           //      myPoint[++myNbExt] = myF.Point(k);
 297           //      myValue[myNbExt] = myF.Value(k);
 298           myPoint[myNbExt] = myF.Point(k);
 299           mySqDist[myNbExt] = myF.SquareDistance(k);
 300           myNbExt++;
 301           // modified by NIZHNY-MKK  Thu Sep 18 14:46:43 2003.END
 302         }
 303       }
 304       // try symmetric point
 305       U *= -1;
 306       MakePreciser(U, P, isMin, anOrtogSection);
 307       E = GetValue(U, myC);
 308       Pe = ProjectPnt (anOrtogSection, myDirection, E),
 309       V = gp_Vec(E,Pe) * gp_Vec(myDirection);
 310       UV(1) = U; UV(2) = V;
 311
 312       aFSR.Perform (myF,UV,UVinf,UVsup);
 313
 314       for (k=1 ; k <= myF.NbExt(); k++) {
 315         if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
 316           // modified by NIZHNY-MKK  Thu Sep 18 14:46:59 2003.BEGIN
 317           //      myPoint[++myNbExt] = myF.Point(k);
 318           //      myValue[myNbExt] = myF.Value(k);
 319           myPoint[myNbExt] = myF.Point(k);
 320           mySqDist[myNbExt] = myF.SquareDistance(k);
 321           myNbExt++;
 322           // modified by NIZHNY-MKK  Thu Sep 18 14:47:04 2003.END
 323         }
 324       }
 325     }
 326   }
 327   myDone = Standard_True;
 328   return;
 329 }
 330
 331 //=============================================================================
 332
 333 Standard_Boolean Extrema_ExtPExtS::IsDone () const { return myDone; }
 334 //=============================================================================
 335
 336 Standard_Integer Extrema_ExtPExtS::NbExt () const
 337 {
 338   if (!IsDone()) { StdFail_NotDone::Raise(); }
 339   if (myIsAnalyticallyComputable)
 340     return myNbExt;
 341   else
 342     return myExtPS.NbExt();
 343 }
 344 //=============================================================================
 345
 346 Standard_Real Extrema_ExtPExtS::SquareDistance (const Standard_Integer N) const
 347 {
 348   if (!IsDone()) { StdFail_NotDone::Raise(); }
 349   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
 350   if (myIsAnalyticallyComputable)
 351     // modified by NIZHNY-MKK  Thu Sep 18 14:48:39 2003.BEGIN
 352     //     return myValue[N];
 353     return mySqDist[N-1];
 354   // modified by NIZHNY-MKK  Thu Sep 18 14:48:42 2003.END
 355   else
 356     return myExtPS.SquareDistance(N);
 357 }
 358 //=============================================================================
 359
 360 Extrema_POnSurf Extrema_ExtPExtS::Point (const Standard_Integer N) const
 361 {
 362   if (!IsDone()) { StdFail_NotDone::Raise(); }
 363   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
 364   if (myIsAnalyticallyComputable) {
 365     // modified by NIZHNY-MKK  Thu Sep 18 14:47:40 2003.BEGIN
 366     //     return myPoint[N];
 367     return myPoint[N-1];
 368   }
 369   // modified by NIZHNY-MKK  Thu Sep 18 14:47:43 2003.END
 370   else
 371     return myExtPS.Point(N);
 372 }
 373 //=============================================================================
 374
 375
 376 static gp_Ax2 GetPosition (const Handle(Adaptor3d_HCurve)& C)
 377 {
 378   switch (C->GetType()) {
 379   case GeomAbs_Line: {
 380     gp_Lin L = C->Line();
 381     gp_Pln Pln = gp_Pln(L.Location(), L.Direction());
 382     //:abv 30.05.02: OCC  - use constructor instead of Set...s() to avoid exception
 383     gp_Ax2 Pos ( Pln.Location(), Pln.Position().Direction(), Pln.Position().XDirection() );
 384 //     Pos.SetAxis(Pln.XAxis());
 385 //     Pos.SetXDirection(Pln.YAxis().Direction());
 386 //     Pos.SetYDirection(Pln.Position().Direction());
 387     return Pos;
 388   }
 389   case GeomAbs_Circle:
 390     return  C->Circle().Position();
 391   case GeomAbs_Ellipse:
 392     return  C->Ellipse().Position();
 393   case GeomAbs_Hyperbola:
 394     return  C->Hyperbola().Position();
 395   case GeomAbs_Parabola:
 396     return  C->Parabola().Position();
 397   default:
 398     return gp_Ax2 ();
 399   }
 400 }
 401 //=============================================================================
 402
 403 static void PerformExtPElC (Extrema_ExtPElC& E,
 404                             const gp_Pnt& P,
 405                             const Handle(Adaptor3d_HCurve)& C,
 406                             const Standard_Real Tol)
 407 {
 408   switch (C->GetType()) {
 409   case GeomAbs_Hyperbola:
 410     E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
 411     return;
 412   case GeomAbs_Line:
 413     E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
 414     return;
 415   case GeomAbs_Circle:
 416     E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
 417     return;
 418   case GeomAbs_Ellipse:
 419     E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
 420     return;
 421   case GeomAbs_Parabola:
 422     E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
 423     return;
 424 #ifndef DEB
 425   default:
 426     return;
 427 #endif
 428   }
 429 }
 430
 431 //=======================================================================
 432 //function : IsCaseAnalyticallyComputable
 433 //purpose  :
 434 //=======================================================================
 435
 436 static Standard_Boolean IsCaseAnalyticallyComputable
 437   (const GeomAbs_CurveType& theType,
 438    const gp_Ax2& theCurvePos,
 439    const gp_Dir& theSurfaceDirection)
 440 {
 441   // check type
 442   switch (theType) {
 443   case GeomAbs_Line:
 444   case GeomAbs_Circle:
 445   case GeomAbs_Ellipse:
 446   case GeomAbs_Hyperbola:
 447   case GeomAbs_Parabola:
 448     break;
 449   default:
 450     return  Standard_False;
 451   }
 452   // check if it is a plane
 453   if (Abs(theCurvePos.Direction() * theSurfaceDirection) <= gp::Resolution())
 454     return Standard_False;
 455   else
 456     return Standard_True;
 457 }
 458 //=======================================================================
 459 //function : GetValue
 460 //purpose  :
 461 //=======================================================================
 462
 463 static gp_Pnt GetValue(const Standard_Real U,
 464                        const Handle(Adaptor3d_HCurve)& C)
 465 {
 466   switch (C->GetType()) {
 467   case GeomAbs_Line:
 468     return ElCLib::Value(U, C->Line());
 469   case GeomAbs_Circle:
 470     return ElCLib::Value(U, C->Circle());
 471   case GeomAbs_Ellipse:
 472     return ElCLib::Value(U, C->Ellipse());
 473   case GeomAbs_Hyperbola:
 474     return ElCLib::Value(U, C->Hyperbola());
 475   case GeomAbs_Parabola:
 476     return ElCLib::Value(U, C->Parabola());
 477   default:
 478     return gp_Pnt();
 479   }
 480 }
 481 //=======================================================================
 482 //function : GetU
 483 //purpose  :
 484 //=======================================================================
 485 //#ifdef DEB
 486 //static Standard_Real GetU(const gp_Vec& vec,
 487 //                        const gp_Pnt& P,
 488 //                        const Handle(Adaptor3d_HCurve)& C)
 489 //{
 490 //  switch (C->GetType()) {
 491 //  case GeomAbs_Line:
 492 //    return ElCLib::Parameter(C->Line().Translated(vec), P);
 493 //  case GeomAbs_Circle:
 494 //    return ElCLib::Parameter(C->Circle().Translated(vec), P);
 495 //  case GeomAbs_Ellipse:
 496 //    return ElCLib::Parameter(C->Ellipse().Translated(vec), P);
 497 //  case GeomAbs_Hyperbola:
 498 //    return ElCLib::Parameter(C->Hyperbola().Translated(vec), P);
 499 //  case GeomAbs_Parabola:
 500 //    return ElCLib::Parameter(C->Parabola().Translated(vec), P);
 501 //  default:
 502 //    return 0;
 503 //  }
 504 //}
 505 //#endif