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