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