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-2012 OPEN CASCADE SAS
   5 //
   6 // The content of this file is subject to the Open CASCADE Technology Public
   7 // License Version 6.5 (the "License"). You may not use the content of this file
   8 // except in compliance with the License. Please obtain a copy of the License
   9 // at http://www.opencascade.org and read it completely before using this file.
  10 //
  11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  13 //
  14 // The Original Code and all software distributed under the License is
  15 // distributed on an "AS IS" basis, without warranty of any kind, and the
  16 // Initial Developer hereby disclaims all such warranties, including without
  17 // limitation, any warranties of merchantability, fitness for a particular
  18 // purpose or non-infringement. Please see the License for the specific terms
  19 // and conditions governing the rights and limitations under the License.
  20
  21
  22
  23 //  modified by Edward AGAPOV (eap) Jan 28 2002 --- DN(), occ143(BUC60654)
  24
  25
  26
  27 #include <Geom2d_OffsetCurve.ixx>
  28 #include <gp.hxx>
  29 #include <Standard_ConstructionError.hxx>
  30 #include <Standard_RangeError.hxx>
  31 #include <Standard_NotImplemented.hxx>
  32 #include <Geom2d_UndefinedDerivative.hxx>
  33 #include <Geom2d_UndefinedValue.hxx>
  34 #include <Geom2d_Line.hxx>
  35 #include <Geom2d_Circle.hxx>
  36 #include <Geom2d_Ellipse.hxx>
  37 #include <Geom2d_Hyperbola.hxx>
  38 #include <Geom2d_Parabola.hxx>
  39 #include <Geom2d_BezierCurve.hxx>
  40 #include <Geom2d_BSplineCurve.hxx>
  41 #include <Geom2d_TrimmedCurve.hxx>
  42 #include <gp_XY.hxx>
  43
  44 typedef Handle(Geom2d_OffsetCurve) Handle(OffsetCurve);
  45 typedef Geom2d_OffsetCurve         OffsetCurve;
  46 typedef Handle(Geom2d_Geometry)    Handle(Geometry);
  47 typedef Handle(Geom2d_Curve)       Handle(Curve);
  48 typedef Geom2d_Curve               Curve;
  49 typedef gp_Dir2d  Dir2d;
  50 typedef gp_Pnt2d  Pnt2d;
  51 typedef gp_Vec2d  Vec2d;
  52 typedef gp_Trsf2d Trsf2d;
  53 typedef gp_XY     XY;
  54
  55
  56
  57 static const int MaxDegree = 9;
  58 //ordre de derivation maximum pour la recherche de la premiere
  59 //derivee non nulle
  60
  61
  62
  63 //=======================================================================
  64 //function : Copy
  65 //purpose  :
  66 //=======================================================================
  67
  68 Handle(Geom2d_Geometry) Geom2d_OffsetCurve::Copy () const
  69 {
  70   Handle(OffsetCurve) C;
  71   C = new OffsetCurve (basisCurve, offsetValue);
  72   return C;
  73 }
  74
  75
  76 //=======================================================================
  77 //function : Geom2d_OffsetCurve
  78 //purpose  :
  79 //=======================================================================
  80
  81 Geom2d_OffsetCurve::Geom2d_OffsetCurve (const Handle(Curve)& C,
  82                                         const Standard_Real Offset)
  83 : offsetValue (Offset)
  84 {
  85   if (C->DynamicType() == STANDARD_TYPE(Geom2d_OffsetCurve)) {
  86     Handle(OffsetCurve) OC = Handle(OffsetCurve)::DownCast(C);
  87     if ((OC->BasisCurve())->Continuity() == GeomAbs_C0)
  88       Standard_ConstructionError::Raise();
  89
  90     basisCurve = Handle(Curve)::DownCast((OC->BasisCurve())->Copy());
  91     offsetValue += OC->Offset();
  92   } else {
  93     if (C->Continuity() == GeomAbs_C0)
  94       Standard_ConstructionError::Raise();
  95
  96     basisCurve = Handle(Curve)::DownCast(C->Copy());
  97   }
  98 }
  99
 100 //=======================================================================
 101 //function : Reverse
 102 //purpose  :
 103 //=======================================================================
 104
 105 void Geom2d_OffsetCurve::Reverse ()
 106 {
 107   basisCurve->Reverse();
 108   offsetValue = -offsetValue;
 109 }
 110
 111 //=======================================================================
 112 //function : ReversedParameter
 113 //purpose  :
 114 //=======================================================================
 115
 116 Standard_Real Geom2d_OffsetCurve::ReversedParameter( const Standard_Real U) const
 117 {
 118   return basisCurve->ReversedParameter( U);
 119 }
 120
 121 //=======================================================================
 122 //function : SetBasisCurve
 123 //purpose  :
 124 //=======================================================================
 125
 126 void Geom2d_OffsetCurve::SetBasisCurve (const Handle(Curve)& C)
 127 {
 128   if (C->Continuity() == GeomAbs_C0)  Standard_ConstructionError::Raise();
 129   basisCurve = Handle(Geom2d_Curve)::DownCast(C->Copy());
 130 }
 131
 132 //=======================================================================
 133 //function : SetOffsetValue
 134 //purpose  :
 135 //=======================================================================
 136
 137 void Geom2d_OffsetCurve::SetOffsetValue (const Standard_Real D) { offsetValue = D; }
 138
 139 //=======================================================================
 140 //function : BasisCurve
 141 //purpose  :
 142 //=======================================================================
 143
 144 Handle(Curve) Geom2d_OffsetCurve::BasisCurve () const
 145 {
 146   return basisCurve;
 147 }
 148
 149 //=======================================================================
 150 //function : Continuity
 151 //purpose  :
 152 //=======================================================================
 153
 154 GeomAbs_Shape Geom2d_OffsetCurve::Continuity () const
 155 {
 156   GeomAbs_Shape OffsetShape=GeomAbs_C0;
 157   switch (basisCurve->Continuity()) {
 158      case GeomAbs_C0 : OffsetShape = GeomAbs_C0;   break;
 159      case GeomAbs_C1 : OffsetShape = GeomAbs_C0;   break;
 160      case GeomAbs_C2 : OffsetShape = GeomAbs_C1;   break;
 161      case GeomAbs_C3 : OffsetShape = GeomAbs_C2;   break;
 162      case GeomAbs_CN : OffsetShape = GeomAbs_CN;   break;
 163      case GeomAbs_G1 : OffsetShape = GeomAbs_G1;   break;
 164      case GeomAbs_G2 : OffsetShape = GeomAbs_G2;   break;
 165   }
 166   return OffsetShape;
 167 }
 168
 169 //=======================================================================
 170 //function : D0
 171 //purpose  :
 172 //=======================================================================
 173
 174 void Geom2d_OffsetCurve::D0 (const Standard_Real   U,
 175                                    Pnt2d& P ) const
 176 {
 177   Vec2d V1;
 178
 179   basisCurve->D1 (U, P, V1);
 180   Standard_Integer Index = 2;
 181   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 182     V1 = basisCurve->DN (U, Index);
 183     Index++;
 184   }
 185   Standard_Real A = V1.Y();
 186   Standard_Real B = - V1.X();
 187   Standard_Real R = Sqrt(A*A + B * B);
 188   if (R <= gp::Resolution())  Geom2d_UndefinedValue::Raise();
 189   A = A * offsetValue/R;
 190   B = B * offsetValue/R;
 191   P.SetCoord(P.X() + A, P.Y() + B);
 192 }
 193
 194 //=======================================================================
 195 //function : D1
 196 //purpose  :
 197 //=======================================================================
 198
 199 void Geom2d_OffsetCurve::D1 (const Standard_Real U, Pnt2d& P, Vec2d& V1) const {
 200
 201    // P(u) = p(u) + Offset * Ndir / R
 202    // with R = || p' ^ Z|| and Ndir = P' ^ Z
 203
 204    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 205
 206   Vec2d V2;
 207   basisCurve->D2 (U, P, V1, V2);
 208   Standard_Integer Index = 2;
 209   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 210     V1 = basisCurve->DN (U, Index);
 211     Index++;
 212   }
 213   if (Index != 2) { V2 = basisCurve->DN (U, Index); }
 214   XY Ndir  (V1.Y(), -V1.X());
 215   XY DNdir (V2.Y(), -V2.X());
 216   Standard_Real R2 = Ndir.SquareModulus();
 217   Standard_Real R  = Sqrt (R2);
 218   Standard_Real R3 = R * R2;
 219   Standard_Real Dr = Ndir.Dot (DNdir);
 220   if (R3 <= gp::Resolution()) {
 221     //We try another computation but the stability is not very good.
 222     if (R2 <= gp::Resolution()) Geom2d_UndefinedDerivative::Raise();
 223     DNdir.Multiply(R);
 224     DNdir.Subtract (Ndir.Multiplied (Dr/R));
 225     DNdir.Multiply (offsetValue/R2);
 226     V1.Add (Vec2d(DNdir));
 227   }
 228   else {
 229     // Same computation as IICURV in EUCLID-IS because the stability is
 230     // better
 231     DNdir.Multiply (offsetValue/R);
 232     DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
 233     V1.Add (Vec2d(DNdir));
 234   }
 235   Ndir.Multiply (offsetValue/R);
 236   Ndir.Add (P.XY());
 237   P.SetXY (Ndir);
 238 }
 239
 240 //=======================================================================
 241 //function : D2
 242 //purpose  :
 243 //=======================================================================
 244
 245 void Geom2d_OffsetCurve::D2 (const Standard_Real U,
 246                                    Pnt2d& P,
 247                                    Vec2d& V1, Vec2d& V2) const
 248 {
 249    // P(u) = p(u) + Offset * Ndir / R
 250    // with R = || p' ^ Z|| and Ndir = P' ^ Z
 251
 252    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 253
 254    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
 255    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 256
 257   Vec2d V3;
 258   basisCurve->D3 (U, P, V1, V2, V3);
 259   Standard_Integer Index = 2;
 260   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 261     V1 = basisCurve->DN (U, Index);
 262     Index++;
 263   }
 264   if (Index != 2) {
 265     V2 = basisCurve->DN (U, Index);
 266     V3 = basisCurve->DN (U, Index + 1);
 267   }
 268   XY Ndir (V1.Y(), -V1.X());
 269   XY DNdir (V2.Y(), -V2.X());
 270   XY D2Ndir (V3.Y(), -V3.X());
 271   Standard_Real R2  = Ndir.SquareModulus();
 272   Standard_Real R   = Sqrt (R2);
 273   Standard_Real R3  = R2 * R;
 274   Standard_Real R4  = R2 * R2;
 275   Standard_Real R5  = R3 * R2;
 276   Standard_Real Dr  = Ndir.Dot (DNdir);
 277   Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
 278   if (R5 <= gp::Resolution()) {
 279     //We try another computation but the stability is not very good
 280     //dixit ISG.
 281     if (R4 <= gp::Resolution()) { Geom2d_UndefinedDerivative::Raise(); }
 282     // V2 = P" (U) :
 283      Standard_Real R4 = R2 * R2;
 284      D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
 285      D2Ndir.Add (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
 286      D2Ndir.Multiply (offsetValue / R);
 287      V2.Add (Vec2d(D2Ndir));
 288      // V1 = P' (U) :
 289      DNdir.Multiply(R);
 290      DNdir.Subtract (Ndir.Multiplied (Dr/R));
 291      DNdir.Multiply (offsetValue/R2);
 292      V1.Add (Vec2d(DNdir));
 293    }
 294    else {
 295      // Same computation as IICURV in EUCLID-IS because the stability is
 296      // better.
 297      // V2 = P" (U) :
 298     D2Ndir.Multiply (offsetValue/R);
 299     D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
 300     D2Ndir.Add (Ndir.Multiplied
 301                      (offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))));
 302     V2.Add (Vec2d(D2Ndir));
 303     // V1 = P' (U)
 304      DNdir.Multiply (offsetValue/R);
 305      DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
 306      V1.Add (Vec2d(DNdir));
 307    }
 308    //P (U) :
 309    Ndir.Multiply (offsetValue/R);
 310    Ndir.Add (P.XY());
 311    P.SetXY (Ndir);
 312 }
 313
 314
 315 //=======================================================================
 316 //function : D3
 317 //purpose  :
 318 //=======================================================================
 319
 320 void Geom2d_OffsetCurve::D3 (const Standard_Real U,
 321                                    Pnt2d& P,
 322                                    Vec2d& V1, Vec2d& V2, Vec2d& V3) const {
 323
 324
 325    // P(u) = p(u) + Offset * Ndir / R
 326    // with R = || p' ^ Z|| and Ndir = P' ^ Z
 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
 339
 340    basisCurve->D3 (U, P, V1, V2, V3);
 341    Vec2d V4 = basisCurve->DN (U, 4);
 342    Standard_Integer Index = 2;
 343    while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 344      V1 = basisCurve->DN (U, Index);
 345      Index++;
 346    }
 347    if (Index != 2) {
 348      V2 = basisCurve->DN (U, Index);
 349      V3 = basisCurve->DN (U, Index + 1);
 350      V4 = basisCurve->DN (U, Index + 2);
 351    }
 352    XY Ndir   (V1.Y(), -V1.X());
 353    XY DNdir  (V2.Y(), -V2.X());
 354    XY D2Ndir (V3.Y(), -V3.X());
 355    XY D3Ndir (V4.Y(), -V4.X());
 356    Standard_Real R2  = Ndir.SquareModulus();
 357    Standard_Real R   = Sqrt (R2);
 358    Standard_Real R3  = R2 * R;
 359    Standard_Real R4  = R2 * R2;
 360    Standard_Real R5  = R3 * R2;
 361    Standard_Real R6  = R3 * R3;
 362    Standard_Real R7  = R5 * R2;
 363    Standard_Real Dr  = Ndir.Dot (DNdir);
 364    Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
 365    Standard_Real D3r = Ndir.Dot (D3Ndir) + 3.0 * DNdir.Dot (D2Ndir);
 366    if (R7 <= gp::Resolution()) {
 367      //We try another computation but the stability is not very good
 368      //dixit ISG.
 369      if (R6 <= gp::Resolution()) Geom2d_UndefinedDerivative::Raise();
 370      // V3 = P"' (U) :
 371      D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * offsetValue * Dr / R2));
 372      D3Ndir.Subtract (
 373      (DNdir.Multiplied ((3.0 * offsetValue) * ((D2r/R2) + (Dr*Dr)/R4))));
 374      D3Ndir.Add (Ndir.Multiplied (
 375      (offsetValue * (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r))
 376      ));
 377      D3Ndir.Multiply (offsetValue/R);
 378      V3.Add (Vec2d(D3Ndir));
 379      // V2 = P" (U) :
 380      Standard_Real R4 = R2 * R2;
 381      D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
 382      D2Ndir.Subtract (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
 383      D2Ndir.Multiply (offsetValue / R);
 384      V2.Add (Vec2d(D2Ndir));
 385      // V1 = P' (U) :
 386      DNdir.Multiply(R);
 387      DNdir.Subtract (Ndir.Multiplied (Dr/R));
 388      DNdir.Multiply (offsetValue/R2);
 389      V1.Add (Vec2d(DNdir));
 390    }
 391    else {
 392      // Same computation as IICURV in EUCLID-IS because the stability is
 393      // better.
 394      // V3 = P"' (U) :
 395       D3Ndir.Multiply (offsetValue/R);
 396      D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * offsetValue * Dr / R3));
 397      D3Ndir.Subtract (DNdir.Multiplied (
 398      ((3.0 * offsetValue) * ((D2r/R3) + (Dr*Dr)/R5))) );
 399      D3Ndir.Add (Ndir.Multiplied (
 400      (offsetValue * (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r))
 401      ));
 402      V3.Add (Vec2d(D3Ndir));
 403      // V2 = P" (U) :
 404      D2Ndir.Multiply (offsetValue/R);
 405      D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
 406      D2Ndir.Subtract (Ndir.Multiplied (
 407                       offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))
 408                                       )
 409                      );
 410      V2.Add (Vec2d(D2Ndir));
 411      // V1 = P' (U) :
 412      DNdir.Multiply (offsetValue/R);
 413      DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
 414      V1.Add (Vec2d(DNdir));
 415    }
 416    //P (U) :
 417    Ndir.Multiply (offsetValue/R);
 418    Ndir.Add (P.XY());
 419    P.SetXY (Ndir);
 420 }
 421
 422 //=======================================================================
 423 //function : DN
 424 //purpose  :
 425 //=======================================================================
 426
 427 Vec2d Geom2d_OffsetCurve::DN (const Standard_Real U,
 428                               const Standard_Integer N) const
 429 {
 430   Standard_RangeError_Raise_if (N < 1, "Geom2d_OffsetCurve::DN()");
 431
 432   gp_Vec2d VN, VBidon;
 433   gp_Pnt2d PBidon;
 434   switch (N) {
 435   case 1: D1( U, PBidon, VN); break;
 436   case 2: D2( U, PBidon, VBidon, VN); break;
 437   case 3: D3( U, PBidon, VBidon, VBidon, VN); break;
 438   default:
 439     Standard_NotImplemented::Raise();
 440   }
 441
 442   return VN;
 443 }
 444
 445
 446 //=======================================================================
 447 //function : Value
 448 //purpose  :
 449 //=======================================================================
 450
 451 void Geom2d_OffsetCurve::Value (const Standard_Real U,
 452                                 Pnt2d& P, Pnt2d& Pbasis,
 453                                 Vec2d& V1basis ) const
 454 {
 455
 456   basisCurve->D1 (U, Pbasis, V1basis);
 457   Standard_Integer Index = 2;
 458   while (V1basis.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 459     V1basis = basisCurve->DN (U, Index);
 460     Index++;
 461   }
 462   Standard_Real A = V1basis.Y();
 463   Standard_Real B = - V1basis.X();
 464   Standard_Real R = Sqrt(A*A + B * B);
 465   if (R <= gp::Resolution())  Geom2d_UndefinedValue::Raise();
 466   A = A * offsetValue/R;
 467   B = B * offsetValue/R;
 468   P.SetCoord (A + Pbasis.X(), B + Pbasis.Y());
 469 }
 470
 471
 472 //=======================================================================
 473 //function : D1
 474 //purpose  :
 475 //=======================================================================
 476
 477 void Geom2d_OffsetCurve::D1 (const Standard_Real U,
 478                              Pnt2d& P, Pnt2d& Pbasis,
 479                              Vec2d& V1, Vec2d& V1basis,
 480                              Vec2d& V2basis ) const
 481 {
 482    // P(u) = p(u) + Offset * Ndir / R
 483    // with R = || p' ^ Z|| and Ndir = P' ^ Z
 484
 485    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 486
 487    basisCurve->D2 (U, Pbasis, V1basis, V2basis);
 488    V1 = V1basis;
 489    Vec2d V2 = V2basis;
 490    Standard_Integer Index = 2;
 491    while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 492      V1 = basisCurve->DN (U, Index);
 493      Index++;
 494    }
 495    if (Index != 2) {
 496      V2 = basisCurve->DN (U, Index);
 497    }
 498    XY Ndir (V1.Y(), -V1.X());
 499    XY DNdir (V2.Y(), -V2.X());
 500    Standard_Real R2 = Ndir.SquareModulus();
 501    Standard_Real R = Sqrt (R2);
 502    Standard_Real R3 = R * R2;
 503    Standard_Real Dr = Ndir.Dot (DNdir);
 504    if (R3 <= gp::Resolution()) {
 505       //We try another computation but the stability is not very good.
 506       if (R2 <= gp::Resolution()) { Geom2d_UndefinedDerivative::Raise(); }
 507       DNdir.Multiply(R);
 508       DNdir.Subtract (Ndir.Multiplied (Dr/R));
 509       DNdir.Multiply (offsetValue / R2);
 510       V1.Add (Vec2d(DNdir));
 511    }
 512    else {
 513       // Same computation as IICURV in EUCLID-IS because the stability is
 514       // better
 515       DNdir.Multiply (offsetValue/R);
 516       DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
 517       V1.Add (Vec2d(DNdir));
 518    }
 519    Ndir.Multiply (offsetValue/R);
 520    Ndir.Add (Pbasis.XY());
 521    P.SetXY (Ndir);
 522 }
 523
 524
 525 //=======================================================================
 526 //function : D2
 527 //purpose  :
 528 //=======================================================================
 529
 530 void Geom2d_OffsetCurve::D2 (const Standard_Real U,
 531                              Pnt2d& P, Pnt2d& Pbasis,
 532                              Vec2d& V1, Vec2d& V2,
 533                              Vec2d& V1basis, Vec2d& V2basis,
 534                              Vec2d& V3basis ) const
 535 {
 536    // P(u) = p(u) + Offset * Ndir / R
 537    // with R = || p' ^ Z|| and Ndir = P' ^ Z
 538
 539    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 540
 541    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
 542    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 543
 544   basisCurve->D3 (U, Pbasis, V1basis, V2basis, V3basis);
 545   Standard_Integer Index = 2;
 546   V1 = V1basis;
 547   V2 = V2basis;
 548   Vec2d V3 = V3basis;
 549   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 550     V1 = basisCurve->DN (U, Index);
 551     Index++;
 552   }
 553   if (Index != 2) {
 554     V2 = basisCurve->DN (U, Index);
 555     V3 = basisCurve->DN (U, Index + 1);
 556   }
 557   XY Ndir (V1.Y(), -V1.X());
 558   XY DNdir (V2.Y(), -V2.X());
 559   XY D2Ndir (V3.Y(), -V3.X());
 560   Standard_Real R2  = Ndir.SquareModulus();
 561   Standard_Real R   = Sqrt (R2);
 562   Standard_Real R3  = R2 * R;
 563   Standard_Real R4  = R2 * R2;
 564   Standard_Real R5  = R3 * R2;
 565   Standard_Real Dr  = Ndir.Dot (DNdir);
 566   Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
 567   if (R5 <= gp::Resolution()) {
 568      //We try another computation but the stability is not very good
 569      //dixit ISG.
 570      if (R4 <= gp::Resolution()) { Geom2d_UndefinedDerivative::Raise(); }
 571      // V2 = P" (U) :
 572      Standard_Real R4 = R2 * R2;
 573      D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
 574      D2Ndir.Subtract (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
 575      D2Ndir.Multiply (offsetValue / R);
 576      V2.Add (Vec2d(D2Ndir));
 577      // V1 = P' (U) :
 578      DNdir.Multiply(R);
 579      DNdir.Subtract (Ndir.Multiplied (Dr/R));
 580      DNdir.Multiply (offsetValue/R2);
 581      V1.Add (Vec2d(DNdir));
 582   }
 583   else {
 584      // Same computation as IICURV in EUCLID-IS because the stability is
 585      // better.
 586      // V2 = P" (U) :
 587      D2Ndir.Multiply (offsetValue/R);
 588      D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
 589      D2Ndir.Subtract (Ndir.Multiplied (
 590                       offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))
 591                                       )
 592                      );
 593      V2.Add (Vec2d(D2Ndir));
 594      // V1 = P' (U) :
 595      DNdir.Multiply (offsetValue/R);
 596      DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
 597      V1.Add (Vec2d(DNdir));
 598   }
 599   //P (U) :
 600   Ndir.Multiply (offsetValue/R);
 601   Ndir.Add (Pbasis.XY());
 602   P.SetXY (Ndir);
 603 }
 604
 605 //=======================================================================
 606 //function : FirstParameter
 607 //purpose  :
 608 //=======================================================================
 609
 610 Standard_Real Geom2d_OffsetCurve::FirstParameter () const
 611 {
 612   return basisCurve->FirstParameter();
 613 }
 614
 615 //=======================================================================
 616 //function : LastParameter
 617 //purpose  :
 618 //=======================================================================
 619
 620 Standard_Real Geom2d_OffsetCurve::LastParameter () const
 621 {
 622   return basisCurve->LastParameter();
 623 }
 624
 625
 626 //=======================================================================
 627 //function : Offset
 628 //purpose  :
 629 //=======================================================================
 630
 631 Standard_Real Geom2d_OffsetCurve::Offset () const { return offsetValue; }
 632
 633 //=======================================================================
 634 //function : IsClosed
 635 //purpose  :
 636 //=======================================================================
 637
 638 Standard_Boolean Geom2d_OffsetCurve::IsClosed () const
 639 {
 640   gp_Pnt2d PF, PL;
 641   D0(FirstParameter(),PF);
 642   D0(LastParameter(),PL);
 643   return ( PF.Distance(PL) <= gp::Resolution());
 644 }
 645
 646 //=======================================================================
 647 //function : IsCN
 648 //purpose  :
 649 //=======================================================================
 650
 651 Standard_Boolean Geom2d_OffsetCurve::IsCN (const Standard_Integer N) const
 652 {
 653   Standard_RangeError_Raise_if (N < 0, " " );
 654   return basisCurve->IsCN (N + 1);
 655 }
 656
 657 //=======================================================================
 658 //function : IsPeriodic
 659 //purpose  :
 660 //=======================================================================
 661
 662 Standard_Boolean Geom2d_OffsetCurve::IsPeriodic () const
 663 {
 664   return basisCurve->IsPeriodic();
 665 }
 666
 667 //=======================================================================
 668 //function : Period
 669 //purpose  :
 670 //=======================================================================
 671
 672 Standard_Real Geom2d_OffsetCurve::Period() const
 673 {
 674   return basisCurve->Period();
 675 }
 676
 677 //=======================================================================
 678 //function : Transform
 679 //purpose  :
 680 //=======================================================================
 681
 682 void Geom2d_OffsetCurve::Transform (const Trsf2d& T)
 683 {
 684   basisCurve->Transform (T);
 685   offsetValue *= Abs(T.ScaleFactor());
 686 }
 687
 688
 689 //=======================================================================
 690 //function : TransformedParameter
 691 //purpose  :
 692 //=======================================================================
 693
 694 Standard_Real Geom2d_OffsetCurve::TransformedParameter(const Standard_Real U,
 695                                                         const gp_Trsf2d& T) const
 696 {
 697   return basisCurve->TransformedParameter(U,T);
 698 }
 699
 700 //=======================================================================
 701 //function : ParametricTransformation
 702 //purpose  :
 703 //=======================================================================
 704
 705 Standard_Real Geom2d_OffsetCurve::ParametricTransformation(const gp_Trsf2d& T) const
 706 {
 707   return basisCurve->ParametricTransformation(T);
 708 }