OCC22610 The algorithm GeomAPI_ProjectPointOnSurf produces wrong results
[occt.git] / src / Extrema / Extrema_GenExtPS.cxx
1 // File:        Extrema_GenExtPS.cxx
2 // Created:     Tue Jul 18 08:21:34 1995
3 // Author:      Modelistation
4 //              <model@metrox>
5
6 //  Modified by skv - Thu Sep 30 15:21:07 2004 OCC593
7
8
9 #include <Extrema_GenExtPS.ixx>
10 #include <StdFail_NotDone.hxx>
11 #include <Standard_OutOfRange.hxx>
12 #include <TColStd_Array2OfInteger.hxx>
13 #include <TColStd_Array2OfReal.hxx>
14 #include <TColgp_Array2OfPnt.hxx>
15 #include <math_FunctionSetRoot.hxx>
16 #include <math_Vector.hxx>
17 #include <math_NewtonFunctionSetRoot.hxx>
18 #include <GeomAbs_IsoType.hxx>
19 #include <Bnd_Sphere.hxx>
20 #include <Extrema_HUBTreeOfSphere.hxx>
21 #include <Extrema_ExtFlag.hxx>
22 #include <Bnd_Array1OfSphere.hxx>
23 #include <Bnd_HArray1OfSphere.hxx>
24
25 //IMPLEMENT_HARRAY1(Extrema_HArray1OfSphere)
26
27
28 class Bnd_SphereUBTreeSelector : public Extrema_UBTreeOfSphere::Selector
29 {
30  public:
31
32   Bnd_SphereUBTreeSelector (const Handle(Bnd_HArray1OfSphere)& theSphereArray,
33 Bnd_Sphere& theSol)
34     : myXYZ(0,0,0),
35       mySphereArray(theSphereArray),
36       mySol(theSol)
37   {
38     //myXYZ = gp_Pnt(0, 0, 0);    
39   }
40
41   void DefineCheckPoint( const gp_Pnt& theXYZ )
42   { myXYZ = theXYZ; }
43
44   Bnd_Sphere& Sphere() const
45   { return mySol; }
46
47   virtual Standard_Boolean Reject( const Bnd_Sphere &theBnd ) const = 0;
48   
49   virtual Standard_Boolean Accept(const Standard_Integer& theObj) = 0;
50
51  protected:
52   gp_Pnt                      myXYZ;
53   const Handle(Bnd_HArray1OfSphere)& mySphereArray;
54   Bnd_Sphere&                                            mySol;
55
56 };
57
58 class Bnd_SphereUBTreeSelectorMin : public Bnd_SphereUBTreeSelector
59 {
60 public:
61   Bnd_SphereUBTreeSelectorMin (const Handle(Bnd_HArray1OfSphere)& theSphereArray,
62                 Bnd_Sphere& theSol)
63                 : Bnd_SphereUBTreeSelector(theSphereArray, theSol),
64                   myMinDist(RealLast())
65   {}
66   
67   void SetMinDist( const Standard_Real theMinDist )
68   { myMinDist = theMinDist; }
69
70   Standard_Real MinDist() const
71   { return myMinDist; }
72
73   Standard_Boolean Reject( const Bnd_Sphere &theBnd ) const
74   { 
75     Bnd_SphereUBTreeSelectorMin* me =
76       const_cast<Bnd_SphereUBTreeSelectorMin*>(this);
77     // myMinDist is decreased each time a nearer object is found
78     return theBnd.IsOut( myXYZ.XYZ(), me->myMinDist );
79   }
80
81   Standard_Boolean Accept(const Standard_Integer&);
82
83 private:
84         Standard_Real   myMinDist;
85 };
86
87 Standard_Boolean Bnd_SphereUBTreeSelectorMin::Accept(const Standard_Integer& theInd)
88 {
89   const Bnd_Sphere& aSph = mySphereArray->Value(theInd);
90   Standard_Real aCurDist;
91
92     if ( (aCurDist = aSph.SquareDistance(myXYZ.XYZ())) < mySol.SquareDistance(myXYZ.XYZ()) )
93     {
94       mySol = aSph;
95       if ( aCurDist < myMinDist ) 
96         myMinDist = aCurDist;
97
98       return Standard_True;
99     }
100
101   return Standard_False;
102 }
103
104 class Bnd_SphereUBTreeSelectorMax : public Bnd_SphereUBTreeSelector
105 {
106 public:
107   Bnd_SphereUBTreeSelectorMax (const Handle(Bnd_HArray1OfSphere)& theSphereArray,
108                 Bnd_Sphere& theSol)
109                 : Bnd_SphereUBTreeSelector(theSphereArray, theSol),
110                   myMaxDist(0)
111   {}
112
113   void SetMaxDist( const Standard_Real theMaxDist )
114   { myMaxDist = theMaxDist; }
115
116   Standard_Real MaxDist() const
117   { return myMaxDist; }
118
119   Standard_Boolean Reject( const Bnd_Sphere &theBnd ) const
120   { 
121     Bnd_SphereUBTreeSelectorMax* me =
122       const_cast<Bnd_SphereUBTreeSelectorMax*>(this);
123     // myMaxDist is decreased each time a nearer object is found
124     return theBnd.IsOut( myXYZ.XYZ(), me->myMaxDist );
125   }
126
127   Standard_Boolean Accept(const Standard_Integer&);
128
129 private:
130         Standard_Real   myMaxDist;
131 };
132
133 Standard_Boolean Bnd_SphereUBTreeSelectorMax::Accept(const Standard_Integer& theInd)
134 {
135   const Bnd_Sphere& aSph = mySphereArray->Value(theInd);
136   Standard_Real aCurDist;
137
138     if ( (aCurDist = aSph.SquareDistance(myXYZ.XYZ())) > mySol.SquareDistance(myXYZ.XYZ()) )
139     {
140       mySol = aSph;
141       if ( aCurDist > myMaxDist ) 
142         myMaxDist = aCurDist;
143
144       return Standard_True;
145     }
146
147   return Standard_False;
148 }
149
150
151
152 //=============================================================================
153
154 /*-----------------------------------------------------------------------------
155 Function:
156    Find all extremum distances between point P and surface
157   S using sampling (NbU,NbV).
158
159 Method:
160    The algorithm bases on the hypothesis that sampling is precise enough 
161   pour que, s'il existe N distances extremales entre le point et la surface,
162   alors il existe aussi N extrema entre le point et la grille.
163   So, the algorithm consists in starting from extrema of the grid to find the 
164   extrema of the surface.
165   The extrema are calculated by the algorithm math_FunctionSetRoot with the
166   following arguments:
167   - F: Extrema_FuncExtPS created from P and S,
168   - UV: math_Vector the components which of are parameters of the extremum on the 
169     grid,
170   - Tol: Min(TolU,TolV), (Prov.:math_FunctionSetRoot does not autorize a vector)
171   - UVinf: math_Vector the components which of are lower limits of u and v,
172   - UVsup: math_Vector the components which of are upper limits of u and v.
173
174 Processing:
175   a- Creation of the table of distances (TbDist(0,NbU+1,0,NbV+1)):
176      The table is expanded at will; lines 0 and NbU+1 and
177      columns 0 and NbV+1 are initialized at RealFirst() or RealLast()
178      to simplify the tests carried out at stage b
179      (there is no need to test if the point is on border of the grid).
180   b- Calculation of extrema:
181      First the minimums and then the maximums are found. These 2 procedured 
182      pass in a similar way:
183   b.a- Initialization:
184       - 'borders' of table  TbDist (RealLast() in case of minimums
185         and  RealLast() in case of maximums),
186       - table TbSel(0,NbU+1,0,NbV+1) of selection of points for 
187         calculation of local extremum (0). When a point will selected,
188         it will not be selectable, as well as the ajacent points 
189         (8 at least). The corresponding addresses will be set to 1.
190   b.b- Calculation of minimums (or maximums):
191        All distances from table TbDist are parsed in a loop:
192       - search minimum (or maximum) in the grid,
193       - calculate extremum on the surface,
194       - update table TbSel.
195 -----------------------------------------------------------------------------*/
196
197 static Standard_Boolean IsoIsDeg  (const Adaptor3d_Surface& S,
198                                    const Standard_Real      Param,
199                                    const GeomAbs_IsoType    IT,
200                                    const Standard_Real      TolMin,
201                                    const Standard_Real      TolMax) 
202 {
203     Standard_Real U1=0.,U2=0.,V1=0.,V2=0.,T;
204     Standard_Boolean Along = Standard_True;
205     U1 = S.FirstUParameter();
206     U2 = S.LastUParameter();
207     V1 = S.FirstVParameter();
208     V2 = S.LastVParameter();
209     gp_Vec D1U,D1V;
210     gp_Pnt P;
211     Standard_Real Step,D1NormMax;
212     if (IT == GeomAbs_IsoV) 
213     {
214       Step = (U2 - U1)/10;
215       D1NormMax=0.;
216       for (T=U1;T<=U2;T=T+Step) 
217       {
218         S.D1(T,Param,P,D1U,D1V);
219         D1NormMax=Max(D1NormMax,D1U.Magnitude());
220       }
221
222       if (D1NormMax >TolMax || D1NormMax < TolMin ) 
223            Along = Standard_False;
224     }
225     else 
226     {
227       Step = (V2 - V1)/10;
228       D1NormMax=0.;
229       for (T=V1;T<=V2;T=T+Step) 
230       {
231         S.D1(Param,T,P,D1U,D1V);
232         D1NormMax=Max(D1NormMax,D1V.Magnitude());
233       }
234
235       if (D1NormMax >TolMax || D1NormMax < TolMin ) 
236            Along = Standard_False;
237
238
239     }
240     return Along;
241 }
242 //----------------------------------------------------------
243 Extrema_GenExtPS::Extrema_GenExtPS() 
244 {
245   myDone = Standard_False;
246   myInit = Standard_False;
247   myFlag = Extrema_ExtFlag_MINMAX;
248   myAlgo = Extrema_ExtAlgo_Grad;
249 }
250
251
252 Extrema_GenExtPS::Extrema_GenExtPS (const gp_Pnt&          P,
253                               const Adaptor3d_Surface& S,
254                               const Standard_Integer NbU, 
255                               const Standard_Integer NbV,
256                               const Standard_Real    TolU, 
257                               const Standard_Real    TolV,
258                               const Extrema_ExtFlag F,
259                               const Extrema_ExtAlgo A) 
260   : myF (P,S), myFlag(F), myAlgo(A)
261 {
262   Initialize(S, NbU, NbV, TolU, TolV);
263   Perform(P);
264 }
265
266 Extrema_GenExtPS::Extrema_GenExtPS (const gp_Pnt&          P,
267                               const Adaptor3d_Surface& S,
268                               const Standard_Integer NbU, 
269                               const Standard_Integer NbV,
270                               const Standard_Real    Umin,
271                               const Standard_Real    Usup,
272                               const Standard_Real    Vmin,
273                               const Standard_Real    Vsup,
274                               const Standard_Real    TolU, 
275                               const Standard_Real    TolV,
276                               const Extrema_ExtFlag F,
277                               const Extrema_ExtAlgo A) 
278   : myF (P,S), myFlag(F), myAlgo(A)
279 {
280   Initialize(S, NbU, NbV, Umin, Usup, Vmin, Vsup, TolU, TolV);
281   Perform(P);
282 }
283
284
285 void Extrema_GenExtPS::Initialize(const Adaptor3d_Surface& S,
286                                   const Standard_Integer NbU, 
287                                   const Standard_Integer NbV,
288                                   const Standard_Real    TolU, 
289                                   const Standard_Real    TolV)
290 {
291   myumin = S.FirstUParameter();
292   myusup = S.LastUParameter();
293   myvmin = S.FirstVParameter();
294   myvsup = S.LastVParameter();
295   Initialize(S,NbU,NbV,myumin,myusup,myvmin,myvsup,TolU,TolV);
296 }
297
298
299 void Extrema_GenExtPS::Initialize(const Adaptor3d_Surface& S,
300                                   const Standard_Integer NbU, 
301                                   const Standard_Integer NbV,
302                                   const Standard_Real    Umin,
303                                   const Standard_Real    Usup,
304                                   const Standard_Real    Vmin,
305                                   const Standard_Real    Vsup,
306                                   const Standard_Real    TolU, 
307                                   const Standard_Real    TolV)
308 {
309   myInit = Standard_True;
310   myS = (Adaptor3d_SurfacePtr)&S;
311   myusample = NbU;
312   myvsample = NbV;
313   mytolu = TolU;
314   mytolv = TolV;
315   myumin = Umin;
316   myusup = Usup;
317   myvmin = Vmin;
318   myvsup = Vsup;
319
320   if ((myusample < 2) ||
321       (myvsample < 2)) { Standard_OutOfRange::Raise(); }
322
323   myF.Initialize(S);
324
325   mySphereUBTree.Nullify();
326
327   if(myAlgo == Extrema_ExtAlgo_Grad)
328   {
329           //If flag was changed and extrema not reinitialized Extrema would fail
330     mypoints = new TColgp_HArray2OfPnt(0,myusample+1,0,myvsample+1);
331   Standard_Real PasU = myusup - myumin;
332   Standard_Real PasV = myvsup - myvmin;
333   Standard_Real U0 = PasU / myusample / 100.;
334   Standard_Real V0 = PasV / myvsample / 100.;
335   gp_Pnt P1;
336   PasU = (PasU - U0) / (myusample - 1);
337   PasV = (PasV - V0) / (myvsample - 1);
338   U0 = U0/2. + myumin;
339   V0 = V0/2. + myvmin;
340
341 // Calcul des distances
342
343   Standard_Integer NoU, NoV;
344   Standard_Real U, V;
345   for ( NoU = 1, U = U0; NoU <= myusample; NoU++, U += PasU) {
346     for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV) {
347       P1 = myS->Value(U, V);
348       mypoints->SetValue(NoU,NoV,P1);
349     }
350   }
351   }
352
353   //mypoints = new TColgp_HArray2OfPnt(0,myusample+1,0,myvsample+1);
354
355 /*
356 a- Constitution du tableau des distances (TbDist(0,myusample+1,0,myvsample+1)):
357    ---------------------------------------------------------------
358 */
359
360 // Parametrage de l echantillon
361
362
363 }
364
365 void Extrema_GenExtPS::BuildTree()
366 {
367   // if tree already exists, assume it is already correctly filled
368   if ( ! mySphereUBTree.IsNull() )
369     return;
370
371   Standard_Real PasU = myusup - myumin;
372   Standard_Real PasV = myvsup - myvmin;
373   Standard_Real U0 = PasU / myusample / 100.;
374   Standard_Real V0 = PasV / myvsample / 100.;
375   gp_Pnt P1;
376   PasU = (PasU - U0) / (myusample - 1);
377   PasV = (PasV - V0) / (myvsample - 1);
378   U0 = U0/2. + myumin;
379   V0 = V0/2. + myvmin;
380
381   // Calcul des distances
382   mySphereUBTree = new Extrema_UBTreeOfSphere;
383   Extrema_UBTreeFillerOfSphere aFiller(*mySphereUBTree);
384   Standard_Integer i = 0;
385   Standard_Real U, V;
386   mySphereArray = new Bnd_HArray1OfSphere(0, myusample * myvsample);
387   Standard_Integer NoU, NoV;
388   for ( NoU = 1, U = U0; NoU <= myusample; NoU++, U += PasU) {
389     for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV) {
390       P1 = myS->Value(U, V);
391       Bnd_Sphere aSph(P1.XYZ(), 0/*mytolu < mytolv ? mytolu : mytolv*/, NoU, NoV);
392       aFiller.Add(i, aSph);
393       mySphereArray->SetValue( i, aSph );
394       i++;
395     }
396   }
397   aFiller.Fill();
398 }
399
400 void Extrema_GenExtPS::FindSolution(const gp_Pnt& P, const math_Vector& UV, const Standard_Real PasU, const Standard_Real PasV, const Extrema_ExtFlag f)
401 {
402   math_Vector Tol(1,2);
403   Tol(1) = mytolu;
404   Tol(2) = mytolv;
405
406   math_Vector UVinf(1,2), UVsup(1,2);
407   UVinf(1) = myumin;
408   UVinf(2) = myvmin;
409   UVsup(1) = myusup;
410   UVsup(2) = myvsup;
411
412   math_Vector errors(1,2);
413   math_Vector root(1, 2);
414   Standard_Real eps = 1.e-9;
415   Standard_Integer nbsubsample = 11;
416
417   Standard_Integer aNbMaxIter = 100;
418
419   if (myF.HasDegIso())
420     aNbMaxIter = 150;
421
422   gp_Pnt PStart = myS->Value(UV(1), UV(2));
423   Standard_Real DistStart = P.SquareDistance(PStart);
424   Standard_Real DistSol = DistStart;
425   
426   math_FunctionSetRoot S (myF,UV,Tol,UVinf,UVsup, aNbMaxIter);
427   Standard_Boolean ToResolveOnSubgrid = Standard_False;
428   if (f == Extrema_ExtFlag_MIN)
429   {
430     if(S.IsDone())
431     {
432       root = S.Root();
433       myF.Value(root, errors);
434       gp_Pnt PSol = myS->Value(root(1), root(2));
435       DistSol = P.SquareDistance(PSol);
436       if(Abs(errors(1)) > eps || Abs(errors(2)) > eps || DistStart < DistSol)
437         //try to improve solution on subgrid of sample points
438         ToResolveOnSubgrid = Standard_True;
439     }
440     else
441       ToResolveOnSubgrid = Standard_True;
442     
443     if (ToResolveOnSubgrid)
444     {
445       Standard_Real u1 = Max(UV(1) - PasU, myumin), u2 = Min(UV(1) + PasU, myusup);
446       Standard_Real v1 = Max(UV(2) - PasV, myvmin), v2 = Min(UV(2) + PasV, myvsup);
447       
448       if(u2 - u1 < 2.*PasU) {
449         if(Abs(u1 - myumin) < 1.e-9) {
450           u2 = u1 + 2.*PasU;
451           u2 = Min(u2, myusup);
452         }
453         if(Abs(u2 - myusup) < 1.e-9) {
454           u1 = u2 - 2.*PasU;
455           u1 = Max(u1, myumin);
456         }
457       }
458       
459       if(v2 - v1 < 2.*PasV) {
460         if(Abs(v1 - myvmin) < 1.e-9) {
461           v2 = v1 + 2.*PasV;
462           v2 = Min(v2, myvsup);
463         }
464         if(Abs(v2 - myvsup) < 1.e-9) {
465           v1 = v2 - 2.*PasV;
466           v1 = Max(v1, myvmin);
467         }
468       }
469       
470       Standard_Real du = (u2 - u1)/(nbsubsample-1);
471       Standard_Real dv = (v2 - v1)/(nbsubsample-1);
472       Standard_Real u, v;
473       Standard_Real dist;
474       
475       Standard_Boolean NewSolution = Standard_False;
476       Standard_Integer Nu, Nv;
477       for (Nu = 1, u = u1; Nu < nbsubsample; Nu++, u += du) {
478         for (Nv = 1, v = v1; Nv < nbsubsample; Nv++, v += dv) {
479           gp_Pnt Puv = myS->Value(u, v);
480           dist = P.SquareDistance(Puv);
481           
482           if(dist < DistSol) {
483             UV(1) = u;
484             UV(2) = v;
485             NewSolution = Standard_True;
486             DistSol = dist;
487           }
488         }
489       }
490       
491       if(NewSolution) {
492         //try to precise
493         math_FunctionSetRoot S (myF,UV,Tol,UVinf,UVsup, aNbMaxIter);
494       }
495     } //end of if (ToResolveOnSubgrid)
496   } //end of if (f == Extrema_ExtFlag_MIN)
497   
498   myDone = Standard_True;
499 }
500
501 void Extrema_GenExtPS::SetFlag(const Extrema_ExtFlag F)
502 {
503   myFlag = F;
504 }
505
506 void Extrema_GenExtPS::SetAlgo(const Extrema_ExtAlgo A)
507 {
508   myAlgo = A;
509 }
510
511 void Extrema_GenExtPS::Perform(const gp_Pnt& P) 
512 {  
513   //BuildTree();
514   myDone = Standard_False;
515   myF.SetPoint(P);
516
517   Standard_Real PasU = myusup - myumin;
518   Standard_Real PasV = myvsup - myvmin;
519   Standard_Real U, U0 = PasU / myusample / 100.;
520   Standard_Real V, V0 = PasV / myvsample / 100.;
521   PasU = (PasU - U0) / (myusample - 1);
522   PasV = (PasV - V0) / (myvsample - 1);
523   U0 = U0/2.+myumin;
524   V0 = V0/2.+myvmin;
525
526   //math_Vector Tol(1, 2);
527   math_Vector UV(1,2);
528
529   if(myAlgo == Extrema_ExtAlgo_Grad)
530   {
531     Standard_Integer NoU,NoV;
532
533     TColStd_Array2OfReal TheDist(0, myusample+1, 0, myvsample+1);
534     for ( NoU = 1, U = U0; NoU <= myusample; NoU++, U += PasU) {
535       for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV) {
536         TheDist(NoU, NoV) = P.SquareDistance(mypoints->Value(NoU, NoV));
537       }
538     }
539
540     TColStd_Array2OfInteger TbSel(0,myusample+1,0,myvsample+1);
541     TbSel.Init(0);
542
543     Standard_Real Dist;
544
545     if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX) 
546     {
547       for (NoV = 0; NoV <= myvsample+1; NoV++) {
548         TheDist(0,NoV) = RealLast();
549         TheDist(myusample+1,NoV) = RealLast();
550       }
551       for (NoU = 1; NoU <= myusample; NoU++) {
552         TheDist(NoU,0) = RealLast();
553         TheDist(NoU,myvsample+1) = RealLast();
554       }
555       for (NoU = 1; NoU <= myusample; NoU++) {
556         for (NoV = 1; NoV <= myvsample; NoV++) {
557           if (TbSel(NoU,NoV) == 0) {
558             Dist = TheDist(NoU,NoV);
559             if ((TheDist(NoU-1,NoV-1) >= Dist) &&
560               (TheDist(NoU-1,NoV  ) >= Dist) &&
561               (TheDist(NoU-1,NoV+1) >= Dist) &&
562               (TheDist(NoU  ,NoV-1) >= Dist) &&
563               (TheDist(NoU  ,NoV+1) >= Dist) &&
564               (TheDist(NoU+1,NoV-1) >= Dist) &&
565               (TheDist(NoU+1,NoV  ) >= Dist) &&
566               (TheDist(NoU+1,NoV+1) >= Dist)) {
567                 //Create array of UV vectors to calculate min
568                 UV(1) = U0 + (NoU - 1) * PasU;
569                 UV(2) = V0 + (NoV - 1) * PasV;
570                 FindSolution(P, UV, PasU, PasV, Extrema_ExtFlag_MIN);
571             }
572           }
573         }
574       }
575     }
576     
577     if(myFlag == Extrema_ExtFlag_MAX || myFlag == Extrema_ExtFlag_MINMAX)
578     {
579       for (NoV = 0; NoV <= myvsample+1; NoV++) {
580         TheDist(0,NoV) = RealFirst();
581         TheDist(myusample+1,NoV) = RealFirst();
582       }
583       for (NoU = 1; NoU <= myusample; NoU++) {
584         TheDist(NoU,0) = RealFirst();
585         TheDist(NoU,myvsample+1) = RealFirst();
586       }
587       for (NoU = 1; NoU <= myusample; NoU++) {
588         for (NoV = 1; NoV <= myvsample; NoV++) {
589           if (TbSel(NoU,NoV) == 0) {
590             Dist = TheDist(NoU,NoV);
591             if ((TheDist(NoU-1,NoV-1) <= Dist) &&
592               (TheDist(NoU-1,NoV  ) <= Dist) &&
593               (TheDist(NoU-1,NoV+1) <= Dist) &&
594               (TheDist(NoU  ,NoV-1) <= Dist) &&
595               (TheDist(NoU  ,NoV+1) <= Dist) &&
596               (TheDist(NoU+1,NoV-1) <= Dist) &&
597               (TheDist(NoU+1,NoV  ) <= Dist) &&
598               (TheDist(NoU+1,NoV+1) <= Dist)) {
599                 //Create array of UV vectors to calculate max
600                 UV(1) = U0 + (NoU - 1) * PasU;
601                 UV(2) = V0 + (NoV - 1) * PasV;
602                 FindSolution(P, UV, PasU, PasV, Extrema_ExtFlag_MAX);
603             }
604           }
605         }
606       }
607     }
608   }
609   else
610   {
611     BuildTree();
612     if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX)
613     {
614       Bnd_Sphere aSol = mySphereArray->Value(0);
615       Bnd_SphereUBTreeSelectorMin aSelector(mySphereArray, aSol);
616       //aSelector.SetMaxDist( RealLast() );
617       aSelector.DefineCheckPoint( P );
618       Standard_Integer aNbSel = mySphereUBTree->Select( aSelector );
619       //TODO: check if no solution in binary tree
620       Bnd_Sphere& aSph = aSelector.Sphere();
621
622       UV(1) = U0 + (aSph.U() - 1) * PasU;
623       UV(2) = V0 + (aSph.V() - 1) * PasV;
624
625       FindSolution(P, UV, PasU, PasV, Extrema_ExtFlag_MIN);
626     }
627     if(myFlag == Extrema_ExtFlag_MAX || myFlag == Extrema_ExtFlag_MINMAX)
628     {
629       Bnd_Sphere aSol = mySphereArray->Value(0);
630       Bnd_SphereUBTreeSelectorMax aSelector(mySphereArray, aSol);
631       //aSelector.SetMaxDist( RealLast() );
632       aSelector.DefineCheckPoint( P );
633       Standard_Integer aNbSel = mySphereUBTree->Select( aSelector );
634       //TODO: check if no solution in binary tree
635       Bnd_Sphere& aSph = aSelector.Sphere();
636
637       UV(1) = U0 + (aSph.U() - 1) * PasU;
638       UV(2) = V0 + (aSph.V() - 1) * PasV;
639
640       FindSolution(P, UV, PasU, PasV, Extrema_ExtFlag_MAX);
641     }
642   }
643 }
644 //=============================================================================
645
646 Standard_Boolean Extrema_GenExtPS::IsDone () const { return myDone; }
647 //=============================================================================
648
649 Standard_Integer Extrema_GenExtPS::NbExt () const
650 {
651   if (!IsDone()) { StdFail_NotDone::Raise(); }
652   return myF.NbExt();
653 }
654 //=============================================================================
655
656 Standard_Real Extrema_GenExtPS::SquareDistance (const Standard_Integer N) const
657 {
658   if (!IsDone()) { StdFail_NotDone::Raise(); }
659   return myF.SquareDistance(N);
660 }
661 //=============================================================================
662
663 Extrema_POnSurf Extrema_GenExtPS::Point (const Standard_Integer N) const
664 {
665   if (!IsDone()) { StdFail_NotDone::Raise(); }
666   return myF.Point(N);
667 }
668 //=============================================================================