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         Bsaux->Segment(u1, u2);
 230         Bs = Bsaux;
 231       }
 232       //OCC566(apo)->
 233       Bnd_Box B1;
 234       Standard_Integer k, k1 = Bs->FirstUKnotIndex(), k2 = Bs->LastUKnotIndex(),
 235                        N = Bs->Degree(), NbKnots = Bs->NbKnots();
 236       TColStd_Array1OfReal Knots(1,NbKnots);
 237       Bs->Knots(Knots);
 238       GeomAdaptor_Curve GACurve(Bs);
 239       Standard_Real first = Knots(k1), last;
 240       for(k = k1 + 1; k <= k2; k++){
 241         last = Knots(k);
 242         tol = Max(FillBox(B1,GACurve,first,last,N), tol);
 243         first = last;
 244       }
 245       if (!B1.IsVoid())
 246       {
 247         B1.Enlarge(weakness*tol);
 248         reduceSplineBox(C, B1, B);
 249         B.Enlarge(Tol);
 250       }
 251       //<-OCC566(apo)
 252       break;
 253     }
 254   default:
 255     {
 256       Bnd_Box B1;
 257       static Standard_Integer N = 33;
 258       tol = FillBox(B1,C,U1,U2,N);
 259       B1.Enlarge(weakness*tol);
 260       Standard_Real x, y, z, X, Y, Z;
 261       B1.Get(x, y, z, X, Y, Z);
 262       B.Update(x, y, z, X, Y, Z);
 263       B.Enlarge(Tol);
 264     }
 265   }
 266 }
 267
 268 //=======================================================================
 269 //function : AddOptimal
 270 //purpose  :
 271 //=======================================================================
 272
 273 void BndLib_Add3dCurve::AddOptimal( const Adaptor3d_Curve& C,
 274                                                       const Standard_Real Tol,
 275                                                       Bnd_Box&      B )
 276 {
 277   BndLib_Add3dCurve::AddOptimal(C,
 278                                                   C.FirstParameter(),
 279                                                   C.LastParameter (),
 280                                                   Tol,B);
 281 }
 282
 283 //=======================================================================
 284 //function : AddOptimal
 285 //purpose  :
 286 //=======================================================================
 287
 288 void BndLib_Add3dCurve::AddOptimal( const Adaptor3d_Curve& C,
 289                                                       const Standard_Real U1,
 290                                                       const Standard_Real U2,
 291                                                       const Standard_Real Tol,
 292                                                       Bnd_Box&            B)
 293 {
 294   switch (C.GetType()) {
 295
 296     case GeomAbs_Line:
 297     {
 298       BndLib::Add(C.Line(),U1,U2,Tol,B);
 299       break;
 300     }
 301     case GeomAbs_Circle:
 302     {
 303       BndLib::Add(C.Circle(),U1,U2,Tol,B);
 304       break;
 305     }
 306     case GeomAbs_Ellipse:
 307     {
 308       BndLib::Add(C.Ellipse(),U1,U2,Tol,B);
 309       break;
 310     }
 311     case GeomAbs_Hyperbola:
 312     {
 313       BndLib::Add(C.Hyperbola(),U1,U2,Tol,B);
 314       break;
 315     }
 316     case GeomAbs_Parabola:
 317     {
 318       BndLib::Add(C.Parabola(),U1,U2,Tol,B);
 319       break;
 320     }
 321     default:
 322     {
 323       AddGenCurv(C, U1, U2, Tol, B);
 324     }
 325   }
 326 }
 327
 328 //=======================================================================
 329 //function : AddGenCurv
 330 //purpose  :
 331 //=======================================================================
 332 void BndLib_Add3dCurve::AddGenCurv(const Adaptor3d_Curve& C,
 333                                    const Standard_Real UMin,
 334                                    const Standard_Real UMax,
 335                                    const Standard_Real Tol,
 336                                    Bnd_Box& B)
 337 {
 338   Standard_Integer Nu = NbSamples(C, UMin, UMax);
 339   //
 340   Standard_Real CoordMin[3] = {RealLast(), RealLast(), RealLast()};
 341   Standard_Real CoordMax[3] = {-RealLast(), -RealLast(), -RealLast()};
 342   Standard_Real DeflMax[3] = {-RealLast(), -RealLast(), -RealLast()};
 343   //
 344   gp_Pnt P;
 345   Standard_Integer i, k;
 346   Standard_Real du = (UMax-UMin)/(Nu-1), du2 = du / 2.;
 347   NCollection_Array1<gp_XYZ> aPnts(1, Nu);
 348   Standard_Real u;
 349   for (i = 1, u = UMin; i <= Nu; i++, u += du)
 350   {
 351     C.D0(u,P);
 352     aPnts(i) = P.XYZ();
 353     //
 354     for(k = 0; k < 3; ++k)
 355     {
 356       if(CoordMin[k] > P.Coord(k+1))
 357       {
 358         CoordMin[k] = P.Coord(k+1);
 359       }
 360       if(CoordMax[k] < P.Coord(k+1))
 361       {
 362         CoordMax[k] = P.Coord(k+1);
 363       }
 364     }
 365     //
 366     if(i > 1)
 367     {
 368       gp_XYZ aPm = 0.5 * (aPnts(i-1) + aPnts(i));
 369       C.D0(u - du2, P);
 370       gp_XYZ aD = (P.XYZ() - aPm);
 371       for(k = 0; k < 3; ++k)
 372       {
 373         if(CoordMin[k] > P.Coord(k+1))
 374         {
 375           CoordMin[k] = P.Coord(k+1);
 376         }
 377         if(CoordMax[k] < P.Coord(k+1))
 378         {
 379           CoordMax[k] = P.Coord(k+1);
 380         }
 381         Standard_Real d = Abs(aD.Coord(k+1));
 382         if(DeflMax[k] < d)
 383         {
 384           DeflMax[k] = d;
 385         }
 386       }
 387     }
 388   }
 389   //
 390   //Adjusting minmax
 391   Standard_Real eps = Max(Tol, Precision::Confusion());
 392   for(k = 0; k < 3; ++k)
 393   {
 394     Standard_Real d = DeflMax[k];
 395     if(d <= eps)
 396     {
 397       continue;
 398     }
 399     Standard_Real CMin = CoordMin[k];
 400     Standard_Real CMax = CoordMax[k];
 401     for(i = 1; i <= Nu; ++i)
 402     {
 403       if(aPnts(i).Coord(k+1) - CMin < d)
 404       {
 405         Standard_Real umin, umax;
 406         umin = UMin + Max(0, i-2) * du;
 407         umax = UMin + Min(Nu-1, i) * du;
 408         Standard_Real cmin = AdjustExtr(C, umin, umax,
 409                                         CMin, k + 1, eps, Standard_True);
 410         if(cmin < CMin)
 411         {
 412           CMin = cmin;
 413         }
 414       }
 415       else if(CMax - aPnts(i).Coord(k+1) < d)
 416       {
 417         Standard_Real umin, umax;
 418         umin = UMin + Max(0, i-2) * du;
 419         umax = UMin + Min(Nu-1, i) * du;
 420         Standard_Real cmax = AdjustExtr(C, umin, umax,
 421                                         CMax, k + 1, eps, Standard_False);
 422         if(cmax > CMax)
 423         {
 424           CMax = cmax;
 425         }
 426       }
 427     }
 428     CoordMin[k] = CMin;
 429     CoordMax[k] = CMax;
 430   }
 431
 432   B.Add(gp_Pnt(CoordMin[0], CoordMin[1], CoordMin[2]));
 433   B.Add(gp_Pnt(CoordMax[0], CoordMax[1], CoordMax[2]));
 434   B.Enlarge(eps);
 435 }
 436 //
 437 class CurvMaxMinCoordMVar : public math_MultipleVarFunction
 438 {
 439 public:
 440   CurvMaxMinCoordMVar(const Adaptor3d_Curve& theCurve,
 441                       const Standard_Real UMin,
 442                       const Standard_Real UMax,
 443                       const Standard_Integer CoordIndx,
 444                       const Standard_Real Sign)
 445 : myCurve(theCurve),
 446   myUMin(UMin),
 447   myUMax(UMax),
 448   myCoordIndx(CoordIndx),
 449   mySign(Sign)
 450   {
 451   }
 452
 453   Standard_Boolean Value (const math_Vector& X,
 454                                 Standard_Real& F)
 455   {
 456     if (!CheckInputData(X(1)))
 457     {
 458       return Standard_False;
 459     }
 460     gp_Pnt aP = myCurve.Value(X(1));
 461
 462     F = mySign * aP.Coord(myCoordIndx);
 463
 464     return Standard_True;
 465   }
 466
 467
 468
 469   Standard_Integer NbVariables() const
 470   {
 471     return 1;
 472   }
 473
 474 private:
 475   CurvMaxMinCoordMVar & operator = (const CurvMaxMinCoordMVar & theOther);
 476
 477   Standard_Boolean CheckInputData(Standard_Real theParam)
 478   {
 479     if (theParam < myUMin ||
 480         theParam > myUMax)
 481       return Standard_False;
 482     return Standard_True;
 483   }
 484
 485   const Adaptor3d_Curve& myCurve;
 486   Standard_Real myUMin;
 487   Standard_Real myUMax;
 488   Standard_Integer myCoordIndx;
 489   Standard_Real mySign;
 490 };
 491 //
 492 class CurvMaxMinCoord : public math_Function
 493 {
 494 public:
 495   CurvMaxMinCoord(const Adaptor3d_Curve& theCurve,
 496                   const Standard_Real UMin,
 497                   const Standard_Real UMax,
 498                   const Standard_Integer CoordIndx,
 499                   const Standard_Real Sign)
 500 : myCurve(theCurve),
 501   myUMin(UMin),
 502   myUMax(UMax),
 503   myCoordIndx(CoordIndx),
 504   mySign(Sign)
 505   {
 506   }
 507
 508   Standard_Boolean Value (const Standard_Real X,
 509                                 Standard_Real& F)
 510   {
 511     if (!CheckInputData(X))
 512     {
 513       return Standard_False;
 514     }
 515     gp_Pnt aP = myCurve.Value(X);
 516
 517     F = mySign * aP.Coord(myCoordIndx);
 518
 519     return Standard_True;
 520   }
 521
 522 private:
 523   CurvMaxMinCoord & operator = (const CurvMaxMinCoord & theOther);
 524
 525   Standard_Boolean CheckInputData(Standard_Real theParam)
 526   {
 527     if (theParam < myUMin ||
 528         theParam > myUMax)
 529       return Standard_False;
 530     return Standard_True;
 531   }
 532
 533   const Adaptor3d_Curve& myCurve;
 534   Standard_Real myUMin;
 535   Standard_Real myUMax;
 536   Standard_Integer myCoordIndx;
 537   Standard_Real mySign;
 538 };
 539
 540 //=======================================================================
 541 //function : AdjustExtr
 542 //purpose  :
 543 //=======================================================================
 544
 545 Standard_Real AdjustExtr(const Adaptor3d_Curve& C,
 546                          const Standard_Real UMin,
 547                          const Standard_Real UMax,
 548                          const Standard_Real Extr0,
 549                          const Standard_Integer CoordIndx,
 550                          const Standard_Real Tol,
 551                          const Standard_Boolean IsMin)
 552 {
 553   Standard_Real aSign = IsMin ? 1.:-1.;
 554   Standard_Real extr = aSign * Extr0;
 555   //
 556   Standard_Real uTol = Max(C.Resolution(Tol), Precision::PConfusion());
 557   Standard_Real Du = (C.LastParameter() - C.FirstParameter());
 558   //
 559   Standard_Real reltol = uTol / Max(Abs(UMin), Abs(UMax));
 560   if(UMax - UMin < 0.01 * Du)
 561   {
 562
 563     math_BrentMinimum anOptLoc(reltol, 100, uTol);
 564     CurvMaxMinCoord aFunc(C, UMin, UMax, CoordIndx, aSign);
 565     anOptLoc.Perform(aFunc, UMin, (UMin+UMax)/2., UMax);
 566     if(anOptLoc.IsDone())
 567     {
 568       extr = anOptLoc.Minimum();
 569       return aSign * extr;
 570     }
 571   }
 572   //
 573   Standard_Integer aNbParticles = Max(8, RealToInt(32 * (UMax - UMin) / Du));
 574   Standard_Real maxstep = (UMax - UMin) / (aNbParticles + 1);
 575   math_Vector aT(1,1);
 576   math_Vector aLowBorder(1,1);
 577   math_Vector aUppBorder(1,1);
 578   math_Vector aSteps(1,1);
 579   aLowBorder(1) = UMin;
 580   aUppBorder(1) = UMax;
 581   aSteps(1) = Min(0.1 * Du, maxstep);
 582
 583   CurvMaxMinCoordMVar aFunc(C, UMin, UMax, CoordIndx, aSign);
 584   math_PSO aFinder(&aFunc, aLowBorder, aUppBorder, aSteps, aNbParticles);
 585   aFinder.Perform(aSteps, extr, aT);
 586   //
 587   math_BrentMinimum anOptLoc(reltol, 100, uTol);
 588   CurvMaxMinCoord aFunc1(C, UMin, UMax, CoordIndx, aSign);
 589   anOptLoc.Perform(aFunc1, Max(aT(1) - aSteps(1), UMin), aT(1), Min(aT(1) + aSteps(1), UMax));
 590
 591   if(anOptLoc.IsDone())
 592   {
 593     extr = anOptLoc.Minimum();
 594     return aSign * extr;
 595   }
 596
 597   return aSign * extr;
 598 }
 599
 600 //=======================================================================
 601 //function : NbSamples
 602 //purpose  :
 603 //=======================================================================
 604
 605 Standard_Integer NbSamples(const Adaptor3d_Curve& C,
 606                            const Standard_Real Umin,
 607                            const Standard_Real Umax)
 608 {
 609   Standard_Integer N;
 610   GeomAbs_CurveType Type = C.GetType();
 611   switch (Type) {
 612   case GeomAbs_BezierCurve:
 613     {
 614       N = 2 * C.NbPoles();
 615       //By default parametric range of Bezier curv is [0, 1]
 616       Standard_Real du = Umax - Umin;
 617       if(du < .9)
 618       {
 619         N = RealToInt(du*N) + 1;
 620         N = Max(N, 5);
 621       }
 622       break;
 623     }
 624   case GeomAbs_BSplineCurve:
 625     {
 626       const Handle(Geom_BSplineCurve)& BC = C.BSpline();
 627       N = 2 * (BC->Degree() + 1)*(BC->NbKnots() -1);
 628       Standard_Real umin = BC->FirstParameter(),
 629                     umax = BC->LastParameter();
 630       Standard_Real du = (Umax - Umin) / (umax - umin);
 631       if(du < .9)
 632       {
 633         N = RealToInt(du*N) + 1;
 634         N = Max(N, 5);
 635       }
 636       break;
 637     }
 638   default:
 639     N = 33;
 640   }
 641   return Min(500, N);
 642 }