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