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-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 // 24-Aug-95 : xab removed C1 and C2 test : appeller  D1 et D2
  23 //             avec discernement !
  24 // 19-09-97  : JPI correction derivee seconde
  25
  26
  27 #include <Geom_OffsetCurve.ixx>
  28
  29
  30 #include <gp.hxx>
  31 #include <gp_XYZ.hxx>
  32 #include <Geom_Line.hxx>
  33 #include <Geom_Circle.hxx>
  34 #include <Geom_Ellipse.hxx>
  35 #include <Geom_Hyperbola.hxx>
  36 #include <Geom_Parabola.hxx>
  37 #include <Geom_BezierCurve.hxx>
  38 #include <Geom_TrimmedCurve.hxx>
  39 #include <Geom_BSplineCurve.hxx>
  40 #include <Geom_Geometry.hxx>
  41
  42 #include <Geom_UndefinedDerivative.hxx>
  43 #include <Geom_UndefinedValue.hxx>
  44 #include <Standard_ConstructionError.hxx>
  45 #include <Standard_RangeError.hxx>
  46 #include <Standard_NotImplemented.hxx>
  47
  48 typedef Geom_OffsetCurve         OffsetCurve;
  49 typedef Handle(Geom_OffsetCurve) Handle(OffsetCurve);
  50 typedef Geom_Curve               Curve;
  51 typedef Handle(Geom_Curve)       Handle(Curve);
  52 typedef Handle(Geom_Geometry)    Handle(Geometry);
  53
  54 typedef gp_Dir  Dir;
  55 typedef gp_Pnt  Pnt;
  56 typedef gp_Trsf Trsf;
  57 typedef gp_Vec  Vec;
  58 typedef gp_XYZ  XYZ;
  59
  60
  61
  62
  63
  64 static const int MaxDegree = 9;
  65    //ordre de derivation maximum pour la recherche de la premiere
  66    //derivee non nulle
  67
  68
  69
  70 //=======================================================================
  71 //function : Copy
  72 //purpose  :
  73 //=======================================================================
  74
  75 Handle(Geom_Geometry) Geom_OffsetCurve::Copy () const {
  76
  77  Handle(OffsetCurve) C;
  78  C = new OffsetCurve (basisCurve, offsetValue, direction);
  79  return C;
  80 }
  81
  82
  83
  84 //=======================================================================
  85 //function : Geom_OffsetCurve
  86 //purpose  :
  87 //=======================================================================
  88
  89 Geom_OffsetCurve::Geom_OffsetCurve (const Handle(Curve)& C,
  90                                     const Standard_Real           Offset,
  91                                     const Dir&           V      )
  92  : direction(V), offsetValue(Offset) {
  93
  94   if (C->DynamicType() == STANDARD_TYPE(Geom_OffsetCurve)) {
  95     Handle(OffsetCurve) OC = Handle(OffsetCurve)::DownCast(C);
  96     if ((OC->BasisCurve())->Continuity() == GeomAbs_C0)
  97       Standard_ConstructionError::Raise();
  98
  99     basisCurve = Handle(Curve)::DownCast((OC->BasisCurve())->Copy());
 100     Standard_Real PrevOff = OC->Offset();
 101     gp_Vec V1(OC->Direction());
 102     gp_Vec V2(direction);
 103     gp_Vec Vdir(PrevOff*V1 + offsetValue*V2);
 104
 105     if (Offset >= 0.) {
 106       offsetValue = Vdir.Magnitude();
 107       direction.SetXYZ(Vdir.XYZ());
 108     } else {
 109       offsetValue = -Vdir.Magnitude();
 110       direction.SetXYZ((-Vdir).XYZ());
 111     }
 112   } else {
 113     if (C->Continuity() == GeomAbs_C0) Standard_ConstructionError::Raise();
 114     basisCurve = Handle(Curve)::DownCast(C->Copy()); // DownCast: 10-03-93
 115   }
 116 }
 117
 118
 119 //=======================================================================
 120 //function : Reverse
 121 //purpose  :
 122 //=======================================================================
 123
 124 void Geom_OffsetCurve::Reverse ()
 125 {
 126   basisCurve->Reverse();
 127   offsetValue = -offsetValue;
 128 }
 129
 130
 131 //=======================================================================
 132 //function : ReversedParameter
 133 //purpose  :
 134 //=======================================================================
 135
 136 Standard_Real Geom_OffsetCurve::ReversedParameter( const Standard_Real U) const
 137 {
 138   return basisCurve->ReversedParameter( U);
 139 }
 140
 141 //=======================================================================
 142 //function : Direction
 143 //purpose  :
 144 //=======================================================================
 145
 146 const gp_Dir& Geom_OffsetCurve::Direction () const
 147   { return direction; }
 148
 149 //=======================================================================
 150 //function : SetDirection
 151 //purpose  :
 152 //=======================================================================
 153
 154 void Geom_OffsetCurve::SetDirection (const Dir& V)
 155   { direction = V; }
 156
 157 //=======================================================================
 158 //function : SetOffsetValue
 159 //purpose  :
 160 //=======================================================================
 161
 162 void Geom_OffsetCurve::SetOffsetValue (const Standard_Real D)
 163   { offsetValue = D; }
 164
 165
 166 //=======================================================================
 167 //function : IsPeriodic
 168 //purpose  :
 169 //=======================================================================
 170
 171 Standard_Boolean Geom_OffsetCurve::IsPeriodic () const
 172 {
 173   return basisCurve->IsPeriodic();
 174 }
 175
 176 //=======================================================================
 177 //function : Period
 178 //purpose  :
 179 //=======================================================================
 180
 181 Standard_Real Geom_OffsetCurve::Period () const
 182 {
 183   return basisCurve->Period();
 184 }
 185
 186 //=======================================================================
 187 //function : SetBasisCurve
 188 //purpose  :
 189 //=======================================================================
 190
 191 void Geom_OffsetCurve::SetBasisCurve (const Handle(Curve)& C) {
 192
 193   if (C->Continuity() == GeomAbs_C0) Standard_ConstructionError::Raise();
 194   basisCurve = Handle(Curve)::DownCast(C->Copy()); // DownCast: 10-03-93
 195 }
 196
 197
 198
 199 //=======================================================================
 200 //function : BasisCurve
 201 //purpose  :
 202 //=======================================================================
 203
 204 Handle(Curve) Geom_OffsetCurve::BasisCurve () const
 205 {
 206   return basisCurve;
 207 }
 208
 209
 210 //=======================================================================
 211 //function : Continuity
 212 //purpose  :
 213 //=======================================================================
 214
 215 GeomAbs_Shape Geom_OffsetCurve::Continuity () const {
 216
 217   GeomAbs_Shape OffsetShape=GeomAbs_C0;
 218   switch (basisCurve->Continuity()) {
 219     case GeomAbs_C0 : OffsetShape = GeomAbs_C0;       break;
 220     case GeomAbs_C1 : OffsetShape = GeomAbs_C0;       break;
 221     case GeomAbs_C2 : OffsetShape = GeomAbs_C1;       break;
 222     case GeomAbs_C3 : OffsetShape = GeomAbs_C2;       break;
 223     case GeomAbs_CN : OffsetShape = GeomAbs_CN;       break;
 224     case GeomAbs_G1 : OffsetShape = GeomAbs_G1;       break;
 225     case GeomAbs_G2 : OffsetShape = GeomAbs_G2;       break;
 226   }
 227   return OffsetShape;
 228 }
 229
 230
 231 //=======================================================================
 232 //function : D0
 233 //purpose  :
 234 //=======================================================================
 235
 236 void Geom_OffsetCurve::D0 (const Standard_Real U, Pnt& P) const
 237 {
 238   gp_Pnt PBasis;
 239   gp_Vec VBasis;
 240   D0(U,P,PBasis,VBasis);
 241 }
 242
 243
 244 //=======================================================================
 245 //function : D1
 246 //purpose  :
 247 //=======================================================================
 248
 249 void Geom_OffsetCurve::D1 (const Standard_Real U, Pnt& P, Vec& V1) const
 250 {
 251   gp_Pnt PBasis;
 252   gp_Vec V1Basis,V2Basis;
 253   D1(U,P,PBasis,V1,V1Basis,V2Basis);
 254 }
 255
 256
 257
 258 //=======================================================================
 259 //function : D2
 260 //purpose  :
 261 //=======================================================================
 262
 263 void Geom_OffsetCurve::D2 (const Standard_Real U, Pnt& P, Vec& V1, Vec& V2) const
 264 {
 265   gp_Pnt PBasis;
 266   gp_Vec V1Basis,V2Basis,V3Basis;
 267   D2(U,P,PBasis,V1,V2,V1Basis,V2Basis,V3Basis);
 268 }
 269
 270
 271 //=======================================================================
 272 //function : D3
 273 //purpose  :
 274 //=======================================================================
 275
 276 void Geom_OffsetCurve::D3 (const Standard_Real U, Pnt& P, Vec& V1, Vec& V2, Vec& V3)
 277 const {
 278
 279
 280    // P(u) = p(u) + Offset * Ndir / R
 281    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
 282
 283    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 284
 285    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
 286    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 287
 288    //P"'(u) = p"'(u) + (Offset / R) * (D3Ndir - (3.0 * Dr/R**2) * D2Ndir -
 289    //         (3.0 * D2r / R2) * DNdir + (3.0 * Dr * Dr / R4) * DNdir -
 290    //         (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir +
 291    //         (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir
 292
 293
 294
 295   basisCurve->D3 (U, P, V1, V2, V3);
 296   Vec V4 = basisCurve->DN (U, 4);
 297   Standard_Integer Index = 2;
 298   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 299     V1 = basisCurve->DN (U, Index);
 300     Index++;
 301   }
 302   if (Index != 2) {
 303     V2 = basisCurve->DN (U, Index);
 304     V3 = basisCurve->DN (U, Index + 1);
 305     V4 = basisCurve->DN (U, Index + 2);
 306   }
 307   XYZ OffsetDir = direction.XYZ();
 308   XYZ Ndir      = (V1.XYZ()).Crossed (OffsetDir);
 309   XYZ DNdir     = (V2.XYZ()).Crossed (OffsetDir);
 310   XYZ D2Ndir    = (V3.XYZ()).Crossed (OffsetDir);
 311   XYZ D3Ndir    = (V4.XYZ()).Crossed (OffsetDir);
 312   Standard_Real R2  = Ndir.SquareModulus();
 313   Standard_Real R   = Sqrt (R2);
 314   Standard_Real R3  = R2 * R;
 315   Standard_Real R4  = R2 * R2;
 316   Standard_Real R5  = R3 * R2;
 317   Standard_Real R6  = R3 * R3;
 318   Standard_Real R7  = R5 * R2;
 319   Standard_Real Dr  = Ndir.Dot (DNdir);
 320   Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
 321   Standard_Real D3r = Ndir.Dot (D3Ndir) + 3.0 * DNdir.Dot (D2Ndir);
 322   if (R7 <= gp::Resolution()) {
 323     if (R6 <= gp::Resolution())  Geom_UndefinedDerivative::Raise();
 324     // V3 = P"' (U) :
 325     D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * Dr / R2));
 326     D3Ndir.Subtract (DNdir.Multiplied (3.0 * ((D2r/R2) + (Dr*Dr/R4))));
 327     D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 -
 328                                  15.0*Dr*Dr*Dr/R6 - D3r));
 329     D3Ndir.Multiply (offsetValue/R);
 330     V3.Add (Vec(D3Ndir));
 331     // V2 = P" (U) :
 332     Standard_Real R4 = R2 * R2;
 333     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
 334     D2Ndir.Subtract (Ndir.Multiplied ((3.0 * Dr * Dr / R4) - (D2r / R2)));
 335     D2Ndir.Multiply (offsetValue / R);
 336     V2.Add (Vec(D2Ndir));
 337     // V1 = P' (U) :
 338     DNdir.Multiply(R);
 339     DNdir.Subtract (Ndir.Multiplied (Dr/R));
 340     DNdir.Multiply (offsetValue/R2);
 341     V1.Add (Vec(DNdir));
 342   }
 343   else {
 344     // V3 = P"' (U) :
 345     D3Ndir.Divide (R);
 346     D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * Dr / R3));
 347     D3Ndir.Subtract (DNdir.Multiplied ((3.0 * ((D2r/R3) + (Dr*Dr)/R5))));
 348     D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 -
 349                 15.0*Dr*Dr*Dr/R7 - D3r));
 350     D3Ndir.Multiply (offsetValue);
 351     V3.Add (Vec(D3Ndir));
 352     // V2 = P" (U) :
 353     D2Ndir.Divide (R);
 354     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R3));
 355     D2Ndir.Subtract (Ndir.Multiplied ((3.0 * Dr * Dr / R5) - (D2r / R3)));
 356     D2Ndir.Multiply (offsetValue);
 357     V2.Add (Vec(D2Ndir));
 358     // V1 = P' (U) :
 359     DNdir.Multiply (offsetValue/R);
 360     DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
 361     V1.Add (Vec(DNdir));
 362   }
 363   //P (U) :
 364   Ndir.Multiply (offsetValue/R);
 365   Ndir.Add (P.XYZ());
 366   P.SetXYZ (Ndir);
 367 }
 368
 369
 370 //=======================================================================
 371 //function : DN
 372 //purpose  :
 373 //=======================================================================
 374
 375 Vec Geom_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const {
 376
 377
 378   if (N < 1)                  Standard_RangeError::Raise();
 379   XYZ OffsetDir = direction.XYZ();
 380   Pnt P;
 381   Vec V1, V2, dummy;
 382   if (N == 1) {
 383     basisCurve->D2 (U, P, V1, V2);
 384     Standard_Integer Index = 2;
 385     while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 386       V1 = basisCurve->DN (U, Index);
 387       Index++;
 388     }
 389     if (Index != 2)  V2 = basisCurve->DN (U, Index);
 390     XYZ Ndir  = (V1.XYZ()).Crossed (OffsetDir);
 391     XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
 392     Standard_Real R2 = Ndir.SquareModulus();
 393     Standard_Real R  = Sqrt (R2);
 394     Standard_Real R3 = R * R2;
 395     Standard_Real Dr = Ndir.Dot (DNdir);
 396     if (R3 <= gp::Resolution()) {
 397       if (R2 <= gp::Resolution())  Geom_UndefinedDerivative::Raise();
 398       Ndir.Multiply (Dr/R);
 399       DNdir.Multiply(R);
 400       DNdir.Subtract (Ndir);
 401       DNdir.Multiply (offsetValue/R2);
 402       V1.Add (Vec(DNdir));
 403     }
 404     else {
 405       Ndir.Multiply (offsetValue * Dr / R3);
 406       DNdir.Multiply (offsetValue/R);
 407       DNdir.Subtract (Ndir);
 408       V1.Add (Vec(DNdir));
 409     }
 410     dummy = V1;
 411   }
 412   else if (N == 2) {
 413     Vec V3;
 414     basisCurve->D3 (U, P, V1, V2, V3);
 415     Standard_Integer Index = 2;
 416     while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 417       V1 = basisCurve->DN (U, Index);
 418       Index++;
 419     }
 420     if (Index != 2) {
 421       V2 = basisCurve->DN (U, Index);
 422       V3 = basisCurve->DN (U, Index + 1);
 423     }
 424     XYZ Ndir   = (V1.XYZ()).Crossed (OffsetDir);
 425     XYZ DNdir  = (V2.XYZ()).Crossed (OffsetDir);
 426     XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir);
 427     Standard_Real R2  = Ndir.SquareModulus();
 428     Standard_Real R   = Sqrt (R2);
 429     Standard_Real R3  = R2 * R;   Standard_Real R4 = R2 * R2;  Standard_Real R5 = R3 * R2;
 430     Standard_Real Dr  = Ndir.Dot (DNdir);
 431     Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
 432     if (R5 <= gp::Resolution()) {
 433       if (R4 <= gp::Resolution())  Geom_UndefinedDerivative::Raise();
 434       Ndir.Multiply ((3.0 * Dr * Dr / R4) - (D2r/R2));
 435       DNdir.Multiply (2.0 * Dr / R2);
 436       D2Ndir.Subtract (DNdir);
 437       D2Ndir.Subtract (Ndir);
 438       D2Ndir.Multiply (offsetValue / R);
 439       V2.Add (Vec(D2Ndir));
 440     }
 441     else {
 442       Ndir.Multiply ((3.0 * Dr * Dr / R4) - (D2r / R2));
 443       DNdir.Multiply (2.0 * Dr / R2);
 444       D2Ndir.Divide (R);
 445       D2Ndir.Subtract (DNdir);
 446       D2Ndir.Subtract (Ndir);
 447       D2Ndir.Multiply (offsetValue);
 448       V2.Add (Vec(D2Ndir));
 449     }
 450     dummy = V2;
 451   }
 452   else if (N == 3) {
 453     Vec V3;
 454     basisCurve->D3 (U, P, V1, V2, V3);
 455     Vec V4 = basisCurve->DN (U, 4);
 456     Standard_Integer Index = 2;
 457     while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 458       V1 = basisCurve->DN (U, Index);
 459       Index++;
 460     }
 461     if (Index != 2) {
 462       V2 = basisCurve->DN (U, Index);
 463       V3 = basisCurve->DN (U, Index + 1);
 464       V4 = basisCurve->DN (U, Index + 2);
 465     }
 466     XYZ Ndir = (V1.XYZ()).Crossed (OffsetDir);
 467     XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
 468     XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir);
 469     XYZ D3Ndir = (V4.XYZ()).Crossed (OffsetDir);
 470     Standard_Real R2  = Ndir.SquareModulus();
 471     Standard_Real R   = Sqrt (R2);  Standard_Real R3 = R2 * R;  Standard_Real R4 = R2 * R2;
 472     Standard_Real R5  = R3 * R2;  Standard_Real R6 = R3 * R3;  Standard_Real R7 = R5 * R2;
 473     Standard_Real Dr  = Ndir.Dot (DNdir);
 474     Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
 475     Standard_Real D3r = Ndir.Dot (D3Ndir) + 3.0 * DNdir.Dot (D2Ndir);
 476     if (R7 <= gp::Resolution()) {
 477       if (R6 <= gp::Resolution()) Geom_UndefinedDerivative::Raise();
 478       D2Ndir.Multiply (3.0 * Dr / R2);
 479       DNdir.Multiply (3.0 * ((D2r/R2) + (Dr*Dr)/R4));
 480       Ndir.Multiply (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r);
 481       D3Ndir.Subtract (D2Ndir);
 482       D3Ndir.Subtract (DNdir);
 483       D3Ndir.Add (Ndir);
 484       D3Ndir.Multiply (offsetValue/R);
 485       V3.Add (Vec(D3Ndir));
 486     }
 487     else {
 488       D2Ndir.Multiply (3.0 * Dr / R3);
 489       DNdir.Multiplied (3.0 * ((D2r/R3) + (Dr*Dr/R5)));
 490       Ndir.Multiply (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r);
 491       D3Ndir.Divide (R);
 492       D3Ndir.Subtract (D2Ndir);
 493       D3Ndir.Subtract (DNdir);
 494       D3Ndir.Add (Ndir);
 495       D3Ndir.Multiply (offsetValue);
 496       V3.Add (Vec(D3Ndir));
 497     }
 498     dummy = V3;
 499   }
 500   else { Standard_NotImplemented::Raise();  }
 501
 502   return dummy;
 503 }
 504
 505 //=======================================================================
 506 //function : D0
 507 //purpose  :
 508 //=======================================================================
 509
 510 void  Geom_OffsetCurve::D0(const Standard_Real U,
 511                            gp_Pnt& P,
 512                            gp_Pnt& Pbasis,
 513                            gp_Vec& V1basis)const
 514 {
 515
 516   basisCurve->D1 (U, Pbasis, V1basis);
 517   Standard_Integer Index = 2;
 518   while (V1basis.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 519     V1basis = basisCurve->DN (U, Index);
 520     Index++;
 521   }
 522   XYZ Ndir = (V1basis.XYZ()).Crossed (direction.XYZ());
 523   Standard_Real R = Ndir.Modulus();
 524   if (R <= gp::Resolution())  Geom_UndefinedValue::Raise();
 525   Ndir.Multiply (offsetValue/R);
 526   Ndir.Add (Pbasis.XYZ());
 527   P.SetXYZ(Ndir);
 528 }
 529
 530 //=======================================================================
 531 //function : D1
 532 //purpose  :
 533 //=======================================================================
 534
 535 void Geom_OffsetCurve::D1 ( const Standard_Real U,
 536                             Pnt& P , Pnt& PBasis ,
 537                             Vec& V1, Vec& V1basis, Vec& V2basis) const {
 538
 539    // P(u) = p(u) + Offset * Ndir / R
 540    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
 541
 542    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 543
 544   basisCurve->D2 (U, PBasis, V1basis, V2basis);
 545   V1 = V1basis;
 546   Vec V2 = V2basis;
 547   Standard_Integer Index = 2;
 548   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 549     V1 = basisCurve->DN (U, Index);
 550     Index++;
 551   }
 552   if (Index != 2)  V2 = basisCurve->DN (U, Index);
 553   XYZ OffsetDir = direction.XYZ();
 554   XYZ Ndir  = (V1.XYZ()).Crossed (OffsetDir);
 555   XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
 556   Standard_Real R2 = Ndir.SquareModulus();
 557   Standard_Real R  = Sqrt (R2);
 558   Standard_Real R3 = R * R2;
 559   Standard_Real Dr = Ndir.Dot (DNdir);
 560   if (R3 <= gp::Resolution()) {
 561     //We try another computation but the stability is not very good.
 562     if (R2 <= gp::Resolution()) Geom_UndefinedDerivative::Raise();
 563     DNdir.Multiply(R);
 564     DNdir.Subtract (Ndir.Multiplied (Dr/R));
 565     DNdir.Multiply (offsetValue/R2);
 566     V1.Add (Vec(DNdir));
 567   }
 568   else {
 569     // Same computation as IICURV in EUCLID-IS because the stability is
 570     // better
 571     DNdir.Multiply (offsetValue/R);
 572     DNdir.Subtract (Ndir.Multiplied (offsetValue * Dr/R3));
 573     V1.Add (Vec(DNdir));
 574   }
 575   Ndir.Multiply (offsetValue/R);
 576   Ndir.Add (PBasis.XYZ());
 577   P.SetXYZ (Ndir);
 578 }
 579
 580
 581 //=======================================================================
 582 //function : D2
 583 //purpose  :
 584 //=======================================================================
 585
 586 void Geom_OffsetCurve::D2 (const Standard_Real U,
 587                            Pnt& P      , Pnt& PBasis ,
 588                            Vec& V1     , Vec& V2     ,
 589                            Vec& V1basis, Vec& V2basis, Vec& V3basis) const {
 590
 591    // P(u) = p(u) + Offset * Ndir / R
 592    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
 593
 594    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 595
 596    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
 597    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 598
 599   basisCurve->D3 (U, PBasis, V1basis, V2basis, V3basis);
 600   Standard_Integer Index = 2;
 601   V1     = V1basis;
 602   V2     = V2basis;
 603   Vec V3 = V3basis;
 604   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 605     V1 = basisCurve->DN (U, Index);
 606     Index++;
 607   }
 608   if (Index != 2) {
 609     V2 = basisCurve->DN (U, Index);
 610     V3 = basisCurve->DN (U, Index + 1);
 611   }
 612   XYZ OffsetDir = direction.XYZ();
 613   XYZ Ndir   = (V1.XYZ()).Crossed (OffsetDir);
 614   XYZ DNdir  = (V2.XYZ()).Crossed (OffsetDir);
 615   XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir);
 616   Standard_Real R2    = Ndir.SquareModulus();
 617   Standard_Real R     = Sqrt (R2);
 618   Standard_Real R3    = R2 * R;
 619   Standard_Real R4    = R2 * R2;
 620   Standard_Real R5    = R3 * R2;
 621   Standard_Real Dr    = Ndir.Dot (DNdir);
 622   Standard_Real D2r   = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
 623   if (R5 <= gp::Resolution()) {
 624     //We try another computation but the stability is not very good
 625     //dixit ISG.
 626     if (R4 <= gp::Resolution())  Geom_UndefinedDerivative::Raise();
 627     // V2 = P" (U) :
 628     Standard_Real R4 = R2 * R2;
 629     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
 630     D2Ndir.Add (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
 631     D2Ndir.Multiply (offsetValue / R);
 632     V2.Add (Vec(D2Ndir));
 633     // V1 = P' (U) :
 634     DNdir.Multiply(R);
 635     DNdir.Subtract (Ndir.Multiplied (Dr/R));
 636     DNdir.Multiply (offsetValue/R2);
 637     V1.Add (Vec(DNdir));
 638   }
 639   else {
 640     // Same computation as IICURV in EUCLID-IS because the stability is
 641     // better.
 642     // V2 = P" (U) :
 643     D2Ndir.Multiply (offsetValue/R);
 644     D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
 645     D2Ndir.Add (Ndir.Multiplied (
 646                      offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))
 647                                      )
 648                     );
 649     V2.Add (Vec(D2Ndir));
 650     // V1 = P' (U) :
 651     DNdir.Multiply (offsetValue/R);
 652     DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
 653     V1.Add (Vec(DNdir));
 654   }
 655   //P (U) :
 656   Ndir.Multiply (offsetValue/R);
 657   Ndir.Add (PBasis.XYZ());
 658   P.SetXYZ (Ndir);
 659 }
 660
 661
 662 //=======================================================================
 663 //function : FirstParameter
 664 //purpose  :
 665 //=======================================================================
 666
 667 Standard_Real Geom_OffsetCurve::FirstParameter () const {
 668
 669    return basisCurve->FirstParameter();
 670 }
 671
 672
 673 //=======================================================================
 674 //function : LastParameter
 675 //purpose  :
 676 //=======================================================================
 677
 678 Standard_Real Geom_OffsetCurve::LastParameter () const {
 679
 680    return basisCurve->LastParameter();
 681 }
 682
 683
 684 //=======================================================================
 685 //function : Offset
 686 //purpose  :
 687 //=======================================================================
 688
 689 Standard_Real Geom_OffsetCurve::Offset () const { return offsetValue; }
 690
 691 //=======================================================================
 692 //function : Value
 693 //purpose  :
 694 //=======================================================================
 695
 696 void Geom_OffsetCurve::Value (
 697 const Standard_Real U, Pnt& P, Pnt& PBasis,  Vec& V1basis) const {
 698
 699   if (basisCurve->Continuity() == GeomAbs_C0)  Geom_UndefinedValue::Raise();
 700   basisCurve->D1 (U, PBasis, V1basis);
 701   Standard_Integer Index = 2;
 702   while (V1basis.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
 703     V1basis = basisCurve->DN (U, Index);
 704     Index++;
 705   }
 706   XYZ Ndir = (V1basis.XYZ()).Crossed (direction.XYZ());
 707   Standard_Real R = Ndir.Modulus();
 708   if (R <= gp::Resolution())  Geom_UndefinedValue::Raise();
 709   Ndir.Multiply (offsetValue/R);
 710   Ndir.Add (PBasis.XYZ());
 711   P.SetXYZ (Ndir);
 712 }
 713
 714
 715 //=======================================================================
 716 //function : IsClosed
 717 //purpose  :
 718 //=======================================================================
 719
 720 Standard_Boolean Geom_OffsetCurve::IsClosed () const
 721 {
 722   gp_Pnt PF,PL;
 723   D0(FirstParameter(),PF);
 724   D0(LastParameter(),PL);
 725   return ( PF.Distance(PL) <= gp::Resolution());
 726 }
 727
 728
 729
 730 //=======================================================================
 731 //function : IsCN
 732 //purpose  :
 733 //=======================================================================
 734
 735 Standard_Boolean Geom_OffsetCurve::IsCN (const Standard_Integer N) const {
 736
 737    Standard_RangeError_Raise_if (N < 0, " ");
 738    return basisCurve->IsCN (N + 1);
 739 }
 740
 741
 742 //=======================================================================
 743 //function : Transform
 744 //purpose  :
 745 //=======================================================================
 746
 747 void Geom_OffsetCurve::Transform (const Trsf& T) {
 748
 749   basisCurve->Transform (T);
 750   direction.Transform(T);
 751   offsetValue *= T.ScaleFactor();
 752 }
 753
 754 //=======================================================================
 755 //function : TransformedParameter
 756 //purpose  :
 757 //=======================================================================
 758
 759 Standard_Real Geom_OffsetCurve::TransformedParameter(const Standard_Real U,
 760                                                      const gp_Trsf& T) const
 761 {
 762   return basisCurve->TransformedParameter(U,T);
 763 }
 764
 765 //=======================================================================
 766 //function : ParametricTransformation
 767 //purpose  :
 768 //=======================================================================
 769
 770 Standard_Real Geom_OffsetCurve::ParametricTransformation(const gp_Trsf& T)
 771 const
 772 {
 773   return basisCurve->ParametricTransformation(T);
 774 }