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