src/Geom/Geom_OffsetCurve.cxx

   1 // Created on: 1991-06-25
   2 // Created by: JCV
   3 // Copyright (c) 1991-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   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
  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.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 // 24-Aug-95 : xab removed C1 and C2 test : appeller  D1 et D2
  18 //             avec discernement !
  19 // 19-09-97  : JPI correction derivee seconde
  20
  21
  22 #include <Geom_OffsetCurve.ixx>
  23
  24
  25 #include <gp.hxx>
  26 #include <gp_XYZ.hxx>
  27 #include <Geom_Line.hxx>
  28 #include <Geom_Circle.hxx>
  29 #include <Geom_Ellipse.hxx>
  30 #include <Geom_Hyperbola.hxx>
  31 #include <Geom_Parabola.hxx>
  32 #include <Geom_BezierCurve.hxx>
  33 #include <Geom_TrimmedCurve.hxx>
  34 #include <Geom_BSplineCurve.hxx>
  35 #include <Geom_Geometry.hxx>
  36
  37 #include <Geom_UndefinedDerivative.hxx>
  38 #include <Geom_UndefinedValue.hxx>
  39 #include <Standard_ConstructionError.hxx>
  40 #include <Standard_RangeError.hxx>
  41 #include <Standard_NotImplemented.hxx>
  42 #include <CSLib_Offset.hxx>
  43
  44 typedef Geom_OffsetCurve         OffsetCurve;
  45 typedef Handle(Geom_OffsetCurve) Handle(OffsetCurve);
  46 typedef Geom_Curve               Curve;
  47 typedef Handle(Geom_Curve)       Handle(Curve);
  48 typedef Handle(Geom_Geometry)    Handle(Geometry);
  49
  50 typedef gp_Dir  Dir;
  51 typedef gp_Pnt  Pnt;
  52 typedef gp_Trsf Trsf;
  53 typedef gp_Vec  Vec;
  54 typedef gp_XYZ  XYZ;
  55
  56
  57
  58
  59 //ordre de derivation maximum pour la recherche de la premiere
  60 //derivee non nulle
  61 static const int maxDerivOrder = 3;
  62 static const Standard_Real MinStep   = 1e-7;
  63 static const Standard_Real MyAngularToleranceForG1 = Precision::Angular();
  64
  65
  66 static gp_Vec dummyDerivative; // used as empty value for unused derivatives in AdjustDerivative
  67 // Recalculate derivatives in the singular point
  68 // Returns true if the direction of derivatives is changed
  69 static Standard_Boolean AdjustDerivative(
  70     const Handle(Geom_Curve)& theCurve, Standard_Integer theMaxDerivative, Standard_Real theU, gp_Vec& theD1,
  71     gp_Vec& theD2 = dummyDerivative, gp_Vec& theD3 = dummyDerivative, gp_Vec& theD4 = dummyDerivative);
  72
  73
  74
  75 //=======================================================================
  76 //function : Copy
  77 //purpose  :
  78 //=======================================================================
  79
  80 Handle(Geom_Geometry) Geom_OffsetCurve::Copy () const {
  81
  82  Handle(OffsetCurve) C;
  83  C = new OffsetCurve (basisCurve, offsetValue, direction);
  84  return C;
  85 }
  86
  87
  88
  89 //=======================================================================
  90 //function : Geom_OffsetCurve
  91 //purpose  : Basis curve cannot be an Offset curve or trimmed from
  92 //            offset curve.
  93 //=======================================================================
  94
  95 Geom_OffsetCurve::Geom_OffsetCurve (const Handle(Geom_Curve)& theCurve,
  96                                     const Standard_Real       theOffset,
  97                                     const gp_Dir&             theDir,
  98                                     const Standard_Boolean    isTheNotCheckC0)
  99  : direction(theDir), offsetValue(theOffset)
 100 {
 101   SetBasisCurve (theCurve, isTheNotCheckC0);
 102 }
 103
 104
 105 //=======================================================================
 106 //function : Reverse
 107 //purpose  :
 108 //=======================================================================
 109
 110 void Geom_OffsetCurve::Reverse ()
 111 {
 112   basisCurve->Reverse();
 113   offsetValue = -offsetValue;
 114 }
 115
 116
 117 //=======================================================================
 118 //function : ReversedParameter
 119 //purpose  :
 120 //=======================================================================
 121
 122 Standard_Real Geom_OffsetCurve::ReversedParameter( const Standard_Real U) const
 123 {
 124   return basisCurve->ReversedParameter( U);
 125 }
 126
 127 //=======================================================================
 128 //function : Direction
 129 //purpose  :
 130 //=======================================================================
 131
 132 const gp_Dir& Geom_OffsetCurve::Direction () const
 133   { return direction; }
 134
 135 //=======================================================================
 136 //function : SetDirection
 137 //purpose  :
 138 //=======================================================================
 139
 140 void Geom_OffsetCurve::SetDirection (const Dir& V)
 141   { direction = V; }
 142
 143 //=======================================================================
 144 //function : SetOffsetValue
 145 //purpose  :
 146 //=======================================================================
 147
 148 void Geom_OffsetCurve::SetOffsetValue (const Standard_Real D)
 149   { offsetValue = D; }
 150
 151
 152 //=======================================================================
 153 //function : IsPeriodic
 154 //purpose  :
 155 //=======================================================================
 156
 157 Standard_Boolean Geom_OffsetCurve::IsPeriodic () const
 158 {
 159   return basisCurve->IsPeriodic();
 160 }
 161
 162 //=======================================================================
 163 //function : Period
 164 //purpose  :
 165 //=======================================================================
 166
 167 Standard_Real Geom_OffsetCurve::Period () const
 168 {
 169   return basisCurve->Period();
 170 }
 171
 172 //=======================================================================
 173 //function : SetBasisCurve
 174 //purpose  :
 175 //=======================================================================
 176
 177 void Geom_OffsetCurve::SetBasisCurve (const Handle(Curve)& C,
 178                                       const Standard_Boolean isNotCheckC0)
 179 {
 180   const Standard_Real aUf = C->FirstParameter(),
 181                       aUl = C->LastParameter();
 182   Handle(Curve) aCheckingCurve =  Handle(Curve)::DownCast(C->Copy());
 183   Standard_Boolean isTrimmed = Standard_False;
 184
 185   while(aCheckingCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)) ||
 186         aCheckingCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve)))
 187   {
 188     if (aCheckingCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)))
 189     {
 190       Handle(Geom_TrimmedCurve) aTrimC =
 191                 Handle(Geom_TrimmedCurve)::DownCast(aCheckingCurve);
 192       aCheckingCurve = aTrimC->BasisCurve();
 193       isTrimmed = Standard_True;
 194     }
 195
 196     if (aCheckingCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve)))
 197     {
 198       Handle(Geom_OffsetCurve) aOC =
 199             Handle(Geom_OffsetCurve)::DownCast(aCheckingCurve);
 200       aCheckingCurve = aOC->BasisCurve();
 201       Standard_Real PrevOff = aOC->Offset();
 202       gp_Vec V1(aOC->Direction());
 203       gp_Vec V2(direction);
 204       gp_Vec Vdir(PrevOff*V1 + offsetValue*V2);
 205
 206       if (offsetValue >= 0.)
 207       {
 208         offsetValue = Vdir.Magnitude();
 209         direction.SetXYZ(Vdir.XYZ());
 210       }
 211       else
 212       {
 213         offsetValue = -Vdir.Magnitude();
 214         direction.SetXYZ((-Vdir).XYZ());
 215       }
 216     }
 217   }
 218
 219   myBasisCurveContinuity = aCheckingCurve->Continuity();
 220
 221   Standard_Boolean isC0 = !isNotCheckC0 &&
 222                           (myBasisCurveContinuity == GeomAbs_C0);
 223
 224   // Basis curve must be at least C1
 225   if (isC0 && aCheckingCurve->IsKind(STANDARD_TYPE(Geom_BSplineCurve)))
 226   {
 227     Handle(Geom_BSplineCurve) aBC = Handle(Geom_BSplineCurve)::DownCast(aCheckingCurve);
 228     if(aBC->IsG1(aUf, aUl, MyAngularToleranceForG1))
 229     {
 230       //Checking if basis curve has more smooth (C1, G2 and above) is not done.
 231       //It can be done in case of need.
 232       myBasisCurveContinuity = GeomAbs_G1;
 233       isC0 = Standard_False;
 234     }
 235
 236     // Raise exception if still C0
 237     if (isC0)
 238       Standard_ConstructionError::Raise("Offset on C0 curve");
 239   }
 240   //
 241   if(isTrimmed)
 242   {
 243     basisCurve = new Geom_TrimmedCurve(aCheckingCurve, aUf, aUl);
 244   }
 245   else
 246   {
 247     basisCurve = aCheckingCurve;
 248   }
 249 }
 250
 251
 252
 253 //=======================================================================
 254 //function : BasisCurve
 255 //purpose  :
 256 //=======================================================================
 257
 258 Handle(Curve) Geom_OffsetCurve::BasisCurve () const
 259 {
 260   return basisCurve;
 261 }
 262
 263
 264 //=======================================================================
 265 //function : Continuity
 266 //purpose  :
 267 //=======================================================================
 268
 269 GeomAbs_Shape Geom_OffsetCurve::Continuity () const {
 270
 271   GeomAbs_Shape OffsetShape=GeomAbs_C0;
 272   switch (myBasisCurveContinuity) {
 273     case GeomAbs_C0 : OffsetShape = GeomAbs_C0;       break;
 274     case GeomAbs_C1 : OffsetShape = GeomAbs_C0;       break;
 275     case GeomAbs_C2 : OffsetShape = GeomAbs_C1;       break;
 276     case GeomAbs_C3 : OffsetShape = GeomAbs_C2;       break;
 277     case GeomAbs_CN : OffsetShape = GeomAbs_CN;       break;
 278     case GeomAbs_G1 : OffsetShape = GeomAbs_G1;       break;
 279     case GeomAbs_G2 : OffsetShape = GeomAbs_G2;       break;
 280   }
 281   return OffsetShape;
 282 }
 283
 284
 285 //=======================================================================
 286 //function : D0
 287 //purpose  :
 288 //=======================================================================
 289
 290 void Geom_OffsetCurve::D0 (const Standard_Real U, Pnt& P) const
 291 {
 292   gp_Pnt PBasis;
 293   gp_Vec VBasis;
 294   D0(U,P,PBasis,VBasis);
 295 }
 296
 297
 298 //=======================================================================
 299 //function : D1
 300 //purpose  :
 301 //=======================================================================
 302
 303 void Geom_OffsetCurve::D1 (const Standard_Real U, Pnt& P, Vec& V1) const
 304 {
 305   gp_Pnt PBasis;
 306   gp_Vec V1Basis,V2Basis;
 307   D1(U,P,PBasis,V1,V1Basis,V2Basis);
 308 }
 309
 310
 311
 312 //=======================================================================
 313 //function : D2
 314 //purpose  :
 315 //=======================================================================
 316
 317 void Geom_OffsetCurve::D2 (const Standard_Real U, Pnt& P, Vec& V1, Vec& V2) const
 318 {
 319   gp_Pnt PBasis;
 320   gp_Vec V1Basis,V2Basis,V3Basis;
 321   D2(U,P,PBasis,V1,V2,V1Basis,V2Basis,V3Basis);
 322 }
 323
 324
 325 //=======================================================================
 326 //function : D3
 327 //purpose  :
 328 //=======================================================================
 329
 330 void Geom_OffsetCurve::D3 (const Standard_Real theU, Pnt& theP, Vec& theV1, Vec& theV2, Vec& theV3) const
 331 {
 332    // P(u) = p(u) + Offset * Ndir / R
 333    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
 334
 335    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 336
 337    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
 338    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 339
 340    //P"'(u) = p"'(u) + (Offset / R) * (D3Ndir - (3.0 * Dr/R**2) * D2Ndir -
 341    //         (3.0 * D2r / R2) * DNdir + (3.0 * Dr * Dr / R4) * DNdir -
 342    //         (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir +
 343    //         (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir
 344
 345   Standard_Boolean IsDirectionChange = Standard_False;
 346
 347   basisCurve->D3 (theU, theP, theV1, theV2, theV3);
 348   Vec aV4 = basisCurve->DN (theU, 4);
 349   if(theV1.SquareMagnitude() <= gp::Resolution())
 350     IsDirectionChange = AdjustDerivative(basisCurve, 4, theU, theV1, theV2, theV3, aV4);
 351
 352   CSLib_Offset::D3(theP, theV1, theV2, theV3, aV4, direction, offsetValue,
 353                    IsDirectionChange, theP, theV1, theV2, theV3);
 354 }
 355
 356
 357 //=======================================================================
 358 //function : DN
 359 //purpose  :
 360 //=======================================================================
 361
 362 Vec Geom_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const
 363   {
 364   Standard_RangeError_Raise_if (N < 1, "Exception: "
 365                               "Geom_OffsetCurve::DN(...). N<1.");
 366
 367   gp_Vec VN, Vtemp;
 368   gp_Pnt Ptemp;
 369   switch (N)
 370     {
 371     case 1:
 372       D1( U, Ptemp, VN);
 373       break;
 374     case 2:
 375       D2( U, Ptemp, Vtemp, VN);
 376       break;
 377     case 3:
 378       D3( U, Ptemp, Vtemp, Vtemp, VN);
 379       break;
 380     default:
 381       Standard_NotImplemented::Raise("Exception: "
 382         "Derivative order is greater than 3. Cannot compute of derivative.");
 383   }
 384
 385   return VN;
 386 }
 387
 388 //=======================================================================
 389 //function : D0
 390 //purpose  :
 391 //=======================================================================
 392
 393 void  Geom_OffsetCurve::D0(const Standard_Real theU, gp_Pnt& theP,
 394                            gp_Pnt& thePbasis, gp_Vec& theV1basis)const
 395 {
 396   basisCurve->D1(theU, thePbasis, theV1basis);
 397   Standard_Boolean IsDirectionChange = Standard_False;
 398   if(theV1basis.SquareMagnitude() <= gp::Resolution())
 399     IsDirectionChange = AdjustDerivative(basisCurve, 1, theU, theV1basis);
 400
 401   CSLib_Offset::D0(thePbasis, theV1basis, direction, offsetValue, IsDirectionChange, theP);
 402 }
 403
 404 //=======================================================================
 405 //function : D1
 406 //purpose  :
 407 //=======================================================================
 408
 409 void Geom_OffsetCurve::D1 ( const Standard_Real theU,
 410                             Pnt& theP , Pnt& thePBasis ,
 411                             Vec& theV1, Vec& theV1basis, Vec& theV2basis) const {
 412
 413    // P(u) = p(u) + Offset * Ndir / R
 414    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
 415
 416    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 417
 418   basisCurve->D2 (theU, thePBasis, theV1basis, theV2basis);
 419
 420   Standard_Boolean IsDirectionChange = Standard_False;
 421   if(theV1basis.SquareMagnitude() <= gp::Resolution())
 422     IsDirectionChange = AdjustDerivative(basisCurve, 2, theU, theV1basis, theV2basis);
 423
 424   CSLib_Offset::D1(thePBasis, theV1basis, theV2basis, direction, offsetValue, IsDirectionChange, theP, theV1);
 425 }
 426
 427
 428 //=======================================================================
 429 //function : D2
 430 //purpose  :
 431 //=======================================================================
 432
 433 void Geom_OffsetCurve::D2 (const Standard_Real theU,
 434                            Pnt& theP, Pnt& thePBasis,
 435                            Vec& theV1, Vec& theV2,
 436                            Vec& theV1basis, Vec& theV2basis, Vec& theV3basis) const
 437 {
 438    // P(u) = p(u) + Offset * Ndir / R
 439    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
 440
 441    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 442
 443    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
 444    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 445
 446   Standard_Boolean IsDirectionChange = Standard_False;
 447
 448   basisCurve->D3 (theU, thePBasis, theV1basis, theV2basis, theV3basis);
 449
 450   if(theV1basis.SquareMagnitude() <= gp::Resolution())
 451     IsDirectionChange = AdjustDerivative(basisCurve, 3, theU, theV1basis, theV2basis, theV3basis);
 452
 453   CSLib_Offset::D2(thePBasis, theV1basis, theV2basis, theV3basis, direction, offsetValue,
 454                    IsDirectionChange, theP, theV1, theV2);
 455 }
 456
 457
 458 //=======================================================================
 459 //function : FirstParameter
 460 //purpose  :
 461 //=======================================================================
 462
 463 Standard_Real Geom_OffsetCurve::FirstParameter () const {
 464
 465    return basisCurve->FirstParameter();
 466 }
 467
 468
 469 //=======================================================================
 470 //function : LastParameter
 471 //purpose  :
 472 //=======================================================================
 473
 474 Standard_Real Geom_OffsetCurve::LastParameter () const {
 475
 476    return basisCurve->LastParameter();
 477 }
 478
 479
 480 //=======================================================================
 481 //function : Offset
 482 //purpose  :
 483 //=======================================================================
 484
 485 Standard_Real Geom_OffsetCurve::Offset () const { return offsetValue; }
 486
 487 //=======================================================================
 488 //function : Value
 489 //purpose  :
 490 //=======================================================================
 491
 492 void Geom_OffsetCurve::Value (const Standard_Real theU, Pnt& theP,
 493                               Pnt& thePbasis,  Vec& theV1basis) const
 494 {
 495   if (myBasisCurveContinuity == GeomAbs_C0)
 496     Geom_UndefinedValue::Raise("Exception: Basis curve is C0 continuity!");
 497
 498   basisCurve->D1(theU, thePbasis, theV1basis);
 499   D0(theU,theP);
 500 }
 501
 502
 503 //=======================================================================
 504 //function : IsClosed
 505 //purpose  :
 506 //=======================================================================
 507
 508 Standard_Boolean Geom_OffsetCurve::IsClosed () const
 509 {
 510   gp_Pnt PF,PL;
 511   D0(FirstParameter(),PF);
 512   D0(LastParameter(),PL);
 513   return ( PF.Distance(PL) <= gp::Resolution());
 514 }
 515
 516
 517
 518 //=======================================================================
 519 //function : IsCN
 520 //purpose  :
 521 //=======================================================================
 522
 523 Standard_Boolean Geom_OffsetCurve::IsCN (const Standard_Integer N) const {
 524
 525    Standard_RangeError_Raise_if (N < 0, " ");
 526    return basisCurve->IsCN (N + 1);
 527 }
 528
 529
 530 //=======================================================================
 531 //function : Transform
 532 //purpose  :
 533 //=======================================================================
 534
 535 void Geom_OffsetCurve::Transform (const Trsf& T) {
 536
 537   basisCurve->Transform (T);
 538   direction.Transform(T);
 539   offsetValue *= T.ScaleFactor();
 540 }
 541
 542 //=======================================================================
 543 //function : TransformedParameter
 544 //purpose  :
 545 //=======================================================================
 546
 547 Standard_Real Geom_OffsetCurve::TransformedParameter(const Standard_Real U,
 548                                                      const gp_Trsf& T) const
 549 {
 550   return basisCurve->TransformedParameter(U,T);
 551 }
 552
 553 //=======================================================================
 554 //function : ParametricTransformation
 555 //purpose  :
 556 //=======================================================================
 557
 558 Standard_Real Geom_OffsetCurve::ParametricTransformation(const gp_Trsf& T)
 559 const
 560 {
 561   return basisCurve->ParametricTransformation(T);
 562 }
 563
 564 //=======================================================================
 565 //function : GetBasisCurveContinuity
 566 //purpose  :
 567 //=======================================================================
 568 GeomAbs_Shape Geom_OffsetCurve::GetBasisCurveContinuity() const
 569 {
 570   return myBasisCurveContinuity;
 571 }
 572
 573
 574 // ============= Auxiliary functions ===================
 575 Standard_Boolean AdjustDerivative(const Handle(Geom_Curve)& theCurve, Standard_Integer theMaxDerivative,
 576                                   Standard_Real theU, gp_Vec& theD1, gp_Vec& theD2,
 577                                   gp_Vec& theD3, gp_Vec& theD4)
 578 {
 579   static const Standard_Real aTol = gp::Resolution();
 580
 581   Standard_Boolean IsDirectionChange = Standard_False;
 582   const Standard_Real anUinfium   = theCurve->FirstParameter();
 583   const Standard_Real anUsupremum = theCurve->LastParameter();
 584
 585   const Standard_Real DivisionFactor = 1.e-3;
 586   Standard_Real du;
 587   if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst()))
 588     du = 0.0;
 589   else
 590     du = anUsupremum - anUinfium;
 591
 592   const Standard_Real aDelta = Max(du * DivisionFactor, MinStep);
 593
 594   //Derivative is approximated by Taylor-series
 595   Standard_Integer anIndex = 1; //Derivative order
 596   gp_Vec V;
 597
 598   do
 599   {
 600     V =  theCurve->DN(theU, ++anIndex);
 601   }
 602   while((V.SquareMagnitude() <= aTol) && anIndex < maxDerivOrder);
 603
 604   Standard_Real u;
 605
 606   if(theU-anUinfium < aDelta)
 607     u = theU+aDelta;
 608   else
 609     u = theU-aDelta;
 610
 611   gp_Pnt P1, P2;
 612   theCurve->D0(Min(theU, u), P1);
 613   theCurve->D0(Max(theU, u), P2);
 614
 615   gp_Vec V1(P1, P2);
 616   IsDirectionChange = V.Dot(V1) < 0.0;
 617   Standard_Real aSign = IsDirectionChange ? -1.0 : 1.0;
 618
 619   theD1 = V * aSign;
 620   gp_Vec* aDeriv[3] = {&theD2, &theD3, &theD4};
 621   for (Standard_Integer i = 1; i < theMaxDerivative; i++)
 622     *(aDeriv[i-1]) = theCurve->DN(theU, anIndex + i) * aSign;
 623
 624   return IsDirectionChange;
 625 }