src/Extrema/Extrema_GenExtPS.cxx

   1 // Created on: 1995-07-18
   2 // Created by: Modelistation
   3 // Copyright (c) 1995-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 //  Modified by skv - Thu Sep 30 15:21:07 2004 OCC593
  23
  24
  25 #include <Extrema_GenExtPS.ixx>
  26 #include <StdFail_NotDone.hxx>
  27 #include <Standard_OutOfRange.hxx>
  28 #include <TColStd_Array2OfInteger.hxx>
  29 #include <TColStd_Array2OfReal.hxx>
  30 #include <TColgp_Array2OfPnt.hxx>
  31 #include <math_FunctionSetRoot.hxx>
  32 #include <math_Vector.hxx>
  33 #include <math_NewtonFunctionSetRoot.hxx>
  34 #include <GeomAbs_IsoType.hxx>
  35 #include <Bnd_Sphere.hxx>
  36 #include <Extrema_HUBTreeOfSphere.hxx>
  37 #include <Extrema_ExtFlag.hxx>
  38 #include <Bnd_Array1OfSphere.hxx>
  39 #include <Bnd_HArray1OfSphere.hxx>
  40 #include <Precision.hxx>
  41 #include <Geom_OffsetSurface.hxx>
  42 #include <Geom_RectangularTrimmedSurface.hxx>
  43 #include <Geom_BSplineSurface.hxx>
  44 #include <Geom_BezierSurface.hxx>
  45 #include <Adaptor3d_HSurface.hxx>
  46 #include <Adaptor3d_HCurve.hxx>
  47 #include <Adaptor3d_Curve.hxx>
  48 #include <Geom_BSplineCurve.hxx>
  49 #include <Geom_BezierCurve.hxx>
  50 //IMPLEMENT_HARRAY1(Extrema_HArray1OfSphere)
  51
  52
  53 class Bnd_SphereUBTreeSelector : public Extrema_UBTreeOfSphere::Selector
  54 {
  55  public:
  56
  57   Bnd_SphereUBTreeSelector (const Handle(Bnd_HArray1OfSphere)& theSphereArray,
  58 Bnd_Sphere& theSol)
  59     : myXYZ(0,0,0),
  60       mySphereArray(theSphereArray),
  61       mySol(theSol)
  62   {
  63     //myXYZ = gp_Pnt(0, 0, 0);
  64   }
  65
  66   void DefineCheckPoint( const gp_Pnt& theXYZ )
  67   { myXYZ = theXYZ; }
  68
  69   Bnd_Sphere& Sphere() const
  70   { return mySol; }
  71
  72   virtual Standard_Boolean Reject( const Bnd_Sphere &theBnd ) const = 0;
  73
  74   virtual Standard_Boolean Accept(const Standard_Integer& theObj) = 0;
  75
  76  protected:
  77   gp_Pnt                      myXYZ;
  78   const Handle(Bnd_HArray1OfSphere)& mySphereArray;
  79   Bnd_Sphere&                                            mySol;
  80
  81 };
  82
  83 class Bnd_SphereUBTreeSelectorMin : public Bnd_SphereUBTreeSelector
  84 {
  85 public:
  86   Bnd_SphereUBTreeSelectorMin (const Handle(Bnd_HArray1OfSphere)& theSphereArray,
  87                 Bnd_Sphere& theSol)
  88                 : Bnd_SphereUBTreeSelector(theSphereArray, theSol),
  89                   myMinDist(RealLast())
  90   {}
  91
  92   void SetMinDist( const Standard_Real theMinDist )
  93   { myMinDist = theMinDist; }
  94
  95   Standard_Real MinDist() const
  96   { return myMinDist; }
  97
  98   Standard_Boolean Reject( const Bnd_Sphere &theBnd ) const
  99   {
 100     Bnd_SphereUBTreeSelectorMin* me =
 101       const_cast<Bnd_SphereUBTreeSelectorMin*>(this);
 102     // myMinDist is decreased each time a nearer object is found
 103     return theBnd.IsOut( myXYZ.XYZ(), me->myMinDist );
 104   }
 105
 106   Standard_Boolean Accept(const Standard_Integer&);
 107
 108 private:
 109         Standard_Real   myMinDist;
 110 };
 111
 112 Standard_Boolean Bnd_SphereUBTreeSelectorMin::Accept(const Standard_Integer& theInd)
 113 {
 114   const Bnd_Sphere& aSph = mySphereArray->Value(theInd);
 115   Standard_Real aCurDist;
 116
 117     if ( (aCurDist = aSph.SquareDistance(myXYZ.XYZ())) < mySol.SquareDistance(myXYZ.XYZ()) )
 118     {
 119       mySol = aSph;
 120       if ( aCurDist < myMinDist )
 121         myMinDist = aCurDist;
 122
 123       return Standard_True;
 124     }
 125
 126   return Standard_False;
 127 }
 128
 129 class Bnd_SphereUBTreeSelectorMax : public Bnd_SphereUBTreeSelector
 130 {
 131 public:
 132   Bnd_SphereUBTreeSelectorMax (const Handle(Bnd_HArray1OfSphere)& theSphereArray,
 133                 Bnd_Sphere& theSol)
 134                 : Bnd_SphereUBTreeSelector(theSphereArray, theSol),
 135                   myMaxDist(0)
 136   {}
 137
 138   void SetMaxDist( const Standard_Real theMaxDist )
 139   { myMaxDist = theMaxDist; }
 140
 141   Standard_Real MaxDist() const
 142   { return myMaxDist; }
 143
 144   Standard_Boolean Reject( const Bnd_Sphere &theBnd ) const
 145   {
 146     Bnd_SphereUBTreeSelectorMax* me =
 147       const_cast<Bnd_SphereUBTreeSelectorMax*>(this);
 148     // myMaxDist is decreased each time a nearer object is found
 149     return theBnd.IsOut( myXYZ.XYZ(), me->myMaxDist );
 150   }
 151
 152   Standard_Boolean Accept(const Standard_Integer&);
 153
 154 private:
 155         Standard_Real   myMaxDist;
 156 };
 157
 158 Standard_Boolean Bnd_SphereUBTreeSelectorMax::Accept(const Standard_Integer& theInd)
 159 {
 160   const Bnd_Sphere& aSph = mySphereArray->Value(theInd);
 161   Standard_Real aCurDist;
 162
 163     if ( (aCurDist = aSph.SquareDistance(myXYZ.XYZ())) > mySol.SquareDistance(myXYZ.XYZ()) )
 164     {
 165       mySol = aSph;
 166       if ( aCurDist > myMaxDist )
 167         myMaxDist = aCurDist;
 168
 169       return Standard_True;
 170     }
 171
 172   return Standard_False;
 173 }
 174
 175 /*
 176  * This function computes the point on surface parameters on edge.
 177  * if it coincides with theParam0 or theParam1, it is returned.
 178  */
 179 static Extrema_POnSurfParams ComputeEdgeParameters
 180       (const Standard_Boolean       IsUEdge,
 181        const Extrema_POnSurfParams &theParam0,
 182        const Extrema_POnSurfParams &theParam1,
 183        const Adaptor3d_SurfacePtr  &theSurf,
 184        const gp_Pnt                &thePoint,
 185        const Standard_Real          theDiffTol)
 186 {
 187   const Standard_Real aSqrDist01 =
 188     theParam0.Value().SquareDistance(theParam1.Value());
 189
 190   if (aSqrDist01 <= theDiffTol) {
 191     // The points are confused. Get the first point and change its type.
 192     return theParam0;
 193   } else {
 194     const Standard_Real aDiffDist =
 195       Abs(theParam0.GetSqrDistance() - theParam1.GetSqrDistance());
 196
 197     if (aDiffDist >= aSqrDist01 - theDiffTol) {
 198       // The shortest distance is one of the nodes.
 199       if (theParam0.GetSqrDistance() > theParam1.GetSqrDistance()) {
 200         // The shortest distance is the point 1.
 201         return theParam1;
 202       } else {
 203         // The shortest distance is the point 0.
 204         return theParam0;
 205       }
 206     } else {
 207       // The shortest distance is inside the edge.
 208       gp_XYZ aPoP(thePoint.XYZ().Subtracted(theParam0.Value().XYZ()));
 209       gp_XYZ aPoP1(theParam1.Value().XYZ().Subtracted(theParam0.Value().XYZ()));
 210       Standard_Real aRatio = aPoP.Dot(aPoP1)/aSqrDist01;
 211       Standard_Real aU[2];
 212       Standard_Real aV[2];
 213
 214       theParam0.Parameter(aU[0], aV[0]);
 215       theParam1.Parameter(aU[1], aV[1]);
 216
 217       Standard_Real aUPar = aU[0];
 218       Standard_Real aVPar = aV[0];
 219
 220       if (IsUEdge) {
 221         aUPar += aRatio*(aU[1] - aU[0]);
 222       } else {
 223         aVPar += aRatio*(aV[1] - aV[0]);
 224       }
 225
 226       Extrema_POnSurfParams aParam(aUPar, aVPar, theSurf->Value(aUPar, aVPar));
 227       Standard_Integer anIndices[2];
 228
 229       theParam0.GetIndices(anIndices[0], anIndices[1]);
 230       aParam.SetElementType(IsUEdge ? Extrema_UIsoEdge : Extrema_VIsoEdge);
 231       aParam.SetSqrDistance(thePoint.SquareDistance(aParam.Value()));
 232       aParam.SetIndices(anIndices[0], anIndices[1]);
 233
 234       return aParam;
 235     }
 236   }
 237 }
 238
 239 //=============================================================================
 240
 241 /*-----------------------------------------------------------------------------
 242 Function:
 243    Find all extremum distances between point P and surface
 244   S using sampling (NbU,NbV).
 245
 246 Method:
 247    The algorithm bases on the hypothesis that sampling is precise enough,
 248   if there exist N extreme distances between the point and the surface,
 249   so there also exist N extrema between the point and the grid.
 250   So, the algorithm consists in starting from extrema of the grid to find the
 251   extrema of the surface.
 252   The extrema are calculated by the algorithm math_FunctionSetRoot with the
 253   following arguments:
 254   - F: Extrema_FuncExtPS created from P and S,
 255   - UV: math_Vector the components which of are parameters of the extremum on the
 256     grid,
 257   - Tol: Min(TolU,TolV), (Prov.:math_FunctionSetRoot does not autorize a vector)
 258   - UVinf: math_Vector the components which of are lower limits of u and v,
 259   - UVsup: math_Vector the components which of are upper limits of u and v.
 260
 261 Processing:
 262   a- Creation of the table of distances (TbDist(0,NbU+1,0,NbV+1)):
 263      The table is expanded at will; lines 0 and NbU+1 and
 264      columns 0 and NbV+1 are initialized at RealFirst() or RealLast()
 265      to simplify the tests carried out at stage b
 266      (there is no need to test if the point is on border of the grid).
 267   b- Calculation of extrema:
 268      First the minimums and then the maximums are found. These 2 procedured
 269      pass in a similar way:
 270   b.a- Initialization:
 271       - 'borders' of table  TbDist (RealLast() in case of minimums
 272         and  RealLast() in case of maximums),
 273       - table TbSel(0,NbU+1,0,NbV+1) of selection of points for
 274         calculation of local extremum (0). When a point will selected,
 275         it will not be selectable, as well as the ajacent points
 276         (8 at least). The corresponding addresses will be set to 1.
 277   b.b- Calculation of minimums (or maximums):
 278        All distances from table TbDist are parsed in a loop:
 279       - search minimum (or maximum) in the grid,
 280       - calculate extremum on the surface,
 281       - update table TbSel.
 282 -----------------------------------------------------------------------------*/
 283
 284 static Standard_Boolean IsoIsDeg  (const Adaptor3d_Surface& S,
 285                                    const Standard_Real      Param,
 286                                    const GeomAbs_IsoType    IT,
 287                                    const Standard_Real      TolMin,
 288                                    const Standard_Real      TolMax)
 289 {
 290     Standard_Real U1=0.,U2=0.,V1=0.,V2=0.,T;
 291     Standard_Boolean Along = Standard_True;
 292     U1 = S.FirstUParameter();
 293     U2 = S.LastUParameter();
 294     V1 = S.FirstVParameter();
 295     V2 = S.LastVParameter();
 296     gp_Vec D1U,D1V;
 297     gp_Pnt P;
 298     Standard_Real Step,D1NormMax;
 299     if (IT == GeomAbs_IsoV)
 300     {
 301       Step = (U2 - U1)/10;
 302       D1NormMax=0.;
 303       for (T=U1;T<=U2;T=T+Step)
 304       {
 305         S.D1(T,Param,P,D1U,D1V);
 306         D1NormMax=Max(D1NormMax,D1U.Magnitude());
 307       }
 308
 309       if (D1NormMax >TolMax || D1NormMax < TolMin )
 310            Along = Standard_False;
 311     }
 312     else
 313     {
 314       Step = (V2 - V1)/10;
 315       D1NormMax=0.;
 316       for (T=V1;T<=V2;T=T+Step)
 317       {
 318         S.D1(Param,T,P,D1U,D1V);
 319         D1NormMax=Max(D1NormMax,D1V.Magnitude());
 320       }
 321
 322       if (D1NormMax >TolMax || D1NormMax < TolMin )
 323            Along = Standard_False;
 324
 325
 326     }
 327     return Along;
 328 }
 329 //----------------------------------------------------------
 330 Extrema_GenExtPS::Extrema_GenExtPS()
 331 {
 332   myDone = Standard_False;
 333   myInit = Standard_False;
 334   myFlag = Extrema_ExtFlag_MINMAX;
 335   myAlgo = Extrema_ExtAlgo_Grad;
 336 }
 337
 338
 339 Extrema_GenExtPS::Extrema_GenExtPS (const gp_Pnt&          P,
 340                               const Adaptor3d_Surface& S,
 341                               const Standard_Integer NbU,
 342                               const Standard_Integer NbV,
 343                               const Standard_Real    TolU,
 344                               const Standard_Real    TolV,
 345                               const Extrema_ExtFlag F,
 346                               const Extrema_ExtAlgo A)
 347   : myF (P,S), myFlag(F), myAlgo(A)
 348 {
 349   Initialize(S, NbU, NbV, TolU, TolV);
 350   Perform(P);
 351 }
 352
 353 Extrema_GenExtPS::Extrema_GenExtPS (const gp_Pnt&          P,
 354                               const Adaptor3d_Surface& S,
 355                               const Standard_Integer NbU,
 356                               const Standard_Integer NbV,
 357                               const Standard_Real    Umin,
 358                               const Standard_Real    Usup,
 359                               const Standard_Real    Vmin,
 360                               const Standard_Real    Vsup,
 361                               const Standard_Real    TolU,
 362                               const Standard_Real    TolV,
 363                               const Extrema_ExtFlag F,
 364                               const Extrema_ExtAlgo A)
 365   : myF (P,S), myFlag(F), myAlgo(A)
 366 {
 367   Initialize(S, NbU, NbV, Umin, Usup, Vmin, Vsup, TolU, TolV);
 368   Perform(P);
 369 }
 370
 371
 372 void Extrema_GenExtPS::Initialize(const Adaptor3d_Surface& S,
 373                                   const Standard_Integer NbU,
 374                                   const Standard_Integer NbV,
 375                                   const Standard_Real    TolU,
 376                                   const Standard_Real    TolV)
 377 {
 378   myumin = S.FirstUParameter();
 379   myusup = S.LastUParameter();
 380   myvmin = S.FirstVParameter();
 381   myvsup = S.LastVParameter();
 382   Initialize(S,NbU,NbV,myumin,myusup,myvmin,myvsup,TolU,TolV);
 383 }
 384
 385
 386 void Extrema_GenExtPS::Initialize(const Adaptor3d_Surface& S,
 387                                   const Standard_Integer NbU,
 388                                   const Standard_Integer NbV,
 389                                   const Standard_Real    Umin,
 390                                   const Standard_Real    Usup,
 391                                   const Standard_Real    Vmin,
 392                                   const Standard_Real    Vsup,
 393                                   const Standard_Real    TolU,
 394                                   const Standard_Real    TolV)
 395 {
 396   myS = (Adaptor3d_SurfacePtr)&S;
 397   myusample = NbU;
 398   myvsample = NbV;
 399   mytolu = TolU;
 400   mytolv = TolV;
 401   myumin = Umin;
 402   myusup = Usup;
 403   myvmin = Vmin;
 404   myvsup = Vsup;
 405
 406   if ((myusample < 2) ||
 407       (myvsample < 2)) { Standard_OutOfRange::Raise(); }
 408
 409   myF.Initialize(S);
 410
 411   mySphereUBTree.Nullify();
 412   myUParams.Nullify();
 413   myVParams.Nullify();
 414   myInit = Standard_False;
 415 }
 416
 417 inline static void fillParams(const TColStd_Array1OfReal& theKnots,
 418                               Standard_Integer theDegree,
 419                               Standard_Real theParMin,
 420                               Standard_Real theParMax,
 421                               Handle_TColStd_HArray1OfReal& theParams,
 422                               Standard_Integer theSample)
 423 {
 424   NCollection_Vector<Standard_Real> aParams;
 425   Standard_Integer i = 1;
 426   Standard_Real aPrevPar = theParMin;
 427   aParams.Append(aPrevPar);
 428   //calculation the array of parametric points depending on the knots array variation and degree of given surface
 429   for ( ; i <  theKnots.Length() && theKnots(i) < (theParMax - Precision::PConfusion()); i++ )
 430   {
 431     if (  theKnots(i+1) < theParMin + Precision::PConfusion())
 432       continue;
 433
 434     Standard_Real aStep = (theKnots(i+1) - theKnots(i))/Max(theDegree,2);
 435     Standard_Integer k =1;
 436     for( ; k <= theDegree ; k++)
 437     {
 438       Standard_Real aPar = theKnots(i) + k * aStep;
 439       if(aPar > theParMax - Precision::PConfusion())
 440         break;
 441       if(aPar > aPrevPar + Precision::PConfusion() )
 442       {
 443         aParams.Append(aPar);
 444         aPrevPar = aPar;
 445       }
 446     }
 447
 448   }
 449   aParams.Append(theParMax);
 450   Standard_Integer nbPar = aParams.Length();
 451   //in case of an insufficient number of points the grid will be built later
 452   if (nbPar < theSample)
 453     return;
 454   theParams = new TColStd_HArray1OfReal(1, nbPar );
 455   for( i = 0; i < nbPar; i++)
 456     theParams->SetValue(i+1,aParams(i));
 457 }
 458
 459 void Extrema_GenExtPS::GetGridPoints( const Adaptor3d_Surface& theSurf)
 460 {
 461   //creation parametric points for BSpline and Bezier surfaces
 462   //with taking into account of Degree and NbKnots of BSpline or Bezier geometry
 463   if(theSurf.GetType() == GeomAbs_OffsetSurface)
 464     GetGridPoints(theSurf.BasisSurface()->Surface());
 465   //parametric points for BSpline surfaces
 466   else if( theSurf.GetType() == GeomAbs_BSplineSurface)
 467   {
 468     Handle(Geom_BSplineSurface) aBspl = theSurf.BSpline();
 469     if(!aBspl.IsNull())
 470     {
 471       TColStd_Array1OfReal aUKnots(1, aBspl->NbUKnots());
 472       aBspl->UKnots( aUKnots);
 473       TColStd_Array1OfReal aVKnots(1, aBspl->NbVKnots());
 474       aBspl->VKnots( aVKnots);
 475       fillParams(aUKnots,aBspl->UDegree(),myumin, myusup, myUParams, myusample);
 476       fillParams(aVKnots,aBspl->VDegree(),myvmin, myvsup, myVParams, myvsample);
 477     }
 478   }
 479   //calculation parametric points for Bezier surfaces
 480   else if(theSurf.GetType() == GeomAbs_BezierSurface)
 481   {
 482     Handle(Geom_BezierSurface) aBezier = theSurf.Bezier();
 483     if(aBezier.IsNull())
 484       return;
 485
 486     TColStd_Array1OfReal aUKnots(1,2);
 487     TColStd_Array1OfReal aVKnots(1,2);
 488     aBezier->Bounds(aUKnots(1), aUKnots(2), aVKnots(1), aVKnots(2));
 489     fillParams(aUKnots, aBezier->UDegree(), myumin, myusup, myUParams, myusample);
 490     fillParams(aVKnots, aBezier->VDegree(), myvmin, myvsup, myVParams, myvsample);
 491
 492   }
 493   //creation points for surfaces based on BSpline or Bezier curves
 494   else if(theSurf.GetType() == GeomAbs_SurfaceOfRevolution ||
 495     theSurf.GetType() == GeomAbs_SurfaceOfExtrusion)
 496   {
 497     Handle(TColStd_HArray1OfReal) anArrKnots;
 498     Standard_Integer aDegree = 0;
 499     if(theSurf.BasisCurve()->Curve().GetType() == GeomAbs_BSplineCurve)
 500     {
 501       Handle(Geom_BSplineCurve) aBspl = theSurf.BasisCurve()->Curve().BSpline();
 502       if(!aBspl.IsNull())
 503       {
 504         anArrKnots = new TColStd_HArray1OfReal(1,aBspl->NbKnots());
 505         aBspl->Knots( anArrKnots->ChangeArray1() );
 506         aDegree = aBspl->Degree();
 507
 508       }
 509
 510     }
 511     if(theSurf.BasisCurve()->Curve().GetType() == GeomAbs_BezierCurve)
 512     {
 513       Handle(Geom_BezierCurve) aBez = theSurf.BasisCurve()->Curve().Bezier();
 514       if(!aBez.IsNull())
 515       {
 516         anArrKnots = new TColStd_HArray1OfReal(1,2);
 517         anArrKnots->SetValue(1, aBez->FirstParameter());
 518         anArrKnots->SetValue(2, aBez->LastParameter());
 519         aDegree = aBez->Degree();
 520
 521       }
 522     }
 523     if(anArrKnots.IsNull())
 524       return;
 525     if( theSurf.GetType() == GeomAbs_SurfaceOfRevolution )
 526       fillParams( anArrKnots->Array1(), aDegree, myvmin, myvsup, myVParams, myvsample);
 527     else
 528       fillParams( anArrKnots->Array1(), aDegree, myumin, myusup, myUParams, myusample);
 529   }
 530   //update the number of points in sample
 531   if(!myUParams.IsNull())
 532     myusample = myUParams->Length();
 533   if( !myVParams.IsNull())
 534     myvsample = myVParams->Length();
 535
 536 }
 537
 538 void Extrema_GenExtPS::BuildGrid(const gp_Pnt &thePoint)
 539 {
 540   Standard_Integer NoU, NoV;
 541
 542   //if grid was already built skip its creation
 543   if (!myInit) {
 544     //build parametric grid in case of a complex surface geometry (BSpline and Bezier surfaces)
 545     GetGridPoints(*myS);
 546
 547     //build grid in other cases
 548     if( myUParams.IsNull() )
 549     {
 550       Standard_Real PasU = myusup - myumin;
 551       Standard_Real U0 = PasU / myusample / 100.;
 552       PasU = (PasU - U0) / (myusample - 1);
 553       U0 = U0/2. + myumin;
 554       myUParams = new TColStd_HArray1OfReal(1,myusample );
 555       Standard_Integer NoU;
 556       Standard_Real U = U0;
 557       for ( NoU = 1 ; NoU <= myusample; NoU++, U += PasU)
 558         myUParams->SetValue(NoU, U);
 559     }
 560
 561     if( myVParams.IsNull())
 562     {
 563       Standard_Real PasV = myvsup - myvmin;
 564       Standard_Real V0 = PasV / myvsample / 100.;
 565       PasV = (PasV - V0) / (myvsample - 1);
 566       V0 = V0/2. + myvmin;
 567
 568       myVParams = new TColStd_HArray1OfReal(1,myvsample );
 569       Standard_Integer  NoV;
 570       Standard_Real V = V0;
 571
 572       for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV)
 573         myVParams->SetValue(NoV, V);
 574     }
 575
 576     //If flag was changed and extrema not reinitialized Extrema would fail
 577     myPoints = new Extrema_HArray2OfPOnSurfParams
 578       (0, myusample + 1, 0, myvsample + 1);
 579     // Calculation of distances
 580
 581     for ( NoU = 1 ; NoU <= myusample; NoU++ ) {
 582       for ( NoV = 1 ; NoV <= myvsample; NoV++) {
 583         gp_Pnt aP1 = myS->Value(myUParams->Value(NoU), myVParams->Value(NoV));
 584         Extrema_POnSurfParams aParam
 585           (myUParams->Value(NoU), myVParams->Value(NoV), aP1);
 586
 587         aParam.SetElementType(Extrema_Node);
 588         aParam.SetIndices(NoU, NoV);
 589         myPoints->SetValue(NoU, NoV, aParam);
 590       }
 591     }
 592
 593     // Fill boundary with negative square distance.
 594     // It is used for computation of Maximum.
 595     for (NoV = 0; NoV <= myvsample + 1; NoV++) {
 596       myPoints->ChangeValue(0, NoV).SetSqrDistance(-1.);
 597       myPoints->ChangeValue(myusample + 1, NoV).SetSqrDistance(-1.);
 598     }
 599
 600     for (NoU = 1; NoU <= myusample; NoU++) {
 601       myPoints->ChangeValue(NoU, 0).SetSqrDistance(-1.);
 602       myPoints->ChangeValue(NoU, myvsample + 1).SetSqrDistance(-1.);
 603     }
 604
 605     myInit = Standard_True;
 606   }
 607
 608   // Compute distances to mesh.
 609   // Step 1. Compute distances to nodes.
 610   for ( NoU = 1 ; NoU <= myusample; NoU++ ) {
 611     for ( NoV = 1 ; NoV <= myvsample; NoV++) {
 612       Extrema_POnSurfParams &aParam = myPoints->ChangeValue(NoU, NoV);
 613
 614       aParam.SetSqrDistance(thePoint.SquareDistance(aParam.Value()));
 615     }
 616   }
 617
 618   // For search of minimum compute distances to mesh.
 619   if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX)
 620   {
 621     // This is the tolerance of difference of squared values.
 622     // No need to set it too small.
 623     const Standard_Real aDiffTol = mytolu + mytolv;
 624
 625     // Step 2. Compute distances to edges.
 626     // Assume UEdge(i, j) = { Point(i, j); Point(i + 1, j    ) }
 627     // Assume VEdge(i, j) = { Point(i, j); Point(i,     j + 1) }
 628     Handle(Extrema_HArray2OfPOnSurfParams) aUEdgePntParams =
 629           new Extrema_HArray2OfPOnSurfParams(1, myusample - 1, 1, myvsample);
 630     Handle(Extrema_HArray2OfPOnSurfParams) aVEdgePntParams =
 631           new Extrema_HArray2OfPOnSurfParams(1, myusample, 1, myvsample - 1);
 632
 633     for ( NoU = 1 ; NoU <= myusample; NoU++ ) {
 634       for ( NoV = 1 ; NoV <= myvsample; NoV++) {
 635         const Extrema_POnSurfParams &aParam0 = myPoints->Value(NoU, NoV);
 636
 637         if (NoU < myusample) {
 638           // Compute parameters to UEdge.
 639           const Extrema_POnSurfParams &aParam1 = myPoints->Value(NoU + 1, NoV);
 640           Extrema_POnSurfParams aUEdgeParam = ComputeEdgeParameters
 641               (Standard_True, aParam0, aParam1, myS, thePoint, aDiffTol);
 642
 643           aUEdgePntParams->SetValue(NoU, NoV, aUEdgeParam);
 644         }
 645
 646         if (NoV < myvsample) {
 647           // Compute parameters to VEdge.
 648           const Extrema_POnSurfParams &aParam1 = myPoints->Value(NoU, NoV + 1);
 649           Extrema_POnSurfParams aVEdgeParam = ComputeEdgeParameters
 650               (Standard_False, aParam0, aParam1, myS, thePoint, aDiffTol);
 651
 652           aVEdgePntParams->SetValue(NoU, NoV, aVEdgeParam);
 653         }
 654       }
 655     }
 656
 657     // Step 3. Compute distances to faces.
 658     // Assume myFacePntParams(i, j) =
 659     // { Point(i, j); Point(i + 1, j); Point(i + 1, j + 1); Point(i, j + 1) }
 660     //   Or
 661     // { UEdge(i, j); VEdge(i + 1, j); UEdge(i, j + 1); VEdge(i, j) }
 662     myFacePntParams =
 663           new Extrema_HArray2OfPOnSurfParams(0, myusample, 0, myvsample);
 664     Standard_Real aSqrDist01;
 665     Standard_Real aDiffDist;
 666     Standard_Boolean isOut;
 667
 668     for ( NoU = 1 ; NoU < myusample; NoU++ ) {
 669       for ( NoV = 1 ; NoV < myvsample; NoV++) {
 670         const Extrema_POnSurfParams &aUE0 = aUEdgePntParams->Value(NoU, NoV);
 671         const Extrema_POnSurfParams &aUE1 = aUEdgePntParams->Value(NoU, NoV+1);
 672         const Extrema_POnSurfParams &aVE0 = aVEdgePntParams->Value(NoU, NoV);
 673         const Extrema_POnSurfParams &aVE1 = aVEdgePntParams->Value(NoU+1, NoV);
 674
 675         aSqrDist01 = aUE0.Value().SquareDistance(aUE1.Value());
 676         aDiffDist = Abs(aUE0.GetSqrDistance() - aUE1.GetSqrDistance());
 677         isOut = Standard_False;
 678
 679         if (aDiffDist >= aSqrDist01 - aDiffTol) {
 680           // The projection is outside the face.
 681           isOut = Standard_True;
 682         } else {
 683           aSqrDist01 = aVE0.Value().SquareDistance(aVE1.Value());
 684           aDiffDist = Abs(aVE0.GetSqrDistance() - aVE1.GetSqrDistance());
 685
 686           if (aDiffDist >= aSqrDist01 - aDiffTol) {
 687             // The projection is outside the face.
 688             isOut = Standard_True;
 689           }
 690         }
 691
 692         if (isOut) {
 693           // Get the closest point on an edge.
 694           const Extrema_POnSurfParams &aUEMin =
 695             aUE0.GetSqrDistance() < aUE1.GetSqrDistance() ? aUE0 : aUE1;
 696           const Extrema_POnSurfParams &aVEMin =
 697             aVE0.GetSqrDistance() < aVE1.GetSqrDistance() ? aVE0 : aVE1;
 698           const Extrema_POnSurfParams &aEMin =
 699             aUEMin.GetSqrDistance() < aVEMin.GetSqrDistance() ? aUEMin : aVEMin;
 700
 701           myFacePntParams->SetValue(NoU, NoV, aEMin);
 702         } else {
 703           // Find closest point inside the face.
 704           Standard_Real aU[2];
 705           Standard_Real aV[2];
 706           Standard_Real aUPar;
 707           Standard_Real aVPar;
 708
 709           // Compute U parameter.
 710           aUE0.Parameter(aU[0], aV[0]);
 711           aUE1.Parameter(aU[1], aV[1]);
 712           aUPar = 0.5*(aU[0] + aU[1]);
 713           // Compute V parameter.
 714           aVE0.Parameter(aU[0], aV[0]);
 715           aVE1.Parameter(aU[1], aV[1]);
 716           aVPar = 0.5*(aV[0] + aV[1]);
 717
 718           Extrema_POnSurfParams aParam(aUPar, aVPar, myS->Value(aUPar, aVPar));
 719
 720           aParam.SetElementType(Extrema_Face);
 721           aParam.SetSqrDistance(thePoint.SquareDistance(aParam.Value()));
 722           aParam.SetIndices(NoU, NoV);
 723           myFacePntParams->SetValue(NoU, NoV, aParam);
 724         }
 725       }
 726     }
 727
 728     // Fill boundary with RealLast square distance.
 729     for (NoV = 0; NoV <= myvsample; NoV++) {
 730       myFacePntParams->ChangeValue(0, NoV).SetSqrDistance(RealLast());
 731       myFacePntParams->ChangeValue(myusample, NoV).SetSqrDistance(RealLast());
 732     }
 733
 734     for (NoU = 1; NoU < myusample; NoU++) {
 735       myFacePntParams->ChangeValue(NoU, 0).SetSqrDistance(RealLast());
 736       myFacePntParams->ChangeValue(NoU, myvsample).SetSqrDistance(RealLast());
 737     }
 738   }
 739 }
 740
 741 /*
 742 a- Constitution of the table of distances (TbDist(0,myusample+1,0,myvsample+1)):
 743    ---------------------------------------------------------------
 744 */
 745
 746 // Parameterisation of the sample
 747
 748 void Extrema_GenExtPS::BuildTree()
 749 {
 750   // if tree already exists, assume it is already correctly filled
 751   if ( ! mySphereUBTree.IsNull() )
 752     return;
 753
 754   Standard_Real PasU = myusup - myumin;
 755   Standard_Real PasV = myvsup - myvmin;
 756   Standard_Real U0 = PasU / myusample / 100.;
 757   Standard_Real V0 = PasV / myvsample / 100.;
 758   gp_Pnt P1;
 759   PasU = (PasU - U0) / (myusample - 1);
 760   PasV = (PasV - V0) / (myvsample - 1);
 761   U0 = U0/2. + myumin;
 762   V0 = V0/2. + myvmin;
 763
 764   //build grid of parametric points
 765   myUParams = new TColStd_HArray1OfReal(1,myusample );
 766   myVParams = new TColStd_HArray1OfReal(1,myvsample );
 767   Standard_Integer NoU, NoV;
 768   Standard_Real U = U0, V = V0;
 769   for ( NoU = 1 ; NoU <= myusample; NoU++, U += PasU)
 770     myUParams->SetValue(NoU, U);
 771   for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV)
 772     myVParams->SetValue(NoV, V);
 773
 774   // Calculation of distances
 775   mySphereUBTree = new Extrema_UBTreeOfSphere;
 776   Extrema_UBTreeFillerOfSphere aFiller(*mySphereUBTree);
 777   Standard_Integer i = 0;
 778
 779   mySphereArray = new Bnd_HArray1OfSphere(0, myusample * myvsample);
 780
 781   for ( NoU = 1; NoU <= myusample; NoU++ ) {
 782     for ( NoV = 1; NoV <= myvsample; NoV++) {
 783       P1 = myS->Value(myUParams->Value(NoU), myVParams->Value(NoV));
 784       Bnd_Sphere aSph(P1.XYZ(), 0/*mytolu < mytolv ? mytolu : mytolv*/, NoU, NoV);
 785       aFiller.Add(i, aSph);
 786       mySphereArray->SetValue( i, aSph );
 787       i++;
 788     }
 789   }
 790   aFiller.Fill();
 791 }
 792
 793 void Extrema_GenExtPS::FindSolution(const gp_Pnt& /*P*/,
 794                                     const Extrema_POnSurfParams &theParams)
 795 {
 796   math_Vector Tol(1,2);
 797   Tol(1) = mytolu;
 798   Tol(2) = mytolv;
 799
 800   math_Vector UV(1, 2);
 801
 802   theParams.Parameter(UV(1), UV(2));
 803
 804   math_Vector UVinf(1,2), UVsup(1,2);
 805   UVinf(1) = myumin;
 806   UVinf(2) = myvmin;
 807   UVsup(1) = myusup;
 808   UVsup(2) = myvsup;
 809
 810   math_Vector errors(1,2);
 811   math_Vector root(1, 2);
 812
 813   Standard_Integer aNbMaxIter = 100;
 814
 815   gp_Pnt PStart = theParams.Value();
 816
 817   math_FunctionSetRoot S (myF,UV,Tol,UVinf,UVsup, aNbMaxIter);
 818
 819   myDone = Standard_True;
 820 }
 821
 822 void Extrema_GenExtPS::SetFlag(const Extrema_ExtFlag F)
 823 {
 824   myFlag = F;
 825 }
 826
 827 void Extrema_GenExtPS::SetAlgo(const Extrema_ExtAlgo A)
 828 {
 829   if(myAlgo != A)
 830      myInit = Standard_False;
 831   myAlgo = A;
 832
 833 }
 834
 835 void Extrema_GenExtPS::Perform(const gp_Pnt& P)
 836 {
 837   myDone = Standard_False;
 838   myF.SetPoint(P);
 839
 840   if(myAlgo == Extrema_ExtAlgo_Grad)
 841   {
 842     BuildGrid(P);
 843     Standard_Integer NoU,NoV;
 844
 845     if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX)
 846     {
 847       Extrema_ElementType anElemType;
 848       Standard_Integer iU;
 849       Standard_Integer iV;
 850       Standard_Integer iU2;
 851       Standard_Integer iV2;
 852       Standard_Boolean isMin;
 853       Standard_Integer i;
 854
 855       for (NoU = 1; NoU < myusample; NoU++) {
 856         for (NoV = 1; NoV < myvsample; NoV++) {
 857           const Extrema_POnSurfParams &aParam =
 858             myFacePntParams->Value(NoU, NoV);
 859
 860           isMin = Standard_False;
 861           anElemType = aParam.GetElementType();
 862
 863           if (anElemType == Extrema_Face) {
 864             isMin = Standard_True;
 865           } else {
 866             // Check if it is a boundary edge or corner vertex.
 867             aParam.GetIndices(iU, iV);
 868
 869             if (anElemType == Extrema_UIsoEdge) {
 870               isMin = (iV == 1 || iV == myvsample);
 871             } else if (anElemType == Extrema_VIsoEdge) {
 872               isMin = (iU == 1 || iU == myusample);
 873             } else if (anElemType == Extrema_Node) {
 874               isMin = (iU == 1 || iU == myusample) &&
 875                       (iV == 1 || iV == myvsample);
 876             }
 877
 878             if (!isMin) {
 879               // This is a middle element.
 880               if (anElemType == Extrema_UIsoEdge ||
 881                 (anElemType == Extrema_Node && (iU == 1 || iU == myusample))) {
 882                 // Check the down face.
 883                 const Extrema_POnSurfParams &aDownParam =
 884                     myFacePntParams->Value(NoU, NoV - 1);
 885
 886                 if (aDownParam.GetElementType() == anElemType) {
 887                   aDownParam.GetIndices(iU2, iV2);
 888                   isMin = (iU == iU2 && iV == iV2);
 889                 }
 890               } else if (anElemType == Extrema_VIsoEdge ||
 891                 (anElemType == Extrema_Node && (iV == 1 || iV == myvsample))) {
 892                 // Check the right face.
 893                 const Extrema_POnSurfParams &aRightParam =
 894                     myFacePntParams->Value(NoU - 1, NoV);
 895
 896                 if (aRightParam.GetElementType() == anElemType) {
 897                   aRightParam.GetIndices(iU2, iV2);
 898                   isMin = (iU == iU2 && iV == iV2);
 899                 }
 900               } else if (iU == NoU && iV == NoV) {
 901                 // Check the lower-left node. For this purpose it is necessary
 902                 // to check lower-left, lower and left faces.
 903                 isMin = Standard_True;
 904
 905                 const Extrema_POnSurfParams *anOtherParam[3] =
 906                 { &myFacePntParams->Value(NoU, NoV - 1),     // Down
 907                   &myFacePntParams->Value(NoU - 1, NoV - 1), // Lower-left
 908                   &myFacePntParams->Value(NoU - 1, NoV) };   // Left
 909
 910                 for (i = 0; i < 3 && isMin; i++) {
 911                   if (anOtherParam[i]->GetElementType() == Extrema_Node) {
 912                     anOtherParam[i]->GetIndices(iU2, iV2);
 913                     isMin = (iU == iU2 && iV == iV2);
 914                   } else {
 915                     isMin = Standard_False;
 916                   }
 917                 }
 918               }
 919             }
 920           }
 921
 922           if (isMin) {
 923             FindSolution(P, aParam);
 924           }
 925         }
 926       }
 927     }
 928
 929     if(myFlag == Extrema_ExtFlag_MAX || myFlag == Extrema_ExtFlag_MINMAX)
 930     {
 931       Standard_Real Dist;
 932
 933       for (NoU = 1; NoU <= myusample; NoU++) {
 934         for (NoV = 1; NoV <= myvsample; NoV++) {
 935           Dist = myPoints->Value(NoU, NoV).GetSqrDistance();
 936           if ((myPoints->Value(NoU-1,NoV-1).GetSqrDistance() <= Dist) &&
 937               (myPoints->Value(NoU-1,NoV  ).GetSqrDistance() <= Dist) &&
 938               (myPoints->Value(NoU-1,NoV+1).GetSqrDistance() <= Dist) &&
 939               (myPoints->Value(NoU  ,NoV-1).GetSqrDistance() <= Dist) &&
 940               (myPoints->Value(NoU  ,NoV+1).GetSqrDistance() <= Dist) &&
 941               (myPoints->Value(NoU+1,NoV-1).GetSqrDistance() <= Dist) &&
 942               (myPoints->Value(NoU+1,NoV  ).GetSqrDistance() <= Dist) &&
 943               (myPoints->Value(NoU+1,NoV+1).GetSqrDistance() <= Dist)) {
 944             // Find maximum.
 945             FindSolution(P, myPoints->Value(NoU, NoV));
 946           }
 947         }
 948       }
 949     }
 950   }
 951   else
 952   {
 953     BuildTree();
 954     if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX)
 955     {
 956       Bnd_Sphere aSol = mySphereArray->Value(0);
 957       Bnd_SphereUBTreeSelectorMin aSelector(mySphereArray, aSol);
 958       //aSelector.SetMaxDist( RealLast() );
 959       aSelector.DefineCheckPoint( P );
 960       mySphereUBTree->Select( aSelector );
 961       //TODO: check if no solution in binary tree
 962       Bnd_Sphere& aSph = aSelector.Sphere();
 963       Standard_Real aU = myUParams->Value(aSph.U());
 964       Standard_Real aV = myVParams->Value(aSph.V());
 965       Extrema_POnSurfParams aParams(aU, aV, myS->Value(aU, aV));
 966
 967       aParams.SetSqrDistance(P.SquareDistance(aParams.Value()));
 968       aParams.SetIndices(aSph.U(), aSph.V());
 969       FindSolution(P, aParams);
 970     }
 971     if(myFlag == Extrema_ExtFlag_MAX || myFlag == Extrema_ExtFlag_MINMAX)
 972     {
 973       Bnd_Sphere aSol = mySphereArray->Value(0);
 974       Bnd_SphereUBTreeSelectorMax aSelector(mySphereArray, aSol);
 975       //aSelector.SetMaxDist( RealLast() );
 976       aSelector.DefineCheckPoint( P );
 977       mySphereUBTree->Select( aSelector );
 978       //TODO: check if no solution in binary tree
 979       Bnd_Sphere& aSph = aSelector.Sphere();
 980       Standard_Real aU = myUParams->Value(aSph.U());
 981       Standard_Real aV = myVParams->Value(aSph.V());
 982       Extrema_POnSurfParams aParams(aU, aV, myS->Value(aU, aV));
 983       aParams.SetSqrDistance(P.SquareDistance(aParams.Value()));
 984       aParams.SetIndices(aSph.U(), aSph.V());
 985
 986       FindSolution(P, aParams);
 987     }
 988   }
 989 }
 990 //=============================================================================
 991
 992 Standard_Boolean Extrema_GenExtPS::IsDone () const { return myDone; }
 993 //=============================================================================
 994
 995 Standard_Integer Extrema_GenExtPS::NbExt () const
 996 {
 997   if (!IsDone()) { StdFail_NotDone::Raise(); }
 998   return myF.NbExt();
 999 }
1000 //=============================================================================
1001
1002 Standard_Real Extrema_GenExtPS::SquareDistance (const Standard_Integer N) const
1003 {
1004   if (!IsDone()) { StdFail_NotDone::Raise(); }
1005   return myF.SquareDistance(N);
1006 }
1007 //=============================================================================
1008
1009 Extrema_POnSurf Extrema_GenExtPS::Point (const Standard_Integer N) const
1010 {
1011   if (!IsDone()) { StdFail_NotDone::Raise(); }
1012   return myF.Point(N);
1013 }
1014 //=============================================================================