0029151: GCC 7.1 warnings "this statement may fall through" [-Wimplicit-fallthrough=]
[occt.git] / src / BndLib / BndLib_AddSurface.cxx
CommitLineData
b311480e 1// Created on: 1995-07-24
2// Created by: Modelistation
3// Copyright (c) 1995-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 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
973c2be1 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.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
17// Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503
18
7fd59977 19#include <Adaptor3d_HSurface.hxx>
42cf5bc1 20#include <Adaptor3d_Surface.hxx>
21#include <Bnd_Box.hxx>
7fd59977 22#include <BndLib.hxx>
42cf5bc1 23#include <BndLib_AddSurface.hxx>
9bf0740b 24#include <BSplCLib.hxx>
7fd59977 25#include <ElSLib.hxx>
3ba87fdb 26#include <ElCLib.hxx>
7fd59977 27#include <Geom_BezierSurface.hxx>
42cf5bc1 28#include <Geom_BSplineSurface.hxx>
29#include <GeomAbs_SurfaceType.hxx>
30#include <gp_Pln.hxx>
31#include <gp_Pnt.hxx>
3ba87fdb 32#include <gp_Cylinder.hxx>
33#include <gp_Cone.hxx>
34#include <gp_Lin.hxx>
7fd59977 35#include <Precision.hxx>
42cf5bc1 36#include <TColgp_Array2OfPnt.hxx>
37#include <TColStd_Array1OfInteger.hxx>
38#include <TColStd_Array1OfReal.hxx>
3ba87fdb 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//
45static Standard_Integer NbUSamples(const Adaptor3d_Surface& S,
46 const Standard_Real Umin,
47 const Standard_Real Umax);
48//
49static Standard_Integer NbVSamples(const Adaptor3d_Surface& S,
50 const Standard_Real Vmin,
51 const Standard_Real Vmax);
52//
53static 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);
7fd59977 62
63//=======================================================================
64//function : Add
65//purpose :
66//=======================================================================
7fd59977 67void 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
83static 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
110static 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
133static 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
168static 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();
7fd59977 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}
7fd59977 208
9bf0740b 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).
215void 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{
447c7e54 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
9bf0740b 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
7fd59977 247//=======================================================================
248//function : Add
249//purpose :
250//=======================================================================
7fd59977 251void 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) {
9bf0740b 286
7fd59977 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() &&
9bf0740b 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);
7fd59977 319 else
9bf0740b 320 BndLib::Add(S.Sphere(),UMin,UMax,VMin,VMax,Tol,B);
7fd59977 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 {
9bf0740b 334 Standard_Boolean isUseConvexHullAlgorithm = Standard_True;
7fd59977 335 Standard_Real PTol = Precision::Parametric(Precision::Confusion());
9bf0740b 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;
7fd59977 461 }
462 }
b1811c1d 463 Standard_FALLTHROUGH
7fd59977 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++){
9bf0740b 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 }
7fd59977 476 }
477 B.Enlarge(Tol);
478 }
479 }
480}
3ba87fdb 481//----- Methods for AddOptimal ---------------------------------------
482
483//=======================================================================
484//function : AddOptimal
485//purpose :
486//=======================================================================
487void 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
503void 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//=======================================================================
568void 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//
712class SurfMaxMinCoord : public math_MultipleVarFunction
713{
714public:
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
803private:
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
831Standard_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
897Standard_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
941Standard_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