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