0031875: Draw Harness, ViewerTest - command vaspects -mostContinuity lacks g1 and...
[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 : myumin(0.0),
222   myusup(0.0),
223   myvmin(0.0),
224   myvsup(0.0),
225   myusample(0),
226   myvsample(0),
227   mytolu(0.0),
228   mytolv(0.0),
229   myS(NULL)
230 {
231   myDone = Standard_False;
232   myInit = Standard_False;
233   myFlag = Extrema_ExtFlag_MINMAX;
234   myAlgo = Extrema_ExtAlgo_Grad;
235 }
236
237
238 Extrema_GenExtPS::Extrema_GenExtPS (const gp_Pnt&          P,
239                                     const Adaptor3d_Surface& S,
240                                     const Standard_Integer NbU, 
241                                     const Standard_Integer NbV,
242                                     const Standard_Real    TolU, 
243                                     const Standard_Real    TolV,
244                                     const Extrema_ExtFlag F,
245                                     const Extrema_ExtAlgo A) 
246 : myF (P,S),
247   myFlag(F),
248   myAlgo(A)
249 {
250   Initialize(S, NbU, NbV, TolU, TolV);
251   Perform(P);
252 }
253
254 Extrema_GenExtPS::Extrema_GenExtPS (const gp_Pnt&          P,
255                                     const Adaptor3d_Surface& S,
256                                     const Standard_Integer NbU, 
257                                     const Standard_Integer NbV,
258                                     const Standard_Real    Umin,
259                                     const Standard_Real    Usup,
260                                     const Standard_Real    Vmin,
261                                     const Standard_Real    Vsup,
262                                     const Standard_Real    TolU, 
263                                     const Standard_Real    TolV,
264                                     const Extrema_ExtFlag F,
265                                     const Extrema_ExtAlgo A) 
266 : myF (P,S),
267   myFlag(F),
268   myAlgo(A)
269 {
270   Initialize(S, NbU, NbV, Umin, Usup, Vmin, Vsup, TolU, TolV);
271   Perform(P);
272 }
273
274
275 void Extrema_GenExtPS::Initialize(const Adaptor3d_Surface& S,
276                                   const Standard_Integer NbU, 
277                                   const Standard_Integer NbV,
278                                   const Standard_Real    TolU, 
279                                   const Standard_Real    TolV)
280 {
281   myumin = S.FirstUParameter();
282   myusup = S.LastUParameter();
283   myvmin = S.FirstVParameter();
284   myvsup = S.LastVParameter();
285   Initialize(S,NbU,NbV,myumin,myusup,myvmin,myvsup,TolU,TolV);
286 }
287
288
289 void Extrema_GenExtPS::Initialize(const Adaptor3d_Surface& S,
290                                   const Standard_Integer NbU, 
291                                   const Standard_Integer NbV,
292                                   const Standard_Real    Umin,
293                                   const Standard_Real    Usup,
294                                   const Standard_Real    Vmin,
295                                   const Standard_Real    Vsup,
296                                   const Standard_Real    TolU, 
297                                   const Standard_Real    TolV)
298 {
299   myS = (Adaptor3d_SurfacePtr)&S;
300   myusample = NbU;
301   myvsample = NbV;
302   mytolu = TolU;
303   mytolv = TolV;
304   myumin = Umin;
305   myusup = Usup;
306   myvmin = Vmin;
307   myvsup = Vsup;
308
309   if ((myusample < 2) ||
310       (myvsample < 2)) { throw Standard_OutOfRange(); }
311
312   myF.Initialize(S);
313
314   mySphereUBTree.Nullify();
315   myUParams.Nullify();
316   myVParams.Nullify();
317   myInit = Standard_False;
318 }
319
320 inline static void fillParams(const TColStd_Array1OfReal& theKnots,
321                               Standard_Integer theDegree,
322                               Standard_Real theParMin,
323                               Standard_Real theParMax,
324                               Handle(TColStd_HArray1OfReal)& theParams,
325                               Standard_Integer theSample)
326 {
327   NCollection_Vector<Standard_Real> aParams;
328   Standard_Integer i = 1;
329   Standard_Real aPrevPar = theParMin;
330   aParams.Append(aPrevPar);
331   //calculation the array of parametric points depending on the knots array variation and degree of given surface
332   for ( ; i <  theKnots.Length() && theKnots(i) < (theParMax - Precision::PConfusion()); i++ )
333   {
334     if (  theKnots(i+1) < theParMin + Precision::PConfusion())
335       continue;
336
337     Standard_Real aStep = (theKnots(i+1) - theKnots(i))/Max(theDegree,2);
338     Standard_Integer k =1;
339     for( ; k <= theDegree ; k++)
340     {
341       Standard_Real aPar = theKnots(i) + k * aStep;
342       if(aPar > theParMax - Precision::PConfusion())
343         break;
344       if(aPar > aPrevPar + Precision::PConfusion() )
345       {
346         aParams.Append(aPar);
347         aPrevPar = aPar;
348       }
349     }
350
351   }
352   aParams.Append(theParMax);
353   Standard_Integer nbPar = aParams.Length();
354   //in case of an insufficient number of points the grid will be built later 
355   if (nbPar < theSample)
356     return;
357   theParams = new TColStd_HArray1OfReal(1, nbPar );
358   for( i = 0; i < nbPar; i++)
359     theParams->SetValue(i+1,aParams(i));
360 }
361
362 void Extrema_GenExtPS::GetGridPoints( const Adaptor3d_Surface& theSurf)
363 {
364   //creation parametric points for BSpline and Bezier surfaces
365   //with taking into account of Degree and NbKnots of BSpline or Bezier geometry
366   if(theSurf.GetType() == GeomAbs_OffsetSurface)
367     GetGridPoints(theSurf.BasisSurface()->Surface());
368   //parametric points for BSpline surfaces
369   else if( theSurf.GetType() == GeomAbs_BSplineSurface) 
370   {
371     Handle(Geom_BSplineSurface) aBspl = theSurf.BSpline();
372     if(!aBspl.IsNull())
373     {
374       TColStd_Array1OfReal aUKnots(1, aBspl->NbUKnots());
375       aBspl->UKnots( aUKnots);
376       TColStd_Array1OfReal aVKnots(1, aBspl->NbVKnots());
377       aBspl->VKnots( aVKnots);
378       fillParams(aUKnots,aBspl->UDegree(),myumin, myusup, myUParams, myusample);
379       fillParams(aVKnots,aBspl->VDegree(),myvmin, myvsup, myVParams, myvsample);
380     }
381   }
382   //calculation parametric points for Bezier surfaces
383   else if(theSurf.GetType() == GeomAbs_BezierSurface)
384   {
385     Handle(Geom_BezierSurface) aBezier = theSurf.Bezier();
386     if(aBezier.IsNull())
387       return;
388
389     TColStd_Array1OfReal aUKnots(1,2);
390     TColStd_Array1OfReal aVKnots(1,2);
391     aBezier->Bounds(aUKnots(1), aUKnots(2), aVKnots(1), aVKnots(2));
392     fillParams(aUKnots, aBezier->UDegree(), myumin, myusup, myUParams, myusample);
393     fillParams(aVKnots, aBezier->VDegree(), myvmin, myvsup, myVParams, myvsample);
394
395   }
396   //creation points for surfaces based on BSpline or Bezier curves
397   else if(theSurf.GetType() == GeomAbs_SurfaceOfRevolution || 
398     theSurf.GetType() == GeomAbs_SurfaceOfExtrusion)
399   {
400     Handle(TColStd_HArray1OfReal) anArrKnots;
401     Standard_Integer aDegree = 0;
402     if(theSurf.BasisCurve()->Curve().GetType() == GeomAbs_BSplineCurve)
403     {
404       Handle(Geom_BSplineCurve) aBspl = theSurf.BasisCurve()->Curve().BSpline();
405       if(!aBspl.IsNull())
406       {
407         anArrKnots = new TColStd_HArray1OfReal(1,aBspl->NbKnots());
408         aBspl->Knots( anArrKnots->ChangeArray1() );
409         aDegree = aBspl->Degree();
410         
411       }
412
413     }
414     if(theSurf.BasisCurve()->Curve().GetType() == GeomAbs_BezierCurve)
415     {
416       Handle(Geom_BezierCurve) aBez = theSurf.BasisCurve()->Curve().Bezier();
417       if(!aBez.IsNull())
418       {
419         anArrKnots = new TColStd_HArray1OfReal(1,2);
420         anArrKnots->SetValue(1, aBez->FirstParameter());
421         anArrKnots->SetValue(2, aBez->LastParameter());
422         aDegree = aBez->Degree();
423         
424       }
425     }
426     if(anArrKnots.IsNull())
427       return;
428     if( theSurf.GetType() == GeomAbs_SurfaceOfRevolution )
429       fillParams( anArrKnots->Array1(), aDegree, myvmin, myvsup, myVParams, myvsample);
430     else
431       fillParams( anArrKnots->Array1(), aDegree, myumin, myusup, myUParams, myusample);
432   }
433   //update the number of points in sample
434   if(!myUParams.IsNull())
435     myusample = myUParams->Length();
436   if( !myVParams.IsNull())
437     myvsample = myVParams->Length();
438   
439 }
440
441 /*
442  * This function computes the point on surface parameters on edge.
443  * if it coincides with theParam0 or theParam1, it is returned.
444  */
445 const Extrema_POnSurfParams& Extrema_GenExtPS::ComputeEdgeParameters
446       (const Standard_Boolean       IsUEdge,
447        const Extrema_POnSurfParams &theParam0,
448        const Extrema_POnSurfParams &theParam1,
449        const gp_Pnt                &thePoint,
450        const Standard_Real          theDiffTol)
451 {
452   const Standard_Real aSqrDist01 =
453     theParam0.Value().SquareDistance(theParam1.Value());
454
455   if (aSqrDist01 <= theDiffTol)
456   {
457     // The points are confused. Get the first point and change its type.
458     return theParam0;
459   }
460   else
461   {
462     const Standard_Real aDiffDist =
463       Abs(theParam0.GetSqrDistance() - theParam1.GetSqrDistance());
464
465     if (aDiffDist >= aSqrDist01 - theDiffTol)
466     {
467       // The shortest distance is one of the nodes.
468       if (theParam0.GetSqrDistance() > theParam1.GetSqrDistance())
469       {
470         // The shortest distance is the point 1.
471         return theParam1;
472       }
473       else
474       {
475         // The shortest distance is the point 0.
476         return theParam0;
477       }
478     }
479     else
480     {
481       // The shortest distance is inside the edge.
482       gp_XYZ aPoP(thePoint.XYZ().Subtracted(theParam0.Value().XYZ()));
483       gp_XYZ aPoP1(theParam1.Value().XYZ().Subtracted(theParam0.Value().XYZ()));
484       Standard_Real aRatio = aPoP.Dot(aPoP1)/aSqrDist01;
485       Standard_Real aU[2];
486       Standard_Real aV[2];
487
488       theParam0.Parameter(aU[0], aV[0]);
489       theParam1.Parameter(aU[1], aV[1]);
490
491       Standard_Real aUPar = aU[0];
492       Standard_Real aVPar = aV[0];
493
494       if (IsUEdge)
495       {
496         aUPar += aRatio*(aU[1] - aU[0]);
497       }
498       else
499       {
500         aVPar += aRatio*(aV[1] - aV[0]);
501       }
502
503       myGridParam.SetParameters(aUPar, aVPar, myS->Value(aUPar, aVPar));
504       Standard_Integer anIndices[2];
505
506       theParam0.GetIndices(anIndices[0], anIndices[1]);
507       myGridParam.SetElementType(IsUEdge ? Extrema_UIsoEdge : Extrema_VIsoEdge);
508       myGridParam.SetSqrDistance(thePoint.SquareDistance(myGridParam.Value()));
509       myGridParam.SetIndices(anIndices[0], anIndices[1]);
510       return myGridParam;
511     }
512   }
513 }
514
515 void Extrema_GenExtPS::BuildGrid(const gp_Pnt &thePoint)
516 {
517   Standard_Integer NoU, NoV;
518
519   //if grid was already built skip its creation
520   if (!myInit) {
521     //build parametric grid in case of a complex surface geometry (BSpline and Bezier surfaces)
522     GetGridPoints(*myS);
523     
524     //build grid in other cases 
525     if( myUParams.IsNull() )
526     {
527       Standard_Real PasU = myusup - myumin;
528       Standard_Real U0 = PasU / myusample / 100.;
529       PasU = (PasU - U0) / (myusample - 1);
530       U0 = U0/2. + myumin;
531       myUParams = new TColStd_HArray1OfReal(1,myusample );
532       Standard_Real U = U0;
533       for ( NoU = 1 ; NoU <= myusample; NoU++, U += PasU) 
534         myUParams->SetValue(NoU, U);
535     }
536   
537     if( myVParams.IsNull())
538     {
539       Standard_Real PasV = myvsup - myvmin;
540       Standard_Real V0 = PasV / myvsample / 100.;
541       PasV = (PasV - V0) / (myvsample - 1);
542       V0 = V0/2. + myvmin;
543       
544       myVParams = new TColStd_HArray1OfReal(1,myvsample );
545       Standard_Real V = V0;
546      
547       for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV)
548         myVParams->SetValue(NoV, V);
549     }
550
551     //If flag was changed and extrema not reinitialized Extrema would fail
552     myPoints = new Extrema_HArray2OfPOnSurfParams
553       (0, myusample + 1, 0, myvsample + 1);
554     // Calculation of distances
555   
556     for ( NoU = 1 ; NoU <= myusample; NoU++ ) {
557       for ( NoV = 1 ; NoV <= myvsample; NoV++) {
558         gp_Pnt aP1 = myS->Value(myUParams->Value(NoU), myVParams->Value(NoV));
559         Extrema_POnSurfParams aParam
560           (myUParams->Value(NoU), myVParams->Value(NoV), aP1);
561
562         aParam.SetElementType(Extrema_Node);
563         aParam.SetIndices(NoU, NoV);
564         myPoints->SetValue(NoU, NoV, aParam);
565       }
566     }
567
568     myFacePntParams = 
569       new Extrema_HArray2OfPOnSurfParams(0, myusample, 0, myvsample);
570
571     myUEdgePntParams =
572       new Extrema_HArray2OfPOnSurfParams(1, myusample - 1, 1, myvsample);
573     myVEdgePntParams =
574       new Extrema_HArray2OfPOnSurfParams(1, myusample, 1, myvsample - 1);
575
576     // Fill boundary with negative square distance.
577     // It is used for computation of Maximum.
578     for (NoV = 0; NoV <= myvsample + 1; NoV++) {
579       myPoints->ChangeValue(0, NoV).SetSqrDistance(-1.);
580       myPoints->ChangeValue(myusample + 1, NoV).SetSqrDistance(-1.);
581     }
582
583     for (NoU = 1; NoU <= myusample; NoU++) {
584       myPoints->ChangeValue(NoU, 0).SetSqrDistance(-1.);
585       myPoints->ChangeValue(NoU, myvsample + 1).SetSqrDistance(-1.);
586     }
587     
588     myInit = Standard_True;
589   }
590
591   // Compute distances to mesh.
592   // Step 1. Compute distances to nodes.
593   for ( NoU = 1 ; NoU <= myusample; NoU++ ) {
594     for ( NoV = 1 ; NoV <= myvsample; NoV++) {
595       Extrema_POnSurfParams &aParam = myPoints->ChangeValue(NoU, NoV);
596
597       aParam.SetSqrDistance(thePoint.SquareDistance(aParam.Value()));
598     }
599   }
600
601   // For search of minimum compute distances to mesh.
602   if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX) 
603   {
604     // This is the tolerance of difference of squared values.
605     // No need to set it too small.
606     const Standard_Real aDiffTol = mytolu + mytolv;
607
608     // Step 2. Compute distances to edges.
609     // Assume UEdge(i, j) = { Point(i, j); Point(i + 1, j    ) }
610     // Assume VEdge(i, j) = { Point(i, j); Point(i,     j + 1) }
611     for ( NoU = 1 ; NoU <= myusample; NoU++ ) 
612     {
613       for ( NoV = 1 ; NoV <= myvsample; NoV++)
614       {
615         const Extrema_POnSurfParams &aParam0 = myPoints->Value(NoU, NoV);
616
617         if (NoU < myusample)
618         {
619           // Compute parameters to UEdge.
620           const Extrema_POnSurfParams &aParam1 = myPoints->Value(NoU + 1, NoV);
621           const Extrema_POnSurfParams &anEdgeParam = ComputeEdgeParameters(Standard_True, aParam0, aParam1, thePoint, aDiffTol);
622
623           myUEdgePntParams->SetValue(NoU, NoV, anEdgeParam);
624         }
625
626         if (NoV < myvsample)
627         {
628           // Compute parameters to VEdge.
629           const Extrema_POnSurfParams &aParam1 = myPoints->Value(NoU, NoV + 1);
630           const Extrema_POnSurfParams &anEdgeParam = ComputeEdgeParameters(Standard_False, aParam0, aParam1, thePoint, aDiffTol);
631
632           myVEdgePntParams->SetValue(NoU, NoV, anEdgeParam);
633         }
634       }
635     }
636
637     // Step 3. Compute distances to faces.
638     // Assume myFacePntParams(i, j) =
639     // { Point(i, j); Point(i + 1, j); Point(i + 1, j + 1); Point(i, j + 1) }
640     //   Or
641     // { UEdge(i, j); VEdge(i + 1, j); UEdge(i, j + 1); VEdge(i, j) }
642     Standard_Real aSqrDist01;
643     Standard_Real aDiffDist;
644     Standard_Boolean isOut;
645
646     for ( NoU = 1 ; NoU < myusample; NoU++ ) {
647       for ( NoV = 1 ; NoV < myvsample; NoV++) {
648         const Extrema_POnSurfParams &aUE0 = myUEdgePntParams->Value(NoU, NoV);
649         const Extrema_POnSurfParams &aUE1 = myUEdgePntParams->Value(NoU, NoV+1);
650         const Extrema_POnSurfParams &aVE0 = myVEdgePntParams->Value(NoU, NoV);
651         const Extrema_POnSurfParams &aVE1 = myVEdgePntParams->Value(NoU+1, NoV);
652
653         aSqrDist01 = aUE0.Value().SquareDistance(aUE1.Value());
654         aDiffDist = Abs(aUE0.GetSqrDistance() - aUE1.GetSqrDistance());
655         isOut = Standard_False;
656
657         if (aDiffDist >= aSqrDist01 - aDiffTol) {
658           // The projection is outside the face.
659           isOut = Standard_True;
660         } else {
661           aSqrDist01 = aVE0.Value().SquareDistance(aVE1.Value());
662           aDiffDist = Abs(aVE0.GetSqrDistance() - aVE1.GetSqrDistance());
663
664           if (aDiffDist >= aSqrDist01 - aDiffTol) {
665             // The projection is outside the face.
666             isOut = Standard_True;
667           }
668         }
669
670         if (isOut) {
671           // Get the closest point on an edge.
672           const Extrema_POnSurfParams &aUEMin =
673             aUE0.GetSqrDistance() < aUE1.GetSqrDistance() ? aUE0 : aUE1;
674           const Extrema_POnSurfParams &aVEMin =
675             aVE0.GetSqrDistance() < aVE1.GetSqrDistance() ? aVE0 : aVE1;
676           const Extrema_POnSurfParams &aEMin =
677             aUEMin.GetSqrDistance() < aVEMin.GetSqrDistance() ? aUEMin : aVEMin;
678
679           myFacePntParams->SetValue(NoU, NoV, aEMin);
680         } else {
681           // Find closest point inside the face.
682           Standard_Real aU[2];
683           Standard_Real aV[2];
684           Standard_Real aUPar;
685           Standard_Real aVPar;
686
687           // Compute U parameter.
688           aUE0.Parameter(aU[0], aV[0]);
689           aUE1.Parameter(aU[1], aV[1]);
690           aUPar = 0.5*(aU[0] + aU[1]);
691           // Compute V parameter.
692           aVE0.Parameter(aU[0], aV[0]);
693           aVE1.Parameter(aU[1], aV[1]);
694           aVPar = 0.5*(aV[0] + aV[1]);
695
696           Extrema_POnSurfParams aParam(aUPar, aVPar, myS->Value(aUPar, aVPar));
697
698           aParam.SetElementType(Extrema_Face);
699           aParam.SetSqrDistance(thePoint.SquareDistance(aParam.Value()));
700           aParam.SetIndices(NoU, NoV);
701           myFacePntParams->SetValue(NoU, NoV, aParam);
702         }
703       }
704     }
705
706     // Fill boundary with RealLast square distance.
707     for (NoV = 0; NoV <= myvsample; NoV++) {
708       myFacePntParams->ChangeValue(0, NoV).SetSqrDistance(RealLast());
709       myFacePntParams->ChangeValue(myusample, NoV).SetSqrDistance(RealLast());
710     }
711
712     for (NoU = 1; NoU < myusample; NoU++) {
713       myFacePntParams->ChangeValue(NoU, 0).SetSqrDistance(RealLast());
714       myFacePntParams->ChangeValue(NoU, myvsample).SetSqrDistance(RealLast());
715     }
716   }
717 }
718
719 // Parametrization of the sample
720 void Extrema_GenExtPS::BuildTree()
721 {
722   // if tree already exists, assume it is already correctly filled
723   if ( ! mySphereUBTree.IsNull() )
724     return;
725
726    if (myS->GetType() == GeomAbs_BSplineSurface) {
727      Handle(Geom_BSplineSurface) aBspl = myS->BSpline();
728      Standard_Integer aUValue = aBspl->UDegree() * aBspl->NbUKnots();
729      Standard_Integer aVValue = aBspl->VDegree() * aBspl->NbVKnots();
730      if (aUValue > myusample)
731        myusample = aUValue;
732      if (aVValue > myvsample)
733        myvsample = aVValue;
734    }
735
736   Standard_Real PasU = myusup - myumin;
737   Standard_Real PasV = myvsup - myvmin;
738   Standard_Real U0 = PasU / myusample / 100.;
739   Standard_Real V0 = PasV / myvsample / 100.;
740   gp_Pnt P1;
741   PasU = (PasU - U0) / (myusample - 1);
742   PasV = (PasV - V0) / (myvsample - 1);
743   U0 = U0/2. + myumin;
744   V0 = V0/2. + myvmin;
745
746   //build grid of parametric points 
747   myUParams = new TColStd_HArray1OfReal(1,myusample );
748   myVParams = new TColStd_HArray1OfReal(1,myvsample );
749   Standard_Integer NoU, NoV;
750   Standard_Real U = U0, V = V0;
751   for ( NoU = 1 ; NoU <= myusample; NoU++, U += PasU) 
752     myUParams->SetValue(NoU, U);
753   for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV)
754     myVParams->SetValue(NoV, V);
755
756   // Calculation of distances
757   mySphereUBTree = new Extrema_UBTreeOfSphere;
758   Extrema_UBTreeFillerOfSphere aFiller(*mySphereUBTree);
759   Standard_Integer i = 0;
760   
761   mySphereArray = new Bnd_HArray1OfSphere(0, myusample * myvsample);
762  
763   for ( NoU = 1; NoU <= myusample; NoU++ ) {
764     for ( NoV = 1; NoV <= myvsample; NoV++) {
765       P1 = myS->Value(myUParams->Value(NoU), myVParams->Value(NoV));
766       Bnd_Sphere aSph(P1.XYZ(), 0/*mytolu < mytolv ? mytolu : mytolv*/, NoU, NoV);
767       aFiller.Add(i, aSph);
768       mySphereArray->SetValue( i, aSph );
769       i++;
770     }
771   }
772   aFiller.Fill();
773 }
774
775 void Extrema_GenExtPS::FindSolution(const gp_Pnt& /*P*/, 
776                                     const Extrema_POnSurfParams &theParams)
777 {
778   math_Vector Tol(1,2);
779   Tol(1) = mytolu;
780   Tol(2) = mytolv;
781
782   math_Vector UV(1, 2);
783   theParams.Parameter(UV(1), UV(2));
784
785   math_Vector UVinf(1,2), UVsup(1,2);
786   UVinf(1) = myumin;
787   UVinf(2) = myvmin;
788   UVsup(1) = myusup;
789   UVsup(2) = myvsup;
790
791   math_FunctionSetRoot S(myF, Tol);
792   S.Perform(myF, UV, UVinf, UVsup);
793
794   myDone = Standard_True;
795 }
796
797 void Extrema_GenExtPS::SetFlag(const Extrema_ExtFlag F)
798 {
799   myFlag = F;
800 }
801
802 void Extrema_GenExtPS::SetAlgo(const Extrema_ExtAlgo A)
803 {
804   if(myAlgo != A)
805      myInit = Standard_False;
806   myAlgo = A;
807  
808 }
809
810 void Extrema_GenExtPS::Perform(const gp_Pnt& P) 
811 {  
812   myDone = Standard_False;
813   myF.SetPoint(P);
814   
815   if(myAlgo == Extrema_ExtAlgo_Grad)
816   {
817     BuildGrid(P);
818     Standard_Integer NoU,NoV;
819
820     if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX) 
821     {
822       Extrema_ElementType anElemType;
823       Standard_Integer iU;
824       Standard_Integer iV;
825       Standard_Integer iU2;
826       Standard_Integer iV2;
827       Standard_Boolean isMin;
828       Standard_Integer i;
829
830       for (NoU = 1; NoU < myusample; NoU++) {
831         for (NoV = 1; NoV < myvsample; NoV++) {
832           const Extrema_POnSurfParams &aParam =
833             myFacePntParams->Value(NoU, NoV);
834
835           isMin = Standard_False;
836           anElemType = aParam.GetElementType();
837
838           if (anElemType == Extrema_Face) {
839             isMin = Standard_True;
840           } else {
841             // Check if it is a boundary edge or corner vertex.
842             aParam.GetIndices(iU, iV);
843
844             if (anElemType == Extrema_UIsoEdge) {
845               isMin = (iV == 1 || iV == myvsample);
846             } else if (anElemType == Extrema_VIsoEdge) {
847               isMin = (iU == 1 || iU == myusample);
848             } else if (anElemType == Extrema_Node) {
849               isMin = (iU == 1 || iU == myusample) &&
850                       (iV == 1 || iV == myvsample);
851             }
852
853             if (!isMin) {
854               // This is a middle element.
855               if (anElemType == Extrema_UIsoEdge ||
856                 (anElemType == Extrema_Node && (iU == 1 || iU == myusample))) {
857                 // Check the down face.
858                 const Extrema_POnSurfParams &aDownParam =
859                     myFacePntParams->Value(NoU, NoV - 1);
860
861                 if (aDownParam.GetElementType() == anElemType) {
862                   aDownParam.GetIndices(iU2, iV2);
863                   isMin = (iU == iU2 && iV == iV2);
864                 }
865               } else if (anElemType == Extrema_VIsoEdge ||
866                 (anElemType == Extrema_Node && (iV == 1 || iV == myvsample))) {
867                 // Check the right face.
868                 const Extrema_POnSurfParams &aRightParam =
869                     myFacePntParams->Value(NoU - 1, NoV);
870
871                 if (aRightParam.GetElementType() == anElemType) {
872                   aRightParam.GetIndices(iU2, iV2);
873                   isMin = (iU == iU2 && iV == iV2);
874                 }
875               } else if (iU == NoU && iV == NoV) {
876                 // Check the lower-left node. For this purpose it is necessary
877                 // to check lower-left, lower and left faces.
878                 isMin = Standard_True;
879
880                 const Extrema_POnSurfParams *anOtherParam[3] =
881                 { &myFacePntParams->Value(NoU, NoV - 1),     // Down
882                   &myFacePntParams->Value(NoU - 1, NoV - 1), // Lower-left
883                   &myFacePntParams->Value(NoU - 1, NoV) };   // Left
884
885                 for (i = 0; i < 3 && isMin; i++) {
886                   if (anOtherParam[i]->GetElementType() == Extrema_Node) {
887                     anOtherParam[i]->GetIndices(iU2, iV2);
888                     isMin = (iU == iU2 && iV == iV2);
889                   } else {
890                     isMin = Standard_False;
891                   }
892                 }
893               }
894             }
895           }
896
897           if (isMin) {
898             FindSolution(P, aParam);
899           }
900         }
901       }
902     }
903     
904     if(myFlag == Extrema_ExtFlag_MAX || myFlag == Extrema_ExtFlag_MINMAX)
905     {
906       Standard_Real Dist;
907
908       for (NoU = 1; NoU <= myusample; NoU++)
909       {
910         for (NoV = 1; NoV <= myvsample; NoV++)
911         {
912           const Extrema_POnSurfParams &aParamMain = myPoints->Value(NoU, NoV);
913           const Extrema_POnSurfParams &aParam1 = myPoints->Value(NoU - 1, NoV - 1);
914           const Extrema_POnSurfParams &aParam2 = myPoints->Value(NoU - 1, NoV);
915           const Extrema_POnSurfParams &aParam3 = myPoints->Value(NoU - 1, NoV + 1);
916           const Extrema_POnSurfParams &aParam4 = myPoints->Value(NoU, NoV - 1);
917           const Extrema_POnSurfParams &aParam5 = myPoints->Value(NoU, NoV + 1);
918           const Extrema_POnSurfParams &aParam6 = myPoints->Value(NoU + 1, NoV - 1);
919           const Extrema_POnSurfParams &aParam7 = myPoints->Value(NoU + 1, NoV);
920           const Extrema_POnSurfParams &aParam8 = myPoints->Value(NoU + 1, NoV + 1);
921
922           Dist = aParamMain.GetSqrDistance();
923
924           if ((aParam1.GetSqrDistance() <= Dist) &&
925               (aParam2.GetSqrDistance() <= Dist) &&
926               (aParam3.GetSqrDistance() <= Dist) &&
927               (aParam4.GetSqrDistance() <= Dist) &&
928               (aParam5.GetSqrDistance() <= Dist) &&
929               (aParam6.GetSqrDistance() <= Dist) &&
930               (aParam7.GetSqrDistance() <= Dist) &&
931               (aParam8.GetSqrDistance() <= Dist))
932           {
933             // Find maximum.
934             FindSolution(P, myPoints->Value(NoU, NoV));
935           }
936         }
937       }
938     }
939   }
940   else
941   {
942     BuildTree();
943     if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX)
944     {
945       Bnd_Sphere aSol = mySphereArray->Value(0);
946       Bnd_SphereUBTreeSelectorMin aSelector(mySphereArray, aSol);
947       //aSelector.SetMaxDist( RealLast() );
948       aSelector.DefineCheckPoint( P );
949       mySphereUBTree->Select( aSelector );
950       //TODO: check if no solution in binary tree
951       Bnd_Sphere& aSph = aSelector.Sphere();
952       Standard_Real aU = myUParams->Value(aSph.U());
953       Standard_Real aV = myVParams->Value(aSph.V());
954       Extrema_POnSurfParams aParams(aU, aV, myS->Value(aU, aV));
955
956       aParams.SetSqrDistance(P.SquareDistance(aParams.Value()));
957       aParams.SetIndices(aSph.U(), aSph.V());
958       FindSolution(P, aParams);
959     }
960     if(myFlag == Extrema_ExtFlag_MAX || myFlag == Extrema_ExtFlag_MINMAX)
961     {
962       Bnd_Sphere aSol = mySphereArray->Value(0);
963       Bnd_SphereUBTreeSelectorMax aSelector(mySphereArray, aSol);
964       //aSelector.SetMaxDist( RealLast() );
965       aSelector.DefineCheckPoint( P );
966       mySphereUBTree->Select( aSelector );
967       //TODO: check if no solution in binary tree
968       Bnd_Sphere& aSph = aSelector.Sphere();
969       Standard_Real aU = myUParams->Value(aSph.U());
970       Standard_Real aV = myVParams->Value(aSph.V());
971       Extrema_POnSurfParams aParams(aU, aV, myS->Value(aU, aV));
972       aParams.SetSqrDistance(P.SquareDistance(aParams.Value()));
973       aParams.SetIndices(aSph.U(), aSph.V());
974
975       FindSolution(P, aParams);
976     }
977   }
978 }
979 //=============================================================================
980
981 Standard_Boolean Extrema_GenExtPS::IsDone () const { return myDone; }
982 //=============================================================================
983
984 Standard_Integer Extrema_GenExtPS::NbExt () const
985 {
986   if (!IsDone()) { throw StdFail_NotDone(); }
987   return myF.NbExt();
988 }
989 //=============================================================================
990
991 Standard_Real Extrema_GenExtPS::SquareDistance (const Standard_Integer N) const
992 {
993   if ((N < 1) || (N > NbExt()))
994   {
995     throw Standard_OutOfRange();
996   }
997
998   return myF.SquareDistance(N);
999 }
1000 //=============================================================================
1001
1002 const Extrema_POnSurf& Extrema_GenExtPS::Point (const Standard_Integer N) const
1003 {
1004   if ((N < 1) || (N > NbExt()))
1005   {
1006     throw Standard_OutOfRange();
1007   }
1008
1009   return myF.Point(N);
1010 }
1011 //=============================================================================