0031642: Visualization - crash in Graphic3d_Structure::SetVisual() on redisplaying...
[occt.git] / src / BndLib / BndLib_Add3dCurve.cxx
1 // Copyright (c) 1996-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15
16 #include <Adaptor3d_Curve.hxx>
17 #include <Bnd_Box.hxx>
18 #include <BndLib.hxx>
19 #include <BndLib_Add3dCurve.hxx>
20 #include <ElCLib.hxx>
21 #include <Geom_BezierCurve.hxx>
22 #include <Geom_BSplineCurve.hxx>
23 #include <GeomAbs_CurveType.hxx>
24 #include <GeomAdaptor_Curve.hxx>
25 #include <gp_Circ.hxx>
26 #include <gp_Elips.hxx>
27 #include <gp_Hypr.hxx>
28 #include <gp_Lin.hxx>
29 #include <gp_Parab.hxx>
30 #include <gp_Pnt.hxx>
31 #include <Precision.hxx>
32 #include <TColgp_Array1OfPnt.hxx>
33 #include <TColStd_Array1OfInteger.hxx>
34 #include <TColStd_Array1OfReal.hxx>
35 #include <math_MultipleVarFunction.hxx>
36 #include <math_Function.hxx>
37 #include <math_PSO.hxx>
38 #include <math_BrentMinimum.hxx>
39 //
40 static Standard_Integer NbSamples(const Adaptor3d_Curve& C, 
41                                    const Standard_Real Umin,
42                                    const Standard_Real Umax);
43
44 static Standard_Real  AdjustExtr(const Adaptor3d_Curve& C, 
45                                  const Standard_Real UMin,
46                                                    const Standard_Real UMax,
47                                  const Standard_Real Extr0,
48                                  const Standard_Integer CoordIndx,                                 
49                                  const Standard_Real Tol, 
50                                  const Standard_Boolean IsMin);
51
52
53 //=======================================================================
54 //function : reduceSplineBox
55 //purpose  : This method intended to reduce box in case of 
56 //           bezier and bspline curve.
57 //=======================================================================
58 static void reduceSplineBox(const Adaptor3d_Curve& theCurve,
59                             const Bnd_Box& theOrigBox,
60                             Bnd_Box & theReducedBox)
61 {
62   // Guaranteed bounding box based on poles of bspline.
63   Bnd_Box aPolesBox;
64   Standard_Real aPolesXMin, aPolesYMin, aPolesZMin,
65                 aPolesXMax, aPolesYMax, aPolesZMax;
66
67   if (theCurve.GetType() == GeomAbs_BSplineCurve)
68   {
69     Handle(Geom_BSplineCurve) aC = theCurve.BSpline();
70     const TColgp_Array1OfPnt& aPoles     = aC->Poles();
71
72     for(Standard_Integer anIdx  = aPoles.Lower();
73         anIdx <= aPoles.Upper();
74         ++anIdx)
75     {
76       aPolesBox.Add(aPoles.Value(anIdx));
77     }
78   }
79   if (theCurve.GetType() == GeomAbs_BezierCurve)
80   {
81     Handle(Geom_BezierCurve) aC = theCurve.Bezier();
82     const TColgp_Array1OfPnt& aPoles     = aC->Poles();
83
84     for(Standard_Integer anIdx  = aPoles.Lower();
85         anIdx <= aPoles.Upper();
86         ++anIdx)
87     {
88       aPolesBox.Add(aPoles.Value(anIdx));
89     }
90   }
91
92   aPolesBox.Get(aPolesXMin, aPolesYMin, aPolesZMin,
93                 aPolesXMax, aPolesYMax, aPolesZMax);
94
95   Standard_Real x, y, z, X, Y, Z;
96   theOrigBox.Get(x, y, z, X, Y, Z);
97
98   // Left bound.
99   if (aPolesXMin > x)
100     x = aPolesXMin;
101   if (aPolesYMin > y)
102     y = aPolesYMin;
103   if (aPolesZMin > z)
104     z = aPolesZMin;
105
106   // Right bound.
107   if (aPolesXMax < X)
108     X = aPolesXMax;
109   if (aPolesYMax < Y)
110     Y = aPolesYMax;
111   if (aPolesZMax < Z)
112     Z = aPolesZMax;
113
114   theReducedBox.Update(x, y, z, X, Y, Z);
115 }
116
117 //=======================================================================
118 //function : Add
119 //purpose  : 
120 //=======================================================================
121 void BndLib_Add3dCurve::Add( const Adaptor3d_Curve& C,
122                            const Standard_Real Tol,
123                                  Bnd_Box&      B )
124 {
125   BndLib_Add3dCurve::Add(C,
126                          C.FirstParameter(),
127                          C.LastParameter (),
128                          Tol,B);
129 }
130
131 //OCC566(apo)->
132 static Standard_Real FillBox(Bnd_Box& B, const Adaptor3d_Curve& C, 
133                              const Standard_Real first, const Standard_Real last, 
134                              const Standard_Integer N)
135 {
136   gp_Pnt P1, P2, P3;
137   C.D0(first,P1);  B.Add(P1);
138   Standard_Real p = first, dp = last-first, tol= 0.;
139   if(Abs(dp) > Precision::PConfusion()){
140     Standard_Integer i;
141     dp /= 2*N; 
142     for(i = 1; i <= N; i++){
143       p += dp;  C.D0(p,P2);  B.Add(P2);
144       p += dp;  C.D0(p,P3);  B.Add(P3);
145       gp_Pnt Pc((P1.XYZ()+P3.XYZ())/2.0);
146       tol = Max(tol,Pc.Distance(P2));
147       P1 = P3;
148     }
149   }else{
150     C.D0(first,P1);  B.Add(P1);
151     C.D0(last,P3);  B.Add(P3);
152     tol = 0.;
153   }
154   return tol;
155 }
156 //<-OCC566(apo)
157 //=======================================================================
158 //function : Add
159 //purpose  : 
160 //=======================================================================
161
162 void BndLib_Add3dCurve::Add( const Adaptor3d_Curve& C,
163                            const Standard_Real U1,
164                            const Standard_Real U2,
165                            const Standard_Real Tol,
166                                  Bnd_Box&      B )
167 {
168   static Standard_Real weakness = 1.5;  //OCC566(apo)
169   Standard_Real tol = 0.0;
170   switch (C.GetType()) {
171
172   case GeomAbs_Line: 
173     {
174       BndLib::Add(C.Line(),U1,U2,Tol,B);
175       break;
176     }
177   case GeomAbs_Circle: 
178     {
179       BndLib::Add(C.Circle(),U1,U2,Tol,B);
180       break;
181     }
182   case GeomAbs_Ellipse: 
183     {
184       BndLib::Add(C.Ellipse(),U1,U2,Tol,B);
185       break;
186     }
187   case GeomAbs_Hyperbola: 
188     {
189       BndLib::Add(C.Hyperbola(),U1,U2,Tol,B);
190       break;
191     }
192   case GeomAbs_Parabola: 
193     {
194       BndLib::Add(C.Parabola(),U1,U2,Tol,B);
195       break;
196     }
197   case GeomAbs_BezierCurve: 
198     {
199       Handle(Geom_BezierCurve) Bz = C.Bezier();
200       Standard_Integer N = Bz->Degree();
201       GeomAdaptor_Curve GACurve(Bz);
202       Bnd_Box B1;
203       tol = FillBox(B1,GACurve,U1,U2,N);
204       B1.Enlarge(weakness*tol);
205       reduceSplineBox(C, B1, B);
206       B.Enlarge(Tol);
207       break;
208     }
209   case GeomAbs_BSplineCurve: 
210     {
211       Handle(Geom_BSplineCurve) Bs = C.BSpline();
212       if(Abs(Bs->FirstParameter() - U1) > Precision::Parametric(Tol)||
213          Abs(Bs->LastParameter()  - U2) > Precision::Parametric(Tol)) {
214
215         Handle(Geom_Geometry) G = Bs->Copy();
216         Handle(Geom_BSplineCurve) Bsaux (Handle(Geom_BSplineCurve)::DownCast (G));
217         Standard_Real u1 = U1, u2 = U2;
218         //// modified by jgv, 24.10.01 for BUC61031 ////
219         if (Bsaux->IsPeriodic())
220           ElCLib::AdjustPeriodic( Bsaux->FirstParameter(), Bsaux->LastParameter(), Precision::PConfusion(), u1, u2 );
221         else {
222           ////////////////////////////////////////////////
223           //  modified by NIZHNY-EAP Fri Dec  3 14:29:14 1999 ___BEGIN___
224           // To avoid exeption in Segment
225           if(Bsaux->FirstParameter() > U1) u1 = Bsaux->FirstParameter();
226           if(Bsaux->LastParameter()  < U2 ) u2  = Bsaux->LastParameter();
227           //  modified by NIZHNY-EAP Fri Dec  3 14:29:18 1999 ___END___
228         }
229         Standard_Real aSegmentTol = Precision::PConfusion();
230         if (Abs(u2 - u1) < aSegmentTol)
231           aSegmentTol = Abs(u2 - u1) * 0.01;
232         Bsaux->Segment(u1, u2, aSegmentTol);
233         Bs = Bsaux;
234       }
235       //OCC566(apo)->
236       Bnd_Box B1;
237       Standard_Integer k, k1 = Bs->FirstUKnotIndex(), k2 = Bs->LastUKnotIndex(),
238                        N = Bs->Degree(), NbKnots = Bs->NbKnots();
239       TColStd_Array1OfReal Knots(1,NbKnots);
240       Bs->Knots(Knots);
241       GeomAdaptor_Curve GACurve(Bs);
242       Standard_Real first = Knots(k1), last;
243       for(k = k1 + 1; k <= k2; k++){
244         last = Knots(k); 
245         tol = Max(FillBox(B1,GACurve,first,last,N), tol);
246         first = last;
247       }
248       if (!B1.IsVoid())
249       {
250         B1.Enlarge(weakness*tol);
251         reduceSplineBox(C, B1, B);
252         B.Enlarge(Tol);
253       }
254       //<-OCC566(apo)
255       break;
256     }
257   default:
258     {
259       Bnd_Box B1;
260       static Standard_Integer N = 33;
261       tol = FillBox(B1,C,U1,U2,N);
262       B1.Enlarge(weakness*tol);
263       Standard_Real x, y, z, X, Y, Z;
264       B1.Get(x, y, z, X, Y, Z);
265       B.Update(x, y, z, X, Y, Z);
266       B.Enlarge(Tol);
267     }
268   }
269 }
270
271 //=======================================================================
272 //function : AddOptimal
273 //purpose  : 
274 //=======================================================================
275
276 void BndLib_Add3dCurve::AddOptimal( const Adaptor3d_Curve& C,
277                                                       const Standard_Real Tol,
278                                                       Bnd_Box&      B )
279 {
280   BndLib_Add3dCurve::AddOptimal(C,
281                                                   C.FirstParameter(),
282                                                   C.LastParameter (),
283                                                   Tol,B);
284 }
285
286 //=======================================================================
287 //function : AddOptimal
288 //purpose  : 
289 //=======================================================================
290
291 void BndLib_Add3dCurve::AddOptimal( const Adaptor3d_Curve& C,
292                                                       const Standard_Real U1,
293                                                       const Standard_Real U2,
294                                                       const Standard_Real Tol,
295                                                       Bnd_Box&            B)
296 {
297   switch (C.GetType()) {
298
299     case GeomAbs_Line: 
300     {
301       BndLib::Add(C.Line(),U1,U2,Tol,B);
302       break;
303     }
304     case GeomAbs_Circle: 
305     {
306       BndLib::Add(C.Circle(),U1,U2,Tol,B);
307       break;
308     }
309     case GeomAbs_Ellipse: 
310     {
311       BndLib::Add(C.Ellipse(),U1,U2,Tol,B);
312       break;
313     }
314     case GeomAbs_Hyperbola: 
315     {
316       BndLib::Add(C.Hyperbola(),U1,U2,Tol,B);
317       break;
318     }
319     case GeomAbs_Parabola: 
320     {
321       BndLib::Add(C.Parabola(),U1,U2,Tol,B);
322       break;
323     }
324     default:
325     {
326       AddGenCurv(C, U1, U2, Tol, B);
327     }
328   }
329 }
330
331 //=======================================================================
332 //function : AddGenCurv
333 //purpose  : 
334 //=======================================================================
335 void BndLib_Add3dCurve::AddGenCurv(const Adaptor3d_Curve& C, 
336                                    const Standard_Real UMin,
337                                    const Standard_Real UMax,
338                                    const Standard_Real Tol,
339                                    Bnd_Box& B)
340 {
341   Standard_Integer Nu = NbSamples(C, UMin, UMax);
342   //
343   Standard_Real CoordMin[3] = {RealLast(), RealLast(), RealLast()}; 
344   Standard_Real CoordMax[3] = {-RealLast(), -RealLast(), -RealLast()};
345   Standard_Real DeflMax[3] = {-RealLast(), -RealLast(), -RealLast()};
346   //
347   gp_Pnt P;
348   Standard_Integer i, k;
349   Standard_Real du = (UMax-UMin)/(Nu-1), du2 = du / 2.;
350   NCollection_Array1<gp_XYZ> aPnts(1, Nu);
351   Standard_Real u;
352   for (i = 1, u = UMin; i <= Nu; i++, u += du)
353   {
354     C.D0(u,P);
355     aPnts(i) = P.XYZ();
356     //
357     for(k = 0; k < 3; ++k)
358     {
359       if(CoordMin[k] > P.Coord(k+1))
360       {
361         CoordMin[k] = P.Coord(k+1);
362       }
363       if(CoordMax[k] < P.Coord(k+1))
364       {
365         CoordMax[k] = P.Coord(k+1);
366       }
367     }
368     //
369     if(i > 1)
370     {
371       gp_XYZ aPm = 0.5 * (aPnts(i-1) + aPnts(i));
372       C.D0(u - du2, P);
373       gp_XYZ aD = (P.XYZ() - aPm);
374       for(k = 0; k < 3; ++k)
375       {
376         if(CoordMin[k] > P.Coord(k+1))
377         {
378           CoordMin[k] = P.Coord(k+1);
379         }
380         if(CoordMax[k] < P.Coord(k+1))
381         {
382           CoordMax[k] = P.Coord(k+1);
383         }
384         Standard_Real d = Abs(aD.Coord(k+1));
385         if(DeflMax[k] < d)
386         {
387           DeflMax[k] = d;
388         }
389       }
390     }
391   }
392   //
393   //Adjusting minmax 
394   Standard_Real eps = Max(Tol, Precision::Confusion());
395   for(k = 0; k < 3; ++k)
396   {
397     Standard_Real d = DeflMax[k];
398     if(d <= eps)
399     {
400       continue;
401     }
402     Standard_Real CMin = CoordMin[k];
403     Standard_Real CMax = CoordMax[k];
404     for(i = 1; i <= Nu; ++i)
405     {
406       if(aPnts(i).Coord(k+1) - CMin < d)
407       {
408         Standard_Real umin, umax;
409         umin = UMin + Max(0, i-2) * du;
410         umax = UMin + Min(Nu-1, i) * du;
411         Standard_Real cmin = AdjustExtr(C, umin, umax,
412                                         CMin, k + 1, eps, Standard_True);
413         if(cmin < CMin)
414         {
415           CMin = cmin;
416         }
417       }
418       else if(CMax - aPnts(i).Coord(k+1) < d)
419       {
420         Standard_Real umin, umax;
421         umin = UMin + Max(0, i-2) * du;
422         umax = UMin + Min(Nu-1, i) * du;
423         Standard_Real cmax = AdjustExtr(C, umin, umax,
424                                         CMax, k + 1, eps, Standard_False);
425         if(cmax > CMax)
426         {
427           CMax = cmax;
428         }
429       }
430     }
431     CoordMin[k] = CMin;
432     CoordMax[k] = CMax;
433   }
434
435   B.Add(gp_Pnt(CoordMin[0], CoordMin[1], CoordMin[2]));
436   B.Add(gp_Pnt(CoordMax[0], CoordMax[1], CoordMax[2]));
437   B.Enlarge(eps);
438 }
439 //
440 class CurvMaxMinCoordMVar : public math_MultipleVarFunction
441 {
442 public:
443   CurvMaxMinCoordMVar(const Adaptor3d_Curve& theCurve, 
444                       const Standard_Real UMin,
445                       const Standard_Real UMax,
446                       const Standard_Integer CoordIndx,                                 
447                       const Standard_Real Sign)
448 : myCurve(theCurve),
449   myUMin(UMin),
450   myUMax(UMax),
451   myCoordIndx(CoordIndx),
452   mySign(Sign)
453   {
454   }
455
456   Standard_Boolean Value (const math_Vector& X,
457                                 Standard_Real& F)
458   {
459     if (!CheckInputData(X(1)))
460     {
461       return Standard_False;
462     }
463     gp_Pnt aP = myCurve.Value(X(1));
464
465     F = mySign * aP.Coord(myCoordIndx);
466
467     return Standard_True;
468   }
469
470   
471
472   Standard_Integer NbVariables() const
473   {
474     return 1;
475   }
476
477 private:
478   CurvMaxMinCoordMVar & operator = (const CurvMaxMinCoordMVar & theOther);
479
480   Standard_Boolean CheckInputData(Standard_Real theParam)
481   {
482     if (theParam < myUMin || 
483         theParam > myUMax)
484       return Standard_False;
485     return Standard_True;
486   }
487
488   const Adaptor3d_Curve& myCurve;
489   Standard_Real myUMin;
490   Standard_Real myUMax;
491   Standard_Integer myCoordIndx;
492   Standard_Real mySign;
493 };
494 //
495 class CurvMaxMinCoord : public math_Function
496 {
497 public:
498   CurvMaxMinCoord(const Adaptor3d_Curve& theCurve, 
499                   const Standard_Real UMin,
500                   const Standard_Real UMax,
501                   const Standard_Integer CoordIndx,                                 
502                   const Standard_Real Sign)
503 : myCurve(theCurve),
504   myUMin(UMin),
505   myUMax(UMax),
506   myCoordIndx(CoordIndx),
507   mySign(Sign)
508   {
509   }
510
511   Standard_Boolean Value (const Standard_Real X,
512                                 Standard_Real& F)
513   {
514     if (!CheckInputData(X))
515     {
516       return Standard_False;
517     }
518     gp_Pnt aP = myCurve.Value(X);
519
520     F = mySign * aP.Coord(myCoordIndx);
521
522     return Standard_True;
523   }
524
525 private:
526   CurvMaxMinCoord & operator = (const CurvMaxMinCoord & theOther);
527
528   Standard_Boolean CheckInputData(Standard_Real theParam)
529   {
530     if (theParam < myUMin || 
531         theParam > myUMax)
532       return Standard_False;
533     return Standard_True;
534   }
535
536   const Adaptor3d_Curve& myCurve;
537   Standard_Real myUMin;
538   Standard_Real myUMax;
539   Standard_Integer myCoordIndx;
540   Standard_Real mySign;
541 };
542
543 //=======================================================================
544 //function : AdjustExtr
545 //purpose  : 
546 //=======================================================================
547
548 Standard_Real AdjustExtr(const Adaptor3d_Curve& C, 
549                          const Standard_Real UMin,
550                          const Standard_Real UMax,
551                          const Standard_Real Extr0,
552                          const Standard_Integer CoordIndx,                                 
553                          const Standard_Real Tol, 
554                          const Standard_Boolean IsMin)
555 {
556   Standard_Real aSign = IsMin ? 1.:-1.;
557   Standard_Real extr = aSign * Extr0;
558   //
559   Standard_Real uTol = Max(C.Resolution(Tol), Precision::PConfusion());
560   Standard_Real Du = (C.LastParameter() - C.FirstParameter());
561   //
562   Standard_Real reltol = uTol / Max(Abs(UMin), Abs(UMax));
563   if(UMax - UMin < 0.01 * Du)
564   {
565
566     math_BrentMinimum anOptLoc(reltol, 100, uTol);
567     CurvMaxMinCoord aFunc(C, UMin, UMax, CoordIndx, aSign);
568     anOptLoc.Perform(aFunc, UMin, (UMin+UMax)/2., UMax);
569     if(anOptLoc.IsDone())
570     {
571       extr = anOptLoc.Minimum();
572       return aSign * extr;
573     }
574   }
575   //
576   Standard_Integer aNbParticles = Max(8, RealToInt(32 * (UMax - UMin) / Du));
577   Standard_Real maxstep = (UMax - UMin) / (aNbParticles + 1);
578   math_Vector aT(1,1);
579   math_Vector aLowBorder(1,1);
580   math_Vector aUppBorder(1,1);
581   math_Vector aSteps(1,1);
582   aLowBorder(1) = UMin;
583   aUppBorder(1) = UMax;
584   aSteps(1) = Min(0.1 * Du, maxstep);
585
586   CurvMaxMinCoordMVar aFunc(C, UMin, UMax, CoordIndx, aSign);
587   math_PSO aFinder(&aFunc, aLowBorder, aUppBorder, aSteps, aNbParticles); 
588   aFinder.Perform(aSteps, extr, aT);
589   //
590   math_BrentMinimum anOptLoc(reltol, 100, uTol);
591   CurvMaxMinCoord aFunc1(C, UMin, UMax, CoordIndx, aSign);
592   anOptLoc.Perform(aFunc1, Max(aT(1) - aSteps(1), UMin), aT(1), Min(aT(1) + aSteps(1), UMax));
593
594   if(anOptLoc.IsDone())
595   {
596     extr = anOptLoc.Minimum();
597     return aSign * extr;
598   }
599
600   return aSign * extr;
601 }
602
603 //=======================================================================
604 //function : NbSamples
605 //purpose  : 
606 //=======================================================================
607
608 Standard_Integer NbSamples(const Adaptor3d_Curve& C, 
609                            const Standard_Real Umin,
610                            const Standard_Real Umax) 
611 {
612   Standard_Integer N;
613   GeomAbs_CurveType Type = C.GetType();
614   switch (Type) {
615   case GeomAbs_BezierCurve: 
616     {
617       N = 2 * C.NbPoles();
618       //By default parametric range of Bezier curv is [0, 1]
619       Standard_Real du = Umax - Umin;
620       if(du < .9)
621       {
622         N = RealToInt(du*N) + 1;
623         N = Max(N, 5);
624       }
625       break;
626     }
627   case GeomAbs_BSplineCurve: 
628     {
629       const Handle(Geom_BSplineCurve)& BC = C.BSpline();
630       N = 2 * (BC->Degree() + 1)*(BC->NbKnots() -1);
631       Standard_Real umin = BC->FirstParameter(), 
632                     umax = BC->LastParameter();
633       Standard_Real du = (Umax - Umin) / (umax - umin);
634       if(du < .9)
635       {
636         N = RealToInt(du*N) + 1;
637         N = Max(N, 5);
638       }
639       break;
640     }
641   default:
642     N = 33;
643   }
644   return Min(500, N);
645 }