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