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