0029151: GCC 7.1 warnings "this statement may fall through" [-Wimplicit-fallthrough=]
[occt.git] / src / BndLib / BndLib_AddSurface.cxx
1 // Created on: 1995-07-24
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 - Fri Aug 27 12:29:04 2004 OCC6503
18
19 #include <Adaptor3d_HSurface.hxx>
20 #include <Adaptor3d_Surface.hxx>
21 #include <Bnd_Box.hxx>
22 #include <BndLib.hxx>
23 #include <BndLib_AddSurface.hxx>
24 #include <BSplCLib.hxx>
25 #include <ElSLib.hxx>
26 #include <ElCLib.hxx>
27 #include <Geom_BezierSurface.hxx>
28 #include <Geom_BSplineSurface.hxx>
29 #include <GeomAbs_SurfaceType.hxx>
30 #include <gp_Pln.hxx>
31 #include <gp_Pnt.hxx>
32 #include <gp_Cylinder.hxx>
33 #include <gp_Cone.hxx>
34 #include <gp_Lin.hxx>
35 #include <Precision.hxx>
36 #include <TColgp_Array2OfPnt.hxx>
37 #include <TColStd_Array1OfInteger.hxx>
38 #include <TColStd_Array1OfReal.hxx>
39 #include <BndLib_Add3dCurve.hxx>
40 #include <math_MultipleVarFunction.hxx>
41 #include <math_PSO.hxx>
42 #include <math_Matrix.hxx>
43 #include <math_Powell.hxx>
44 //
45 static Standard_Integer NbUSamples(const Adaptor3d_Surface& S, 
46                                    const Standard_Real Umin,
47                                    const Standard_Real Umax);
48 //
49 static Standard_Integer NbVSamples(const Adaptor3d_Surface& S, 
50                                    const Standard_Real Vmin,
51                                    const Standard_Real Vmax);
52 //
53 static Standard_Real  AdjustExtr(const Adaptor3d_Surface& S, 
54                                  const Standard_Real UMin,
55                                                    const Standard_Real UMax,
56                                                    const Standard_Real VMin,
57                                                    const Standard_Real VMax,
58                                  const Standard_Real Extr0,
59                                  const Standard_Integer CoordIndx,                                 
60                                  const Standard_Real Tol, 
61                                  const Standard_Boolean IsMin);
62
63 //=======================================================================
64 //function : Add
65 //purpose  : 
66 //=======================================================================
67 void BndLib_AddSurface::Add(const Adaptor3d_Surface& S,
68                             const Standard_Real Tol,
69                             Bnd_Box& B ) 
70 {
71
72   BndLib_AddSurface::Add(S,
73                          S.FirstUParameter(),
74                          S.LastUParameter (),
75                          S.FirstVParameter(),
76                          S.LastVParameter (),Tol,B);
77 }
78 //=======================================================================
79 //function : NbUSamples
80 //purpose  : 
81 //=======================================================================
82
83 static Standard_Integer NbUSamples(const Adaptor3d_Surface& S) 
84 {
85   Standard_Integer N;
86   GeomAbs_SurfaceType Type = S.GetType();
87   switch (Type) {
88   case GeomAbs_BezierSurface: 
89     {
90       N = 2*S.NbUPoles();
91       break;
92     }
93   case GeomAbs_BSplineSurface: 
94     {
95       const Handle(Geom_BSplineSurface)& BS = S.BSpline();
96       N = 2*(BS->UDegree() + 1)*(BS->NbUKnots() -1);
97       break;
98     }
99   default:
100     N = 33;
101   }
102   return Min (50,N);
103 }
104
105 //=======================================================================
106 //function : NbVSamples
107 //purpose  : 
108 //=======================================================================
109
110 static Standard_Integer NbVSamples(const Adaptor3d_Surface& S) 
111 {
112   Standard_Integer N;
113   GeomAbs_SurfaceType Type = S.GetType();
114   switch (Type) {
115   case GeomAbs_BezierSurface:
116     {
117       N = 2*S.NbVPoles();
118       break;
119     }
120   case GeomAbs_BSplineSurface:
121     {
122       const Handle(Geom_BSplineSurface)& BS = S.BSpline();
123       N = 2*(BS->VDegree() + 1)*(BS->NbVKnots() - 1) ;
124       break;
125     }
126   default:
127     N = 33;
128   }
129   return Min(50,N);
130 }
131
132 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 Begin
133 static gp_Pnt BaryCenter(const gp_Pln        &aPlane,
134                          const Standard_Real  aUMin,
135                          const Standard_Real  aUMax,
136                          const Standard_Real  aVMin,
137                          const Standard_Real  aVMax)
138 {
139   Standard_Real aU, aV;
140   Standard_Boolean isU1Inf = Precision::IsInfinite(aUMin);
141   Standard_Boolean isU2Inf = Precision::IsInfinite(aUMax);
142   Standard_Boolean isV1Inf = Precision::IsInfinite(aVMin);
143   Standard_Boolean isV2Inf = Precision::IsInfinite(aVMax);
144
145   if (isU1Inf && isU2Inf)
146     aU = 0;
147   else if (isU1Inf)
148     aU = aUMax - 10.;
149   else if (isU2Inf)
150     aU = aUMin + 10.;
151   else
152     aU = (aUMin + aUMax)/2.;
153     
154   if (isV1Inf && isV2Inf)
155     aV = 0;
156   else if (isV1Inf)
157     aV = aVMax - 10.;
158   else if (isV2Inf)
159     aV = aVMin + 10.;
160   else
161     aV = (aVMin + aVMax)/2.;
162
163   gp_Pnt aCenter = ElSLib::Value(aU, aV, aPlane);
164
165   return aCenter;
166 }
167
168 static void TreatInfinitePlane(const gp_Pln        &aPlane,
169                                const Standard_Real  aUMin,
170                                const Standard_Real  aUMax,
171                                const Standard_Real  aVMin,
172                                const Standard_Real  aVMax,
173                                const Standard_Real  aTol,
174                                      Bnd_Box       &aB)
175 {
176   // Get 3 coordinate axes of the plane.
177   const gp_Dir        &aNorm        = aPlane.Axis().Direction();
178   const Standard_Real  anAngularTol = RealEpsilon();
179
180   // Get location of the plane as its barycenter
181   gp_Pnt aLocation = BaryCenter(aPlane, aUMin, aUMax, aVMin, aVMax);
182
183   if (aNorm.IsParallel(gp::DX(), anAngularTol)) {
184     aB.Add(aLocation);
185     aB.OpenYmin();
186     aB.OpenYmax();
187     aB.OpenZmin();
188     aB.OpenZmax();
189   } else if (aNorm.IsParallel(gp::DY(), anAngularTol)) {
190     aB.Add(aLocation);
191     aB.OpenXmin();
192     aB.OpenXmax();
193     aB.OpenZmin();
194     aB.OpenZmax();
195   } else if (aNorm.IsParallel(gp::DZ(), anAngularTol)) {
196     aB.Add(aLocation);
197     aB.OpenXmin();
198     aB.OpenXmax();
199     aB.OpenYmin();
200     aB.OpenYmax();
201   } else {
202     aB.SetWhole();
203     return;
204   }
205
206   aB.Enlarge(aTol);
207 }
208
209 // Compute start and finish indexes used in convex hull.
210 // theMinIdx - minimum poles index, that can be used.
211 // theMaxIdx - maximum poles index, that can be used.
212 // theShiftCoeff - shift between flatknots array and poles array.
213 // This vaule should be equal to 1 in case of non periodic BSpline,
214 // and (degree + 1) - mults(the lowest index).
215 void ComputePolesIndexes(const TColStd_Array1OfReal &theFlatKnots,
216                          const Standard_Integer theDegree,
217                          const Standard_Real theMin,
218                          const Standard_Real theMax,
219                          const Standard_Integer theMinIdx,
220                          const Standard_Integer theMaxIdx,
221                          const Standard_Integer theShiftCoeff,
222                          Standard_Integer &theOutMinIdx,
223                          Standard_Integer &theOutMaxIdx)
224 {
225   // Set initial values for the result indexes to handle situation when requested parameter space
226   // is slightly greater than B-spline parameter space.
227   theOutMinIdx = theFlatKnots.Lower();
228   theOutMaxIdx = theFlatKnots.Upper();
229
230   // Compute first and last used flat knots.
231   for(Standard_Integer aKnotIdx = theFlatKnots.Lower();
232       aKnotIdx < theFlatKnots.Upper();
233       aKnotIdx++)
234   {
235     if (theFlatKnots(aKnotIdx) <= theMin)
236       theOutMinIdx = aKnotIdx;
237
238     if (theFlatKnots(theFlatKnots.Upper() - aKnotIdx + theFlatKnots.Lower()) >= theMax)
239       theOutMaxIdx = theFlatKnots.Upper() - aKnotIdx + theFlatKnots.Lower();
240   }
241
242   theOutMinIdx = Max(theOutMinIdx - 2 * theDegree + 2 - theShiftCoeff, theMinIdx);
243   theOutMaxIdx = Min(theOutMaxIdx - 2 + theDegree + 1 - theShiftCoeff, theMaxIdx);
244 }
245
246 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 End
247 //=======================================================================
248 //function : Add
249 //purpose  : 
250 //=======================================================================
251 void BndLib_AddSurface::Add(const Adaptor3d_Surface& S,
252                             const Standard_Real UMin,
253                             const Standard_Real UMax,
254                             const Standard_Real VMin,
255                             const Standard_Real VMax,
256                             const Standard_Real Tol,
257                             Bnd_Box& B ) 
258 {
259   GeomAbs_SurfaceType Type = S.GetType(); // skv OCC6503
260
261   if (Precision::IsInfinite(VMin) ||
262       Precision::IsInfinite(VMax) ||
263       Precision::IsInfinite(UMin) ||
264       Precision::IsInfinite(UMax)   ) {
265 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 Begin
266 //     B.SetWhole();
267 //     return;
268     switch (Type) {
269     case GeomAbs_Plane: 
270       {
271         TreatInfinitePlane(S.Plane(), UMin, UMax, VMin, VMax, Tol, B);
272         return;
273       }
274     default: 
275       {
276         B.SetWhole();
277         return;
278       }
279     }
280 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 End
281   }
282
283 //   GeomAbs_SurfaceType Type = S.GetType(); // skv OCC6503
284
285   switch (Type) {
286
287   case GeomAbs_Plane: 
288     {
289       gp_Pln Plan = S.Plane();
290       B.Add(ElSLib::Value(UMin,VMin,Plan)); 
291       B.Add(ElSLib::Value(UMin,VMax,Plan)); 
292       B.Add(ElSLib::Value(UMax,VMin,Plan)); 
293       B.Add(ElSLib::Value(UMax,VMax,Plan)); 
294       B.Enlarge(Tol);
295       break;
296     }
297   case GeomAbs_Cylinder: 
298     {
299       BndLib::Add(S.Cylinder(),UMin,UMax,VMin,VMax,Tol,B);
300       break;
301     }
302   case GeomAbs_Cone: 
303     {
304       BndLib::Add(S.Cone(),UMin,UMax,VMin,VMax,Tol,B);
305       break;
306     }
307   case GeomAbs_Torus: 
308     {
309       BndLib::Add(S.Torus(),UMin,UMax,VMin,VMax,Tol,B);
310       break;
311     }
312   case GeomAbs_Sphere: 
313     {
314       if (Abs(UMin) < Precision::Angular() &&
315           Abs(UMax - 2.*M_PI) < Precision::Angular() &&
316           Abs(VMin + M_PI/2.) < Precision::Angular() &&
317           Abs(VMax - M_PI/2.) < Precision::Angular()) // a whole sphere
318         BndLib::Add(S.Sphere(),Tol,B);
319       else
320         BndLib::Add(S.Sphere(),UMin,UMax,VMin,VMax,Tol,B);
321       break;
322     }
323   case GeomAbs_OffsetSurface: 
324     {
325       Handle(Adaptor3d_HSurface) HS = S.BasisSurface();
326       Add (HS->Surface(),UMin,UMax,VMin,VMax,Tol,B);
327       B.Enlarge(S.OffsetValue());
328       B.Enlarge(Tol);
329       break;
330     } 
331   case GeomAbs_BezierSurface:
332   case GeomAbs_BSplineSurface: 
333     {
334       Standard_Boolean isUseConvexHullAlgorithm = Standard_True;
335       Standard_Real PTol = Precision::Parametric(Precision::Confusion());
336       // Borders of underlying geometry.
337       Standard_Real anUMinParam = UMin, anUMaxParam = UMax,// BSpline case.
338                      aVMinParam = VMin,  aVMaxParam = VMax;
339       if (Type == GeomAbs_BezierSurface)
340       {
341         // Bezier surface:
342         // All of poles used for any parameter,
343         // thats why in case of trimmed parameters handled by grid algorithm.
344
345         if (Abs(UMin-S.FirstUParameter()) > PTol ||
346             Abs(VMin-S.FirstVParameter()) > PTol ||
347             Abs(UMax-S.LastUParameter ()) > PTol ||
348             Abs(VMax-S.LastVParameter ()) > PTol )
349         {
350           // Borders not equal to topology borders.
351           isUseConvexHullAlgorithm = Standard_False;
352         }
353       }
354       else
355       {
356         // BSpline:
357         // If Umin, Vmin, Umax, Vmax lies inside geometry bounds then:
358         // use convex hull algorithm,
359         // if Umin, VMin, Umax, Vmax lies outside then:
360         // use grid algorithm on analytic continuation (default case).
361         S.BSpline()->Bounds(anUMinParam, anUMaxParam, aVMinParam, aVMaxParam);
362
363         if ( (UMin - anUMinParam) < -PTol ||
364              (VMin -  aVMinParam) < -PTol ||
365              (UMax - anUMaxParam) >  PTol ||
366              (VMax -  aVMaxParam) >  PTol )
367         {
368           // Out of geometry borders.
369           isUseConvexHullAlgorithm = Standard_False;
370         }
371       }
372
373       if (isUseConvexHullAlgorithm)
374       {
375           TColgp_Array2OfPnt Tp(1,S.NbUPoles(),1,S.NbVPoles());
376           Standard_Integer UMinIdx = 0, UMaxIdx = 0;
377           Standard_Integer VMinIdx = 0, VMaxIdx = 0;
378           if (Type == GeomAbs_BezierSurface)
379           {
380             S.Bezier()->Poles(Tp);
381
382             UMinIdx = Tp.LowerRow();
383             UMaxIdx = Tp.UpperRow();
384             VMinIdx = Tp.LowerCol();
385             VMaxIdx = Tp.UpperCol();
386           }
387           else
388           {
389             S.BSpline()->Poles(Tp);
390
391             UMinIdx = Tp.LowerRow();
392             UMaxIdx = Tp.UpperRow();
393             VMinIdx = Tp.LowerCol();
394             VMaxIdx = Tp.UpperCol();
395
396             if (UMin > anUMinParam ||
397                 UMax < anUMaxParam)
398             {
399               Standard_Integer anUFlatKnotsCount = S.BSpline()->NbUPoles() + S.BSpline()->UDegree() + 1;
400               Standard_Integer aShift = 1;
401
402               if (S.BSpline()->IsUPeriodic())
403               {
404                 TColStd_Array1OfInteger aMults(1, S.BSpline()->NbUKnots());
405                 S.BSpline()->UMultiplicities(aMults);
406                 anUFlatKnotsCount = BSplCLib::KnotSequenceLength(aMults, S.BSpline()->UDegree(), Standard_True);
407
408                 aShift = S.BSpline()->UDegree() + 1 - S.BSpline()->UMultiplicity(1);
409               }
410
411               TColStd_Array1OfReal anUFlatKnots(1, anUFlatKnotsCount);
412               S.BSpline()->UKnotSequence(anUFlatKnots);
413
414               ComputePolesIndexes(anUFlatKnots,
415                                   S.BSpline()->UDegree(),
416                                   UMin, UMax,
417                                   UMinIdx, UMaxIdx,  // Min and Max Indexes
418                                   aShift,
419                                   UMinIdx, UMaxIdx); // the Output indexes
420             }
421
422             if (VMin > aVMinParam ||
423                 VMax < aVMaxParam)
424             {
425               Standard_Integer anVFlatKnotsCount = S.BSpline()->NbVPoles() + S.BSpline()->VDegree() + 1;
426               Standard_Integer aShift = 1;
427
428               if (S.BSpline()->IsVPeriodic())
429               {
430                 TColStd_Array1OfInteger aMults(1, S.BSpline()->NbVKnots());
431                 S.BSpline()->VMultiplicities(aMults);
432                 anVFlatKnotsCount = BSplCLib::KnotSequenceLength(aMults, S.BSpline()->VDegree(), Standard_True);
433
434                 aShift = S.BSpline()->VDegree() + 1 - S.BSpline()->VMultiplicity(1);
435               }
436
437               TColStd_Array1OfReal anVFlatKnots(1, anVFlatKnotsCount);
438               S.BSpline()->VKnotSequence(anVFlatKnots);
439
440               ComputePolesIndexes(anVFlatKnots,
441                                   S.BSpline()->VDegree(),
442                                   VMin, VMax,
443                                   VMinIdx, VMaxIdx,  // Min and Max Indexes
444                                   aShift,
445                                   VMinIdx, VMaxIdx); // the Output indexes
446             }
447
448           }
449
450           // Use poles to build convex hull.
451           for (Standard_Integer i = UMinIdx; i <= UMaxIdx; i++)
452           {
453             for (Standard_Integer j = VMinIdx; j <= VMaxIdx; j++)
454             {
455               B.Add(Tp(i,j));
456             }
457           }
458
459           B.Enlarge(Tol);
460           break;
461       }
462     }
463     Standard_FALLTHROUGH
464   default: 
465     {
466       Standard_Integer Nu = NbUSamples(S);
467       Standard_Integer Nv = NbVSamples(S);
468       gp_Pnt P;
469       for (Standard_Integer i =1 ;i<=Nu;i++){
470         Standard_Real U = UMin + ((UMax-UMin)*(i-1)/(Nu-1));
471         for (Standard_Integer j=1 ;j<=Nv;j++){
472           Standard_Real V = VMin + ((VMax-VMin)*(j-1)/(Nv-1));
473           S.D0(U,V,P);
474           B.Add(P);
475         }
476       }
477       B.Enlarge(Tol);
478     }
479   }
480 }
481 //----- Methods for AddOptimal ---------------------------------------
482
483 //=======================================================================
484 //function : AddOptimal
485 //purpose  : 
486 //=======================================================================
487 void BndLib_AddSurface::AddOptimal(const Adaptor3d_Surface& S,
488                                                      const Standard_Real Tol,
489                                                      Bnd_Box& B ) 
490 {
491
492   BndLib_AddSurface::AddOptimal(S,
493                                                   S.FirstUParameter(),
494                                                   S.LastUParameter (),
495                                                   S.FirstVParameter(),
496                                                   S.LastVParameter (),Tol,B);
497 }
498 //=======================================================================
499 //function : AddOptimal
500 //purpose  : 
501 //=======================================================================
502
503 void BndLib_AddSurface::AddOptimal(const Adaptor3d_Surface& S,
504                             const Standard_Real UMin,
505                             const Standard_Real UMax,
506                             const Standard_Real VMin,
507                             const Standard_Real VMax,
508                             const Standard_Real Tol,
509                             Bnd_Box& B ) 
510 {
511   GeomAbs_SurfaceType Type = S.GetType(); 
512
513   if (Precision::IsInfinite(VMin) ||
514       Precision::IsInfinite(VMax) ||
515       Precision::IsInfinite(UMin) ||
516       Precision::IsInfinite(UMax)   ) {
517     switch (Type) {
518       case GeomAbs_Plane: 
519       {
520         TreatInfinitePlane(S.Plane(), UMin, UMax, VMin, VMax, Tol, B);
521         return;
522       }
523       default: 
524       {
525               B.SetWhole();
526               return;
527       }
528     }
529   }
530
531   switch (Type) {
532     
533     case GeomAbs_Plane: 
534     {
535       gp_Pln Plan = S.Plane();
536       B.Add(ElSLib::Value(UMin,VMin,Plan)); 
537       B.Add(ElSLib::Value(UMin,VMax,Plan)); 
538       B.Add(ElSLib::Value(UMax,VMin,Plan)); 
539       B.Add(ElSLib::Value(UMax,VMax,Plan)); 
540       B.Enlarge(Tol);
541       break;
542     }
543     case GeomAbs_Cylinder: 
544     {
545       BndLib::Add(S.Cylinder(), UMin, UMax, VMin, VMax, Tol, B);
546       break;
547     }
548     case GeomAbs_Cone: 
549     {
550       BndLib::Add(S.Cone(), UMin, UMax, VMin, VMax, Tol, B);
551       break;
552     }
553     case GeomAbs_Sphere: 
554     {
555       BndLib::Add(S.Sphere(), UMin, UMax, VMin, VMax, Tol, B); 
556       break;
557     }
558     default: 
559     {
560       AddGenSurf(S, UMin, UMax, VMin, VMax, Tol, B);
561     }
562   }
563 }
564 //=======================================================================
565 //function : AddGenSurf
566 //purpose  : 
567 //=======================================================================
568 void BndLib_AddSurface::AddGenSurf(const Adaptor3d_Surface& S, 
569                                    const Standard_Real UMin,
570                                    const Standard_Real UMax,
571                                    const Standard_Real VMin,
572                                    const Standard_Real VMax,
573                                    const Standard_Real Tol,
574                                    Bnd_Box& B)
575 {
576   Standard_Integer Nu = NbUSamples(S, UMin, UMax);
577   Standard_Integer Nv = NbVSamples(S, VMin, VMax);
578   //
579   Standard_Real CoordMin[3] = {RealLast(), RealLast(), RealLast()}; 
580   Standard_Real CoordMax[3] = {-RealLast(), -RealLast(), -RealLast()};
581   Standard_Real DeflMax[3] = {-RealLast(), -RealLast(), -RealLast()};
582   //
583   //
584   Standard_Real du = (UMax-UMin)/(Nu-1), du2 = du / 2.;
585   Standard_Real dv = (VMax-VMin)/(Nv-1), dv2 = dv / 2.;
586   NCollection_Array2<gp_XYZ> aPnts(1, Nu, 1, Nv);
587   Standard_Real u, v;
588   Standard_Integer i, j, k;
589   gp_Pnt P;
590   for (i = 1, u = UMin; i <= Nu; i++, u += du){
591     for (j = 1, v = VMin;j <= Nv; j++, v += dv){
592       S.D0(u,v,P);
593       aPnts(i, j) = P.XYZ();
594       //
595       for(k = 0; k < 3; ++k)
596       {
597         if(CoordMin[k] > P.Coord(k+1))
598         {
599           CoordMin[k] = P.Coord(k+1);
600         }
601         if(CoordMax[k] < P.Coord(k+1))
602         {
603           CoordMax[k] = P.Coord(k+1);
604         }
605       }
606       //
607       if(i > 1)
608       {
609         gp_XYZ aPm = 0.5 * (aPnts(i-1,j) + aPnts(i, j));
610         S.D0(u - du2, v, P);
611         gp_XYZ aD = (P.XYZ() - aPm);
612         for(k = 0; k < 3; ++k)
613         {
614           if(CoordMin[k] > P.Coord(k+1))
615           {
616             CoordMin[k] = P.Coord(k+1);
617           }
618           if(CoordMax[k] < P.Coord(k+1))
619           {
620             CoordMax[k] = P.Coord(k+1);
621           }
622           Standard_Real d = Abs(aD.Coord(k+1));
623           if(DeflMax[k] < d)
624           {
625             DeflMax[k] = d;
626           }
627         }
628       }
629       if(j > 1)
630       {
631         gp_XYZ aPm = 0.5 * (aPnts(i,j-1) + aPnts(i, j));
632         S.D0(u , v - dv2, P);
633         gp_XYZ aD = (P.XYZ() - aPm);
634         for(k = 0; k < 3; ++k)
635         {
636           if(CoordMin[k] > P.Coord(k+1))
637           {
638             CoordMin[k] = P.Coord(k+1);
639           }
640           if(CoordMax[k] < P.Coord(k+1))
641           {
642             CoordMax[k] = P.Coord(k+1);
643           }
644           Standard_Real d = Abs(aD.Coord(k+1));
645           if(DeflMax[k] < d)
646           {
647             DeflMax[k] = d;
648           }
649         }
650       }
651     }
652   }
653   //
654   //Adjusting minmax 
655   Standard_Real eps = Max(Tol, Precision::Confusion()); 
656   for(k = 0; k < 3; ++k)
657   {
658     Standard_Real d = DeflMax[k];
659     if(d <= eps)
660     {
661       continue;
662     }
663
664     Standard_Real CMin = CoordMin[k];
665     Standard_Real CMax = CoordMax[k];
666     for(i = 1; i <= Nu; ++i)
667     {
668       for(j = 1; j <= Nv; ++j)
669       {
670         if(aPnts(i,j).Coord(k+1) - CMin < d)
671         {
672           Standard_Real umin, umax, vmin, vmax;
673           umin = UMin + Max(0, i-2) * du;
674           umax = UMin + Min(Nu-1, i) * du;
675           vmin = VMin + Max(0, j-2) * dv;
676           vmax = VMin + Min(Nv-1, j) * dv;
677           Standard_Real cmin = AdjustExtr(S, umin, umax, vmin, vmax,
678                                           CMin, k + 1, eps, Standard_True);
679           if(cmin < CMin)
680           {
681             CMin = cmin;
682           }
683         }
684         else if(CMax - aPnts(i,j).Coord(k+1) < d)
685         {
686           Standard_Real umin, umax, vmin, vmax;
687           umin = UMin + Max(0, i-2) * du;
688           umax = UMin + Min(Nu-1, i) * du;
689           vmin = VMin + Max(0, j-2) * dv;
690           vmax = VMin + Min(Nv-1, j) * dv;
691           Standard_Real cmax = AdjustExtr(S, umin, umax, vmin, vmax,
692                                           CMax, k + 1, eps, Standard_False);
693           if(cmax > CMax)
694           {
695             CMax = cmax;
696           }
697         }
698       }
699     }
700     CoordMin[k] = CMin;
701     CoordMax[k] = CMax;
702
703   }
704
705   B.Add(gp_Pnt(CoordMin[0], CoordMin[1], CoordMin[2]));
706   B.Add(gp_Pnt(CoordMax[0], CoordMax[1], CoordMax[2]));
707   B.Enlarge(eps);
708 }
709 //
710
711 //
712 class SurfMaxMinCoord : public math_MultipleVarFunction
713 {
714 public:
715   SurfMaxMinCoord(const Adaptor3d_Surface& theSurf, 
716               const Standard_Real UMin,
717               const Standard_Real UMax,
718               const Standard_Real VMin,
719               const Standard_Real VMax,
720               const Standard_Integer CoordIndx,                                 
721               const Standard_Real Sign)
722 : mySurf(theSurf),
723   myUMin(UMin),
724   myUMax(UMax),
725   myVMin(VMin),
726   myVMax(VMax),
727   myCoordIndx(CoordIndx),
728   mySign(Sign)
729   {
730     math_Vector X(1,2);
731     X(1) = UMin;
732     X(2) = (VMin + VMax) / 2.;
733     Standard_Real F1, F2;
734     Value(X, F1);
735     X(1) = UMax;
736     Value(X, F2);
737     Standard_Real DU = Abs((F2 - F1) / (UMax - UMin));
738     X(1) = (UMin + UMax) / 2.;
739     X(2) = VMin;
740     Value(X, F1);
741     X(2) = VMax;
742     Value(X, F2);
743     Standard_Real DV = Abs((F2 - F1) / (VMax - VMin));
744     myPenalty = 10. * Max(DU, DV);
745     myPenalty = Max(myPenalty, 1.);
746   }
747
748   Standard_Boolean Value (const math_Vector& X,
749                                   Standard_Real& F)
750   {
751     if (CheckInputData(X))
752     {
753       gp_Pnt aP = mySurf.Value(X(1), X(2));
754       F = mySign * aP.Coord(myCoordIndx);
755     }
756     else
757     {
758       Standard_Real UPen = 0., VPen = 0., u0, v0;
759       if(X(1) < myUMin)
760       {
761         UPen = myPenalty * (myUMin - X(1));
762         u0 = myUMin;
763       }
764       else if(X(1) > myUMax)
765       {
766         UPen = myPenalty * (X(1) - myUMax);
767         u0 = myUMax;
768       }
769       else
770       {
771         u0 = X(1);
772       }
773       //
774       if(X(2) < myVMin)
775       {
776         VPen = myPenalty * (myVMin - X(2));
777         v0 = myVMin;
778       }
779       else if(X(2) > myVMax)
780       {
781         VPen = myPenalty * (X(2) - myVMax);
782         v0 = myVMax;
783       }
784       else
785       {
786         v0 = X(2);
787       }
788       //
789       gp_Pnt aP = mySurf.Value(u0, v0);
790       F = mySign * aP.Coord(myCoordIndx) + UPen + VPen;
791     }
792
793     return Standard_True;
794   }
795
796   
797
798   Standard_Integer NbVariables() const
799   {
800     return 2;
801   }
802
803 private:
804   SurfMaxMinCoord & operator = (const SurfMaxMinCoord & theOther);
805
806   Standard_Boolean CheckInputData(const math_Vector theParams)
807   {
808     if (theParams(1) < myUMin || 
809         theParams(1) > myUMax || 
810         theParams(2) < myVMin || 
811         theParams(2) > myVMax)
812       return Standard_False;
813     return Standard_True;
814   }
815
816   const Adaptor3d_Surface& mySurf;
817   Standard_Real myUMin;
818   Standard_Real myUMax;
819   Standard_Real myVMin;
820   Standard_Real myVMax;
821   Standard_Integer myCoordIndx;
822   Standard_Real mySign;
823   Standard_Real myPenalty;
824 };
825
826 //=======================================================================
827 //function : AdjustExtr
828 //purpose  : 
829 //=======================================================================
830
831 Standard_Real AdjustExtr(const Adaptor3d_Surface& S, 
832                          const Standard_Real UMin,
833                          const Standard_Real UMax,
834                          const Standard_Real VMin,
835                          const Standard_Real VMax,
836                          const Standard_Real Extr0,
837                          const Standard_Integer CoordIndx,                                 
838                          const Standard_Real Tol, 
839                          const Standard_Boolean IsMin)
840 {
841   Standard_Real aSign = IsMin ? 1.:-1.;
842   Standard_Real extr = aSign * Extr0;
843   Standard_Real relTol = 2.*Tol;
844   if(Abs(extr) > Tol)
845   {
846     relTol /= Abs(extr);
847   }
848   Standard_Real Du = (S.LastUParameter() - S.FirstUParameter());
849   Standard_Real Dv = (S.LastVParameter() - S.FirstVParameter());
850   //
851   math_Vector aT(1,2);
852   math_Vector aLowBorder(1,2);
853   math_Vector aUppBorder(1,2);
854   math_Vector aSteps(1,2);
855   aLowBorder(1) = UMin;
856   aUppBorder(1) = UMax;
857   aLowBorder(2) = VMin;
858   aUppBorder(2) = VMax;
859
860   Standard_Integer aNbU = Max(8, RealToInt(32 * (UMax - UMin) / Du));
861   Standard_Integer aNbV = Max(8, RealToInt(32 * (VMax - VMin) / Dv));
862   Standard_Integer aNbParticles = aNbU * aNbV;
863   Standard_Real aMaxUStep = (UMax - UMin) / (aNbU + 1);
864   aSteps(1) = Min(0.1 * Du, aMaxUStep);  
865   Standard_Real aMaxVStep = (VMax - VMin) / (aNbV + 1);
866   aSteps(2) = Min(0.1 * Dv, aMaxVStep);
867  
868   SurfMaxMinCoord aFunc(S, UMin, UMax, VMin, VMax, CoordIndx, aSign);
869   math_PSO aFinder(&aFunc, aLowBorder, aUppBorder, aSteps, aNbParticles); 
870   aFinder.Perform(aSteps, extr, aT);
871   
872   //Refinement of extremal value
873   math_Matrix aDir(1, 2, 1, 2, 0.0);
874   aDir(1, 1) = 1.;
875   aDir(2, 1) = 0.;
876   aDir(1, 2) = 0.;
877   aDir(2, 2) = 1.;
878
879   Standard_Integer aNbIter = 200;
880   math_Powell powell(aFunc, relTol, aNbIter, Tol);
881   powell.Perform(aFunc, aT, aDir);
882
883   if (powell.IsDone())
884   {
885     powell.Location(aT);
886     extr = powell.Minimum();
887   }
888
889   return aSign * extr;
890 }
891
892 //=======================================================================
893 //function : NbUSamples
894 //purpose  : 
895 //=======================================================================
896
897 Standard_Integer NbUSamples(const Adaptor3d_Surface& S, 
898                             const Standard_Real Umin,
899                             const Standard_Real Umax) 
900 {
901   Standard_Integer N;
902   GeomAbs_SurfaceType Type = S.GetType();
903   switch (Type) {
904   case GeomAbs_BezierSurface: 
905     {
906       N = 2*S.NbUPoles();
907       //By default parametric range of Bezier surf is [0, 1] [0, 1]
908       Standard_Real du = Umax - Umin;
909       if(du < .9)
910       {
911         N = RealToInt(du*N) + 1;
912         N = Max(N, 5);
913       }
914       break;
915     }
916   case GeomAbs_BSplineSurface: 
917     {
918       const Handle(Geom_BSplineSurface)& BS = S.BSpline();
919       N = 2*(BS->UDegree() + 1)*(BS->NbUKnots() -1);
920       Standard_Real umin, umax, vmin, vmax;
921       BS->Bounds(umin, umax, vmin, vmax);
922       Standard_Real du = (Umax - Umin) / (umax - umin);
923       if(du < .9)
924       {
925         N = RealToInt(du*N) + 1;
926         N = Max(N, 5);
927       }
928       break;
929     }
930   default:
931     N = 33;
932   }
933   return Min (50,N);
934 }
935
936 //=======================================================================
937 //function : NbVSamples
938 //purpose  : 
939 //=======================================================================
940
941 Standard_Integer NbVSamples(const Adaptor3d_Surface& S, 
942                             const Standard_Real Vmin,
943                             const Standard_Real Vmax) 
944 {
945   Standard_Integer N;
946   GeomAbs_SurfaceType Type = S.GetType();
947   switch (Type) {
948   case GeomAbs_BezierSurface:
949     {
950       N = 2*S.NbVPoles();
951       //By default parametric range of Bezier surf is [0, 1] [0, 1]
952       Standard_Real dv = Vmax - Vmin;
953       if(dv < .9)
954       {
955         N = RealToInt(dv*N) + 1;
956         N = Max(N, 5);
957       }
958       break;
959     }
960   case GeomAbs_BSplineSurface:
961     {
962       const Handle(Geom_BSplineSurface)& BS = S.BSpline();
963       N = 2*(BS->VDegree() + 1)*(BS->NbVKnots() - 1) ;
964      Standard_Real umin, umax, vmin, vmax;
965       BS->Bounds(umin, umax, vmin, vmax);
966       Standard_Real dv = (Vmax - Vmin) / (vmax - vmin);
967       if(dv < .9)
968       {
969         N = RealToInt(dv*N) + 1;
970         N = Max(N, 5);
971       }
972       break;
973     }
974   default:
975     N = 33;
976   }
977   return Min(50,N);
978 }
979