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