src/Extrema/Extrema_ExtPExtS.cxx

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