0024682: Move out B-spline cache from curves and surfaces to dedicated classes BSplCL...
[occt.git] / src / Geom2d / Geom2d_OffsetCurve.cxx
index 07126a7..79d86ba 100644 (file)
@@ -23,6 +23,7 @@
 #include <Standard_ConstructionError.hxx>
 #include <Standard_RangeError.hxx>
 #include <Standard_NotImplemented.hxx>
+#include <CSLib_Offset.hxx>
 #include <Geom2d_UndefinedDerivative.hxx>
 #include <Geom2d_UndefinedValue.hxx>
 #include <Geom2d_Line.hxx>
@@ -54,6 +55,14 @@ static const int maxDerivOrder = 3;
 static const Standard_Real MinStep   = 1e-7;
 static const Standard_Real MyAngularToleranceForG1 = Precision::Angular();
 
+static gp_Vec2d dummyDerivative; // used as empty value for unused derivatives in AdjustDerivative
+
+// Recalculate derivatives in the singular point
+// Returns true if the direction of derivatives is changed
+static Standard_Boolean AdjustDerivative(const Handle(Geom2d_Curve)& theCurve, Standard_Integer theMaxDerivative,
+                                         Standard_Real theU, gp_Vec2d& theD1, gp_Vec2d& theD2 = dummyDerivative,
+                                         gp_Vec2d& theD3 = dummyDerivative, gp_Vec2d& theD4 = dummyDerivative);
+
 //=======================================================================
 //function : Copy
 //purpose  : 
@@ -211,163 +220,38 @@ GeomAbs_Shape Geom2d_OffsetCurve::Continuity () const
 //purpose  : 
 //=======================================================================
 
-void Geom2d_OffsetCurve::D0 (const Standard_Real   theU,
-                                  Pnt2d& theP ) const 
+void Geom2d_OffsetCurve::D0 (const Standard_Real theU,
+                                   Pnt2d&        theP) const
   {
-  const Standard_Real aTol = gp::Resolution();
-
   Vec2d vD1;
-
   basisCurve->D1 (theU, theP, vD1);
-  Standard_Real Ndu = vD1.Magnitude();
 
-  if(Ndu <= aTol)
-    {
-    const Standard_Real anUinfium   = basisCurve->FirstParameter();
-    const Standard_Real anUsupremum = basisCurve->LastParameter();
-
-    const Standard_Real DivisionFactor = 1.e-3;
-    Standard_Real du;
-    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) 
-      du = 0.0;
-    else
-      du = anUsupremum-anUinfium;
-      
-    const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
-//Derivative is approximated by Taylor-series
-    
-    Standard_Integer anIndex = 1; //Derivative order
-    Vec2d V;
-    
-    do
-      {
-      V =  basisCurve->DN(theU,++anIndex);
-      Ndu = V.Magnitude();
-      }
-    while((Ndu <= aTol) && anIndex < maxDerivOrder);
-    
-    Standard_Real u;
-    
-    if(theU-anUinfium < aDelta)
-      u = theU+aDelta;
-    else
-      u = theU-aDelta;
-    
-    Pnt2d P1, P2;
-    basisCurve->D0(Min(theU, u),P1);
-    basisCurve->D0(Max(theU, u),P2);
-    
-    Vec2d V1(P1,P2);
-    Standard_Real aDirFactor = V.Dot(V1);
-    
-    if(aDirFactor < 0.0)
-      vD1 = -V;
-    else
-      vD1 = V;
-
-    Ndu = vD1.Magnitude();
-    }//if(Ndu <= aTol)
-
-  if (Ndu <= aTol)
-    Geom2d_UndefinedValue::Raise("Exception: Undefined normal vector "
-          "because tangent vector has zero-magnitude!");
-  
-  Standard_Real A = vD1.Y();
-  Standard_Real B = - vD1.X();
-  A = A * offsetValue/Ndu;
-  B = B * offsetValue/Ndu;
-  theP.SetCoord(theP.X() + A, theP.Y() + B);
-  }
+  Standard_Boolean IsDirectionChange = Standard_False;
+  if(vD1.SquareMagnitude() <= gp::Resolution())
+    IsDirectionChange = AdjustDerivative(basisCurve, 1, theU, vD1);
+
+  CSLib_Offset::D0(theP, vD1, offsetValue, IsDirectionChange, theP);
+}
 
 //=======================================================================
 //function : D1
 //purpose  : 
 //=======================================================================
-void Geom2d_OffsetCurve::D1 (const Standard_Real theU, Pnt2d& P, Vec2d& theV1) const
-  {
+void Geom2d_OffsetCurve::D1 (const Standard_Real theU, Pnt2d& theP, Vec2d& theV1) const
+{
    // P(u) = p(u) + Offset * Ndir / R
    // with R = || p' ^ Z|| and Ndir = P' ^ Z
 
    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 
-  const Standard_Real aTol = gp::Resolution();
-
   Vec2d V2;
-  basisCurve->D2 (theU, P, theV1, V2);
-  
-  if(theV1.Magnitude() <= aTol)
-    {
-    const Standard_Real anUinfium   = basisCurve->FirstParameter();
-    const Standard_Real anUsupremum = basisCurve->LastParameter();
-
-    const Standard_Real DivisionFactor = 1.e-3;
-    Standard_Real du;
-    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) 
-      du = 0.0;
-    else
-      du = anUsupremum-anUinfium;
-      
-    const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
-//Derivative is approximated by Taylor-series
-    
-    Standard_Integer anIndex = 1; //Derivative order
-    Vec2d V;
-    
-    do
-      {
-      V =  basisCurve->DN(theU,++anIndex);
-      }
-    while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder);
-    
-    Standard_Real u;
-    
-    if(theU-anUinfium < aDelta)
-      u = theU+aDelta;
-    else
-      u = theU-aDelta;
-    
-    Pnt2d P1, P2;
-    basisCurve->D0(Min(theU, u),P1);
-    basisCurve->D0(Max(theU, u),P2);
-    
-    Vec2d V1(P1,P2);
-    Standard_Real aDirFactor = V.Dot(V1);
-    
-    if(aDirFactor < 0.0)
-      {
-      theV1 = -V;
-      V2 = - basisCurve->DN (theU, anIndex+1);
-      }
-    else
-      {
-      theV1 = V;
-      V2 = basisCurve->DN (theU, anIndex+1);
-      }
-    }//if(theV1.Magnitude() <= aTol)
-
-  XY Ndir  (theV1.Y(), -theV1.X());
-  XY DNdir (V2.Y(), -V2.X());
-  Standard_Real R2 = Ndir.SquareModulus();
-  Standard_Real R  = Sqrt (R2);
-  Standard_Real R3 = R * R2;
-  Standard_Real Dr = Ndir.Dot (DNdir);
-  if (R3 <= gp::Resolution()) {
-    //We try another computation but the stability is not very good.
-    if (R2 <= gp::Resolution()) Geom2d_UndefinedDerivative::Raise();
-    DNdir.Multiply(R);
-    DNdir.Subtract (Ndir.Multiplied (Dr/R));
-    DNdir.Multiply (offsetValue/R2);
-    theV1.Add (Vec2d(DNdir));
-  }
-  else {
-    // Same computation as IICURV in EUCLID-IS because the stability is
-    // better
-    DNdir.Multiply (offsetValue/R);
-    DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));        
-    theV1.Add (Vec2d(DNdir));
-  }
+  basisCurve->D2 (theU, theP, theV1, V2);
 
-  D0(theU, P);
+  Standard_Boolean IsDirectionChange = Standard_False;
+  if(theV1.SquareMagnitude() <= gp::Resolution())
+    IsDirectionChange = AdjustDerivative(basisCurve, 2, theU, theV1, V2);
+
+  CSLib_Offset::D1(theP, theV1, V2, offsetValue, IsDirectionChange, theP, theV1);
 }
 
 //=======================================================================
@@ -376,8 +260,8 @@ void Geom2d_OffsetCurve::D1 (const Standard_Real theU, Pnt2d& P, Vec2d& theV1) c
 //=======================================================================
 
 void Geom2d_OffsetCurve::D2 (const Standard_Real theU, 
-                                  Pnt2d& P, 
-                                  Vec2d& theV1, Vec2d& V2) const 
+                                   Pnt2d& theP, 
+                                   Vec2d& theV1, Vec2d& theV2) const 
 {
    // P(u) = p(u) + Offset * Ndir / R
    // with R = || p' ^ Z|| and Ndir = P' ^ Z
@@ -388,124 +272,13 @@ void Geom2d_OffsetCurve::D2 (const Standard_Real theU,
    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 
   Vec2d V3;
-  basisCurve->D3 (theU, P, theV1, V2, V3);
-
-  const Standard_Real aTol = gp::Resolution();
+  basisCurve->D3 (theU, theP, theV1, theV2, V3);
 
   Standard_Boolean IsDirectionChange = Standard_False;
+  if(theV1.SquareMagnitude() <= gp::Resolution())
+    IsDirectionChange = AdjustDerivative(basisCurve, 3, theU, theV1, theV2, V3);
 
-  if(theV1.Magnitude() <= aTol)
-    {
-    const Standard_Real anUinfium   = basisCurve->FirstParameter();
-    const Standard_Real anUsupremum = basisCurve->LastParameter();
-
-    const Standard_Real DivisionFactor = 1.e-3;
-    Standard_Real du;
-    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) 
-      du = 0.0;
-    else
-      du = anUsupremum-anUinfium;
-      
-    const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
-//Derivative is approximated by Taylor-series
-    
-    Standard_Integer anIndex = 1; //Derivative order
-    Vec2d V;
-    
-    do
-      {
-      V =  basisCurve->DN(theU,++anIndex);
-      }
-    while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder);
-    
-    Standard_Real u;
-    
-    if(theU-anUinfium < aDelta)
-      u = theU+aDelta;
-    else
-      u = theU-aDelta;
-    
-    Pnt2d P1, P2;
-    basisCurve->D0(Min(theU, u),P1);
-    basisCurve->D0(Max(theU, u),P2);
-    
-    Vec2d V1(P1,P2);
-    Standard_Real aDirFactor = V.Dot(V1);
-    
-    if(aDirFactor < 0.0)
-      {
-      theV1 = -V;
-      V2 = -basisCurve->DN (theU, anIndex+1);
-      V3 = -basisCurve->DN (theU, anIndex + 2);
-
-      IsDirectionChange = Standard_True;
-      }
-    else
-      {
-      theV1 = V;
-      V2 = basisCurve->DN (theU, anIndex+1);
-      V3 = basisCurve->DN (theU, anIndex + 2);
-      }
-    }//if(V1.Magnitude() <= aTol)
-
-  XY Ndir (theV1.Y(), -theV1.X());
-  XY DNdir (V2.Y(), -V2.X());
-  XY D2Ndir (V3.Y(), -V3.X());
-  Standard_Real R2  = Ndir.SquareModulus();
-  Standard_Real R   = Sqrt (R2);
-  Standard_Real R3  = R2 * R;
-  Standard_Real R4  = R2 * R2;
-  Standard_Real R5  = R3 * R2;
-  Standard_Real Dr  = Ndir.Dot (DNdir);
-  Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
-  if (R5 <= gp::Resolution())
-    {
-    //We try another computation but the stability is not very good
-    //dixit ISG.
-    if (R4 <= gp::Resolution())
-      {
-      Geom2d_UndefinedDerivative::Raise();
-      }
-    // V2 = P" (U) :
-    Standard_Real R4 = R2 * R2;
-    D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
-    D2Ndir.Add (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
-    D2Ndir.Multiply (offsetValue / R);
-
-    if(IsDirectionChange)
-      V2=-V2;
-
-    V2.Add (Vec2d(D2Ndir));
-
-    // V1 = P' (U) :
-    DNdir.Multiply(R);
-    DNdir.Subtract (Ndir.Multiplied (Dr/R));
-    DNdir.Multiply (offsetValue/R2);
-    theV1.Add (Vec2d(DNdir));
-    }
-  else
-    {
-    // Same computation as IICURV in EUCLID-IS because the stability is
-    // better.
-    // V2 = P" (U) :
-    D2Ndir.Multiply (offsetValue/R);
-    D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
-    D2Ndir.Add (Ndir.Multiplied 
-                    (offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))));
-
-    if(IsDirectionChange)
-      V2=-V2;
-
-    V2.Add (Vec2d(D2Ndir));
-
-    // V1 = P' (U) 
-    DNdir.Multiply (offsetValue/R);
-    DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));        
-    theV1.Add (Vec2d(DNdir));
-    }
-
-  //P (U) :
-  D0(theU, P);
+  CSLib_Offset::D2(theP, theV1, theV2, V3, offsetValue, IsDirectionChange, theP, theV1, theV2);
 }
 
 
@@ -515,9 +288,9 @@ void Geom2d_OffsetCurve::D2 (const Standard_Real theU,
 //=======================================================================
 
 void Geom2d_OffsetCurve::D3 (const Standard_Real theU, 
-                                   Pnt2d& P, 
-                                   Vec2d& theV1, Vec2d& V2, Vec2d& V3) const {
-
+                                   Pnt2d& theP, 
+                                   Vec2d& theV1, Vec2d& theV2, Vec2d& theV3) const
+{
 
    // P(u) = p(u) + Offset * Ndir / R
    // with R = || p' ^ Z|| and Ndir = P' ^ Z
@@ -532,149 +305,16 @@ void Geom2d_OffsetCurve::D3 (const Standard_Real theU,
    //         (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir +
    //         (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir
 
-  const Standard_Real aTol = gp::Resolution();
+  basisCurve->D3 (theU, theP, theV1, theV2, theV3);
+  Vec2d V4 = basisCurve->DN (theU, 4);
 
   Standard_Boolean IsDirectionChange = Standard_False;
+  if(theV1.SquareMagnitude() <= gp::Resolution())
+    IsDirectionChange = AdjustDerivative(basisCurve, 4, theU, theV1, theV2, theV3, V4);
 
-  basisCurve->D3 (theU, P, theV1, V2, V3);
-  Vec2d V4 = basisCurve->DN (theU, 4);
-
-  if(theV1.Magnitude() <= aTol)
-    {
-    const Standard_Real anUinfium   = basisCurve->FirstParameter();
-    const Standard_Real anUsupremum = basisCurve->LastParameter();
-
-    const Standard_Real DivisionFactor = 1.e-3;
-    Standard_Real du;
-    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) 
-      du = 0.0;
-    else
-      du = anUsupremum-anUinfium;
-      
-    const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
-//Derivative is approximated by Taylor-series
-    
-    Standard_Integer anIndex = 1; //Derivative order
-    Vec2d V;
-    
-    do
-      {
-      V =  basisCurve->DN(theU,++anIndex);
-      }
-    while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder);
-    
-    Standard_Real u;
-    
-    if(theU-anUinfium < aDelta)
-      u = theU+aDelta;
-    else
-      u = theU-aDelta;
-    
-    Pnt2d P1, P2;
-    basisCurve->D0(Min(theU, u),P1);
-    basisCurve->D0(Max(theU, u),P2);
-    
-    Vec2d V1(P1,P2);
-    Standard_Real aDirFactor = V.Dot(V1);
-    
-    if(aDirFactor < 0.0)
-      {
-      theV1 = -V;
-      V2 = -basisCurve->DN (theU, anIndex + 1);
-      V3 = -basisCurve->DN (theU, anIndex + 2);
-      V4 = -basisCurve->DN (theU, anIndex + 3);
-
-      IsDirectionChange = Standard_True;
-      }
-    else
-      {
-      theV1 = V;
-      V2 = basisCurve->DN (theU, anIndex + 1);
-      V3 = basisCurve->DN (theU, anIndex + 2);
-      V4 = basisCurve->DN (theU, anIndex + 3);
-      }
-    }//if(V1.Magnitude() <= aTol)
-
-  XY Ndir   (theV1.Y(), -theV1.X());
-  XY DNdir  (V2.Y(), -V2.X());
-  XY D2Ndir (V3.Y(), -V3.X());
-  XY D3Ndir (V4.Y(), -V4.X());
-  Standard_Real R2  = Ndir.SquareModulus();
-  Standard_Real R   = Sqrt (R2);
-  Standard_Real R3  = R2 * R;
-  Standard_Real R4  = R2 * R2;
-  Standard_Real R5  = R3 * R2;
-  Standard_Real R6  = R3 * R3;
-  Standard_Real R7  = R5 * R2;
-  Standard_Real Dr  = Ndir.Dot (DNdir);
-  Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
-  Standard_Real D3r = Ndir.Dot (D3Ndir) + 3.0 * DNdir.Dot (D2Ndir);
-
-  if (R7 <= gp::Resolution()) 
-    {
-    //We try another computation but the stability is not very good
-    //dixit ISG.
-
-    if (R6 <= gp::Resolution())
-      Geom2d_UndefinedDerivative::Raise();
-
-    // V3 = P"' (U) :
-    D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * offsetValue * Dr / R2));
-    D3Ndir.Subtract (
-      (DNdir.Multiplied ((3.0 * offsetValue) * ((D2r/R2) + (Dr*Dr)/R4))));
-    D3Ndir.Add (Ndir.Multiplied (
-      (offsetValue * (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r))));
-    D3Ndir.Multiply (offsetValue/R);
-
-    if(IsDirectionChange)
-      V3=-V3;
-
-    V3.Add (Vec2d(D3Ndir));
-
-
-    // V2 = P" (U) :
-    Standard_Real R4 = R2 * R2;
-    D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
-    D2Ndir.Subtract (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
-    D2Ndir.Multiply (offsetValue / R);
-    V2.Add (Vec2d(D2Ndir));
-    // V1 = P' (U) :
-    DNdir.Multiply(R);
-    DNdir.Subtract (Ndir.Multiplied (Dr/R));
-    DNdir.Multiply (offsetValue/R2);
-    theV1.Add (Vec2d(DNdir));
-    }
-  else
-    {
-    // Same computation as IICURV in EUCLID-IS because the stability is
-    // better.
-    // V3 = P"' (U) :
-    D3Ndir.Multiply (offsetValue/R);
-    D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * offsetValue * Dr / R3));
-    D3Ndir.Subtract (DNdir.Multiplied (
-      ((3.0 * offsetValue) * ((D2r/R3) + (Dr*Dr)/R5))) );
-    D3Ndir.Add (Ndir.Multiplied (
-      (offsetValue * (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r))));
-
-    if(IsDirectionChange)
-      V3=-V3;
-
-    V3.Add (Vec2d(D3Ndir));
-
-    // V2 = P" (U) :
-    D2Ndir.Multiply (offsetValue/R);
-    D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
-    D2Ndir.Subtract (Ndir.Multiplied (
-      offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))));
-    V2.Add (Vec2d(D2Ndir));
-    // V1 = P' (U) :
-    DNdir.Multiply (offsetValue/R);
-    DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));        
-    theV1.Add (Vec2d(DNdir));
-    }
-  //P (U) :
-  D0(theU, P);
-  }
+  CSLib_Offset::D3(theP, theV1, theV2, theV3, V4, offsetValue, IsDirectionChange,
+                   theP, theV1, theV2, theV3);
+}
 
 //=======================================================================
 //function : DN
@@ -709,10 +349,10 @@ Standard_RangeError_Raise_if (N < 1, "Exception: Geom2d_OffsetCurve::DN(). N<1."
 void Geom2d_OffsetCurve::Value (const Standard_Real theU, 
                                 Pnt2d& theP, Pnt2d& thePbasis,
                                 Vec2d& theV1basis ) const 
-  {
+{
   basisCurve->D1(theU, thePbasis, theV1basis);
   D0(theU,theP);
-  }
+}
 
 
 //=======================================================================
@@ -741,30 +381,8 @@ void Geom2d_OffsetCurve::D1 (const Standard_Real U,
    if (Index != 2) {
      V2 = basisCurve->DN (U, Index);
    }
-   XY Ndir (V1.Y(), -V1.X());
-   XY DNdir (V2.Y(), -V2.X());
-   Standard_Real R2 = Ndir.SquareModulus();
-   Standard_Real R = Sqrt (R2);
-   Standard_Real R3 = R * R2;
-   Standard_Real Dr = Ndir.Dot (DNdir);
-   if (R3 <= gp::Resolution()) {
-      //We try another computation but the stability is not very good.
-      if (R2 <= gp::Resolution()) { Geom2d_UndefinedDerivative::Raise(); }
-      DNdir.Multiply(R);
-      DNdir.Subtract (Ndir.Multiplied (Dr/R));
-      DNdir.Multiply (offsetValue / R2);
-      V1.Add (Vec2d(DNdir));
-   }
-   else {
-      // Same computation as IICURV in EUCLID-IS because the stability is
-      // better
-      DNdir.Multiply (offsetValue/R);
-      DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));        
-      V1.Add (Vec2d(DNdir));
-   }
-   Ndir.Multiply (offsetValue/R);
-   Ndir.Add (Pbasis.XY());
-   P.SetXY (Ndir);
+
+   CSLib_Offset::D1(P, V1, V2, offsetValue, Standard_False, P, V1);
 }
 
 
@@ -800,52 +418,8 @@ void Geom2d_OffsetCurve::D2 (const Standard_Real U,
     V2 = basisCurve->DN (U, Index);
     V3 = basisCurve->DN (U, Index + 1);
   }
-  XY Ndir (V1.Y(), -V1.X());
-  XY DNdir (V2.Y(), -V2.X());
-  XY D2Ndir (V3.Y(), -V3.X());
-  Standard_Real R2  = Ndir.SquareModulus();
-  Standard_Real R   = Sqrt (R2);
-  Standard_Real R3  = R2 * R;
-  Standard_Real R4  = R2 * R2;
-  Standard_Real R5  = R3 * R2;
-  Standard_Real Dr  = Ndir.Dot (DNdir);
-  Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
-  if (R5 <= gp::Resolution()) {
-     //We try another computation but the stability is not very good
-     //dixit ISG.
-     if (R4 <= gp::Resolution()) { Geom2d_UndefinedDerivative::Raise(); }
-     // V2 = P" (U) :
-     Standard_Real R4 = R2 * R2;
-     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
-     D2Ndir.Subtract (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
-     D2Ndir.Multiply (offsetValue / R);
-     V2.Add (Vec2d(D2Ndir));
-     // V1 = P' (U) :
-     DNdir.Multiply(R);
-     DNdir.Subtract (Ndir.Multiplied (Dr/R));
-     DNdir.Multiply (offsetValue/R2);
-     V1.Add (Vec2d(DNdir));
-  }
-  else {
-     // Same computation as IICURV in EUCLID-IS because the stability is
-     // better.
-     // V2 = P" (U) :
-     D2Ndir.Multiply (offsetValue/R);
-     D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
-     D2Ndir.Subtract (Ndir.Multiplied (
-                      offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))
-                                      )
-                     );
-     V2.Add (Vec2d(D2Ndir));
-     // V1 = P' (U) :
-     DNdir.Multiply (offsetValue/R);
-     DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));        
-     V1.Add (Vec2d(DNdir));
-  }
-  //P (U) :
-  Ndir.Multiply (offsetValue/R);
-  Ndir.Add (Pbasis.XY());
-  P.SetXY (Ndir);
+
+  CSLib_Offset::D2(P, V1, V2, V3, offsetValue, Standard_False, P, V1, V2);
 }
 
 //=======================================================================
@@ -961,3 +535,57 @@ GeomAbs_Shape Geom2d_OffsetCurve::GetBasisCurveContinuity() const
 {
   return myBasisCurveContinuity;
 }
+
+
+// ============= Auxiliary functions ===================
+Standard_Boolean AdjustDerivative(const Handle(Geom2d_Curve)& theCurve, Standard_Integer theMaxDerivative,
+                                  Standard_Real theU, gp_Vec2d& theD1, gp_Vec2d& theD2,
+                                  gp_Vec2d& theD3, gp_Vec2d& theD4)
+{
+  static const Standard_Real aTol = gp::Resolution();
+
+  Standard_Boolean IsDirectionChange = Standard_False;
+  const Standard_Real anUinfium   = theCurve->FirstParameter();
+  const Standard_Real anUsupremum = theCurve->LastParameter();
+
+  const Standard_Real DivisionFactor = 1.e-3;
+  Standard_Real du;
+  if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) 
+    du = 0.0;
+  else
+    du = anUsupremum - anUinfium;
+
+  const Standard_Real aDelta = Max(du * DivisionFactor, MinStep);
+
+  //Derivative is approximated by Taylor-series
+  Standard_Integer anIndex = 1; //Derivative order
+  Vec2d V;
+
+  do
+  {
+    V =  theCurve->DN(theU, ++anIndex);
+  }
+  while((V.SquareMagnitude() <= aTol) && anIndex < maxDerivOrder);
+
+  Standard_Real u;
+
+  if(theU-anUinfium < aDelta)
+    u = theU+aDelta;
+  else
+    u = theU-aDelta;
+
+  Pnt2d P1, P2;
+  theCurve->D0(Min(theU, u),P1);
+  theCurve->D0(Max(theU, u),P2);
+
+  Vec2d V1(P1,P2);
+  IsDirectionChange = V.Dot(V1) < 0.0;
+  Standard_Real aSign = IsDirectionChange ? -1.0 : 1.0;
+
+  theD1 = V * aSign;
+  gp_Vec2d* aDeriv[3] = {&theD2, &theD3, &theD4};
+  for (Standard_Integer i = 1; i < theMaxDerivative; i++)
+    *(aDeriv[i-1]) = theCurve->DN(theU, anIndex + i) * aSign;
+
+  return IsDirectionChange;
+}