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