0024377: [Regression-like] OCC 6.7.0 beta contaminates log with more unnecessary...
[occt.git] / src / Extrema / Extrema_FuncExtCC.gxx
index cc3cbfc..19317b8 100755 (executable)
 //                  GetStateNumber()
 
 /*-----------------------------------------------------------------------------
- Fonctions permettant de rechercher une distance extremale entre 2 courbes
+Fonctions permettant de rechercher une distance extremale entre 2 courbes
 C1 et C2 (en partant de points approches C1(u0) et C2(v0)).
- Cette classe herite de math_FunctionSetWithDerivatives et est utilisee par
+Cette classe herite de math_FunctionSetWithDerivatives et est utilisee par
 l'algorithme math_FunctionSetRoot.
- Si on note Du et Dv, les derivees en u et v, les 2 fonctions a annuler sont:
-  { F1(u,v) = (C2(v)-C1(u)).Du(u)/||Du|| }
-  { F2(u,v) = (C2(v)-C1(u)).Dv(v)||Dv|| }
- Si on note Duu et Dvv, les derivees secondes de C1 et C2, les derivees de F1
+Si on note Du et Dv, les derivees en u et v, les 2 fonctions a annuler sont:
+{ F1(u,v) = (C2(v)-C1(u)).Du(u)/||Du|| }
+{ F2(u,v) = (C2(v)-C1(u)).Dv(v)||Dv|| }
+Si on note Duu et Dvv, les derivees secondes de C1 et C2, les derivees de F1
 et F2 sont egales a:
-  { Duf1(u,v) = -||Du|| + C1C2.Duu/||Du||- F1(u,v)*Duu*Du/||Du||**2
-  { Dvf1(u,v) = Dv.Du/||Du|| 
-  { Duf2(u,v) = -Du.Dv/||Dv||
-  { Dvf2(u,v) = ||Dv|| + C2C1.Dvv/||Dv||- F2(u,v)*Dv*Dvv/||Dv||**2
+{ Duf1(u,v) = -||Du|| + C1C2.Duu/||Du||- F1(u,v)*Duu*Du/||Du||**2
+{ Dvf1(u,v) = Dv.Du/||Du|| 
+{ Duf2(u,v) = -Du.Dv/||Dv||
+{ Dvf2(u,v) = ||Dv|| + C2C1.Dvv/||Dv||- F2(u,v)*Dv*Dvv/||Dv||**2
 
 ----------------------------------------------------------------------------*/
+
+#include <Precision.hxx>
+
+
+static const Standard_Real MinTol    = 1.e-20;
+static const Standard_Real TolFactor = 1.e-12;
+static const Standard_Real MinStep   = 1e-7;
+
+static const Standard_Integer MaxOrder = 3;
+
+
+
 //=============================================================================
+Standard_Real Extrema_FuncExtCC::SearchOfTolerance(const Standard_Address C)
+{
+  const Standard_Integer NPoint = 10;
+  Standard_Real aStartParam, anEndParam;
+
+  if(C==myC1)
+  {
+    aStartParam = myUinfium;
+    anEndParam = myUsupremum;
+  }
+  else if(C==myC2)
+  {
+    aStartParam = myVinfium;
+    anEndParam = myVsupremum;
+  }
+  else
+  {
+    //Warning: No curve for tolerance computing!
+    return MinTol;
+  }
+
+  const Standard_Real aStep = (anEndParam - aStartParam)/(Standard_Real)NPoint;
+
+  Standard_Integer aNum = 0;
+  Standard_Real aMax = -Precision::Infinite();  //Maximum value of 1st derivative
+  //(it is computed with using NPoint point)
+
+  do
+  {
+    Standard_Real u = aStartParam + aNum*aStep;    //parameter for every point
+    if(u > anEndParam)
+      u = anEndParam;
 
-#define Tol  1.e-20
-#define delta 1.e-9
+    Pnt Ptemp;  //empty point (is not used below)
+    Vec VDer;   // 1st derivative vector
+    Tool1::D1(*((Curve1*)C), u, Ptemp, VDer);
+    Standard_Real vm = VDer.Magnitude();
+    if(vm > aMax)
+      aMax = vm;
+  }
+  while(++aNum < NPoint+1);
+
+  return Max(aMax*TolFactor,MinTol);
+}
 
+//=============================================================================
 Extrema_FuncExtCC::Extrema_FuncExtCC(const Standard_Real thetol) : myC1 (0), myC2 (0), myTol (thetol)
 {
+  math_Vector V1(1,2), V2(1,2);
+  V1(1) = 0.0;
+  V2(1) = 0.0;
+  V1(2) = 0.0;
+  V2(2) = 0.0;
+  SubIntervalInitialize(V1, V2);
+  myMaxDerivOrderC1 = 0;
+  myTolC1=MinTol;
+  myMaxDerivOrderC2 = 0;
+  myTolC2=MinTol;
 }
 
+//=============================================================================
 Extrema_FuncExtCC::Extrema_FuncExtCC (const Curve1& C1,
-                                     const Curve2& C2,
-                                     const Standard_Real thetol) :
-                                       myC1 ((Standard_Address)&C1), myC2 ((Standard_Address)&C2),
-                                       myTol (thetol)
+                                      const Curve2& C2,
+                                      const Standard_Real thetol) :
+myC1 ((Standard_Address)&C1), myC2 ((Standard_Address)&C2),
+myTol (thetol)
 {
+  math_Vector V1(1,2), V2(1,2);
+  V1(1) = Tool1::FirstParameter(*((Curve1*)myC1));
+  V2(1) = Tool1::LastParameter(*((Curve1*)myC1));
+  V1(2) = Tool2::FirstParameter(*((Curve2*)myC2));
+  V2(2) = Tool2::LastParameter(*((Curve2*)myC2));
+  SubIntervalInitialize(V1, V2);
+
+  switch(Tool1::GetType(*((Curve1*)myC1)))
+  {
+  case GeomAbs_BezierCurve:
+  case GeomAbs_BSplineCurve:
+  case GeomAbs_OtherCurve:
+    myMaxDerivOrderC1 = MaxOrder;
+    myTolC1 = SearchOfTolerance((Standard_Address)&C1);
+    break;
+  default:
+    myMaxDerivOrderC1 = 0;
+    myTolC1=MinTol;
+    break;
+  }
+
+  switch(Tool2::GetType(*((Curve2*)myC2)))
+  {
+  case GeomAbs_BezierCurve:
+  case GeomAbs_BSplineCurve:
+  case GeomAbs_OtherCurve:
+    myMaxDerivOrderC2 = MaxOrder;
+    myTolC2 = SearchOfTolerance((Standard_Address)&C2);
+    break;
+  default:
+    myMaxDerivOrderC2 = 0;
+    myTolC2=MinTol;
+    break;
+  }
 }
 
-Standard_Boolean Extrema_FuncExtCC::Value (const math_Vector& UV, 
-                                          math_Vector&       F)
+//=============================================================================
+void Extrema_FuncExtCC::SetCurve (const Standard_Integer theRank, const Curve1& C)
 {
-  
-  Vec Du, Dv;
+  Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_FuncExtCC::SetCurve()")
 
+    if (theRank == 1)
+    {
+      myC1 = (Standard_Address)&C;
+      switch(/*Tool1::GetType(*((Curve1*)myC1))*/ C.GetType())
+      {
+      case GeomAbs_BezierCurve:
+      case GeomAbs_BSplineCurve:
+      case GeomAbs_OtherCurve:
+        myMaxDerivOrderC1 = MaxOrder;
+        myTolC1 = SearchOfTolerance((Standard_Address)&C);
+        break;
+      default:
+        myMaxDerivOrderC1 = 0;
+        myTolC1=MinTol;
+        break;
+      }
+    }
+    else if (theRank == 2)
+    {
+      myC2 = (Standard_Address)&C;
+      switch(/*Tool2::GetType(*((Curve2*)myC2))*/C.GetType())
+      {
+      case GeomAbs_BezierCurve:
+      case GeomAbs_BSplineCurve:
+      case GeomAbs_OtherCurve:
+        myMaxDerivOrderC2 = MaxOrder;
+        myTolC2 = SearchOfTolerance((Standard_Address)&C);
+        break;
+      default:
+        myMaxDerivOrderC2 = 0;
+        myTolC2=MinTol;
+        break;
+      }
+    }
+}
+
+//=============================================================================
+Standard_Boolean Extrema_FuncExtCC::Value (const math_Vector& UV, math_Vector& F)
+{
   myU = UV(1);
   myV = UV(2);
-  Tool1::D1(*((Curve1*)myC1), myU,myP1,Du);
-  Tool2::D1(*((Curve2*)myC2), myV,myP2,Dv);
+  Tool1::D1(*((Curve1*)myC1), myU,myP1,myDu);
+  Tool2::D1(*((Curve2*)myC2), myV,myP2,myDv);
+
   Vec P1P2 (myP1,myP2);
 
-  Standard_Real Ndu = Du.Magnitude();
-  if (Ndu <= Tol) {
-    Pnt P1, P2;
-    P1 = Tool1::Value(*((Curve1*)myC1), myU-delta);
-    P2 = Tool1::Value(*((Curve1*)myC1), myU+delta);
-    Vec V(P1,P2);
-    Du = V;
-    Ndu = Du.Magnitude();
-    if (Ndu <= Tol) {
-      return Standard_False;
-    }
+  Standard_Real Ndu = myDu.Magnitude();
+
+
+  if(myMaxDerivOrderC1 != 0)
+  {
+    if (Ndu <= myTolC1)
+    {
+      //Derivative is approximated by Taylor-series
+      const Standard_Real DivisionFactor = 1.e-3;
+      Standard_Real du;
+      if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
+        du = 0.0;
+      else
+        du = myUsupremum-myUinfium;
+
+      const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
+
+      Standard_Integer n = 1; //Derivative order
+      Vec V;
+      Standard_Boolean IsDeriveFound;
+
+      do
+      {
+        V = Tool1::DN(*((Curve1*)myC1),myU,++n);
+        Ndu = V.Magnitude();
+        IsDeriveFound = (Ndu > myTolC1);
+      }
+      while(!IsDeriveFound && n < myMaxDerivOrderC1);
+
+      if(IsDeriveFound)
+      {
+        Standard_Real u;
+
+        if(myU-myUinfium < aDelta)
+          u = myU+aDelta;
+        else
+          u = myU-aDelta;
+
+        Pnt P1, P2;
+        Tool1::D0(*((Curve1*)myC1),Min(myU, u),P1);
+        Tool1::D0(*((Curve1*)myC1),Max(myU, u),P2);
+
+        Vec V1(P1,P2);
+        Standard_Real aDirFactor = V.Dot(V1);
+
+        if(aDirFactor < 0.0)
+          myDu = -V;
+        else
+          myDu = V;
+      }//if(IsDeriveFound)
+      else
+      {
+        //Derivative is approximated by three points
+
+        Pnt Ptemp; //(0,0,0)-coordinate
+        Pnt P1, P2, P3;
+        Standard_Boolean IsParameterGrown;
+
+        if(myU-myUinfium < 2*aDelta)
+        {
+          Tool1::D0(*((Curve1*)myC1),myU,P1);
+          Tool1::D0(*((Curve1*)myC1),myU+aDelta,P2);
+          Tool1::D0(*((Curve1*)myC1),myU+2*aDelta,P3);
+          IsParameterGrown = Standard_True;
+        }
+        else
+        {
+          Tool1::D0(*((Curve1*)myC1),myU-2*aDelta,P1);
+          Tool1::D0(*((Curve1*)myC1),myU-aDelta,P2);
+          Tool1::D0(*((Curve1*)myC1),myU,P3);
+          IsParameterGrown = Standard_False;
+        }
+
+        Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
+
+        if(IsParameterGrown)
+          myDu=-3*V1+4*V2-V3;
+        else
+          myDu=V1-4*V2+3*V3;
+      }//else of if(IsDeriveFound)
+      Ndu = myDu.Magnitude();      
+    }//if (Ndu <= myTolC1) condition
+  }//if(myMaxDerivOrder != 0)
+
+  if (Ndu <= MinTol)
+  {
+    //Warning: 1st derivative of C1 is equal to zero!
+    return Standard_False;
   }
 
-  Standard_Real Ndv = Dv.Magnitude();
-  if (Ndv <= Tol) { 
-    // Traitement des singularite, on approche la Tangente
-    // par une corde
-    Pnt P1, P2;
-    P1 = Tool2::Value(*((Curve2*)myC2), myV-delta);
-    P2 = Tool2::Value(*((Curve2*)myC2), myV+delta);
-    Vec V(P1,P2);
-    Dv = V;
-    Ndv = Dv.Magnitude();
-    if (Ndv <= Tol) {
-      return Standard_False;
-    }
+  Standard_Real Ndv = myDv.Magnitude();
+
+  if(myMaxDerivOrderC2 != 0)
+  {
+    if (Ndv <= myTolC2)
+    { 
+      const Standard_Real DivisionFactor = 1.e-3;
+      Standard_Real dv;
+      if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
+        dv = 0.0;
+      else
+        dv = myVsupremum-myVinfium;
+
+      const Standard_Real aDelta = Max(dv*DivisionFactor,MinStep);
+
+      //Derivative is approximated by Taylor-series
+
+      Standard_Integer n = 1; //Derivative order
+      Vec V;
+      Standard_Boolean IsDeriveFound;
+
+      do
+      {
+        V = Tool2::DN(*((Curve2*)myC2),myV,++n);
+        Ndv = V.Magnitude();
+        IsDeriveFound = (Ndv > myTolC2);
+      }
+      while(!IsDeriveFound && n < myMaxDerivOrderC2);
+
+      if(IsDeriveFound)
+      {
+        Standard_Real v;
+
+        if(myV-myVinfium < aDelta)
+          v = myV+aDelta;
+        else
+          v = myV-aDelta;
+
+        Pnt P1, P2;
+        Tool2::D0(*((Curve2*)myC2),Min(myV, v),P1);
+        Tool2::D0(*((Curve2*)myC2),Max(myV, v),P2);
+
+        Vec V1(P1,P2);
+        Standard_Real aDirFactor = V.Dot(V1);
+
+        if(aDirFactor < 0.0)
+          myDv = -V;
+        else
+          myDv = V;
+      }//if(IsDeriveFound)
+      else
+      {
+        //Derivative is approximated by three points
+
+        Pnt Ptemp; //(0,0,0)-coordinate
+        Pnt P1, P2, P3;
+        Standard_Boolean IsParameterGrown;
+
+        if(myV-myVinfium < 2*aDelta)
+        {
+          Tool2::D0(*((Curve2*)myC2),myV,P1);
+          Tool2::D0(*((Curve2*)myC2),myV+aDelta,P2);
+          Tool2::D0(*((Curve2*)myC2),myV+2*aDelta,P3);
+          IsParameterGrown = Standard_True;
+        }
+        else
+        {
+          Tool2::D0(*((Curve2*)myC2),myV-2*aDelta,P1);
+          Tool2::D0(*((Curve2*)myC2),myV-aDelta,P2);
+          Tool2::D0(*((Curve2*)myC2),myV,P3);
+          IsParameterGrown = Standard_False;
+        }
+
+        Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
+
+        if(IsParameterGrown)
+          myDv=-3*V1+4*V2-V3;
+        else
+          myDv=V1-4*V2+3*V3;
+      }//else of if(IsDeriveFound)
+
+      Ndv = myDv.Magnitude();
+    }//if (Ndv <= myTolC2)
+  }//if(myMaxDerivOrder != 0)
+
+  if (Ndv <= MinTol)
+  {
+    //1st derivative of C2 is equal to zero!
+    return Standard_False;
   }
 
-  F(1) = P1P2.Dot(Du)/Ndu;
-  F(2) = P1P2.Dot(Dv)/Ndv;
+  F(1) = P1P2.Dot(myDu)/Ndu;
+  F(2) = P1P2.Dot(myDv)/Ndv;
   return Standard_True;
 }
 //=============================================================================
 
 Standard_Boolean Extrema_FuncExtCC::Derivatives (const math_Vector& UV, 
-                                                math_Matrix& Df)
+                                                 math_Matrix& Df)
 {
   math_Vector F(1,2);
   return Values(UV,F,Df);
@@ -110,72 +401,296 @@ Standard_Boolean Extrema_FuncExtCC::Derivatives (const math_Vector& UV,
 //=============================================================================
 
 Standard_Boolean Extrema_FuncExtCC::Values (const math_Vector& UV, 
-                                           math_Vector& F,
-                                           math_Matrix& Df)
+                                            math_Vector& F,
+                                            math_Matrix& Df)
 {
   myU = UV(1);
   myV = UV(2);
+
+  if(Value(UV, F) == Standard_False) //Computes F, myDu, myDv
+  {
+    //Warning: No function value found!
+    return Standard_False;
+  }//if(Value(UV, F) == Standard_False)
+
   Vec Du, Dv, Duu, Dvv;
   Tool1::D2(*((Curve1*)myC1), myU,myP1,Du,Duu);
   Tool2::D2(*((Curve2*)myC2), myV,myP2,Dv,Dvv);
 
+  //Calling of "Value(...)" function change class member values. 
+  //After running it is necessary to return to previous values.  
+  const Standard_Real myU_old  = myU,    myV_old = myV;
+  const Pnt           myP1_old = myP1,  myP2_old = myP2;
+  const Vec           myDu_old = myDu,  myDv_old = myDv;
+
+
+  //Attention:  aDelta value must be greater than same value for "Value(...)"
+  //            function to avoid of points' collisions.
+
+  const Standard_Real DivisionFactor = 0.01;
+
+  Standard_Real du;
+  if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
+    du = 0.0;
+  else
+    du = myUsupremum-myUinfium;
+
+  const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep);
+
+  Standard_Real dv;
+  if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
+    dv = 0.0;
+  else
+    dv = myVsupremum-myVinfium;
+
+  const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep); 
+
   Vec P1P2 (myP1,myP2);
 
-  Standard_Real Ndu = Du.Magnitude();
-  if (Ndu <= Tol) {
-      Pnt P1, P2;
-    Vec V1;
-    Tool1::D1(*((Curve1*)myC1),myU+delta, P2, Duu);
-    Tool1::D1(*((Curve1*)myC1),myU-delta, P1, V1);
-    Vec V(P1,P2);
-    Du = V;
-    Duu -= V1;
-    Ndu = Du.Magnitude();
-    if (Ndu <= Tol) {
-      return Standard_False;
-    }  
-  }
+  if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
+  {
+    //Derivative is approximated by three points
+
+    math_Vector FF1(1,2), FF2(1,2), FF3(1,2);
+    Standard_Real F1, F2, F3;
+
+    /////////////////////////// Search of DF1_u derivative (begin) ///////////////////
+    if(myU-myUinfium < 2*aDeltaU)
+    {
+      F1=F(1);
+      math_Vector UV2(1,2), UV3(1,2);
+      UV2(1)=myU+aDeltaU;
+      UV2(2)=myV;
+      UV3(1)=myU+2*aDeltaU;
+      UV3(2)=myV;
+      if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
+      {
+        //There are many points close to singularity points and which have zero-derivative.
+        //Try to decrease aDelta variable's value.
+        return Standard_False;
+      }
+
+      F2 = FF2(1);
+      F3 = FF3(1);
+
+      Df(1,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU);
+    }//if(myU-myUinfium < 2*aDeltaU)
+    else
+    {
+      F3 = F(1);
+      math_Vector UV2(1,2), UV1(1,2);
+      UV2(1)=myU-aDeltaU;
+      UV2(2)=myV;
+      UV1(1)=myU-2*aDeltaU;
+      UV1(2)=myV;
+
+      if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
+      {
+        //There are many points close to singularity points and which have zero-derivative.
+        //Try to decrease aDelta variable's value.
+        return Standard_False;
+      }
+
+      F1 = FF1(1);
+      F2 = FF2(1);
+
+      Df(1,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU);
+    }//else of if(myU-myUinfium < 2*aDeltaU) condition
+    /////////////////////////// Search of DF1_u derivative (end) ///////////////////
+
+    //Return to previous values
+    myU = myU_old;
+    myV = myV_old;
+
+    /////////////////////////// Search of DF1_v derivative (begin) ///////////////////
+    if(myV-myVinfium < 2*aDeltaV)
+    {
+      F1=F(1);
+      math_Vector UV2(1,2), UV3(1,2);
+      UV2(1)=myU;
+      UV2(2)=myV+aDeltaV;
+      UV3(1)=myU;
+      UV3(2)=myV+2*aDeltaV;
+
+      if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
+      {
+        //There are many points close to singularity points and which have zero-derivative.
+        //Try to decrease aDelta variable's value.
+        return Standard_False;
+      }
+      F2 = FF2(1);
+      F3 = FF3(1);
+
+      Df(1,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV);
+    }//if(myV-myVinfium < 2*aDeltaV)
+    else
+    {
+      F3 = F(1);
+      math_Vector UV2(1,2), UV1(1,2);
+      UV2(1)=myU;
+      UV2(2)=myV-aDeltaV;
+      UV1(1)=myU;
+      UV1(2)=myV-2*aDeltaV;
+      if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
+      {
+        //There are many points close to singularity points and which have zero-derivative.
+        //Try to decrease aDelta variable's value.
+        return Standard_False;
+      }
+
+      F1 = FF1(1);
+      F2 = FF2(1);
+
+      Df(1,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV);
+    }//else of if(myV-myVinfium < 2*aDeltaV)
+    /////////////////////////// Search of DF1_v derivative (end) ///////////////////
+
+    //Return to previous values
+    myU = myU_old;
+    myV = myV_old;
+    myP1 = myP1_old,  myP2 = myP2_old;
+    myDu = myDu_old,  myDv = myDv_old;
+  }//if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
+  else
+  {
+    const Standard_Real Ndu = myDu.Magnitude();
+    Df(1,1) = - Ndu + (P1P2.Dot(Duu)/Ndu) - F(1)*(myDu.Dot(Duu)/(Ndu*Ndu));
+    Df(1,2) = myDv.Dot(myDu)/Ndu;
+  }//else of if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
+
+  if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
+  {
+    //Derivative is approximated by three points
+
+    math_Vector FF1(1,2), FF2(1,2), FF3(1,2);
+    Standard_Real F1, F2, F3;
+
+    /////////////////////////// Search of DF2_v derivative (begin) ///////////////////
+    if(myV-myVinfium < 2*aDeltaV)
+    {
+      F1=F(2);
+      math_Vector UV2(1,2), UV3(1,2);
+      UV2(1)=myU;
+      UV2(2)=myV+aDeltaV;
+      UV3(1)=myU;
+      UV3(2)=myV+2*aDeltaV;
+
+      if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
+      {
+        //There are many points close to singularity points and which have zero-derivative.
+        //Try to decrease aDelta variable's value.
+        return Standard_False;
+      }
+
+      F2 = FF2(2);
+      F3 = FF3(2);
+
+      Df(2,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV);
+
+    }//if(myV-myVinfium < 2*aDeltaV)
+    else
+    {
+      F3 = F(2);
+      math_Vector UV2(1,2), UV1(1,2);
+      UV2(1)=myU;
+      UV2(2)=myV-aDeltaV;
+      UV1(1)=myU;
+      UV1(2)=myV-2*aDeltaV;
+
+      if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
+      {
+        //There are many points close to singularity points and which have zero-derivative.
+        //Try to decrease aDelta variable's value.
+        return Standard_False;
+      }
+
+      F1 = FF1(2);
+      F2 = FF2(2);
+
+      Df(2,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV);
+    }//else of if(myV-myVinfium < 2*aDeltaV)
+    /////////////////////////// Search of DF2_v derivative (end) ///////////////////
+
+    //Return to previous values
+    myU = myU_old;
+    myV = myV_old;
+
+    /////////////////////////// Search of DF2_u derivative (begin) ///////////////////
+    if(myU-myUinfium < 2*aDeltaU)
+    {
+      F1=F(2);
+      math_Vector UV2(1,2), UV3(1,2);
+      UV2(1)=myU+aDeltaU;
+      UV2(2)=myV;
+      UV3(1)=myU+2*aDeltaU;
+      UV3(2)=myV;
+      if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
+      {
+        //There are many points close to singularity points and which have zero-derivative.
+        //Try to decrease aDelta variable's value.
+        return Standard_False;
+      }
+
+      F2 = FF2(2);
+      F3 = FF3(2);
+
+      Df(2,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU);
+    }//if(myU-myUinfium < 2*aDelta)
+    else
+    {
+      F3 = F(2);
+      math_Vector UV2(1,2), UV1(1,2);
+      UV2(1)=myU-aDeltaU;
+      UV2(2)=myV;
+      UV1(1)=myU-2*aDeltaU;
+      UV1(2)=myV;
+
+      if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
+      {
+        //There are many points close to singularity points
+        //and which have zero-derivative.
+        //Try to decrease aDelta variable's value.
+        return Standard_False;
+      }
+
+      F1 = FF1(2);
+      F2 = FF2(2);
+
+      Df(2,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU);
+    }//else of if(myU-myUinfium < 2*aDeltaU)
+    /////////////////////////// Search of DF2_u derivative (end) ///////////////////
+
+    //Return to previous values
+    myU = myU_old;
+    myV = myV_old;
+    myP1 = myP1_old;
+    myP2 = myP2_old;
+    myDu = myDu_old;
+    myDv = myDv_old;
+  }//if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
+  else
+  {
+    Standard_Real Ndv = myDv.Magnitude();
+    Df(2,2) = Ndv + (P1P2.Dot(Dvv)/Ndv) - F(2)*(myDv.Dot(Dvv)/(Ndv*Ndv));
+    Df(2,1) = -myDu.Dot(myDv)/Ndv;    
+  }//else of if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
 
-  Standard_Real Ndv = Dv.Magnitude();
- if (Ndv <= Tol) {
-    Pnt P1, P2;
-    Vec V1;
-    Tool2::D1(*((Curve2*)myC2),myV+delta, P2, Dvv);
-    Tool2::D1(*((Curve2*)myC2),myV-delta, P1, V1);
-    Vec V(P1,P2);
-    Dv = V;
-    Dvv -= V1;
-    Ndv = Dv.Magnitude();
-    if (Ndv <= Tol) {
-      return Standard_False;
-    }  
-  }
-  
-  F(1) = P1P2.Dot(Du)/Ndu;
-  F(2) = P1P2.Dot(Dv)/Ndv;
-  
-  Df(1,1) = - Ndu + (P1P2.Dot(Duu)/Ndu) - F(1)*(Du.Dot(Duu)/(Ndu*Ndu));
-  Df(1,2) = Dv.Dot(Du)/Ndu;
-  Df(2,1) = -Du.Dot(Dv)/Ndv;
-  Df(2,2) = Ndv + (P1P2.Dot(Dvv)/Ndv) - F(2)*(Dv.Dot(Dvv)/(Ndv*Ndv));
   return Standard_True;
 
-}
+}//end of function
 //=============================================================================
 
 Standard_Integer Extrema_FuncExtCC::GetStateNumber ()
 {
-  Vec Du, Dv;
-  Pnt P1, P2;
-  Tool1::D1(*((Curve1*)myC1), myU, P1, Du);
-  Tool2::D1(*((Curve2*)myC2), myV, P2, Dv);
-  Vec P1P2 (P1, P2);
+  Vec Du (myDu), Dv (myDv);
+  Vec P1P2 (myP1, myP2);
+
   Standard_Real mod = Du.Magnitude();
-  if(mod > Tol) {
+  if(mod > myTolC1) {
     Du /= mod;
   }
   mod = Dv.Magnitude();
-  if(mod > Tol) {
+  if(mod > myTolC2) {
     Dv /= mod;
   }
 
@@ -189,17 +704,20 @@ Standard_Integer Extrema_FuncExtCC::GetStateNumber ()
 //=============================================================================
 
 void Extrema_FuncExtCC::Points (const Standard_Integer N,
-                               POnC& P1,
-                               POnC& P2) const
+                                POnC& P1,
+                                POnC& P2) const
 {
   P1 = myPoints.Value(2*N-1);
   P2 = myPoints.Value(2*N);
 }
 //=============================================================================
 
-
-
-
-
-
-
+void Extrema_FuncExtCC::SubIntervalInitialize(const math_Vector& theInfBound,
+                                              const math_Vector& theSupBound)
+{
+  myUinfium = theInfBound(1);
+  myUsupremum = theSupBound(1);
+  myVinfium = theInfBound(2);
+  myVsupremum = theSupBound(2);
+}
+//=============================================================================