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 Extrema_GenExtPS::Extrema_GenExtPS()
 285 {
 286   myDone = Standard_False;
 287   myInit = Standard_False;
 288   myFlag = Extrema_ExtFlag_MINMAX;
 289   myAlgo = Extrema_ExtAlgo_Grad;
 290 }
 291
 292
 293 Extrema_GenExtPS::Extrema_GenExtPS (const gp_Pnt&          P,
 294                               const Adaptor3d_Surface& S,
 295                               const Standard_Integer NbU,
 296                               const Standard_Integer NbV,
 297                               const Standard_Real    TolU,
 298                               const Standard_Real    TolV,
 299                               const Extrema_ExtFlag F,
 300                               const Extrema_ExtAlgo A)
 301   : myF (P,S), myFlag(F), myAlgo(A)
 302 {
 303   Initialize(S, NbU, NbV, TolU, TolV);
 304   Perform(P);
 305 }
 306
 307 Extrema_GenExtPS::Extrema_GenExtPS (const gp_Pnt&          P,
 308                               const Adaptor3d_Surface& S,
 309                               const Standard_Integer NbU,
 310                               const Standard_Integer NbV,
 311                               const Standard_Real    Umin,
 312                               const Standard_Real    Usup,
 313                               const Standard_Real    Vmin,
 314                               const Standard_Real    Vsup,
 315                               const Standard_Real    TolU,
 316                               const Standard_Real    TolV,
 317                               const Extrema_ExtFlag F,
 318                               const Extrema_ExtAlgo A)
 319   : myF (P,S), myFlag(F), myAlgo(A)
 320 {
 321   Initialize(S, NbU, NbV, Umin, Usup, Vmin, Vsup, TolU, TolV);
 322   Perform(P);
 323 }
 324
 325
 326 void Extrema_GenExtPS::Initialize(const Adaptor3d_Surface& S,
 327                                   const Standard_Integer NbU,
 328                                   const Standard_Integer NbV,
 329                                   const Standard_Real    TolU,
 330                                   const Standard_Real    TolV)
 331 {
 332   myumin = S.FirstUParameter();
 333   myusup = S.LastUParameter();
 334   myvmin = S.FirstVParameter();
 335   myvsup = S.LastVParameter();
 336   Initialize(S,NbU,NbV,myumin,myusup,myvmin,myvsup,TolU,TolV);
 337 }
 338
 339
 340 void Extrema_GenExtPS::Initialize(const Adaptor3d_Surface& S,
 341                                   const Standard_Integer NbU,
 342                                   const Standard_Integer NbV,
 343                                   const Standard_Real    Umin,
 344                                   const Standard_Real    Usup,
 345                                   const Standard_Real    Vmin,
 346                                   const Standard_Real    Vsup,
 347                                   const Standard_Real    TolU,
 348                                   const Standard_Real    TolV)
 349 {
 350   myS = (Adaptor3d_SurfacePtr)&S;
 351   myusample = NbU;
 352   myvsample = NbV;
 353   mytolu = TolU;
 354   mytolv = TolV;
 355   myumin = Umin;
 356   myusup = Usup;
 357   myvmin = Vmin;
 358   myvsup = Vsup;
 359
 360   if ((myusample < 2) ||
 361       (myvsample < 2)) { Standard_OutOfRange::Raise(); }
 362
 363   myF.Initialize(S);
 364
 365   mySphereUBTree.Nullify();
 366   myUParams.Nullify();
 367   myVParams.Nullify();
 368   myInit = Standard_False;
 369 }
 370
 371 inline static void fillParams(const TColStd_Array1OfReal& theKnots,
 372                               Standard_Integer theDegree,
 373                               Standard_Real theParMin,
 374                               Standard_Real theParMax,
 375                               Handle_TColStd_HArray1OfReal& theParams,
 376                               Standard_Integer theSample)
 377 {
 378   NCollection_Vector<Standard_Real> aParams;
 379   Standard_Integer i = 1;
 380   Standard_Real aPrevPar = theParMin;
 381   aParams.Append(aPrevPar);
 382   //calculation the array of parametric points depending on the knots array variation and degree of given surface
 383   for ( ; i <  theKnots.Length() && theKnots(i) < (theParMax - Precision::PConfusion()); i++ )
 384   {
 385     if (  theKnots(i+1) < theParMin + Precision::PConfusion())
 386       continue;
 387
 388     Standard_Real aStep = (theKnots(i+1) - theKnots(i))/Max(theDegree,2);
 389     Standard_Integer k =1;
 390     for( ; k <= theDegree ; k++)
 391     {
 392       Standard_Real aPar = theKnots(i) + k * aStep;
 393       if(aPar > theParMax - Precision::PConfusion())
 394         break;
 395       if(aPar > aPrevPar + Precision::PConfusion() )
 396       {
 397         aParams.Append(aPar);
 398         aPrevPar = aPar;
 399       }
 400     }
 401
 402   }
 403   aParams.Append(theParMax);
 404   Standard_Integer nbPar = aParams.Length();
 405   //in case of an insufficient number of points the grid will be built later
 406   if (nbPar < theSample)
 407     return;
 408   theParams = new TColStd_HArray1OfReal(1, nbPar );
 409   for( i = 0; i < nbPar; i++)
 410     theParams->SetValue(i+1,aParams(i));
 411 }
 412
 413 void Extrema_GenExtPS::GetGridPoints( const Adaptor3d_Surface& theSurf)
 414 {
 415   //creation parametric points for BSpline and Bezier surfaces
 416   //with taking into account of Degree and NbKnots of BSpline or Bezier geometry
 417   if(theSurf.GetType() == GeomAbs_OffsetSurface)
 418     GetGridPoints(theSurf.BasisSurface()->Surface());
 419   //parametric points for BSpline surfaces
 420   else if( theSurf.GetType() == GeomAbs_BSplineSurface)
 421   {
 422     Handle(Geom_BSplineSurface) aBspl = theSurf.BSpline();
 423     if(!aBspl.IsNull())
 424     {
 425       TColStd_Array1OfReal aUKnots(1, aBspl->NbUKnots());
 426       aBspl->UKnots( aUKnots);
 427       TColStd_Array1OfReal aVKnots(1, aBspl->NbVKnots());
 428       aBspl->VKnots( aVKnots);
 429       fillParams(aUKnots,aBspl->UDegree(),myumin, myusup, myUParams, myusample);
 430       fillParams(aVKnots,aBspl->VDegree(),myvmin, myvsup, myVParams, myvsample);
 431     }
 432   }
 433   //calculation parametric points for Bezier surfaces
 434   else if(theSurf.GetType() == GeomAbs_BezierSurface)
 435   {
 436     Handle(Geom_BezierSurface) aBezier = theSurf.Bezier();
 437     if(aBezier.IsNull())
 438       return;
 439
 440     TColStd_Array1OfReal aUKnots(1,2);
 441     TColStd_Array1OfReal aVKnots(1,2);
 442     aBezier->Bounds(aUKnots(1), aUKnots(2), aVKnots(1), aVKnots(2));
 443     fillParams(aUKnots, aBezier->UDegree(), myumin, myusup, myUParams, myusample);
 444     fillParams(aVKnots, aBezier->VDegree(), myvmin, myvsup, myVParams, myvsample);
 445
 446   }
 447   //creation points for surfaces based on BSpline or Bezier curves
 448   else if(theSurf.GetType() == GeomAbs_SurfaceOfRevolution ||
 449     theSurf.GetType() == GeomAbs_SurfaceOfExtrusion)
 450   {
 451     Handle(TColStd_HArray1OfReal) anArrKnots;
 452     Standard_Integer aDegree = 0;
 453     if(theSurf.BasisCurve()->Curve().GetType() == GeomAbs_BSplineCurve)
 454     {
 455       Handle(Geom_BSplineCurve) aBspl = theSurf.BasisCurve()->Curve().BSpline();
 456       if(!aBspl.IsNull())
 457       {
 458         anArrKnots = new TColStd_HArray1OfReal(1,aBspl->NbKnots());
 459         aBspl->Knots( anArrKnots->ChangeArray1() );
 460         aDegree = aBspl->Degree();
 461
 462       }
 463
 464     }
 465     if(theSurf.BasisCurve()->Curve().GetType() == GeomAbs_BezierCurve)
 466     {
 467       Handle(Geom_BezierCurve) aBez = theSurf.BasisCurve()->Curve().Bezier();
 468       if(!aBez.IsNull())
 469       {
 470         anArrKnots = new TColStd_HArray1OfReal(1,2);
 471         anArrKnots->SetValue(1, aBez->FirstParameter());
 472         anArrKnots->SetValue(2, aBez->LastParameter());
 473         aDegree = aBez->Degree();
 474
 475       }
 476     }
 477     if(anArrKnots.IsNull())
 478       return;
 479     if( theSurf.GetType() == GeomAbs_SurfaceOfRevolution )
 480       fillParams( anArrKnots->Array1(), aDegree, myvmin, myvsup, myVParams, myvsample);
 481     else
 482       fillParams( anArrKnots->Array1(), aDegree, myumin, myusup, myUParams, myusample);
 483   }
 484   //update the number of points in sample
 485   if(!myUParams.IsNull())
 486     myusample = myUParams->Length();
 487   if( !myVParams.IsNull())
 488     myvsample = myVParams->Length();
 489
 490 }
 491
 492 void Extrema_GenExtPS::BuildGrid(const gp_Pnt &thePoint)
 493 {
 494   Standard_Integer NoU, NoV;
 495
 496   //if grid was already built skip its creation
 497   if (!myInit) {
 498     //build parametric grid in case of a complex surface geometry (BSpline and Bezier surfaces)
 499     GetGridPoints(*myS);
 500
 501     //build grid in other cases
 502     if( myUParams.IsNull() )
 503     {
 504       Standard_Real PasU = myusup - myumin;
 505       Standard_Real U0 = PasU / myusample / 100.;
 506       PasU = (PasU - U0) / (myusample - 1);
 507       U0 = U0/2. + myumin;
 508       myUParams = new TColStd_HArray1OfReal(1,myusample );
 509       Standard_Integer NoU;
 510       Standard_Real U = U0;
 511       for ( NoU = 1 ; NoU <= myusample; NoU++, U += PasU)
 512         myUParams->SetValue(NoU, U);
 513     }
 514
 515     if( myVParams.IsNull())
 516     {
 517       Standard_Real PasV = myvsup - myvmin;
 518       Standard_Real V0 = PasV / myvsample / 100.;
 519       PasV = (PasV - V0) / (myvsample - 1);
 520       V0 = V0/2. + myvmin;
 521
 522       myVParams = new TColStd_HArray1OfReal(1,myvsample );
 523       Standard_Integer  NoV;
 524       Standard_Real V = V0;
 525
 526       for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV)
 527         myVParams->SetValue(NoV, V);
 528     }
 529
 530     //If flag was changed and extrema not reinitialized Extrema would fail
 531     myPoints = new Extrema_HArray2OfPOnSurfParams
 532       (0, myusample + 1, 0, myvsample + 1);
 533     // Calculation of distances
 534
 535     for ( NoU = 1 ; NoU <= myusample; NoU++ ) {
 536       for ( NoV = 1 ; NoV <= myvsample; NoV++) {
 537         gp_Pnt aP1 = myS->Value(myUParams->Value(NoU), myVParams->Value(NoV));
 538         Extrema_POnSurfParams aParam
 539           (myUParams->Value(NoU), myVParams->Value(NoV), aP1);
 540
 541         aParam.SetElementType(Extrema_Node);
 542         aParam.SetIndices(NoU, NoV);
 543         myPoints->SetValue(NoU, NoV, aParam);
 544       }
 545     }
 546
 547     // Fill boundary with negative square distance.
 548     // It is used for computation of Maximum.
 549     for (NoV = 0; NoV <= myvsample + 1; NoV++) {
 550       myPoints->ChangeValue(0, NoV).SetSqrDistance(-1.);
 551       myPoints->ChangeValue(myusample + 1, NoV).SetSqrDistance(-1.);
 552     }
 553
 554     for (NoU = 1; NoU <= myusample; NoU++) {
 555       myPoints->ChangeValue(NoU, 0).SetSqrDistance(-1.);
 556       myPoints->ChangeValue(NoU, myvsample + 1).SetSqrDistance(-1.);
 557     }
 558
 559     myInit = Standard_True;
 560   }
 561
 562   // Compute distances to mesh.
 563   // Step 1. Compute distances to nodes.
 564   for ( NoU = 1 ; NoU <= myusample; NoU++ ) {
 565     for ( NoV = 1 ; NoV <= myvsample; NoV++) {
 566       Extrema_POnSurfParams &aParam = myPoints->ChangeValue(NoU, NoV);
 567
 568       aParam.SetSqrDistance(thePoint.SquareDistance(aParam.Value()));
 569     }
 570   }
 571
 572   // For search of minimum compute distances to mesh.
 573   if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX)
 574   {
 575     // This is the tolerance of difference of squared values.
 576     // No need to set it too small.
 577     const Standard_Real aDiffTol = mytolu + mytolv;
 578
 579     // Step 2. Compute distances to edges.
 580     // Assume UEdge(i, j) = { Point(i, j); Point(i + 1, j    ) }
 581     // Assume VEdge(i, j) = { Point(i, j); Point(i,     j + 1) }
 582     Handle(Extrema_HArray2OfPOnSurfParams) aUEdgePntParams =
 583           new Extrema_HArray2OfPOnSurfParams(1, myusample - 1, 1, myvsample);
 584     Handle(Extrema_HArray2OfPOnSurfParams) aVEdgePntParams =
 585           new Extrema_HArray2OfPOnSurfParams(1, myusample, 1, myvsample - 1);
 586
 587     for ( NoU = 1 ; NoU <= myusample; NoU++ ) {
 588       for ( NoV = 1 ; NoV <= myvsample; NoV++) {
 589         const Extrema_POnSurfParams &aParam0 = myPoints->Value(NoU, NoV);
 590
 591         if (NoU < myusample) {
 592           // Compute parameters to UEdge.
 593           const Extrema_POnSurfParams &aParam1 = myPoints->Value(NoU + 1, NoV);
 594           Extrema_POnSurfParams aUEdgeParam = ComputeEdgeParameters
 595               (Standard_True, aParam0, aParam1, myS, thePoint, aDiffTol);
 596
 597           aUEdgePntParams->SetValue(NoU, NoV, aUEdgeParam);
 598         }
 599
 600         if (NoV < myvsample) {
 601           // Compute parameters to VEdge.
 602           const Extrema_POnSurfParams &aParam1 = myPoints->Value(NoU, NoV + 1);
 603           Extrema_POnSurfParams aVEdgeParam = ComputeEdgeParameters
 604               (Standard_False, aParam0, aParam1, myS, thePoint, aDiffTol);
 605
 606           aVEdgePntParams->SetValue(NoU, NoV, aVEdgeParam);
 607         }
 608       }
 609     }
 610
 611     // Step 3. Compute distances to faces.
 612     // Assume myFacePntParams(i, j) =
 613     // { Point(i, j); Point(i + 1, j); Point(i + 1, j + 1); Point(i, j + 1) }
 614     //   Or
 615     // { UEdge(i, j); VEdge(i + 1, j); UEdge(i, j + 1); VEdge(i, j) }
 616     myFacePntParams =
 617           new Extrema_HArray2OfPOnSurfParams(0, myusample, 0, myvsample);
 618     Standard_Real aSqrDist01;
 619     Standard_Real aDiffDist;
 620     Standard_Boolean isOut;
 621
 622     for ( NoU = 1 ; NoU < myusample; NoU++ ) {
 623       for ( NoV = 1 ; NoV < myvsample; NoV++) {
 624         const Extrema_POnSurfParams &aUE0 = aUEdgePntParams->Value(NoU, NoV);
 625         const Extrema_POnSurfParams &aUE1 = aUEdgePntParams->Value(NoU, NoV+1);
 626         const Extrema_POnSurfParams &aVE0 = aVEdgePntParams->Value(NoU, NoV);
 627         const Extrema_POnSurfParams &aVE1 = aVEdgePntParams->Value(NoU+1, NoV);
 628
 629         aSqrDist01 = aUE0.Value().SquareDistance(aUE1.Value());
 630         aDiffDist = Abs(aUE0.GetSqrDistance() - aUE1.GetSqrDistance());
 631         isOut = Standard_False;
 632
 633         if (aDiffDist >= aSqrDist01 - aDiffTol) {
 634           // The projection is outside the face.
 635           isOut = Standard_True;
 636         } else {
 637           aSqrDist01 = aVE0.Value().SquareDistance(aVE1.Value());
 638           aDiffDist = Abs(aVE0.GetSqrDistance() - aVE1.GetSqrDistance());
 639
 640           if (aDiffDist >= aSqrDist01 - aDiffTol) {
 641             // The projection is outside the face.
 642             isOut = Standard_True;
 643           }
 644         }
 645
 646         if (isOut) {
 647           // Get the closest point on an edge.
 648           const Extrema_POnSurfParams &aUEMin =
 649             aUE0.GetSqrDistance() < aUE1.GetSqrDistance() ? aUE0 : aUE1;
 650           const Extrema_POnSurfParams &aVEMin =
 651             aVE0.GetSqrDistance() < aVE1.GetSqrDistance() ? aVE0 : aVE1;
 652           const Extrema_POnSurfParams &aEMin =
 653             aUEMin.GetSqrDistance() < aVEMin.GetSqrDistance() ? aUEMin : aVEMin;
 654
 655           myFacePntParams->SetValue(NoU, NoV, aEMin);
 656         } else {
 657           // Find closest point inside the face.
 658           Standard_Real aU[2];
 659           Standard_Real aV[2];
 660           Standard_Real aUPar;
 661           Standard_Real aVPar;
 662
 663           // Compute U parameter.
 664           aUE0.Parameter(aU[0], aV[0]);
 665           aUE1.Parameter(aU[1], aV[1]);
 666           aUPar = 0.5*(aU[0] + aU[1]);
 667           // Compute V parameter.
 668           aVE0.Parameter(aU[0], aV[0]);
 669           aVE1.Parameter(aU[1], aV[1]);
 670           aVPar = 0.5*(aV[0] + aV[1]);
 671
 672           Extrema_POnSurfParams aParam(aUPar, aVPar, myS->Value(aUPar, aVPar));
 673
 674           aParam.SetElementType(Extrema_Face);
 675           aParam.SetSqrDistance(thePoint.SquareDistance(aParam.Value()));
 676           aParam.SetIndices(NoU, NoV);
 677           myFacePntParams->SetValue(NoU, NoV, aParam);
 678         }
 679       }
 680     }
 681
 682     // Fill boundary with RealLast square distance.
 683     for (NoV = 0; NoV <= myvsample; NoV++) {
 684       myFacePntParams->ChangeValue(0, NoV).SetSqrDistance(RealLast());
 685       myFacePntParams->ChangeValue(myusample, NoV).SetSqrDistance(RealLast());
 686     }
 687
 688     for (NoU = 1; NoU < myusample; NoU++) {
 689       myFacePntParams->ChangeValue(NoU, 0).SetSqrDistance(RealLast());
 690       myFacePntParams->ChangeValue(NoU, myvsample).SetSqrDistance(RealLast());
 691     }
 692   }
 693 }
 694
 695 /*
 696 a- Constitution of the table of distances (TbDist(0,myusample+1,0,myvsample+1)):
 697    ---------------------------------------------------------------
 698 */
 699
 700 // Parameterisation of the sample
 701
 702 void Extrema_GenExtPS::BuildTree()
 703 {
 704   // if tree already exists, assume it is already correctly filled
 705   if ( ! mySphereUBTree.IsNull() )
 706     return;
 707
 708   Standard_Real PasU = myusup - myumin;
 709   Standard_Real PasV = myvsup - myvmin;
 710   Standard_Real U0 = PasU / myusample / 100.;
 711   Standard_Real V0 = PasV / myvsample / 100.;
 712   gp_Pnt P1;
 713   PasU = (PasU - U0) / (myusample - 1);
 714   PasV = (PasV - V0) / (myvsample - 1);
 715   U0 = U0/2. + myumin;
 716   V0 = V0/2. + myvmin;
 717
 718   //build grid of parametric points
 719   myUParams = new TColStd_HArray1OfReal(1,myusample );
 720   myVParams = new TColStd_HArray1OfReal(1,myvsample );
 721   Standard_Integer NoU, NoV;
 722   Standard_Real U = U0, V = V0;
 723   for ( NoU = 1 ; NoU <= myusample; NoU++, U += PasU)
 724     myUParams->SetValue(NoU, U);
 725   for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV)
 726     myVParams->SetValue(NoV, V);
 727
 728   // Calculation of distances
 729   mySphereUBTree = new Extrema_UBTreeOfSphere;
 730   Extrema_UBTreeFillerOfSphere aFiller(*mySphereUBTree);
 731   Standard_Integer i = 0;
 732
 733   mySphereArray = new Bnd_HArray1OfSphere(0, myusample * myvsample);
 734
 735   for ( NoU = 1; NoU <= myusample; NoU++ ) {
 736     for ( NoV = 1; NoV <= myvsample; NoV++) {
 737       P1 = myS->Value(myUParams->Value(NoU), myVParams->Value(NoV));
 738       Bnd_Sphere aSph(P1.XYZ(), 0/*mytolu < mytolv ? mytolu : mytolv*/, NoU, NoV);
 739       aFiller.Add(i, aSph);
 740       mySphereArray->SetValue( i, aSph );
 741       i++;
 742     }
 743   }
 744   aFiller.Fill();
 745 }
 746
 747 void Extrema_GenExtPS::FindSolution(const gp_Pnt& /*P*/,
 748                                     const Extrema_POnSurfParams &theParams)
 749 {
 750   math_Vector Tol(1,2);
 751   Tol(1) = mytolu;
 752   Tol(2) = mytolv;
 753
 754   math_Vector UV(1, 2);
 755
 756   theParams.Parameter(UV(1), UV(2));
 757
 758   math_Vector UVinf(1,2), UVsup(1,2);
 759   UVinf(1) = myumin;
 760   UVinf(2) = myvmin;
 761   UVsup(1) = myusup;
 762   UVsup(2) = myvsup;
 763
 764   math_Vector errors(1,2);
 765   math_Vector root(1, 2);
 766
 767   Standard_Integer aNbMaxIter = 100;
 768
 769   gp_Pnt PStart = theParams.Value();
 770
 771   math_FunctionSetRoot S (myF,UV,Tol,UVinf,UVsup, aNbMaxIter);
 772
 773   myDone = Standard_True;
 774 }
 775
 776 void Extrema_GenExtPS::SetFlag(const Extrema_ExtFlag F)
 777 {
 778   myFlag = F;
 779 }
 780
 781 void Extrema_GenExtPS::SetAlgo(const Extrema_ExtAlgo A)
 782 {
 783   if(myAlgo != A)
 784      myInit = Standard_False;
 785   myAlgo = A;
 786
 787 }
 788
 789 void Extrema_GenExtPS::Perform(const gp_Pnt& P)
 790 {
 791   myDone = Standard_False;
 792   myF.SetPoint(P);
 793
 794   if(myAlgo == Extrema_ExtAlgo_Grad)
 795   {
 796     BuildGrid(P);
 797     Standard_Integer NoU,NoV;
 798
 799     if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX)
 800     {
 801       Extrema_ElementType anElemType;
 802       Standard_Integer iU;
 803       Standard_Integer iV;
 804       Standard_Integer iU2;
 805       Standard_Integer iV2;
 806       Standard_Boolean isMin;
 807       Standard_Integer i;
 808
 809       for (NoU = 1; NoU < myusample; NoU++) {
 810         for (NoV = 1; NoV < myvsample; NoV++) {
 811           const Extrema_POnSurfParams &aParam =
 812             myFacePntParams->Value(NoU, NoV);
 813
 814           isMin = Standard_False;
 815           anElemType = aParam.GetElementType();
 816
 817           if (anElemType == Extrema_Face) {
 818             isMin = Standard_True;
 819           } else {
 820             // Check if it is a boundary edge or corner vertex.
 821             aParam.GetIndices(iU, iV);
 822
 823             if (anElemType == Extrema_UIsoEdge) {
 824               isMin = (iV == 1 || iV == myvsample);
 825             } else if (anElemType == Extrema_VIsoEdge) {
 826               isMin = (iU == 1 || iU == myusample);
 827             } else if (anElemType == Extrema_Node) {
 828               isMin = (iU == 1 || iU == myusample) &&
 829                       (iV == 1 || iV == myvsample);
 830             }
 831
 832             if (!isMin) {
 833               // This is a middle element.
 834               if (anElemType == Extrema_UIsoEdge ||
 835                 (anElemType == Extrema_Node && (iU == 1 || iU == myusample))) {
 836                 // Check the down face.
 837                 const Extrema_POnSurfParams &aDownParam =
 838                     myFacePntParams->Value(NoU, NoV - 1);
 839
 840                 if (aDownParam.GetElementType() == anElemType) {
 841                   aDownParam.GetIndices(iU2, iV2);
 842                   isMin = (iU == iU2 && iV == iV2);
 843                 }
 844               } else if (anElemType == Extrema_VIsoEdge ||
 845                 (anElemType == Extrema_Node && (iV == 1 || iV == myvsample))) {
 846                 // Check the right face.
 847                 const Extrema_POnSurfParams &aRightParam =
 848                     myFacePntParams->Value(NoU - 1, NoV);
 849
 850                 if (aRightParam.GetElementType() == anElemType) {
 851                   aRightParam.GetIndices(iU2, iV2);
 852                   isMin = (iU == iU2 && iV == iV2);
 853                 }
 854               } else if (iU == NoU && iV == NoV) {
 855                 // Check the lower-left node. For this purpose it is necessary
 856                 // to check lower-left, lower and left faces.
 857                 isMin = Standard_True;
 858
 859                 const Extrema_POnSurfParams *anOtherParam[3] =
 860                 { &myFacePntParams->Value(NoU, NoV - 1),     // Down
 861                   &myFacePntParams->Value(NoU - 1, NoV - 1), // Lower-left
 862                   &myFacePntParams->Value(NoU - 1, NoV) };   // Left
 863
 864                 for (i = 0; i < 3 && isMin; i++) {
 865                   if (anOtherParam[i]->GetElementType() == Extrema_Node) {
 866                     anOtherParam[i]->GetIndices(iU2, iV2);
 867                     isMin = (iU == iU2 && iV == iV2);
 868                   } else {
 869                     isMin = Standard_False;
 870                   }
 871                 }
 872               }
 873             }
 874           }
 875
 876           if (isMin) {
 877             FindSolution(P, aParam);
 878           }
 879         }
 880       }
 881     }
 882
 883     if(myFlag == Extrema_ExtFlag_MAX || myFlag == Extrema_ExtFlag_MINMAX)
 884     {
 885       Standard_Real Dist;
 886
 887       for (NoU = 1; NoU <= myusample; NoU++) {
 888         for (NoV = 1; NoV <= myvsample; NoV++) {
 889           Dist = myPoints->Value(NoU, NoV).GetSqrDistance();
 890           if ((myPoints->Value(NoU-1,NoV-1).GetSqrDistance() <= Dist) &&
 891               (myPoints->Value(NoU-1,NoV  ).GetSqrDistance() <= Dist) &&
 892               (myPoints->Value(NoU-1,NoV+1).GetSqrDistance() <= Dist) &&
 893               (myPoints->Value(NoU  ,NoV-1).GetSqrDistance() <= Dist) &&
 894               (myPoints->Value(NoU  ,NoV+1).GetSqrDistance() <= Dist) &&
 895               (myPoints->Value(NoU+1,NoV-1).GetSqrDistance() <= Dist) &&
 896               (myPoints->Value(NoU+1,NoV  ).GetSqrDistance() <= Dist) &&
 897               (myPoints->Value(NoU+1,NoV+1).GetSqrDistance() <= Dist)) {
 898             // Find maximum.
 899             FindSolution(P, myPoints->Value(NoU, NoV));
 900           }
 901         }
 902       }
 903     }
 904   }
 905   else
 906   {
 907     BuildTree();
 908     if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX)
 909     {
 910       Bnd_Sphere aSol = mySphereArray->Value(0);
 911       Bnd_SphereUBTreeSelectorMin aSelector(mySphereArray, aSol);
 912       //aSelector.SetMaxDist( RealLast() );
 913       aSelector.DefineCheckPoint( P );
 914       mySphereUBTree->Select( aSelector );
 915       //TODO: check if no solution in binary tree
 916       Bnd_Sphere& aSph = aSelector.Sphere();
 917       Standard_Real aU = myUParams->Value(aSph.U());
 918       Standard_Real aV = myVParams->Value(aSph.V());
 919       Extrema_POnSurfParams aParams(aU, aV, myS->Value(aU, aV));
 920
 921       aParams.SetSqrDistance(P.SquareDistance(aParams.Value()));
 922       aParams.SetIndices(aSph.U(), aSph.V());
 923       FindSolution(P, aParams);
 924     }
 925     if(myFlag == Extrema_ExtFlag_MAX || myFlag == Extrema_ExtFlag_MINMAX)
 926     {
 927       Bnd_Sphere aSol = mySphereArray->Value(0);
 928       Bnd_SphereUBTreeSelectorMax aSelector(mySphereArray, aSol);
 929       //aSelector.SetMaxDist( RealLast() );
 930       aSelector.DefineCheckPoint( P );
 931       mySphereUBTree->Select( aSelector );
 932       //TODO: check if no solution in binary tree
 933       Bnd_Sphere& aSph = aSelector.Sphere();
 934       Standard_Real aU = myUParams->Value(aSph.U());
 935       Standard_Real aV = myVParams->Value(aSph.V());
 936       Extrema_POnSurfParams aParams(aU, aV, myS->Value(aU, aV));
 937       aParams.SetSqrDistance(P.SquareDistance(aParams.Value()));
 938       aParams.SetIndices(aSph.U(), aSph.V());
 939
 940       FindSolution(P, aParams);
 941     }
 942   }
 943 }
 944 //=============================================================================
 945
 946 Standard_Boolean Extrema_GenExtPS::IsDone () const { return myDone; }
 947 //=============================================================================
 948
 949 Standard_Integer Extrema_GenExtPS::NbExt () const
 950 {
 951   if (!IsDone()) { StdFail_NotDone::Raise(); }
 952   return myF.NbExt();
 953 }
 954 //=============================================================================
 955
 956 Standard_Real Extrema_GenExtPS::SquareDistance (const Standard_Integer N) const
 957 {
 958   if (!IsDone()) { StdFail_NotDone::Raise(); }
 959   return myF.SquareDistance(N);
 960 }
 961 //=============================================================================
 962
 963 Extrema_POnSurf Extrema_GenExtPS::Point (const Standard_Integer N) const
 964 {
 965   if (!IsDone()) { StdFail_NotDone::Raise(); }
 966   return myF.Point(N);
 967 }
 968 //=============================================================================