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