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