0026980: Intersection algorithm spends much system time and system memory
[occt.git] / src / IntWalk / IntWalk_PWalking.cxx
index 539b1ee..34c895c 100644 (file)
 // Alternatively, this file may be used under the terms of Open CASCADE
 // commercial license or contractual agreement.
 
-#include <IntWalk_PWalking.ixx>
-
-#include <IntWalk_StatusDeflection.hxx>
-
-#include <TColgp_Array1OfPnt.hxx>
-#include <TColStd_Array1OfReal.hxx>
-
-#include <IntImp_ComputeTangence.hxx>
 
 #include <Adaptor3d_HSurface.hxx>
 #include <Adaptor3d_HSurfaceTool.hxx>
-
-#include <Precision.hxx>
-
-#include <math_FunctionSetRoot.hxx>
+#include <Extrema_GenLocateExtPS.hxx>
 #include <Geom_Surface.hxx>
-
-#include <Standard_Failure.hxx>
+#include <gp_Dir.hxx>
+#include <gp_Pnt.hxx>
 #include <gp_Pnt2d.hxx>
+#include <IntImp_ComputeTangence.hxx>
+#include <IntSurf_LineOn2S.hxx>
+#include <IntSurf_PntOn2S.hxx>
+#include <IntWalk_PWalking.hxx>
+#include <IntWalk_StatusDeflection.hxx>
+#include <math_FunctionSetRoot.hxx>
+#include <Precision.hxx>
+#include <Standard_Failure.hxx>
+#include <Standard_OutOfRange.hxx>
+#include <StdFail_NotDone.hxx>
+#include <TColgp_Array1OfPnt.hxx>
+#include <TColStd_Array1OfReal.hxx>
 
 //==================================================================================
-// function : IntWalk_PWalking::IntWalk_PWalking
-// purpose  :
-// estimate of max step : To avoid abrupt changes 
-// during change of isos 
+// function : ComputePasInit
+// purpose  : estimate of max step : To avoid abrupt changes during change of isos 
 //==================================================================================
-void ComputePasInit(Standard_Real *pasuv,
-                    Standard_Real Um1,Standard_Real UM1,
-                    Standard_Real Vm1,Standard_Real VM1,
-                    Standard_Real Um2,Standard_Real UM2,
-                    Standard_Real Vm2,Standard_Real VM2,
-                    Standard_Real _Um1,Standard_Real _UM1,
-                    Standard_Real _Vm1,Standard_Real _VM1,
-                    Standard_Real _Um2,Standard_Real _UM2,
-                    Standard_Real _Vm2,Standard_Real _VM2,
-                    const Handle(Adaptor3d_HSurface)& ,
-                    const Handle(Adaptor3d_HSurface)& ,
-                    const Standard_Real Increment) 
-{ 
-  Standard_Real du1=Abs(UM1-Um1);
-  Standard_Real dv1=Abs(VM1-Vm1);
-  Standard_Real du2=Abs(UM2-Um2);
-  Standard_Real dv2=Abs(VM2-Vm2);
-
-  Standard_Real _du1=Abs(_UM1-_Um1);
-  Standard_Real _dv1=Abs(_VM1-_Vm1);
-  Standard_Real _du2=Abs(_UM2-_Um2);
-  Standard_Real _dv2=Abs(_VM2-_Vm2);
+void IntWalk_PWalking::ComputePasInit(const Standard_Real theDeltaU1,
+                                      const Standard_Real theDeltaV1,
+                                      const Standard_Real theDeltaU2,
+                                      const Standard_Real theDeltaV2)
+{
+  const Standard_Real aRangePart = 0.01;
+  const Standard_Real Increment = 2.0*pasMax;
+  const Handle(Adaptor3d_HSurface)& 
+          Caro1 = myIntersectionOn2S.Function().AuxillarSurface1();
+  const Handle(Adaptor3d_HSurface)& 
+          Caro2 = myIntersectionOn2S.Function().AuxillarSurface2();
+
+  const Standard_Real aDeltaU1=Abs(UM1-Um1);
+  const Standard_Real aDeltaV1=Abs(VM1-Vm1);
+  const Standard_Real aDeltaU2=Abs(UM2-Um2);
+  const Standard_Real aDeltaV2=Abs(VM2-Vm2);
 
   //-- limit the reduction of uv box estimate to 0.01 natural box
-  //--  du1 : On box of Inter
-  //-- _du1 : On parametric space
-  if(_du1<1e50 && du1<0.01*_du1) du1=0.01*_du1;
-  if(_dv1<1e50 && dv1<0.01*_dv1) dv1=0.01*_dv1;
-  if(_du2<1e50 && du2<0.01*_du2) du2=0.01*_du2;
-  if(_dv2<1e50 && dv2<0.01*_dv2) dv2=0.01*_dv2;
-
-  pasuv[0]=Increment*du1;
-  pasuv[1]=Increment*dv1;
-  pasuv[2]=Increment*du2;
-  pasuv[3]=Increment*dv2;
+  //--  theDeltaU1 : On box of Inter
+  //-- aDeltaU1 : On parametric space
+  if(!Precision::IsInfinite(aDeltaU1))
+    pasuv[0]=Max(Increment*Max(theDeltaU1, aRangePart*aDeltaU1), pasuv[0]);
+  else
+    pasuv[0]=Max(Increment*theDeltaU1, pasuv[0]);
+
+  if(!Precision::IsInfinite(aDeltaV1))
+    pasuv[1]=Max(Increment*Max(theDeltaV1, aRangePart*aDeltaV1), pasuv[1]);
+  else
+    pasuv[1]=Max(Increment*theDeltaV1, pasuv[1]);
+
+  if(!Precision::IsInfinite(aDeltaU2))
+    pasuv[2]=Max(Increment*Max(theDeltaU2, aRangePart*aDeltaU2), pasuv[2]);
+  else
+    pasuv[2]=Max(Increment*theDeltaU2, pasuv[2]);
+
+  if(!Precision::IsInfinite(aDeltaV2))
+    pasuv[3]=Max(Increment*Max(theDeltaV2, aRangePart*aDeltaV2), pasuv[3]);
+  else
+    pasuv[3]=Max(Increment*theDeltaV2, pasuv[3]);
+
+  const Standard_Real ResoU1tol = Adaptor3d_HSurfaceTool::UResolution(Caro1, tolconf);
+  const Standard_Real ResoV1tol = Adaptor3d_HSurfaceTool::VResolution(Caro1, tolconf);
+  const Standard_Real ResoU2tol = Adaptor3d_HSurfaceTool::UResolution(Caro2, tolconf);
+  const Standard_Real ResoV2tol = Adaptor3d_HSurfaceTool::VResolution(Caro2, tolconf);
+
+  myStepMin[0] = Max(myStepMin[0], 2.0*ResoU1tol);
+  myStepMin[1] = Max(myStepMin[1], 2.0*ResoV1tol);
+  myStepMin[2] = Max(myStepMin[2], 2.0*ResoU2tol);
+  myStepMin[3] = Max(myStepMin[3], 2.0*ResoV2tol);
+
+  for(Standard_Integer i  = 0; i < 4; i++)
+  {
+    pasuv[i]=Max(myStepMin[i], pasuv[i]);
+  }
 }
 
 //=======================================================================
@@ -251,6 +269,7 @@ done(Standard_True),
 close(Standard_False),
 fleche(Deflection),
 tolconf(Epsilon),
+myTolTang(TolTangency),
 sensCheminement(1),
 myIntersectionOn2S(Caro1,Caro2,TolTangency),
 STATIC_BLOCAGE_SUR_PAS_TROP_GRAND(0),
@@ -361,6 +380,11 @@ STATIC_PRECEDENT_INFLEXION(0)
     }
   }
 
+  myStepMin[0] = 100.0*ResoU1;
+  myStepMin[1] = 100.0*ResoV1;
+  myStepMin[2] = 100.0*ResoU2;
+  myStepMin[3] = 100.0*ResoV2;
+
   //-- ComputePasInit(pasuv,Um1,UM1,Vm1,VM1,Um2,UM2,Vm2,VM2,Caro1,Caro2);
 
   for (Standard_Integer i = 0; i<=3;i++) {
@@ -391,7 +415,8 @@ done(Standard_True),
 close(Standard_False),
 fleche(Deflection),
 tolconf(Epsilon),
-sensCheminement(1),       
+myTolTang(TolTangency),
+sensCheminement(1),
 myIntersectionOn2S(Caro1,Caro2,TolTangency),
 STATIC_BLOCAGE_SUR_PAS_TROP_GRAND(0),
 STATIC_PRECEDENT_INFLEXION(0)
@@ -527,6 +552,12 @@ STATIC_PRECEDENT_INFLEXION(0)
   if(ResoV1>0.0001*pasuv[1]) ResoV1=0.00001*pasuv[1];
   if(ResoU2>0.0001*pasuv[2]) ResoU2=0.00001*pasuv[2];
   if(ResoV2>0.0001*pasuv[3]) ResoV2=0.00001*pasuv[3];
+  
+  myStepMin[0] = 100.0*ResoU1;
+  myStepMin[1] = 100.0*ResoV1;
+  myStepMin[2] = 100.0*ResoU2;
+  myStepMin[3] = 100.0*ResoV2;
+
   //
   TColStd_Array1OfReal Par(1,4);
   Par(1) = U1;
@@ -575,6 +606,96 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep)
 {
   Perform(ParDep,Um1,Vm1,Um2,Vm2,UM1,VM1,UM2,VM2);
 }
+
+//=======================================================================
+//function : SQDistPointSurface
+//purpose  : Returns square distance between thePnt and theSurf.
+//            (theU0, theV0) is initial point for extrema
+//=======================================================================
+static Standard_Real SQDistPointSurface(const gp_Pnt &thePnt,
+                                        const Adaptor3d_Surface& theSurf,
+                                        const Standard_Real theU0,
+                                        const Standard_Real theV0)
+{
+  const Extrema_GenLocateExtPS aExtPS(thePnt, theSurf, theU0, theV0,
+                      Precision::PConfusion(), Precision::PConfusion());
+  if(!aExtPS.IsDone())
+    return RealLast();
+  
+  return aExtPS.SquareDistance();
+}
+
+//==================================================================================
+// function : IsTangentExtCheck
+// purpose  : Additional check if the surfaces are tangent.
+//            Checks if any point in one surface lie in another surface
+//            (with given tolerance)
+//==================================================================================
+static Standard_Boolean IsTangentExtCheck(const Handle(Adaptor3d_HSurface)& theSurf1,
+                                          const Handle(Adaptor3d_HSurface)& theSurf2,
+                                          const Standard_Real theU10,
+                                          const Standard_Real theV10,
+                                          const Standard_Real theU20,
+                                          const Standard_Real theV20,
+                                          const Standard_Real theToler,
+                                          const Standard_Real theArrStep[])
+{
+  {
+    gp_Pnt aPt;
+    gp_Vec aDu1, aDv1, aDu2, aDv2;
+    theSurf1->D1(theU10, theV10, aPt, aDu1, aDv1);
+    theSurf2->D1(theU20, theV20, aPt, aDu2, aDv2);
+
+    const gp_Vec  aN1(aDu1.Crossed(aDv1)),
+                  aN2(aDu2.Crossed(aDv2));
+    const Standard_Real aDP = aN1.Dot(aN2),
+                        aSQ1 = aN1.SquareMagnitude(),
+                        aSQ2 = aN2.SquareMagnitude();
+
+    if((aSQ1 < RealSmall()) || (aSQ2 < RealSmall()))
+      return Standard_True; //Tangent
+
+    if(aDP*aDP < 0.9998*aSQ1*aSQ2)
+    {//cos(ang N1<->N2) < 0.9999
+      return Standard_False; //Not tangent
+    }
+  }
+
+  //For two faces (2^2 = 4)
+  const Standard_Real aSQToler = 4.0*theToler*theToler;
+  const Standard_Integer aNbItems = 4;
+  const Standard_Real aParUS1[aNbItems] = { theU10 + theArrStep[0],
+                                            theU10 - theArrStep[0],
+                                            theU10, theU10};
+  const Standard_Real aParVS1[aNbItems] = { theV10, theV10,
+                                            theV10 + theArrStep[1],
+                                            theV10 - theArrStep[1]};
+  const Standard_Real aParUS2[aNbItems] = { theU20 + theArrStep[2],
+                                            theU20 - theArrStep[2],
+                                            theU20, theU20};
+  const Standard_Real aParVS2[aNbItems] = { theV20, theV20,
+                                            theV20 + theArrStep[3],
+                                            theV20 - theArrStep[3]};
+
+  for(Standard_Integer i = 0; i < aNbItems; i++)
+  {
+    gp_Pnt aP(theSurf1->Value(aParUS1[i], aParVS1[i]));
+    const Standard_Real aSqDist = SQDistPointSurface(aP, theSurf2->Surface(), theU20, theV20);
+    if(aSqDist > aSQToler)
+      return Standard_False;
+  }
+
+  for(Standard_Integer i = 0; i < aNbItems; i++)
+  {
+    gp_Pnt aP(theSurf2->Value(aParUS2[i], aParVS2[i]));
+    const Standard_Real aSqDist = SQDistPointSurface(aP, theSurf1->Surface(), theU10, theV10);
+    if(aSqDist > aSQToler)
+      return Standard_False;
+  }
+
+  return Standard_True;
+}
+
 //==================================================================================
 // function : Perform
 // purpose  : 
@@ -613,22 +734,8 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
   const Standard_Real ULast2  = Adaptor3d_HSurfaceTool::LastUParameter (Caro2);
   const Standard_Real VLast2  = Adaptor3d_HSurfaceTool::LastVParameter (Caro2);
   //
-  ComputePasInit(pasuv,u1min,u1max,v1min,v1max,u2min,u2max,v2min,v2max,
-    Um1,UM1,Vm1,VM1,Um2,UM2,Vm2,VM2,Caro1,Caro2,pasMax+pasMax);
-  //
-  if(pasuv[0]<100.0*ResoU1) {
-    pasuv[0]=100.0*ResoU1; 
-  }
-  if(pasuv[1]<100.0*ResoV1) {
-    pasuv[1]=100.0*ResoV1; 
-  }
-  if(pasuv[2]<100.0*ResoU2) {
-    pasuv[2]=100.0*ResoU2;
-  }
-  if(pasuv[3]<100.0*ResoV2) {
-    pasuv[3]=100.0*ResoV2;
-  }
-  //
+  ComputePasInit(u1max - u1min,v1max - v1min,u2max - u2min,v2max - v2min);
+
   for (Standard_Integer i=0; i<4; ++i)
   {
     if(pasuv[i]>10)
@@ -701,6 +808,10 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
   Standard_Boolean bTestFirstPoint = Standard_True;
 
   previousPoint.Parameters(Param(1),Param(2),Param(3),Param(4));
+
+  if(IsTangentExtCheck(Caro1, Caro2, Param(1), Param(2), Param(3), Param(4), myTolTang, pasuv))
+    return;
+
   AddAPoint(line,previousPoint);
   //
   IntWalk_StatusDeflection Status = IntWalk_OK;
@@ -710,6 +821,11 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
   Standard_Integer LevelOfPointConfondu = 0; 
   Standard_Integer LevelOfIterWithoutAppend = -1;
   //
+
+  const Standard_Real aTol[4] = { Epsilon(u1max - u1min),
+                                  Epsilon(v1max - v1min),
+                                  Epsilon(u2max - u2min),
+                                  Epsilon(v2max - v2min)};
   Arrive = Standard_False;
   while(!Arrive) //010
   {
@@ -792,8 +908,23 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
 
       if (myIntersectionOn2S.IsDone() && !myIntersectionOn2S.IsEmpty())
       {
+        //If we go along any surface boundary then it is possible 
+        //to find "outboundaried" point.
+        //Nevertheless, if this deflection is quite small, we will be
+        //able to adjust this point to the boundary.
+
         Standard_Real aNewPnt[4], anAbsParamDist[4];
         myIntersectionOn2S.Point().Parameters(aNewPnt[0], aNewPnt[1], aNewPnt[2], aNewPnt[3]);
+        const Standard_Real aParMin[4] = {u1min, v1min, u2min, v2min};
+        const Standard_Real aParMax[4] = {u1max, v1max, u2max, v2max};
+
+        for(Standard_Integer i = 0; i < 4; i++)
+        {
+          if(Abs(aNewPnt[i] - aParMin[i]) < aTol[i])
+            aNewPnt[i] = aParMin[i];
+          else if(Abs(aNewPnt[i] - aParMax[i]) < aTol[i])
+            aNewPnt[i] = aParMax[i];
+        }
 
         if (aNewPnt[0] < u1min || aNewPnt[0] > u1max ||
             aNewPnt[1] < v1min || aNewPnt[1] > v1max ||
@@ -803,6 +934,11 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
           break; // Out of borders, handle this later.
         }
 
+        myIntersectionOn2S.ChangePoint().SetValue(aNewPnt[0],
+                                                  aNewPnt[1],
+                                                  aNewPnt[2],
+                                                  aNewPnt[3]);
+
         anAbsParamDist[0] = Abs(Param(1) - dP1 - aNewPnt[0]);
         anAbsParamDist[1] = Abs(Param(2) - dP2 - aNewPnt[1]);
         anAbsParamDist[2] = Abs(Param(3) - dP3 - aNewPnt[2]);
@@ -885,7 +1021,7 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
             LevelOfEmptyInmyIntersectionOn2S=0;
             if(LevelOfIterWithoutAppend < 10)
             {
-              Status = TestDeflection();
+              Status = TestDeflection(ChoixIso);
             }                  
             else
             {
@@ -1003,28 +1139,18 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
 
             if(LevelOfIterWithoutAppend > 5)
             {
-              if(pasSav[0]<pasInit[0])
+              for (Standard_Integer i = 0; i < 4; i++)
               {
-                pasInit[0]-=(pasInit[0]-pasSav[0])*0.25;
-                LevelOfIterWithoutAppend=0;
-              }
-
-              if(pasSav[1]<pasInit[1])
-              {
-                pasInit[1]-=(pasInit[1]-pasSav[1])*0.25;
-                LevelOfIterWithoutAppend=0;
-              }
+                if (pasSav[i] > pasInit[i])
+                  continue;
 
-              if(pasSav[2]<pasInit[2])
-              {
-                pasInit[2]-=(pasInit[2]-pasSav[2])*0.25;
-                LevelOfIterWithoutAppend=0;
-              }
+                const Standard_Real aDelta = (pasInit[i]-pasSav[i])*0.25;
 
-              if(pasSav[3]<pasInit[3])
-              {
-                pasInit[3]-=(pasInit[3]-pasSav[3])*0.25;
-                LevelOfIterWithoutAppend=0;
+                if(aDelta > Epsilon(pasInit[i]))
+                {
+                  pasInit[i] -= aDelta;
+                  LevelOfIterWithoutAppend=0;
+                }
               }
             }
 
@@ -1042,28 +1168,13 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
               {
                 pastroppetit=Standard_True;
 
-                if(pasuv[0]<pasInit[0])
-                {
-                  pasuv[0]+=(pasInit[0]-pasuv[0])*0.25;
-                  pastroppetit=Standard_False;
-                }
-
-                if(pasuv[1]<pasInit[1])
-                {
-                  pasuv[1]+=(pasInit[1]-pasuv[1])*0.25;
-                  pastroppetit=Standard_False;
-                }
-
-                if(pasuv[2]<pasInit[2])
+                for(Standard_Integer i = 0; i < 4; i++)
                 {
-                  pasuv[2]+=(pasInit[2]-pasuv[2])*0.25;
-                  pastroppetit=Standard_False; 
-                }
-
-                if(pasuv[3]<pasInit[3])
-                {
-                  pasuv[3]+=(pasInit[3]-pasuv[3])*0.25;
-                  pastroppetit=Standard_False;
+                  if(pasuv[i]<pasInit[i])
+                  {
+                    pasuv[i]+=(pasInit[i]-pasuv[i])*0.25;
+                    pastroppetit=Standard_False;
+                  }
                 }
 
                 if(pastroppetit)
@@ -1087,6 +1198,32 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
 
             break;
           }
+        case IntWalk_StepTooSmall:
+          {
+            Standard_Boolean hasStepBeenIncreased = Standard_False;
+
+            for(Standard_Integer i = 0; i < 4; i++)
+            {
+              const Standard_Real aNewStep = Min(1.5*pasuv[i], pasInit[i]);
+              if(aNewStep > pasuv[i])
+              {
+                pasuv[i] = aNewStep;
+                hasStepBeenIncreased = Standard_True;
+              }
+            }
+
+            if(hasStepBeenIncreased)
+            {
+              Param(1)=SvParam[0];
+              Param(2)=SvParam[1];
+              Param(3)=SvParam[2];
+              Param(4)=SvParam[3];
+
+              LevelOfIterWithoutAppend = 0;
+
+              break;
+            }
+          }
         case IntWalk_OK:
         case IntWalk_ArretSurPoint://006
           {
@@ -1175,6 +1312,7 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
 
                     if(RejectIndex >= RejectIndexMAX)
                     {
+                      Arrive = Standard_True;
                       break;
                     }
 
@@ -1262,6 +1400,7 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
 
                       if(RejectIndex >= RejectIndexMAX)
                       {
+                        Arrive = Standard_True;
                         break;
                       }
 
@@ -1337,8 +1476,11 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
                           gp_Vec2d V2onS1(gp_Pnt2d(pU1, pV1), gp_Pnt2d(u1, v1));
                           gp_Vec2d V1onS2(gp_Pnt2d(ppU2, ppV2), gp_Pnt2d(pU2, pV2));
                           gp_Vec2d V2onS2(gp_Pnt2d(pU2, pV2), gp_Pnt2d(u2, v2));
-                          if (V1onS1 * V2onS1 < 0. ||
-                              V1onS2 * V2onS2 < 0.)
+
+                          const Standard_Real aDot1 = V1onS1 * V2onS1;
+                          const Standard_Real aDot2 = V1onS2 * V2onS2;
+
+                          if ((aDot1 < 0.0) || (aDot2 < 0.0))
                           {
                             Arrive = Standard_True;
                             break;
@@ -1433,7 +1575,58 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
                             Arrive = Standard_True;
                             break;
                           }
-                        }
+
+                          {
+                            //Check singularity.
+                            //I.e. check if we are walking along direction, which does not
+                            //result in comming to any point (i.e. derivative 
+                            //3D-intersection curve along this direction is equal to 0).
+                            //A sphere with direction {dU=1, dV=0} from point
+                            //(U=0, V=M_PI/2) can be considered as example for
+                            //this case (we cannot find another 3D-point if we go thus).
+
+                            //Direction chosen along 1st and 2nd surface correspondingly
+                            const gp_Vec2d  aDirS1(prevPntOnS1, curPntOnS1),
+                                            aDirS2(prevPntOnS2, curPntOnS2);
+
+                            gp_Pnt aPtemp;
+                            gp_Vec aDuS1, aDvS1, aDuS2, aDvS2;
+
+                            myIntersectionOn2S.Function().AuxillarSurface1()->
+                                  D1(curPntOnS1.X(), curPntOnS1.Y(), aPtemp, aDuS1, aDvS1);
+                            myIntersectionOn2S.Function().AuxillarSurface2()->
+                                  D1(curPntOnS2.X(), curPntOnS2.Y(), aPtemp, aDuS2, aDvS2);
+
+                            //Derivative WLine along (it is vector-function indeed)
+                            //directions chosen
+                            //(https://en.wikipedia.org/wiki/Directional_derivative#Variation_using_only_direction_of_vector).
+                            //F1 - on the 1st surface, F2 - on the 2nd surface.
+                            //x, y, z - coordinates of derivative vector.
+                            const Standard_Real aF1x =  aDuS1.X()*aDirS1.X() + 
+                                                        aDvS1.X()*aDirS1.Y();
+                            const Standard_Real aF1y =  aDuS1.Y()*aDirS1.X() +
+                                                        aDvS1.Y()*aDirS1.Y();
+                            const Standard_Real aF1z =  aDuS1.Z()*aDirS1.X() +
+                                                        aDvS1.Z()*aDirS1.Y();
+                            const Standard_Real aF2x =  aDuS2.X()*aDirS2.X() +
+                                                        aDvS2.X()*aDirS2.Y();
+                            const Standard_Real aF2y =  aDuS2.Y()*aDirS2.X() +
+                                                        aDvS2.Y()*aDirS2.Y();
+                            const Standard_Real aF2z =  aDuS2.Z()*aDirS2.X() +
+                                                        aDvS2.Z()*aDirS2.Y();
+
+                            const Standard_Real aF1 = aF1x*aF1x + aF1y*aF1y + aF1z*aF1z;
+                            const Standard_Real aF2 = aF2x*aF2x + aF2y*aF2y + aF2z*aF2z;
+
+                            if((aF1 < gp::Resolution()) && (aF2 < gp::Resolution()))
+                            {
+                              //All derivative are equal to 0. Therefore, there is
+                              //no point in going along direction chosen.
+                              Arrive = Standard_True;
+                              break;
+                            }
+                          }
+                        }//if (previoustg) cond.
 
                         ////////////////////////////////////////
                         AddAPoint(line,previousPoint);
@@ -1441,6 +1634,7 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
 
                         if(RejectIndex >= RejectIndexMAX)
                         {
+                          Arrive = Standard_True;
                           break;
                         }
 
@@ -1581,7 +1775,7 @@ Standard_Boolean IntWalk_PWalking::ExtendLineInCommonZone(const IntImp_ConstIsop
         return bOutOfTangentZone;
       }
 
-      Status = TestDeflection();
+      Status = TestDeflection(ChoixIso);
 
       if(Status == IntWalk_OK) {
 
@@ -1917,37 +2111,20 @@ DistanceMinimizeByGradient( const Handle(Adaptor3d_HSurface)& theASurf1,
   const Standard_Real aTol = 1.0e-14;
   Handle(Geom_Surface) aS1, aS2;
 
-  switch(theASurf1->GetType())
-  {
-  case GeomAbs_BezierSurface:
-    aS1 = theASurf1->Surface().Bezier();
-    break;
-  case GeomAbs_BSplineSurface:
-    aS1 = theASurf1->Surface().BSpline();
-    break;
-  default:
-    return Standard_True;
-  }
-
-  switch(theASurf2->GetType())
-  {
-  case GeomAbs_BezierSurface:
-    aS2 = theASurf2->Surface().Bezier();
-    break;
-  case GeomAbs_BSplineSurface:
-    aS2 = theASurf2->Surface().BSpline();
-    break;
-  default:
-    return Standard_True;
-  }
+  if (theASurf1->GetType() != GeomAbs_BezierSurface &&
+      theASurf1->GetType() != GeomAbs_BSplineSurface)
+      return Standard_True;
+  if (theASurf2->GetType() != GeomAbs_BezierSurface &&
+      theASurf2->GetType() != GeomAbs_BSplineSurface)
+      return Standard_True;
 
   Standard_Boolean aStatus = Standard_False;
 
   gp_Pnt aP1, aP2;
   gp_Vec aD1u, aD1v, aD2U, aD2V;
 
-  aS1->D1(theU1, theV1, aP1, aD1u, aD1v);
-  aS2->D1(theU2, theV2, aP2, aD2U, aD2V);
+  theASurf1->D1(theU1, theV1, aP1, aD1u, aD1v);
+  theASurf2->D1(theU2, theV2, aP2, aD2U, aD2V);
 
   Standard_Real aSQDistPrev = aP1.SquareDistance(aP2);
 
@@ -1984,8 +2161,8 @@ DistanceMinimizeByGradient( const Handle(Adaptor3d_HSurface)& theASurf1,
 
     gp_Pnt aPt1, aPt2;
 
-    aS1->D1(aPARu, aPARv, aPt1, aD1u, aD1v);
-    aS2->D1(aParU, aParV, aPt2, aD2U, aD2V);
+    theASurf1->D1(aPARu, aPARv, aPt1, aD1u, aD1v);
+    theASurf2->D1(aParU, aParV, aPt2, aD2U, aD2V);
 
     Standard_Real aSQDist = aPt1.SquareDistance(aPt2);
 
@@ -2009,14 +2186,14 @@ DistanceMinimizeByGradient( const Handle(Adaptor3d_HSurface)& theASurf1,
       }
       else
       {
-        aS1->D1(theU1, theV1, aPt1, aD1u, aD1v);
-        aS2->D1(theU2, theV2, aPt2, aD2U, aD2V);
-
-        gp_Vec aP12(aPt1, aPt2);
-        aGradFu = -aP12.Dot(aD1u);
-        aGradFv = -aP12.Dot(aD1v);
-        aGradFU = aP12.Dot(aD2U);
-        aGradFV = aP12.Dot(aD2V);
+        theASurf1->D1(theU1, theV1, aPt1, aD1u, aD1v);
+        theASurf2->D1(theU2, theV2, aPt2, aD2U, aD2V);
+
+        gp_Vec aPt12(aPt1, aPt2);
+        aGradFu = -aPt12.Dot(aD1u);
+        aGradFv = -aPt12.Dot(aD1v);
+        aGradFU = aPt12.Dot(aD2U);
+        aGradFV = aPt12.Dot(aD2V);
         aSTEPuv = theStep0U1V1;
         aStepUV = theStep0U2V2;
       }
@@ -2209,7 +2386,6 @@ PutToBoundary(const Handle(Adaptor3d_HSurface)& theASurf1,
   IsParallel(line, Standard_True, aTol, isU1parallel, isV1parallel);
   IsParallel(line, Standard_False, aTol, isU2parallel, isV2parallel);
 
-  const Standard_Integer aNbPnts = line->NbPoints();
   Standard_Real u1, v1, u2, v2;
   line->Value(1).Parameters(u1, v1, u2, v2);
   Standard_Real aDelta = 0.0;
@@ -2297,6 +2473,7 @@ PutToBoundary(const Handle(Adaptor3d_HSurface)& theASurf1,
       v1, u2, v2, Standard_True);
   }
 
+  const Standard_Integer aNbPnts = line->NbPoints();
   isNeedAdding = Standard_False;
   line->Value(aNbPnts).Parameters(u1, v1, u2, v2);
 
@@ -2622,7 +2799,7 @@ namespace {
   static const Standard_Real d = 7.0;
 }
 
-IntWalk_StatusDeflection  IntWalk_PWalking::TestDeflection()
+IntWalk_StatusDeflection  IntWalk_PWalking::TestDeflection(const IntImp_ConstIsoparametric choixIso)
 
 // test if vector is observed by calculating an increase of vector 
 //     or the previous point and its tangent, the new calculated point and its  
@@ -2643,8 +2820,11 @@ IntWalk_StatusDeflection  IntWalk_PWalking::TestDeflection()
   }
 
   IntWalk_StatusDeflection Status = IntWalk_OK;
-  Standard_Real FlecheCourante ,Ratio;
+  Standard_Real FlecheCourante , Ratio = 1.0;
 
+  // Caro1 and Caro2
+  const Handle(Adaptor3d_HSurface)& Caro1 = myIntersectionOn2S.Function().AuxillarSurface1();
+  const Handle(Adaptor3d_HSurface)& Caro2 = myIntersectionOn2S.Function().AuxillarSurface2();
 
   const IntSurf_PntOn2S& CurrentPoint = myIntersectionOn2S.Point(); 
   //==================================================================================
@@ -2656,10 +2836,12 @@ IntWalk_StatusDeflection  IntWalk_PWalking::TestDeflection()
 
   const gp_Dir& TgCourante = myIntersectionOn2S.Direction();
 
+  const Standard_Real aCosBetweenTangent = TgCourante.Dot(previousd);
+
   //==================================================================================
   //=========   R i s k   o f    i n f l e x i o n   p o i n t  ============
   //==================================================================================  
-  if (TgCourante.Dot(previousd)<0) {
+  if (aCosBetweenTangent < 0) {
     //------------------------------------------------------------
     //-- Risk of inflexion point : Divide the step by 2
     //-- Initialize STATIC_PRECEDENT_INFLEXION so that 
@@ -2677,7 +2859,6 @@ IntWalk_StatusDeflection  IntWalk_PWalking::TestDeflection()
     else 
       return IntWalk_PasTropGrand;
   }
-
   else {
     if(STATIC_PRECEDENT_INFLEXION  > 0) { 
       STATIC_PRECEDENT_INFLEXION -- ;
@@ -2689,15 +2870,70 @@ IntWalk_StatusDeflection  IntWalk_PWalking::TestDeflection()
   //=========  D e t e c t    c o n f u s e d    P o in t s       ===========
   //==================================================================================
 
-  Standard_Real Dist = previousPoint.Value().
+  const Standard_Real aSqDist = previousPoint.Value().
     SquareDistance(CurrentPoint.Value());
 
 
-  if (Dist < tolconf*tolconf ) { 
-    pasuv[0] = Max(5.*ResoU1,Min(1.5*pasuv[0],pasInit[0]));
-    pasuv[1] = Max(5.*ResoV1,Min(1.5*pasuv[1],pasInit[1]));
-    pasuv[2] = Max(5.*ResoU2,Min(1.5*pasuv[2],pasInit[2]));
-    pasuv[3] = Max(5.*ResoV2,Min(1.5*pasuv[3],pasInit[3]));
+  if (aSqDist < tolconf*tolconf) { 
+    pasInit[0] = Max(pasInit[0], 5.0*ResoU1);
+    pasInit[1] = Max(pasInit[1], 5.0*ResoV1);
+    pasInit[2] = Max(pasInit[2], 5.0*ResoU2);
+    pasInit[3] = Max(pasInit[3], 5.0*ResoV2);
+
+    for(Standard_Integer i = 0; i < 4; i++)
+    {
+      pasuv[i] = Max(pasuv[i], Min(1.5*pasuv[i], pasInit[i]));
+    }
+    //Compute local resolution: for OCC26717
+    if (Abs(pasuv[choixIso] - pasInit[choixIso]) <= Precision::Confusion())
+    {
+      Standard_Real CurU, CurV;
+      if (choixIso == IntImp_UIsoparametricOnCaro1 ||
+          choixIso == IntImp_VIsoparametricOnCaro1)
+        previousPoint.ParametersOnS1(CurU, CurV);
+      else
+        previousPoint.ParametersOnS2(CurU, CurV);
+      gp_Pnt CurPnt = (choixIso == IntImp_UIsoparametricOnCaro1 ||
+                       choixIso == IntImp_VIsoparametricOnCaro1)?
+        Adaptor3d_HSurfaceTool::Value(Caro1, CurU, CurV) :
+        Adaptor3d_HSurfaceTool::Value(Caro2, CurU, CurV);
+      gp_Pnt OffsetPnt;
+      switch(choixIso)
+      {
+      case IntImp_UIsoparametricOnCaro1:
+        OffsetPnt =
+          Adaptor3d_HSurfaceTool::Value(Caro1,
+                                        CurU + sensCheminement*pasuv[0],
+                                        CurV);
+        break;
+      case IntImp_VIsoparametricOnCaro1:
+        OffsetPnt =
+          Adaptor3d_HSurfaceTool::Value(Caro1,
+                                        CurU,
+                                        CurV + sensCheminement*pasuv[1]);
+        break;
+      case IntImp_UIsoparametricOnCaro2:
+        OffsetPnt =
+          Adaptor3d_HSurfaceTool::Value(Caro2,
+                                        CurU + sensCheminement*pasuv[2],
+                                        CurV);
+        break;
+      case IntImp_VIsoparametricOnCaro2:
+        OffsetPnt =
+          Adaptor3d_HSurfaceTool::Value(Caro2,
+                                        CurU,
+                                        CurV + sensCheminement*pasuv[3]);
+        break;
+      default:break;
+      }
+      Standard_Real RefDist = CurPnt.Distance(OffsetPnt);
+      Standard_Real LocalResol = 0.;
+      if (RefDist > gp::Resolution())
+        LocalResol = pasuv[choixIso] * tolconf / RefDist;
+      if (pasuv[choixIso] <= LocalResol)
+        pasuv[choixIso] = pasInit[choixIso] = 2*LocalResol;
+    }
+    ////////////////////////////////////////
     Status = IntWalk_PointConfondu;
   }
 
@@ -2816,7 +3052,7 @@ IntWalk_StatusDeflection  IntWalk_PWalking::TestDeflection()
   //-- Estimate of the vector           --
   //---------------------------------------
   FlecheCourante =
-    Sqrt(Abs((previousd.XYZ()-TgCourante.XYZ()).SquareModulus()*Dist))/8.;
+    Sqrt(Abs((previousd.XYZ()-TgCourante.XYZ()).SquareModulus()*aSqDist))/8.;
 
   if ( FlecheCourante<= fleche*0.5) {     //-- Current step too small
     if(FlecheCourante>1e-16) { 
@@ -2883,10 +3119,114 @@ IntWalk_StatusDeflection  IntWalk_PWalking::TestDeflection()
       Ratio = 0.75 * (fleche / FlecheCourante);
     }
   }
-  pasuv[0] = Max(5.*ResoU1,Min(Min(Ratio*AbsDu1,pasuv[0]),pasInit[0]));
-  pasuv[1] = Max(5.*ResoV1,Min(Min(Ratio*AbsDv1,pasuv[1]),pasInit[1]));
-  pasuv[2] = Max(5.*ResoU2,Min(Min(Ratio*AbsDu2,pasuv[2]),pasInit[2]));
-  pasuv[3] = Max(5.*ResoV2,Min(Min(Ratio*AbsDv2,pasuv[3]),pasInit[3]));
+
+  if(Status != IntWalk_PointConfondu)
+  {
+    //Here, aCosBetweenTangent >= 0.0 definitely.
+
+    /*
+    Brief algorithm description.
+    We have two (not-coincindent) intersection point (P1 and P2). In every point,
+    vector of tangent (to the intersection curve) is known (vectors T1 and T2).
+    Basing on these data, we create osculating circle.
+
+                                    * - arc of osculating circle
+                                *      * 
+                           P1 x----------x P2
+                             /            \
+                            /              \
+                          Vec(T1)         Vec(T2)
+
+    Let me draw your attention to the following facts:
+      1. Vectors T1 and T2 direct FROM (not TO) points P1 and P2. Therefore,
+        one of previously computed vector should be reversed.
+
+    In this case, the absolute (!) value of the deflection between the arc of
+    the osculating circle and the P1P2 segment can be computed as follows:
+          e = d*(1-sin(B/2))/(2*cos(B/2)),      (1)
+    where d is the length of P1P2 segment, B is the angle between vectors T1 and T2.
+    At that,
+              pi/2 <= B <= pi,
+              cos(B/2) >= 0,
+              sin(B/2) > 0,
+              sin(B) > 0,
+              cos(B) < 0.
+
+    Later, the normal state of algorithm work is (as we apply) 
+              tolconf/2 <= e <= tolconf.
+    In this case, we keep previous step.
+
+    If e < tolconf/2 then the local curvature of the intersection curve is small.
+    As result, the step should be increased.
+
+    If e > tolconf then the step is too big. Therefore, we should decrease one.
+
+    Condition (1) is equivalent to
+              sin(B/2) = 1 - 2/(1+(d/(2*e))^2) = Fs(e),
+              cos(B) = 1 - 2*Fs(e)^2 = Fd(e),
+    where Fs(e)and Fd(e) are some function with parameter "deflection".
+
+    Let mean that Fs(e) is decreasing function. Fd(e) is increasing function,
+    in the range, where Fs(e) > 0.0 (i.e. when e < d/2).
+
+    Now, let substitute required deflection (tolconf or tolconf/2) to z. Then
+    it is necessary to check if e < z or if e > z.
+
+    In this case, it is enough to comapare Fs(e) and Fs(z).
+    At that Fs(e) > 0 because sin(B/2) > 0 always.
+
+    Therefore, if Fs(z) < 0.0 then Fs(e) > Fs(z) ==> e < z definitely.
+    If Fs(z) > 0.0 then we can compare Fs(z)^2 and Fs(e)^2 or, in substance,
+    values Fd(e) and Fd(z). If Fd(e) > Fd(z) then e > z and vice versa.
+    */
+
+    //Fd(e) is already known (Fd(e) == -aCosBetweenTangent)
+
+    const Standard_Real anInvSqAbsArcDeflMax = 0.25*aSqDist/(tolconf*tolconf);
+    const Standard_Real aSinB2Max = 1.0 - 2.0/(1.0 + anInvSqAbsArcDeflMax);
+
+    if(aSinB2Max >= 0.0 && (aCosBetweenTangent <= 2.0 * aSinB2Max * aSinB2Max - 1.0))
+    {//Real deflection is greater or equal than tolconf
+      Status = IntWalk_PasTropGrand;
+    }
+    else
+    {//Real deflection is less than tolconf
+      const Standard_Real anInvSqAbsArcDeflMin = 4.0*anInvSqAbsArcDeflMax;
+      const Standard_Real aSinB2Min = 1.0 - 2.0/(1.0 + anInvSqAbsArcDeflMin);
+
+      if((aSinB2Min < 0.0) || (aCosBetweenTangent >= 2.0 * aSinB2Min * aSinB2Min - 1.0))
+      {//Real deflection is less than tolconf/2.0
+        Status = IntWalk_StepTooSmall;
+      }
+    }
+
+    if(Status == IntWalk_PasTropGrand)
+    {
+      pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
+      return Status;
+    }
+
+    if(Status == IntWalk_StepTooSmall)
+    {
+      pasuv[0] = Max(pasuv[0], AbsDu1);
+      pasuv[1] = Max(pasuv[1], AbsDv1);
+      pasuv[2] = Max(pasuv[2], AbsDu2);
+      pasuv[3] = Max(pasuv[3], AbsDv2);
+
+      pasInit[0] = Max(pasInit[0], AbsDu1);
+      pasInit[1] = Max(pasInit[1], AbsDv1);
+      pasInit[2] = Max(pasInit[2], AbsDu2);
+      pasInit[3] = Max(pasInit[3], AbsDv2);
+
+      return Status;
+    }
+  }
+
+  pasuv[0] = Max(myStepMin[0],Min(Min(Ratio*AbsDu1,pasuv[0]),pasInit[0]));
+  pasuv[1] = Max(myStepMin[1],Min(Min(Ratio*AbsDv1,pasuv[1]),pasInit[1]));
+  pasuv[2] = Max(myStepMin[2],Min(Min(Ratio*AbsDu2,pasuv[2]),pasInit[2]));
+  pasuv[3] = Max(myStepMin[3],Min(Min(Ratio*AbsDv2,pasuv[3]),pasInit[3]));
+
   if(Status == IntWalk_OK) STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=0;
   return Status;
 }