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