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     SetBasisCurve (OC->BasisCurve());
  88     offsetValue += OC->Offset();
  89   }
  90   else {
  91     SetBasisCurve(C);
  92   }
  93 }
  94
  95 //=======================================================================
  96 //function : Reverse
  97 //purpose  :
  98 //=======================================================================
  99
 100 void Geom2d_OffsetCurve::Reverse ()
 101 {
 102   basisCurve->Reverse();
 103   offsetValue = -offsetValue;
 104 }
 105
 106 //=======================================================================
 107 //function : ReversedParameter
 108 //purpose  :
 109 //=======================================================================
 110
 111 Standard_Real Geom2d_OffsetCurve::ReversedParameter( const Standard_Real U) const
 112 {
 113   return basisCurve->ReversedParameter( U);
 114 }
 115
 116 //=======================================================================
 117 //function : SetBasisCurve
 118 //purpose  :
 119 //=======================================================================
 120
 121 void Geom2d_OffsetCurve::SetBasisCurve (const Handle(Curve)& C)
 122 {
 123   Handle(Geom2d_Curve) aBasisCurve = Handle(Geom2d_Curve)::DownCast(C->Copy());
 124
 125   // Basis curve must be at least C1
 126   if (aBasisCurve->Continuity() == GeomAbs_C0)
 127   {
 128     // For B-splines it is sometimes possible to increase continuity by removing
 129     // unnecessarily duplicated knots
 130     if (aBasisCurve->IsKind(STANDARD_TYPE(Geom2d_BSplineCurve)))
 131     {
 132       Handle(Geom2d_BSplineCurve) aBCurve = Handle(Geom2d_BSplineCurve)::DownCast(aBasisCurve);
 133       Standard_Integer degree = aBCurve->Degree();
 134       Standard_Real Toler = 1e-7;
 135       Standard_Integer start = aBCurve->IsPeriodic() ? 1 :  aBCurve->FirstUKnotIndex(),
 136                        finish = aBCurve->IsPeriodic() ? aBCurve->NbKnots() :  aBCurve->LastUKnotIndex();
 137       for (Standard_Integer i = start; i <= finish; i++)
 138       {
 139         Standard_Integer mult = aBCurve->Multiplicity(i);
 140         if ( mult == degree )
 141           aBCurve->RemoveKnot(i,degree - 1, Toler);
 142       }
 143     }
 144
 145     // Raise exception if still C0
 146     if (aBasisCurve->Continuity() == GeomAbs_C0)
 147       Standard_ConstructionError::Raise("Offset on C0 curve");
 148   }
 149
 150   basisCurve = aBasisCurve;
 151 }
 152
 153 //=======================================================================
 154 //function : SetOffsetValue
 155 //purpose  :
 156 //=======================================================================
 157
 158 void Geom2d_OffsetCurve::SetOffsetValue (const Standard_Real D) { offsetValue = D; }
 159
 160 //=======================================================================
 161 //function : BasisCurve
 162 //purpose  :
 163 //=======================================================================
 164
 165 Handle(Curve) Geom2d_OffsetCurve::BasisCurve () const
 166 {
 167   return basisCurve;
 168 }
 169
 170 //=======================================================================
 171 //function : Continuity
 172 //purpose  :
 173 //=======================================================================
 174
 175 GeomAbs_Shape Geom2d_OffsetCurve::Continuity () const
 176 {
 177   GeomAbs_Shape OffsetShape=GeomAbs_C0;
 178   switch (basisCurve->Continuity()) {
 179      case GeomAbs_C0 : OffsetShape = GeomAbs_C0;   break;
 180      case GeomAbs_C1 : OffsetShape = GeomAbs_C0;   break;
 181      case GeomAbs_C2 : OffsetShape = GeomAbs_C1;   break;
 182      case GeomAbs_C3 : OffsetShape = GeomAbs_C2;   break;
 183      case GeomAbs_CN : OffsetShape = GeomAbs_CN;   break;
 184      case GeomAbs_G1 : OffsetShape = GeomAbs_G1;   break;
 185      case GeomAbs_G2 : OffsetShape = GeomAbs_G2;   break;
 186   }
 187   return OffsetShape;
 188 }
 189
 190 //=======================================================================
 191 //function : D0
 192 //purpose  :
 193 //=======================================================================
 194
 195 void Geom2d_OffsetCurve::D0 (const Standard_Real   U,
 196                                    Pnt2d& P ) const
 197 {
 198   Vec2d V1;
 199
 200   basisCurve->D1 (U, P, V1);
 201   Standard_Integer Index = 2;
 202   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 203     V1 = basisCurve->DN (U, Index);
 204     Index++;
 205   }
 206   Standard_Real A = V1.Y();
 207   Standard_Real B = - V1.X();
 208   Standard_Real R = Sqrt(A*A + B * B);
 209   if (R <= gp::Resolution())  Geom2d_UndefinedValue::Raise();
 210   A = A * offsetValue/R;
 211   B = B * offsetValue/R;
 212   P.SetCoord(P.X() + A, P.Y() + B);
 213 }
 214
 215 //=======================================================================
 216 //function : D1
 217 //purpose  :
 218 //=======================================================================
 219
 220 void Geom2d_OffsetCurve::D1 (const Standard_Real U, Pnt2d& P, Vec2d& V1) const {
 221
 222    // P(u) = p(u) + Offset * Ndir / R
 223    // with R = || p' ^ Z|| and Ndir = P' ^ Z
 224
 225    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 226
 227   Vec2d V2;
 228   basisCurve->D2 (U, P, V1, V2);
 229   Standard_Integer Index = 2;
 230   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 231     V1 = basisCurve->DN (U, Index);
 232     Index++;
 233   }
 234   if (Index != 2) { V2 = basisCurve->DN (U, Index); }
 235   XY Ndir  (V1.Y(), -V1.X());
 236   XY DNdir (V2.Y(), -V2.X());
 237   Standard_Real R2 = Ndir.SquareModulus();
 238   Standard_Real R  = Sqrt (R2);
 239   Standard_Real R3 = R * R2;
 240   Standard_Real Dr = Ndir.Dot (DNdir);
 241   if (R3 <= gp::Resolution()) {
 242     //We try another computation but the stability is not very good.
 243     if (R2 <= gp::Resolution()) Geom2d_UndefinedDerivative::Raise();
 244     DNdir.Multiply(R);
 245     DNdir.Subtract (Ndir.Multiplied (Dr/R));
 246     DNdir.Multiply (offsetValue/R2);
 247     V1.Add (Vec2d(DNdir));
 248   }
 249   else {
 250     // Same computation as IICURV in EUCLID-IS because the stability is
 251     // better
 252     DNdir.Multiply (offsetValue/R);
 253     DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
 254     V1.Add (Vec2d(DNdir));
 255   }
 256   Ndir.Multiply (offsetValue/R);
 257   Ndir.Add (P.XY());
 258   P.SetXY (Ndir);
 259 }
 260
 261 //=======================================================================
 262 //function : D2
 263 //purpose  :
 264 //=======================================================================
 265
 266 void Geom2d_OffsetCurve::D2 (const Standard_Real U,
 267                                    Pnt2d& P,
 268                                    Vec2d& V1, Vec2d& V2) const
 269 {
 270    // P(u) = p(u) + Offset * Ndir / R
 271    // with R = || p' ^ Z|| and Ndir = P' ^ Z
 272
 273    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 274
 275    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
 276    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 277
 278   Vec2d V3;
 279   basisCurve->D3 (U, P, V1, V2, V3);
 280   Standard_Integer Index = 2;
 281   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 282     V1 = basisCurve->DN (U, Index);
 283     Index++;
 284   }
 285   if (Index != 2) {
 286     V2 = basisCurve->DN (U, Index);
 287     V3 = basisCurve->DN (U, Index + 1);
 288   }
 289   XY Ndir (V1.Y(), -V1.X());
 290   XY DNdir (V2.Y(), -V2.X());
 291   XY D2Ndir (V3.Y(), -V3.X());
 292   Standard_Real R2  = Ndir.SquareModulus();
 293   Standard_Real R   = Sqrt (R2);
 294   Standard_Real R3  = R2 * R;
 295   Standard_Real R4  = R2 * R2;
 296   Standard_Real R5  = R3 * R2;
 297   Standard_Real Dr  = Ndir.Dot (DNdir);
 298   Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
 299   if (R5 <= gp::Resolution()) {
 300     //We try another computation but the stability is not very good
 301     //dixit ISG.
 302     if (R4 <= gp::Resolution()) { Geom2d_UndefinedDerivative::Raise(); }
 303     // V2 = P" (U) :
 304      Standard_Real R4 = R2 * R2;
 305      D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
 306      D2Ndir.Add (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
 307      D2Ndir.Multiply (offsetValue / R);
 308      V2.Add (Vec2d(D2Ndir));
 309      // V1 = P' (U) :
 310      DNdir.Multiply(R);
 311      DNdir.Subtract (Ndir.Multiplied (Dr/R));
 312      DNdir.Multiply (offsetValue/R2);
 313      V1.Add (Vec2d(DNdir));
 314    }
 315    else {
 316      // Same computation as IICURV in EUCLID-IS because the stability is
 317      // better.
 318      // V2 = P" (U) :
 319     D2Ndir.Multiply (offsetValue/R);
 320     D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
 321     D2Ndir.Add (Ndir.Multiplied
 322                      (offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))));
 323     V2.Add (Vec2d(D2Ndir));
 324     // V1 = P' (U)
 325      DNdir.Multiply (offsetValue/R);
 326      DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
 327      V1.Add (Vec2d(DNdir));
 328    }
 329    //P (U) :
 330    Ndir.Multiply (offsetValue/R);
 331    Ndir.Add (P.XY());
 332    P.SetXY (Ndir);
 333 }
 334
 335
 336 //=======================================================================
 337 //function : D3
 338 //purpose  :
 339 //=======================================================================
 340
 341 void Geom2d_OffsetCurve::D3 (const Standard_Real U,
 342                                    Pnt2d& P,
 343                                    Vec2d& V1, Vec2d& V2, Vec2d& V3) const {
 344
 345
 346    // P(u) = p(u) + Offset * Ndir / R
 347    // with R = || p' ^ Z|| and Ndir = P' ^ Z
 348
 349    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 350
 351    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
 352    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 353
 354    //P"'(u) = p"'(u) + (Offset / R) * (D3Ndir - (3.0 * Dr/R**2 ) * D2Ndir -
 355    //         (3.0 * D2r / R2) * DNdir) + (3.0 * Dr * Dr / R4) * DNdir -
 356    //         (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir +
 357    //         (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir
 358
 359
 360
 361    basisCurve->D3 (U, P, V1, V2, V3);
 362    Vec2d V4 = basisCurve->DN (U, 4);
 363    Standard_Integer Index = 2;
 364    while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 365      V1 = basisCurve->DN (U, Index);
 366      Index++;
 367    }
 368    if (Index != 2) {
 369      V2 = basisCurve->DN (U, Index);
 370      V3 = basisCurve->DN (U, Index + 1);
 371      V4 = basisCurve->DN (U, Index + 2);
 372    }
 373    XY Ndir   (V1.Y(), -V1.X());
 374    XY DNdir  (V2.Y(), -V2.X());
 375    XY D2Ndir (V3.Y(), -V3.X());
 376    XY D3Ndir (V4.Y(), -V4.X());
 377    Standard_Real R2  = Ndir.SquareModulus();
 378    Standard_Real R   = Sqrt (R2);
 379    Standard_Real R3  = R2 * R;
 380    Standard_Real R4  = R2 * R2;
 381    Standard_Real R5  = R3 * R2;
 382    Standard_Real R6  = R3 * R3;
 383    Standard_Real R7  = R5 * R2;
 384    Standard_Real Dr  = Ndir.Dot (DNdir);
 385    Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
 386    Standard_Real D3r = Ndir.Dot (D3Ndir) + 3.0 * DNdir.Dot (D2Ndir);
 387    if (R7 <= gp::Resolution()) {
 388      //We try another computation but the stability is not very good
 389      //dixit ISG.
 390      if (R6 <= gp::Resolution()) Geom2d_UndefinedDerivative::Raise();
 391      // V3 = P"' (U) :
 392      D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * offsetValue * Dr / R2));
 393      D3Ndir.Subtract (
 394      (DNdir.Multiplied ((3.0 * offsetValue) * ((D2r/R2) + (Dr*Dr)/R4))));
 395      D3Ndir.Add (Ndir.Multiplied (
 396      (offsetValue * (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r))
 397      ));
 398      D3Ndir.Multiply (offsetValue/R);
 399      V3.Add (Vec2d(D3Ndir));
 400      // V2 = P" (U) :
 401      Standard_Real R4 = R2 * R2;
 402      D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
 403      D2Ndir.Subtract (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
 404      D2Ndir.Multiply (offsetValue / R);
 405      V2.Add (Vec2d(D2Ndir));
 406      // V1 = P' (U) :
 407      DNdir.Multiply(R);
 408      DNdir.Subtract (Ndir.Multiplied (Dr/R));
 409      DNdir.Multiply (offsetValue/R2);
 410      V1.Add (Vec2d(DNdir));
 411    }
 412    else {
 413      // Same computation as IICURV in EUCLID-IS because the stability is
 414      // better.
 415      // V3 = P"' (U) :
 416       D3Ndir.Multiply (offsetValue/R);
 417      D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * offsetValue * Dr / R3));
 418      D3Ndir.Subtract (DNdir.Multiplied (
 419      ((3.0 * offsetValue) * ((D2r/R3) + (Dr*Dr)/R5))) );
 420      D3Ndir.Add (Ndir.Multiplied (
 421      (offsetValue * (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r))
 422      ));
 423      V3.Add (Vec2d(D3Ndir));
 424      // V2 = P" (U) :
 425      D2Ndir.Multiply (offsetValue/R);
 426      D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
 427      D2Ndir.Subtract (Ndir.Multiplied (
 428                       offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))
 429                                       )
 430                      );
 431      V2.Add (Vec2d(D2Ndir));
 432      // V1 = P' (U) :
 433      DNdir.Multiply (offsetValue/R);
 434      DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
 435      V1.Add (Vec2d(DNdir));
 436    }
 437    //P (U) :
 438    Ndir.Multiply (offsetValue/R);
 439    Ndir.Add (P.XY());
 440    P.SetXY (Ndir);
 441 }
 442
 443 //=======================================================================
 444 //function : DN
 445 //purpose  :
 446 //=======================================================================
 447
 448 Vec2d Geom2d_OffsetCurve::DN (const Standard_Real U,
 449                               const Standard_Integer N) const
 450 {
 451   Standard_RangeError_Raise_if (N < 1, "Geom2d_OffsetCurve::DN()");
 452
 453   gp_Vec2d VN, VBidon;
 454   gp_Pnt2d PBidon;
 455   switch (N) {
 456   case 1: D1( U, PBidon, VN); break;
 457   case 2: D2( U, PBidon, VBidon, VN); break;
 458   case 3: D3( U, PBidon, VBidon, VBidon, VN); break;
 459   default:
 460     Standard_NotImplemented::Raise();
 461   }
 462
 463   return VN;
 464 }
 465
 466
 467 //=======================================================================
 468 //function : Value
 469 //purpose  :
 470 //=======================================================================
 471
 472 void Geom2d_OffsetCurve::Value (const Standard_Real U,
 473                                 Pnt2d& P, Pnt2d& Pbasis,
 474                                 Vec2d& V1basis ) const
 475 {
 476
 477   basisCurve->D1 (U, Pbasis, V1basis);
 478   Standard_Integer Index = 2;
 479   while (V1basis.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 480     V1basis = basisCurve->DN (U, Index);
 481     Index++;
 482   }
 483   Standard_Real A = V1basis.Y();
 484   Standard_Real B = - V1basis.X();
 485   Standard_Real R = Sqrt(A*A + B * B);
 486   if (R <= gp::Resolution())  Geom2d_UndefinedValue::Raise();
 487   A = A * offsetValue/R;
 488   B = B * offsetValue/R;
 489   P.SetCoord (A + Pbasis.X(), B + Pbasis.Y());
 490 }
 491
 492
 493 //=======================================================================
 494 //function : D1
 495 //purpose  :
 496 //=======================================================================
 497
 498 void Geom2d_OffsetCurve::D1 (const Standard_Real U,
 499                              Pnt2d& P, Pnt2d& Pbasis,
 500                              Vec2d& V1, Vec2d& V1basis,
 501                              Vec2d& V2basis ) const
 502 {
 503    // P(u) = p(u) + Offset * Ndir / R
 504    // with R = || p' ^ Z|| and Ndir = P' ^ Z
 505
 506    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 507
 508    basisCurve->D2 (U, Pbasis, V1basis, V2basis);
 509    V1 = V1basis;
 510    Vec2d V2 = V2basis;
 511    Standard_Integer Index = 2;
 512    while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 513      V1 = basisCurve->DN (U, Index);
 514      Index++;
 515    }
 516    if (Index != 2) {
 517      V2 = basisCurve->DN (U, Index);
 518    }
 519    XY Ndir (V1.Y(), -V1.X());
 520    XY DNdir (V2.Y(), -V2.X());
 521    Standard_Real R2 = Ndir.SquareModulus();
 522    Standard_Real R = Sqrt (R2);
 523    Standard_Real R3 = R * R2;
 524    Standard_Real Dr = Ndir.Dot (DNdir);
 525    if (R3 <= gp::Resolution()) {
 526       //We try another computation but the stability is not very good.
 527       if (R2 <= gp::Resolution()) { Geom2d_UndefinedDerivative::Raise(); }
 528       DNdir.Multiply(R);
 529       DNdir.Subtract (Ndir.Multiplied (Dr/R));
 530       DNdir.Multiply (offsetValue / R2);
 531       V1.Add (Vec2d(DNdir));
 532    }
 533    else {
 534       // Same computation as IICURV in EUCLID-IS because the stability is
 535       // better
 536       DNdir.Multiply (offsetValue/R);
 537       DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
 538       V1.Add (Vec2d(DNdir));
 539    }
 540    Ndir.Multiply (offsetValue/R);
 541    Ndir.Add (Pbasis.XY());
 542    P.SetXY (Ndir);
 543 }
 544
 545
 546 //=======================================================================
 547 //function : D2
 548 //purpose  :
 549 //=======================================================================
 550
 551 void Geom2d_OffsetCurve::D2 (const Standard_Real U,
 552                              Pnt2d& P, Pnt2d& Pbasis,
 553                              Vec2d& V1, Vec2d& V2,
 554                              Vec2d& V1basis, Vec2d& V2basis,
 555                              Vec2d& V3basis ) const
 556 {
 557    // P(u) = p(u) + Offset * Ndir / R
 558    // with R = || p' ^ Z|| and Ndir = P' ^ Z
 559
 560    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 561
 562    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
 563    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 564
 565   basisCurve->D3 (U, Pbasis, V1basis, V2basis, V3basis);
 566   Standard_Integer Index = 2;
 567   V1 = V1basis;
 568   V2 = V2basis;
 569   Vec2d V3 = V3basis;
 570   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 571     V1 = basisCurve->DN (U, Index);
 572     Index++;
 573   }
 574   if (Index != 2) {
 575     V2 = basisCurve->DN (U, Index);
 576     V3 = basisCurve->DN (U, Index + 1);
 577   }
 578   XY Ndir (V1.Y(), -V1.X());
 579   XY DNdir (V2.Y(), -V2.X());
 580   XY D2Ndir (V3.Y(), -V3.X());
 581   Standard_Real R2  = Ndir.SquareModulus();
 582   Standard_Real R   = Sqrt (R2);
 583   Standard_Real R3  = R2 * R;
 584   Standard_Real R4  = R2 * R2;
 585   Standard_Real R5  = R3 * R2;
 586   Standard_Real Dr  = Ndir.Dot (DNdir);
 587   Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
 588   if (R5 <= gp::Resolution()) {
 589      //We try another computation but the stability is not very good
 590      //dixit ISG.
 591      if (R4 <= gp::Resolution()) { Geom2d_UndefinedDerivative::Raise(); }
 592      // V2 = P" (U) :
 593      Standard_Real R4 = R2 * R2;
 594      D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
 595      D2Ndir.Subtract (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
 596      D2Ndir.Multiply (offsetValue / R);
 597      V2.Add (Vec2d(D2Ndir));
 598      // V1 = P' (U) :
 599      DNdir.Multiply(R);
 600      DNdir.Subtract (Ndir.Multiplied (Dr/R));
 601      DNdir.Multiply (offsetValue/R2);
 602      V1.Add (Vec2d(DNdir));
 603   }
 604   else {
 605      // Same computation as IICURV in EUCLID-IS because the stability is
 606      // better.
 607      // V2 = P" (U) :
 608      D2Ndir.Multiply (offsetValue/R);
 609      D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
 610      D2Ndir.Subtract (Ndir.Multiplied (
 611                       offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))
 612                                       )
 613                      );
 614      V2.Add (Vec2d(D2Ndir));
 615      // V1 = P' (U) :
 616      DNdir.Multiply (offsetValue/R);
 617      DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
 618      V1.Add (Vec2d(DNdir));
 619   }
 620   //P (U) :
 621   Ndir.Multiply (offsetValue/R);
 622   Ndir.Add (Pbasis.XY());
 623   P.SetXY (Ndir);
 624 }
 625
 626 //=======================================================================
 627 //function : FirstParameter
 628 //purpose  :
 629 //=======================================================================
 630
 631 Standard_Real Geom2d_OffsetCurve::FirstParameter () const
 632 {
 633   return basisCurve->FirstParameter();
 634 }
 635
 636 //=======================================================================
 637 //function : LastParameter
 638 //purpose  :
 639 //=======================================================================
 640
 641 Standard_Real Geom2d_OffsetCurve::LastParameter () const
 642 {
 643   return basisCurve->LastParameter();
 644 }
 645
 646
 647 //=======================================================================
 648 //function : Offset
 649 //purpose  :
 650 //=======================================================================
 651
 652 Standard_Real Geom2d_OffsetCurve::Offset () const { return offsetValue; }
 653
 654 //=======================================================================
 655 //function : IsClosed
 656 //purpose  :
 657 //=======================================================================
 658
 659 Standard_Boolean Geom2d_OffsetCurve::IsClosed () const
 660 {
 661   gp_Pnt2d PF, PL;
 662   D0(FirstParameter(),PF);
 663   D0(LastParameter(),PL);
 664   return ( PF.Distance(PL) <= gp::Resolution());
 665 }
 666
 667 //=======================================================================
 668 //function : IsCN
 669 //purpose  :
 670 //=======================================================================
 671
 672 Standard_Boolean Geom2d_OffsetCurve::IsCN (const Standard_Integer N) const
 673 {
 674   Standard_RangeError_Raise_if (N < 0, " " );
 675   return basisCurve->IsCN (N + 1);
 676 }
 677
 678 //=======================================================================
 679 //function : IsPeriodic
 680 //purpose  :
 681 //=======================================================================
 682
 683 Standard_Boolean Geom2d_OffsetCurve::IsPeriodic () const
 684 {
 685   return basisCurve->IsPeriodic();
 686 }
 687
 688 //=======================================================================
 689 //function : Period
 690 //purpose  :
 691 //=======================================================================
 692
 693 Standard_Real Geom2d_OffsetCurve::Period() const
 694 {
 695   return basisCurve->Period();
 696 }
 697
 698 //=======================================================================
 699 //function : Transform
 700 //purpose  :
 701 //=======================================================================
 702
 703 void Geom2d_OffsetCurve::Transform (const Trsf2d& T)
 704 {
 705   basisCurve->Transform (T);
 706   offsetValue *= Abs(T.ScaleFactor());
 707 }
 708
 709
 710 //=======================================================================
 711 //function : TransformedParameter
 712 //purpose  :
 713 //=======================================================================
 714
 715 Standard_Real Geom2d_OffsetCurve::TransformedParameter(const Standard_Real U,
 716                                                         const gp_Trsf2d& T) const
 717 {
 718   return basisCurve->TransformedParameter(U,T);
 719 }
 720
 721 //=======================================================================
 722 //function : ParametricTransformation
 723 //purpose  :
 724 //=======================================================================
 725
 726 Standard_Real Geom2d_OffsetCurve::ParametricTransformation(const gp_Trsf2d& T) const
 727 {
 728   return basisCurve->ParametricTransformation(T);
 729 }