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
  43 typedef Geom_OffsetCurve         OffsetCurve;
  44 typedef Handle(Geom_OffsetCurve) Handle(OffsetCurve);
  45 typedef Geom_Curve               Curve;
  46 typedef Handle(Geom_Curve)       Handle(Curve);
  47 typedef Handle(Geom_Geometry)    Handle(Geometry);
  48
  49 typedef gp_Dir  Dir;
  50 typedef gp_Pnt  Pnt;
  51 typedef gp_Trsf Trsf;
  52 typedef gp_Vec  Vec;
  53 typedef gp_XYZ  XYZ;
  54
  55
  56
  57
  58 //ordre de derivation maximum pour la recherche de la premiere
  59 //derivee non nulle
  60 static const int maxDerivOrder = 3;
  61 static const Standard_Real MinStep   = 1e-7;
  62 static const Standard_Real MyAngularToleranceForG1 = Precision::Angular();
  63
  64
  65
  66
  67 //=======================================================================
  68 //function : Copy
  69 //purpose  :
  70 //=======================================================================
  71
  72 Handle(Geom_Geometry) Geom_OffsetCurve::Copy () const {
  73
  74  Handle(OffsetCurve) C;
  75  C = new OffsetCurve (basisCurve, offsetValue, direction);
  76  return C;
  77 }
  78
  79
  80
  81 //=======================================================================
  82 //function : Geom_OffsetCurve
  83 //purpose  : Basis curve cannot be an Offset curve or trimmed from
  84 //            offset curve.
  85 //=======================================================================
  86
  87 Geom_OffsetCurve::Geom_OffsetCurve (const Handle(Geom_Curve)& theCurve,
  88                                     const Standard_Real       theOffset,
  89                                     const gp_Dir&             theDir,
  90                                     const Standard_Boolean    isTheNotCheckC0)
  91  : direction(theDir), offsetValue(theOffset)
  92 {
  93   SetBasisCurve (theCurve, isTheNotCheckC0);
  94 }
  95
  96
  97 //=======================================================================
  98 //function : Reverse
  99 //purpose  :
 100 //=======================================================================
 101
 102 void Geom_OffsetCurve::Reverse ()
 103 {
 104   basisCurve->Reverse();
 105   offsetValue = -offsetValue;
 106 }
 107
 108
 109 //=======================================================================
 110 //function : ReversedParameter
 111 //purpose  :
 112 //=======================================================================
 113
 114 Standard_Real Geom_OffsetCurve::ReversedParameter( const Standard_Real U) const
 115 {
 116   return basisCurve->ReversedParameter( U);
 117 }
 118
 119 //=======================================================================
 120 //function : Direction
 121 //purpose  :
 122 //=======================================================================
 123
 124 const gp_Dir& Geom_OffsetCurve::Direction () const
 125   { return direction; }
 126
 127 //=======================================================================
 128 //function : SetDirection
 129 //purpose  :
 130 //=======================================================================
 131
 132 void Geom_OffsetCurve::SetDirection (const Dir& V)
 133   { direction = V; }
 134
 135 //=======================================================================
 136 //function : SetOffsetValue
 137 //purpose  :
 138 //=======================================================================
 139
 140 void Geom_OffsetCurve::SetOffsetValue (const Standard_Real D)
 141   { offsetValue = D; }
 142
 143
 144 //=======================================================================
 145 //function : IsPeriodic
 146 //purpose  :
 147 //=======================================================================
 148
 149 Standard_Boolean Geom_OffsetCurve::IsPeriodic () const
 150 {
 151   return basisCurve->IsPeriodic();
 152 }
 153
 154 //=======================================================================
 155 //function : Period
 156 //purpose  :
 157 //=======================================================================
 158
 159 Standard_Real Geom_OffsetCurve::Period () const
 160 {
 161   return basisCurve->Period();
 162 }
 163
 164 //=======================================================================
 165 //function : SetBasisCurve
 166 //purpose  :
 167 //=======================================================================
 168
 169 void Geom_OffsetCurve::SetBasisCurve (const Handle(Curve)& C,
 170                                       const Standard_Boolean isNotCheckC0)
 171 {
 172   const Standard_Real aUf = C->FirstParameter(),
 173                       aUl = C->LastParameter();
 174   Handle(Curve) aCheckingCurve =  Handle(Curve)::DownCast(C->Copy());
 175   Standard_Boolean isTrimmed = Standard_False;
 176
 177   while(aCheckingCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)) ||
 178         aCheckingCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve)))
 179   {
 180     if (aCheckingCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)))
 181     {
 182       Handle(Geom_TrimmedCurve) aTrimC =
 183                 Handle(Geom_TrimmedCurve)::DownCast(aCheckingCurve);
 184       aCheckingCurve = aTrimC->BasisCurve();
 185       isTrimmed = Standard_True;
 186     }
 187
 188     if (aCheckingCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve)))
 189     {
 190       Handle(Geom_OffsetCurve) aOC =
 191             Handle(Geom_OffsetCurve)::DownCast(aCheckingCurve);
 192       aCheckingCurve = aOC->BasisCurve();
 193       Standard_Real PrevOff = aOC->Offset();
 194       gp_Vec V1(aOC->Direction());
 195       gp_Vec V2(direction);
 196       gp_Vec Vdir(PrevOff*V1 + offsetValue*V2);
 197
 198       if (offsetValue >= 0.)
 199       {
 200         offsetValue = Vdir.Magnitude();
 201         direction.SetXYZ(Vdir.XYZ());
 202       }
 203       else
 204       {
 205         offsetValue = -Vdir.Magnitude();
 206         direction.SetXYZ((-Vdir).XYZ());
 207       }
 208     }
 209   }
 210
 211   myBasisCurveContinuity = aCheckingCurve->Continuity();
 212
 213   Standard_Boolean isC0 = !isNotCheckC0 &&
 214                           (myBasisCurveContinuity == GeomAbs_C0);
 215
 216   // Basis curve must be at least C1
 217   if (isC0 && aCheckingCurve->IsKind(STANDARD_TYPE(Geom_BSplineCurve)))
 218   {
 219     Handle(Geom_BSplineCurve) aBC = Handle(Geom_BSplineCurve)::DownCast(aCheckingCurve);
 220     if(aBC->IsG1(aUf, aUl, MyAngularToleranceForG1))
 221     {
 222       //Checking if basis curve has more smooth (C1, G2 and above) is not done.
 223       //It can be done in case of need.
 224       myBasisCurveContinuity = GeomAbs_G1;
 225       isC0 = Standard_False;
 226     }
 227
 228     // Raise exception if still C0
 229     if (isC0)
 230       Standard_ConstructionError::Raise("Offset on C0 curve");
 231   }
 232   //
 233   if(isTrimmed)
 234   {
 235     basisCurve = new Geom_TrimmedCurve(aCheckingCurve, aUf, aUl);
 236   }
 237   else
 238   {
 239     basisCurve = aCheckingCurve;
 240   }
 241 }
 242
 243
 244
 245 //=======================================================================
 246 //function : BasisCurve
 247 //purpose  :
 248 //=======================================================================
 249
 250 Handle(Curve) Geom_OffsetCurve::BasisCurve () const
 251 {
 252   return basisCurve;
 253 }
 254
 255
 256 //=======================================================================
 257 //function : Continuity
 258 //purpose  :
 259 //=======================================================================
 260
 261 GeomAbs_Shape Geom_OffsetCurve::Continuity () const {
 262
 263   GeomAbs_Shape OffsetShape=GeomAbs_C0;
 264   switch (myBasisCurveContinuity) {
 265     case GeomAbs_C0 : OffsetShape = GeomAbs_C0;       break;
 266     case GeomAbs_C1 : OffsetShape = GeomAbs_C0;       break;
 267     case GeomAbs_C2 : OffsetShape = GeomAbs_C1;       break;
 268     case GeomAbs_C3 : OffsetShape = GeomAbs_C2;       break;
 269     case GeomAbs_CN : OffsetShape = GeomAbs_CN;       break;
 270     case GeomAbs_G1 : OffsetShape = GeomAbs_G1;       break;
 271     case GeomAbs_G2 : OffsetShape = GeomAbs_G2;       break;
 272   }
 273   return OffsetShape;
 274 }
 275
 276
 277 //=======================================================================
 278 //function : D0
 279 //purpose  :
 280 //=======================================================================
 281
 282 void Geom_OffsetCurve::D0 (const Standard_Real U, Pnt& P) const
 283 {
 284   gp_Pnt PBasis;
 285   gp_Vec VBasis;
 286   D0(U,P,PBasis,VBasis);
 287 }
 288
 289
 290 //=======================================================================
 291 //function : D1
 292 //purpose  :
 293 //=======================================================================
 294
 295 void Geom_OffsetCurve::D1 (const Standard_Real U, Pnt& P, Vec& V1) const
 296 {
 297   gp_Pnt PBasis;
 298   gp_Vec V1Basis,V2Basis;
 299   D1(U,P,PBasis,V1,V1Basis,V2Basis);
 300 }
 301
 302
 303
 304 //=======================================================================
 305 //function : D2
 306 //purpose  :
 307 //=======================================================================
 308
 309 void Geom_OffsetCurve::D2 (const Standard_Real U, Pnt& P, Vec& V1, Vec& V2) const
 310 {
 311   gp_Pnt PBasis;
 312   gp_Vec V1Basis,V2Basis,V3Basis;
 313   D2(U,P,PBasis,V1,V2,V1Basis,V2Basis,V3Basis);
 314 }
 315
 316
 317 //=======================================================================
 318 //function : D3
 319 //purpose  :
 320 //=======================================================================
 321
 322 void Geom_OffsetCurve::D3 (const Standard_Real theU, Pnt& P, Vec& theV1, Vec& V2, Vec& V3)
 323 const {
 324
 325
 326    // P(u) = p(u) + Offset * Ndir / R
 327    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
 328
 329    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 330
 331    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
 332    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 333
 334    //P"'(u) = p"'(u) + (Offset / R) * (D3Ndir - (3.0 * Dr/R**2) * D2Ndir -
 335    //         (3.0 * D2r / R2) * DNdir + (3.0 * Dr * Dr / R4) * DNdir -
 336    //         (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir +
 337    //         (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir
 338
 339   const Standard_Real aTol = gp::Resolution();
 340
 341   Standard_Boolean IsDirectionChange = Standard_False;
 342
 343   basisCurve->D3 (theU, P, theV1, V2, V3);
 344   Vec V4 = basisCurve->DN (theU, 4);
 345   if(theV1.Magnitude() <= aTol)
 346     {
 347     const Standard_Real anUinfium   = basisCurve->FirstParameter();
 348     const Standard_Real anUsupremum = basisCurve->LastParameter();
 349
 350     const Standard_Real DivisionFactor = 1.e-3;
 351     Standard_Real du;
 352     if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst()))
 353       du = 0.0;
 354     else
 355       du = anUsupremum-anUinfium;
 356
 357     const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
 358 //Derivative is approximated by Taylor-series
 359
 360     Standard_Integer anIndex = 1; //Derivative order
 361     Vec V;
 362
 363     do
 364       {
 365       V =  basisCurve->DN(theU,++anIndex);
 366       }
 367     while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder);
 368
 369     Standard_Real u;
 370
 371     if(theU-anUinfium < aDelta)
 372       u = theU+aDelta;
 373     else
 374       u = theU-aDelta;
 375
 376     Pnt P1, P2;
 377     basisCurve->D0(Min(theU, u),P1);
 378     basisCurve->D0(Max(theU, u),P2);
 379
 380     Vec V1(P1,P2);
 381     Standard_Real aDirFactor = V.Dot(V1);
 382
 383     if(aDirFactor < 0.0)
 384       {
 385       theV1 = -V;
 386       V2 = -basisCurve->DN (theU, anIndex + 1);
 387       V3 = -basisCurve->DN (theU, anIndex + 2);
 388       V4 = -basisCurve->DN (theU, anIndex + 3);
 389
 390       IsDirectionChange = Standard_True;
 391       }
 392     else
 393       {
 394       theV1 = V;
 395       V2 = basisCurve->DN (theU, anIndex + 1);
 396       V3 = basisCurve->DN (theU, anIndex + 2);
 397       V4 = basisCurve->DN (theU, anIndex + 3);
 398       }
 399     }//if(V1.Magnitude() <= aTol)
 400
 401
 402   XYZ OffsetDir = direction.XYZ();
 403   XYZ Ndir      = (theV1.XYZ()).Crossed (OffsetDir);
 404   XYZ DNdir     = (V2.XYZ()).Crossed (OffsetDir);
 405   XYZ D2Ndir    = (V3.XYZ()).Crossed (OffsetDir);
 406   XYZ D3Ndir    = (V4.XYZ()).Crossed (OffsetDir);
 407   Standard_Real R2  = Ndir.SquareModulus();
 408   Standard_Real R   = Sqrt (R2);
 409   Standard_Real R3  = R2 * R;
 410   Standard_Real R4  = R2 * R2;
 411   Standard_Real R5  = R3 * R2;
 412   Standard_Real R6  = R3 * R3;
 413   Standard_Real R7  = R5 * R2;
 414   Standard_Real Dr  = Ndir.Dot (DNdir);
 415   Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
 416   Standard_Real D3r = Ndir.Dot (D3Ndir) + 3.0 * DNdir.Dot (D2Ndir);
 417   if (R7 <= gp::Resolution()) {
 418     if (R6 <= gp::Resolution())  Geom_UndefinedDerivative::Raise();
 419     // V3 = P"' (U) :
 420     D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * Dr / R2));
 421     D3Ndir.Subtract (DNdir.Multiplied (3.0 * ((D2r/R2) + (Dr*Dr/R4))));
 422     D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 -
 423                                  15.0*Dr*Dr*Dr/R6 - D3r));
 424     D3Ndir.Multiply (offsetValue/R);
 425
 426     if(IsDirectionChange)
 427       V3=-V3;
 428
 429     V3.Add (Vec(D3Ndir));
 430
 431     // V2 = P" (U) :
 432     Standard_Real R4 = R2 * R2;
 433     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
 434     D2Ndir.Subtract (Ndir.Multiplied ((3.0 * Dr * Dr / R4) - (D2r / R2)));
 435     D2Ndir.Multiply (offsetValue / R);
 436     V2.Add (Vec(D2Ndir));
 437     // V1 = P' (U) :
 438     DNdir.Multiply(R);
 439     DNdir.Subtract (Ndir.Multiplied (Dr/R));
 440     DNdir.Multiply (offsetValue/R2);
 441     theV1.Add (Vec(DNdir));
 442   }
 443   else {
 444     // V3 = P"' (U) :
 445     D3Ndir.Divide (R);
 446     D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * Dr / R3));
 447     D3Ndir.Subtract (DNdir.Multiplied ((3.0 * ((D2r/R3) + (Dr*Dr)/R5))));
 448     D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 -
 449                 15.0*Dr*Dr*Dr/R7 - D3r));
 450     D3Ndir.Multiply (offsetValue);
 451
 452     if(IsDirectionChange)
 453       V3=-V3;
 454
 455     V3.Add (Vec(D3Ndir));
 456
 457     // V2 = P" (U) :
 458     D2Ndir.Divide (R);
 459     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R3));
 460     D2Ndir.Subtract (Ndir.Multiplied ((3.0 * Dr * Dr / R5) - (D2r / R3)));
 461     D2Ndir.Multiply (offsetValue);
 462     V2.Add (Vec(D2Ndir));
 463     // V1 = P' (U) :
 464     DNdir.Multiply (offsetValue/R);
 465     DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
 466     theV1.Add (Vec(DNdir));
 467   }
 468   //P (U) :
 469   D0(theU,P);
 470 }
 471
 472
 473 //=======================================================================
 474 //function : DN
 475 //purpose  :
 476 //=======================================================================
 477
 478 Vec Geom_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const
 479   {
 480   Standard_RangeError_Raise_if (N < 1, "Exception: "
 481                               "Geom_OffsetCurve::DN(...). N<1.");
 482
 483   gp_Vec VN, Vtemp;
 484   gp_Pnt Ptemp;
 485   switch (N)
 486     {
 487     case 1:
 488       D1( U, Ptemp, VN);
 489       break;
 490     case 2:
 491       D2( U, Ptemp, Vtemp, VN);
 492       break;
 493     case 3:
 494       D3( U, Ptemp, Vtemp, Vtemp, VN);
 495       break;
 496     default:
 497       Standard_NotImplemented::Raise("Exception: "
 498         "Derivative order is greater than 3. Cannot compute of derivative.");
 499   }
 500
 501   return VN;
 502 }
 503
 504 //=======================================================================
 505 //function : D0
 506 //purpose  :
 507 //=======================================================================
 508
 509 void  Geom_OffsetCurve::D0(const Standard_Real theU, gp_Pnt& theP,
 510                            gp_Pnt& thePbasis, gp_Vec& theV1basis)const
 511   {
 512   const Standard_Real aTol = gp::Resolution();
 513
 514   basisCurve->D1 (theU, thePbasis, theV1basis);
 515   Standard_Real Ndu = theV1basis.Magnitude();
 516
 517   if(Ndu <= aTol)
 518     {
 519     const Standard_Real anUinfium   = basisCurve->FirstParameter();
 520     const Standard_Real anUsupremum = basisCurve->LastParameter();
 521
 522     const Standard_Real DivisionFactor = 1.e-3;
 523     Standard_Real du;
 524     if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst()))
 525       du = 0.0;
 526     else
 527       du = anUsupremum-anUinfium;
 528
 529     const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
 530 //Derivative is approximated by Taylor-series
 531
 532     Standard_Integer anIndex = 1; //Derivative order
 533     gp_Vec V;
 534
 535     do
 536       {
 537       V =  basisCurve->DN(theU,++anIndex);
 538       Ndu = V.Magnitude();
 539       }
 540     while((Ndu <= aTol) && anIndex < maxDerivOrder);
 541
 542     Standard_Real u;
 543
 544     if(theU-anUinfium < aDelta)
 545       u = theU+aDelta;
 546     else
 547       u = theU-aDelta;
 548
 549     gp_Pnt P1, P2;
 550     basisCurve->D0(Min(theU, u),P1);
 551     basisCurve->D0(Max(theU, u),P2);
 552
 553     gp_Vec V1(P1,P2);
 554     Standard_Real aDirFactor = V.Dot(V1);
 555
 556     if(aDirFactor < 0.0)
 557       theV1basis = -V;
 558     else
 559       theV1basis = V;
 560
 561     Ndu = theV1basis.Magnitude();
 562     }//if(Ndu <= aTol)
 563
 564   XYZ Ndir = (theV1basis.XYZ()).Crossed (direction.XYZ());
 565   Standard_Real R = Ndir.Modulus();
 566   if (R <= gp::Resolution())
 567     Geom_UndefinedValue::Raise("Exception: Undefined normal vector "
 568           "because tangent vector has zero-magnitude!");
 569
 570   Ndir.Multiply (offsetValue/R);
 571   Ndir.Add (thePbasis.XYZ());
 572   theP.SetXYZ(Ndir);
 573 }
 574
 575 //=======================================================================
 576 //function : D1
 577 //purpose  :
 578 //=======================================================================
 579
 580 void Geom_OffsetCurve::D1 ( const Standard_Real theU,
 581                             Pnt& P , Pnt& PBasis ,
 582                             Vec& theV1, Vec& V1basis, Vec& V2basis) const {
 583
 584    // P(u) = p(u) + Offset * Ndir / R
 585    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
 586
 587    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 588
 589   const Standard_Real aTol = gp::Resolution();
 590
 591   basisCurve->D2 (theU, PBasis, V1basis, V2basis);
 592   theV1 = V1basis;
 593   Vec V2 = V2basis;
 594
 595   if(theV1.Magnitude() <= aTol)
 596     {
 597     const Standard_Real anUinfium   = basisCurve->FirstParameter();
 598     const Standard_Real anUsupremum = basisCurve->LastParameter();
 599
 600     const Standard_Real DivisionFactor = 1.e-3;
 601     Standard_Real du;
 602     if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst()))
 603       du = 0.0;
 604     else
 605       du = anUsupremum-anUinfium;
 606
 607     const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
 608 //Derivative is approximated by Taylor-series
 609
 610     Standard_Integer anIndex = 1; //Derivative order
 611     Vec V;
 612
 613     do
 614       {
 615       V =  basisCurve->DN(theU,++anIndex);
 616       }
 617     while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder);
 618
 619     Standard_Real u;
 620
 621     if(theU-anUinfium < aDelta)
 622       u = theU+aDelta;
 623     else
 624       u = theU-aDelta;
 625
 626     Pnt P1, P2;
 627     basisCurve->D0(Min(theU, u),P1);
 628     basisCurve->D0(Max(theU, u),P2);
 629
 630     Vec V1(P1,P2);
 631     Standard_Real aDirFactor = V.Dot(V1);
 632
 633     if(aDirFactor < 0.0)
 634       {
 635       theV1 = -V;
 636       V2 = - basisCurve->DN (theU, anIndex+1);
 637       }
 638     else
 639       {
 640       theV1 = V;
 641       V2 = basisCurve->DN (theU, anIndex+1);
 642       }
 643
 644     V2basis = V2;
 645     V1basis = theV1;
 646     }//if(theV1.Magnitude() <= aTol)
 647
 648   XYZ OffsetDir = direction.XYZ();
 649   XYZ Ndir  = (theV1.XYZ()).Crossed (OffsetDir);
 650   XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
 651   Standard_Real R2 = Ndir.SquareModulus();
 652   Standard_Real R  = Sqrt (R2);
 653   Standard_Real R3 = R * R2;
 654   Standard_Real Dr = Ndir.Dot (DNdir);
 655   if (R3 <= gp::Resolution()) {
 656     //We try another computation but the stability is not very good.
 657     if (R2 <= gp::Resolution()) Geom_UndefinedDerivative::Raise();
 658     DNdir.Multiply(R);
 659     DNdir.Subtract (Ndir.Multiplied (Dr/R));
 660     DNdir.Multiply (offsetValue/R2);
 661     theV1.Add (Vec(DNdir));
 662   }
 663   else {
 664     // Same computation as IICURV in EUCLID-IS because the stability is
 665     // better
 666     DNdir.Multiply (offsetValue/R);
 667     DNdir.Subtract (Ndir.Multiplied (offsetValue * Dr/R3));
 668     theV1.Add (Vec(DNdir));
 669   }
 670   D0(theU,P);
 671 }
 672
 673
 674 //=======================================================================
 675 //function : D2
 676 //purpose  :
 677 //=======================================================================
 678
 679 void Geom_OffsetCurve::D2 (const Standard_Real theU,
 680                            Pnt& P      , Pnt& PBasis ,
 681                            Vec& theV1     , Vec& V2     ,
 682                            Vec& V1basis, Vec& V2basis, Vec& V3basis) const {
 683
 684    // P(u) = p(u) + Offset * Ndir / R
 685    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
 686
 687    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 688
 689    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
 690    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 691
 692   const Standard_Real aTol = gp::Resolution();
 693
 694   Standard_Boolean IsDirectionChange = Standard_False;
 695
 696   basisCurve->D3 (theU, PBasis, V1basis, V2basis, V3basis);
 697
 698   theV1  = V1basis;
 699   V2     = V2basis;
 700   Vec V3 = V3basis;
 701
 702   if(theV1.Magnitude() <= aTol)
 703     {
 704     const Standard_Real anUinfium   = basisCurve->FirstParameter();
 705     const Standard_Real anUsupremum = basisCurve->LastParameter();
 706
 707     const Standard_Real DivisionFactor = 1.e-3;
 708     Standard_Real du;
 709     if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst()))
 710       du = 0.0;
 711     else
 712       du = anUsupremum-anUinfium;
 713
 714     const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
 715 //Derivative is approximated by Taylor-series
 716
 717     Standard_Integer anIndex = 1; //Derivative order
 718     Vec V;
 719
 720     do
 721       {
 722       V =  basisCurve->DN(theU,++anIndex);
 723       }
 724     while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder);
 725
 726     Standard_Real u;
 727
 728     if(theU-anUinfium < aDelta)
 729       u = theU+aDelta;
 730     else
 731       u = theU-aDelta;
 732
 733     Pnt P1, P2;
 734     basisCurve->D0(Min(theU, u),P1);
 735     basisCurve->D0(Max(theU, u),P2);
 736
 737     Vec V1(P1,P2);
 738     Standard_Real aDirFactor = V.Dot(V1);
 739
 740     if(aDirFactor < 0.0)
 741       {
 742       theV1 = -V;
 743       V2 = -basisCurve->DN (theU, anIndex+1);
 744       V3 = -basisCurve->DN (theU, anIndex + 2);
 745
 746       IsDirectionChange = Standard_True;
 747       }
 748     else
 749       {
 750       theV1 = V;
 751       V2 = basisCurve->DN (theU, anIndex+1);
 752       V3 = basisCurve->DN (theU, anIndex + 2);
 753       }
 754
 755     V2basis = V2;
 756     V1basis = theV1;
 757     }//if(V1.Magnitude() <= aTol)
 758
 759   XYZ OffsetDir = direction.XYZ();
 760   XYZ Ndir   = (theV1.XYZ()).Crossed (OffsetDir);
 761   XYZ DNdir  = (V2.XYZ()).Crossed (OffsetDir);
 762   XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir);
 763   Standard_Real R2    = Ndir.SquareModulus();
 764   Standard_Real R     = Sqrt (R2);
 765   Standard_Real R3    = R2 * R;
 766   Standard_Real R4    = R2 * R2;
 767   Standard_Real R5    = R3 * R2;
 768   Standard_Real Dr    = Ndir.Dot (DNdir);
 769   Standard_Real D2r   = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
 770
 771   if (R5 <= gp::Resolution()) {
 772     //We try another computation but the stability is not very good
 773     //dixit ISG.
 774     if (R4 <= gp::Resolution())  Geom_UndefinedDerivative::Raise();
 775     // V2 = P" (U) :
 776     Standard_Real R4 = R2 * R2;
 777     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
 778     D2Ndir.Add (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
 779     D2Ndir.Multiply (offsetValue / R);
 780
 781     if(IsDirectionChange)
 782       V2=-V2;
 783
 784     V2.Add (Vec(D2Ndir));
 785
 786     // V1 = P' (U) :
 787     DNdir.Multiply(R);
 788     DNdir.Subtract (Ndir.Multiplied (Dr/R));
 789     DNdir.Multiply (offsetValue/R2);
 790     theV1.Add (Vec(DNdir));
 791   }
 792   else {
 793     // Same computation as IICURV in EUCLID-IS because the stability is
 794     // better.
 795     // V2 = P" (U) :
 796     D2Ndir.Multiply (offsetValue/R);
 797     D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
 798     D2Ndir.Add (Ndir.Multiplied (
 799                      offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))
 800                                      )
 801                     );
 802
 803     if(IsDirectionChange)
 804       V2=-V2;
 805
 806     V2.Add (Vec(D2Ndir));
 807
 808     // V1 = P' (U) :
 809     DNdir.Multiply (offsetValue/R);
 810     DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
 811     theV1.Add (Vec(DNdir));
 812   }
 813   //P (U) :
 814   D0(theU,P);
 815 }
 816
 817
 818 //=======================================================================
 819 //function : FirstParameter
 820 //purpose  :
 821 //=======================================================================
 822
 823 Standard_Real Geom_OffsetCurve::FirstParameter () const {
 824
 825    return basisCurve->FirstParameter();
 826 }
 827
 828
 829 //=======================================================================
 830 //function : LastParameter
 831 //purpose  :
 832 //=======================================================================
 833
 834 Standard_Real Geom_OffsetCurve::LastParameter () const {
 835
 836    return basisCurve->LastParameter();
 837 }
 838
 839
 840 //=======================================================================
 841 //function : Offset
 842 //purpose  :
 843 //=======================================================================
 844
 845 Standard_Real Geom_OffsetCurve::Offset () const { return offsetValue; }
 846
 847 //=======================================================================
 848 //function : Value
 849 //purpose  :
 850 //=======================================================================
 851
 852 void Geom_OffsetCurve::Value (const Standard_Real theU, Pnt& theP,
 853                               Pnt& thePbasis,  Vec& theV1basis) const
 854 {
 855   if (myBasisCurveContinuity == GeomAbs_C0)
 856     Geom_UndefinedValue::Raise("Exception: Basis curve is C0 continuity!");
 857
 858   basisCurve->D1(theU, thePbasis, theV1basis);
 859   D0(theU,theP);
 860 }
 861
 862
 863 //=======================================================================
 864 //function : IsClosed
 865 //purpose  :
 866 //=======================================================================
 867
 868 Standard_Boolean Geom_OffsetCurve::IsClosed () const
 869 {
 870   gp_Pnt PF,PL;
 871   D0(FirstParameter(),PF);
 872   D0(LastParameter(),PL);
 873   return ( PF.Distance(PL) <= gp::Resolution());
 874 }
 875
 876
 877
 878 //=======================================================================
 879 //function : IsCN
 880 //purpose  :
 881 //=======================================================================
 882
 883 Standard_Boolean Geom_OffsetCurve::IsCN (const Standard_Integer N) const {
 884
 885    Standard_RangeError_Raise_if (N < 0, " ");
 886    return basisCurve->IsCN (N + 1);
 887 }
 888
 889
 890 //=======================================================================
 891 //function : Transform
 892 //purpose  :
 893 //=======================================================================
 894
 895 void Geom_OffsetCurve::Transform (const Trsf& T) {
 896
 897   basisCurve->Transform (T);
 898   direction.Transform(T);
 899   offsetValue *= T.ScaleFactor();
 900 }
 901
 902 //=======================================================================
 903 //function : TransformedParameter
 904 //purpose  :
 905 //=======================================================================
 906
 907 Standard_Real Geom_OffsetCurve::TransformedParameter(const Standard_Real U,
 908                                                      const gp_Trsf& T) const
 909 {
 910   return basisCurve->TransformedParameter(U,T);
 911 }
 912
 913 //=======================================================================
 914 //function : ParametricTransformation
 915 //purpose  :
 916 //=======================================================================
 917
 918 Standard_Real Geom_OffsetCurve::ParametricTransformation(const gp_Trsf& T)
 919 const
 920 {
 921   return basisCurve->ParametricTransformation(T);
 922 }
 923
 924 //=======================================================================
 925 //function : GetBasisCurveContinuity
 926 //purpose  :
 927 //=======================================================================
 928 GeomAbs_Shape Geom_OffsetCurve::GetBasisCurveContinuity() const
 929 {
 930   return myBasisCurveContinuity;
 931 }