0027190: IntPatch_ImpPrmIntersection algorithm does not split intersection curve by the seam-edge of the quadric

1. Processing when IntPatch_WLine/IntPatch_RLine goes through the seam edge has been improved in DecomposeResult(...) function (see IntPatch_ImpPrmIntersection.cxx).
2. Incorrect initialization of last point of IntPatch_WLine/IntPatch_RLine has been eliminated. Earlier it was the reason of exception.

Creation of test case for this issue.

Adjusting some test cases according to their new behavior. Namely:

1) tests\bugs\modalg_4\bug825 (bug825_2)
Details are described in issue #25915. In short, new intersection algorithm works better than old (WLine without "jumping"). However, Boolean operation loses degenerated edges of the sphere. Consequently, we get the result with Not-closed face.

2) tests\bugs\modalg_6\bug26684_2
TolReached of intersection curve has become smaller. Consequently, intersection algorithm works better than earlier.

src/IntPatch/IntPatch_ImpPrmIntersection.cxx

@@ -69,6 +69,7 @@
 #include <IntPatch_PointLine.hxx>
 
 #include <Extrema_GenLocateExtPS.hxx>
+#include <math_FunctionSetRoot.hxx>
 
 static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLine,
                                         const Standard_Boolean       IsReversed,
@@ -77,6 +78,7 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
                                         const Handle(Adaptor3d_HSurface)&  theQSurf,
                                         const Handle(Adaptor3d_HSurface)&  theOtherSurf,
                                         const Standard_Real                theArcTol,
+                                        const Standard_Real                theTolTang,
                                         IntPatch_SequenceOfLine&           theLines);
 static 
   void ComputeTangency (const IntPatch_TheSOnBounds& solrst,
@@ -106,6 +108,190 @@ static
                               const Standard_Real theToler2D,
                               const Standard_Real thePeriod);
 
+enum PrePoint_Type
+{
+  PrePoint_NONE,
+  PrePoint_SEAMU,
+  PrePoint_SEAMV,
+  PrePoint_SEAMUV,
+  PrePoint_POLESEAMU,
+  PrePoint_POLE
+};
+
+static PrePoint_Type IsSeamOrPole(const Handle(Adaptor3d_HSurface)& theQSurf,
+                                  const Handle(IntSurf_LineOn2S)& theLine,
+                                  const Standard_Boolean IsReversed,
+                                  const Standard_Integer theRefIndex,
+                                  const Standard_Real theDeltaMax)
+{
+  if((theRefIndex < 1) || (theRefIndex >= theLine->NbPoints()))
+    return PrePoint_NONE;
+
+  //Parameters on Quadric and on parametric for reference point
+  Standard_Real aUQRef, aVQRef, aUPRef, aVPRef;
+  Standard_Real aUQNext, aVQNext, aUPNext, aVPNext;
+
+  if(IsReversed)
+  {
+    theLine->Value(theRefIndex).Parameters  (aUPRef, aVPRef, aUQRef, aVQRef);
+    theLine->Value(theRefIndex+1).Parameters(aUPNext, aVPNext, aUQNext, aVQNext);
+  }
+  else
+  {
+    theLine->Value(theRefIndex).Parameters  (aUQRef, aVQRef, aUPRef, aVPRef);
+    theLine->Value(theRefIndex+1).Parameters(aUQNext, aVQNext, aUPNext, aVPNext);
+  }
+
+  const GeomAbs_SurfaceType aType = theQSurf->GetType();
+
+  const Standard_Real aDeltaU = Abs(aUQRef - aUQNext);
+
+  if((aType != GeomAbs_Torus) && (aDeltaU < theDeltaMax))
+    return PrePoint_NONE;
+
+  switch(aType)
+  {
+  case GeomAbs_Cylinder:
+    return PrePoint_SEAMU;
+
+  case GeomAbs_Torus:
+    {
+      const Standard_Real aDeltaV = Abs(aVQRef - aVQNext);
+
+      if((aDeltaU >= theDeltaMax) && (aDeltaV >= theDeltaMax))
+        return PrePoint_SEAMUV;
+
+      if(aDeltaU >= theDeltaMax)
+        return PrePoint_SEAMU;
+
+      if(aDeltaV >= theDeltaMax)
+        return PrePoint_SEAMV;
+    }
+
+    break;
+  case GeomAbs_Sphere:
+  case GeomAbs_Cone:
+    return PrePoint_POLESEAMU;
+  default:
+    break;
+  }
+
+  return PrePoint_NONE;
+}
+
+// The function for searching intersection point, which 
+// lies in the seam-edge of the quadric definetely.
+class FuncPreciseSeam: public math_FunctionSetWithDerivatives
+{
+public:
+  FuncPreciseSeam(const Handle(Adaptor3d_HSurface)& theQSurf, const Handle(Adaptor3d_HSurface)& thePSurf, const Standard_Boolean isTheUSeam): myQSurf(theQSurf), myPSurf(thePSurf), myIsUSeam(isTheUSeam) {};
+  
+  Standard_EXPORT virtual Standard_Integer NbVariables() const
+  {
+    return 3;
+  };
+
+  Standard_EXPORT virtual Standard_Integer NbEquations() const
+  {
+    return 3;
+  }
+
+  Standard_EXPORT virtual Standard_Boolean Value (const math_Vector& theX, math_Vector& theF)
+  {
+    try
+    {
+      const Standard_Integer anIndX = theX.Lower(), anIndF = theF.Lower();
+      const gp_Pnt aP1(myPSurf->Value(theX(anIndX), theX(anIndX+1)));
+      const gp_Pnt aP2(myIsUSeam? myQSurf->Value(0.0, theX(anIndX+2)) : myQSurf->Value(theX(anIndX+2), 0.0));
+
+      (aP1.XYZ()-aP2.XYZ()).Coord(theF(anIndF), theF(anIndF+1), theF(anIndF+2));
+    }
+    catch(Standard_Failure)
+    {
+      return Standard_False;
+    }
+
+    return Standard_True;
+  };
+
+  Standard_EXPORT virtual Standard_Boolean Derivatives (const math_Vector& theX, math_Matrix& theD)
+  {
+    try
+    {
+      const Standard_Integer anIndX = theX.Lower(), anIndRD = theD.LowerRow(), anIndCD = theD.LowerCol();
+      gp_Pnt aPt;
+      gp_Vec aD1u, aD1v, aD2u, aD2v;
+      myPSurf->D1(theX(anIndX), theX(anIndX+1), aPt, aD1u, aD1v);
+      if(myIsUSeam)
+        myQSurf->D1(0.0, theX(anIndX+2), aPt, aD2u, aD2v);
+      else
+        myQSurf->D1(theX(anIndX+2), 0.0, aPt, aD2u, aD2v);
+
+      // d/dX1
+      aD1u.Coord(theD(anIndRD, anIndCD), theD(anIndRD+1, anIndCD), theD(anIndRD+2, anIndCD));
+
+      // d/dX1
+      aD1v.Coord(theD(anIndRD, anIndCD+1), theD(anIndRD+1, anIndCD+1), theD(anIndRD+2, anIndCD+1));
+
+      // d/dX3
+      if(myIsUSeam)
+        aD2v.Reversed().Coord(theD(anIndRD, anIndCD+2), theD(anIndRD+1, anIndCD+2), theD(anIndRD+2, anIndCD+2));
+      else
+        aD2u.Reversed().Coord(theD(anIndRD, anIndCD+2), theD(anIndRD+1, anIndCD+2), theD(anIndRD+2, anIndCD+2));
+    }
+    catch(Standard_Failure)
+    {
+      return Standard_False;
+    }
+
+    return Standard_True;
+  };
+
+  Standard_EXPORT virtual Standard_Boolean Values (const math_Vector& theX, math_Vector& theF, math_Matrix& theD)
+  {
+    try
+    {
+      const Standard_Integer anIndX = theX.Lower(), anIndF = theF.Lower(), anIndRD = theD.LowerRow(), anIndCD = theD.LowerCol();
+      gp_Pnt aP1, aP2;
+      gp_Vec aD1u, aD1v, aD2u, aD2v;
+      myPSurf->D1(theX(anIndX), theX(anIndX+1), aP1, aD1u, aD1v);
+      if(myIsUSeam)
+        myQSurf->D1(0.0, theX(anIndX+2), aP2, aD2u, aD2v);
+      else
+        myQSurf->D1(theX(anIndX+2), 0.0, aP2, aD2u, aD2v);
+
+      //Value
+      (aP1.XYZ()-aP2.XYZ()).Coord(theF(anIndF), theF(anIndF+1), theF(anIndF+2));
+
+      // d/dX1
+      aD1u.Coord(theD(anIndRD, anIndCD), theD(anIndRD+1, anIndCD), theD(anIndRD+2, anIndCD));
+
+      // d/dX1
+      aD1v.Coord(theD(anIndRD, anIndCD+1), theD(anIndRD+1, anIndCD+1), theD(anIndRD+2, anIndCD+1));
+
+      // d/dX3
+      if(myIsUSeam)
+        aD2v.Reversed().Coord(theD(anIndRD, anIndCD+2), theD(anIndRD+1, anIndCD+2), theD(anIndRD+2, anIndCD+2));
+      else
+        aD2u.Reversed().Coord(theD(anIndRD, anIndCD+2), theD(anIndRD+1, anIndCD+2), theD(anIndRD+2, anIndCD+2));
+    }
+    catch(Standard_Failure)
+    {
+      return Standard_False;
+    }
+
+    return Standard_True;
+  }
+
+protected:
+  FuncPreciseSeam operator=(FuncPreciseSeam&);
+
+private:
+  const Handle(Adaptor3d_HSurface)& myQSurf;
+  const Handle(Adaptor3d_HSurface)& myPSurf;
+  const Standard_Boolean myIsUSeam;
+};
+
 //=======================================================================
 //function : IntPatch_ImpPrmIntersection
 //purpose  : 
@@ -1501,7 +1687,9 @@ void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur
     return;
 
   Standard_Boolean isDecomposeRequired =  (Quad.TypeQuadric() == GeomAbs_Cone) || 
-                                          (Quad.TypeQuadric() == GeomAbs_Sphere);
+                                          (Quad.TypeQuadric() == GeomAbs_Sphere) ||
+                                          (Quad.TypeQuadric() == GeomAbs_Cylinder) ||
+                                          (Quad.TypeQuadric() == GeomAbs_Torus);
 
   if(!isDecomposeRequired)
     return;
@@ -1516,7 +1704,7 @@ void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur
   {
     if(DecomposeResult( Handle(IntPatch_PointLine)::DownCast(slin(i)),
                                         reversed, Quad, PDomain, aQSurf,
-                                        anOtherSurf, TolArc, dslin))
+                                        anOtherSurf, TolArc, aTol3d, dslin))
     {
       isDecompose = Standard_True;
     }
@@ -2318,7 +2506,7 @@ static Handle(IntPatch_WLine) MakeSplitWLine (Handle(IntPatch_WLine)&        WLi
   TPntL.SetParameters(uu1,vv1,uu2,vv2);
   TPntL.SetParameter((Standard_Real)sline->NbPoints());
   wline->AddVertex(TPntL);
-  wline->SetLastPoint(sline->NbPoints());
+  wline->SetLastPoint(wline->NbVertex());
 
   return wline;
 }
@@ -2369,6 +2557,7 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
                                         const Handle(Adaptor3d_HSurface)&  theQSurf, //quadric
                                         const Handle(Adaptor3d_HSurface)&  thePSurf, //parametric
                                         const Standard_Real                theArcTol,
+                                        const Standard_Real                theTolTang,
                                         IntPatch_SequenceOfLine&           theLines)
 {
   if(theLine->ArcType() == IntPatch_Restriction)
@@ -2430,12 +2619,7 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
   // build WLine parts (if any)
   Standard_Boolean flNextLine = Standard_True;
   Standard_Boolean hasBeenDecomposed = Standard_False;
-  enum PrePoint_Type
-  {
-    PrePoint_NONE,
-    PrePoint_SEAM,
-    PrePoint_POLE
-  }PrePointExist = PrePoint_NONE;
+  PrePoint_Type aPrePointExist = PrePoint_NONE;
 
   IntSurf_PntOn2S PrePoint;
   while(flNextLine)
@@ -2443,22 +2627,18 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
     // reset variables
     flNextLine = Standard_False;
     Standard_Boolean isDecomposited = Standard_False;
-    Standard_Real U1 = 0., U2 = 0., V1 = 0., V2 = 0., AnU1 = 0.;
+    Standard_Real U1 = 0., U2 = 0., V1 = 0., V2 = 0.;
 
     Handle(IntSurf_LineOn2S) sline = new IntSurf_LineOn2S();
 
     //if((Lindex-Findex+1) <= 2 )
-    if((aLindex <= aFindex) && (PrePointExist != PrePoint_POLE))
+    if((aLindex <= aFindex) && !aPrePointExist)
     {
       //break of "while(flNextLine)" cycle
       break;
     }
 
-    if (PrePointExist == PrePoint_SEAM)
-    {
-      sline->Add(PrePoint);
-    }
-    else if(PrePointExist == PrePoint_POLE)
+    if(aPrePointExist)
     {
       //The last point of the line is the pole of the quadric.
       //Therefore, Walking-line has been broken in this point.
@@ -2483,12 +2663,19 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
       //Consequently, when the line goes throug the pole, @U_{q}@ can be
       //changed on @\pi /2 @ (but not less).
 
+      //Here, we forbid "jumping" between two neighbor Walking-point
+      //with step greater than pi/4
       const Standard_Real aPeriod = M_PI_2, aHalfPeriod = M_PI_4;
       const IntSurf_PntOn2S& aRefPt = aSSLine->Value(aFindex);
 
-      IntSurf_PntOn2S aFirstPoint = PrePoint;
+      const Standard_Real aURes = theQSurf->UResolution(theArcTol),
+                          aVRes = theQSurf->UResolution(theArcTol);
 
-      if(!aFirstPoint.IsSame(aRefPt, Precision::Confusion()))
+      const Standard_Real aTol2d = (aPrePointExist == PrePoint_POLE) ? 0.0 : 
+              (aPrePointExist == PrePoint_SEAMV)? aVRes :
+              (aPrePointExist == PrePoint_SEAMUV)? Max(aURes, aVRes) : aURes;
+
+      if(!PrePoint.IsSame(aRefPt, Precision::Confusion(), aTol2d))
       {
         Standard_Real aURef = 0.0, aVRef = 0.0;
         Standard_Real aUquad = 0.0, aVquad = 0.0;
@@ -2496,15 +2683,16 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
         //Take parameters on quadric
         if(IsReversed)
         {
-          aFirstPoint.ParametersOnS2(aUquad, aVquad);
+          PrePoint.ParametersOnS2(aUquad, aVquad);
           aRefPt.ParametersOnS2(aURef, aVRef);
         }
         else
         {
-          aFirstPoint.ParametersOnS1(aUquad, aVquad);
+          PrePoint.ParametersOnS1(aUquad, aVquad);
           aRefPt.ParametersOnS1(aURef, aVRef);
         }
 
+        if(theQSurf->IsUPeriodic())
         {
           Standard_Real aDeltaPar = aURef-aUquad;
           const Standard_Real anIncr = Sign(aPeriod, aDeltaPar);
@@ -2515,8 +2703,19 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
           }
         }
 
-        aFirstPoint.SetValue(!IsReversed, aUquad, aVquad);
-        sline->Add(aFirstPoint);
+        if(theQSurf->IsVPeriodic())
+        {
+          Standard_Real aDeltaPar = aVRef-aVquad;
+          const Standard_Real anIncr = Sign(aPeriod, aDeltaPar);
+          while((aDeltaPar > aHalfPeriod) || (aDeltaPar < -aHalfPeriod))
+          {
+            aVquad += anIncr;
+            aDeltaPar = aVRef-aVquad;
+          }
+        }
+
+        PrePoint.SetValue(!IsReversed, aUquad, aVquad);
+        sline->Add(PrePoint);
       }
       else
       {
@@ -2525,24 +2724,15 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
       }
     }
 
-    PrePointExist = PrePoint_NONE;
+    aPrePointExist = PrePoint_NONE;
 
     // analyze other points
     for(Standard_Integer k = aFindex; k <= aLindex; k++)
     {
       if( k == aFindex )
       {
-        if(IsReversed)
-        {
-          aSSLine->Value(k).ParametersOnS2(AnU1,V1);   // S2 - quadric, set U,V by Pnt3D
-        }
-        else
-        {
-          aSSLine->Value(k).ParametersOnS1(AnU1,V1);    // S1 - quadric, set U,V by Pnt3D
-        }
-
-        sline->Add(aSSLine->Value(k));
         PrePoint = aSSLine->Value(k);
+        sline->Add(PrePoint);
         continue;
       }
 
@@ -2555,77 +2745,172 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
         aSSLine->Value(k).ParametersOnS1(U1,V1);    // S1 - quadric, set U,V by Pnt3D
       }
 
-      if(Abs(U1-AnU1) > aDeltaUmax)
+      aPrePointExist = IsSeamOrPole(theQSurf, aSSLine, IsReversed, k-1, aDeltaUmax);
+
+      if(aPrePointExist != PrePoint_NONE)
       {
         aBindex = k;
         isDecomposited = Standard_True;
         ////
-        if (Abs(U1) <= Precision::PConfusion() ||
-            Abs(U1 - 2*M_PI) <= Precision::PConfusion())
+        const Standard_Real aPeriod = M_PI+M_PI, aHalfPeriod = M_PI;
+        const IntSurf_PntOn2S& aRefPt = aSSLine->Value(aBindex-1);
+
+        //Not quadric point
+        Standard_Real aU0 = 0.0, aV0 = 0.0;
+        //Quadric point
+        Standard_Real aUQuadRef = 0.0, aVQuadRef = 0.0;
+
+        if(IsReversed)
         {
-          IntSurf_PntOn2S NewPoint;
-          IntSurf_PntOn2S CurPoint = aSSLine->Value(k);
-          gp_Pnt thePnt = CurPoint.Value();
-          Standard_Real theU1, theV1, theU2, theV2;
-          theU1 = (Abs(U1) <= Precision::PConfusion())? 2*M_PI : 0.;
-          theV1 = V1;
-          NewPoint.SetValue(thePnt);
-          if (!IsReversed)
-          {
-            CurPoint.ParametersOnS2(theU2, theV2);
-            NewPoint.SetValue(theU1, theV1, theU2, theV2);
-          }
-          else
-          {
-            CurPoint.ParametersOnS1(theU2, theV2);
-            NewPoint.SetValue(theU2, theV2, theU1, theV1);
-          }
-          sline->Add(NewPoint);
-        }
-        else if (Abs(AnU1) <= Precision::PConfusion() ||
-                 Abs(AnU1 - 2*M_PI) <= Precision::PConfusion())
-        {
-          //Modify <PrePoint>
-          PrePointExist = PrePoint_SEAM;
-          Standard_Real theU1, theV1;
-          if (!IsReversed)
-          {
-            PrePoint.ParametersOnS1(theU1, theV1);
-            theU1 = (Abs(AnU1) <= Precision::PConfusion())? 2*M_PI : 0.;
-            PrePoint.SetValue(Standard_True, //on first
-                              theU1, theV1);
-          }
-          else
-          {
-            PrePoint.ParametersOnS2(theU1, theV1);
-            theU1 = (Abs(AnU1) <= Precision::PConfusion())? 2*M_PI : 0.;
-            PrePoint.SetValue(Standard_False, //on second
-                              theU1, theV1);
-          }
+          aRefPt.Parameters(aU0, aV0, aUQuadRef, aVQuadRef);
         }
         else
-        {//Check if WLine goes through pole
-          const Standard_Real aTol = Precision::Confusion();
-          const Standard_Real aPeriod = M_PI+M_PI, aHalfPeriod = M_PI;
-          const IntSurf_PntOn2S& aRefPt = aSSLine->Value(aBindex-1);
-          
-          //Not quadric point
-          Standard_Real aU0 = 0.0, aV0 = 0.0;
-          //Quadric point
-          Standard_Real aUQuadRef = 0.0, aVQuadRef = 0.0;
+        {
+          aRefPt.Parameters(aUQuadRef, aVQuadRef, aU0, aV0);
+        }
 
+        if(aPrePointExist == PrePoint_SEAMUV)
+        {
+          aPrePointExist = PrePoint_NONE;
+
+          gp_Pnt aPQuad;
+          Standard_Real aUquad = 0.0;
+          Standard_Real aVquad = 0.0; 
+
+          theQSurf->D0(aUquad, aVquad, aPQuad);
+
+          Extrema_GenLocateExtPS anExtr(aPQuad, thePSurf->Surface(), aU0, aV0,
+                                        Precision::PConfusion(),
+                                        Precision::PConfusion());
+
+          if(!anExtr.IsDone())
+          {
+            break;
+          }
+
+          if(anExtr.SquareDistance() < theTolTang*theTolTang)
+          {
+            anExtr.Point().Parameter(aU0, aV0);
+            gp_Pnt aP0(anExtr.Point().Value());
+
+            IntSurf_PntOn2S aNewPoint;
+            aNewPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), IsReversed, aU0, aV0);
+
+            if(!aNewPoint.IsSame(aRefPt, Precision::Confusion()))
+            {
+              //Adjust found U-paramter to previous point of the Walking-line
+              Standard_Real aDeltaPar = aUQuadRef-aUquad;
+              const Standard_Real anIncrU = Sign(aPeriod, aDeltaPar);
+              while((aDeltaPar > aHalfPeriod) || (aDeltaPar < -aHalfPeriod))
+              {
+                aUquad += anIncrU;
+                aDeltaPar = aUQuadRef-aUquad;
+              }
+
+              //Adjust found V-paramter to previous point of the Walking-line
+              aDeltaPar = aVQuadRef-aVquad;
+              const Standard_Real anIncrV = Sign(aPeriod, aDeltaPar);
+              while((aDeltaPar > aHalfPeriod) || (aDeltaPar < -aHalfPeriod))
+              {
+                aVquad += anIncrV;
+                aDeltaPar = aVQuadRef-aVquad;
+              }
+
+              aNewPoint.SetValue(!IsReversed, aUquad, aVquad);
+              
+              sline->Add(aNewPoint);
+              aPrePointExist = PrePoint_SEAMUV;
+              PrePoint = aNewPoint;
+            }
+          }
+        }
+        else if(aPrePointExist == PrePoint_SEAMV)
+        {//WLine goes through seam
+          aPrePointExist = PrePoint_NONE;
+
+          FuncPreciseSeam aF(theQSurf, thePSurf, Standard_False);
+          math_Vector aTol(1, 3), aStartPoint(1,3);
+          
+          //Parameters on parametric surface
+          Standard_Real aUp = 0.0, aVp = 0.0;
           if(IsReversed)
           {
-            aRefPt.Parameters(aU0, aV0, aUQuadRef, aVQuadRef);
+            aSSLine->Value(k).ParametersOnS1(aUp, aVp);
           }
           else
           {
-            aRefPt.Parameters(aUQuadRef, aVQuadRef, aU0, aV0);
+            aSSLine->Value(k).ParametersOnS2(aUp, aVp);
           }
 
-          //Transforms parametric surface in coordinate-system of the quadric
-          gp_Trsf aTr;
-          aTr.SetTransformation(theQuad.Sphere().Position());
+          aTol(1) = thePSurf->UResolution(theArcTol);
+          aTol(2) = thePSurf->VResolution(theArcTol);
+          aTol(3) = theQSurf->UResolution(theArcTol);
+          aStartPoint(1) = 0.5*(aU0 + aUp);
+          aStartPoint(2) = 0.5*(aV0 + aVp);
+          aStartPoint(3) = 0.5*(aUQuadRef + U1);
+
+          math_FunctionSetRoot aSRF(aF, aTol);
+          aSRF.Perform(aF, aStartPoint);
+
+          if(!aSRF.IsDone())
+          {
+            break;
+          }
+
+          // Now aStartPoint is useless. Therefore, we use it for keeping
+          // new point.
+          aSRF.Root(aStartPoint);
+          
+          //On parametric
+          aU0 = aStartPoint(1);
+          aV0 = aStartPoint(2);
+
+          //On quadric
+          Standard_Real aUquad = aStartPoint(3);
+          Standard_Real aVquad = 0.0; 
+          const gp_Pnt aPQuad(theQSurf->Value(aUquad, aVquad));
+          const gp_Pnt aP0(thePSurf->Value(aU0, aV0));
+
+          {
+            //Adjust found U-paramter to previous point of the Walking-line
+            Standard_Real aDeltaPar = aVQuadRef-aVquad;
+            const Standard_Real anIncr = Sign(aPeriod, aDeltaPar);
+            while((aDeltaPar > aHalfPeriod) || (aDeltaPar < -aHalfPeriod))
+            {
+              aVquad += anIncr;
+              aDeltaPar = aVQuadRef-aVquad;
+            }
+          }
+
+          IntSurf_PntOn2S aNewPoint;
+          if(IsReversed)
+            aNewPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aU0, aV0, aUquad, aVquad);
+          else
+            aNewPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aUquad, aVquad, aU0, aV0);
+
+          if(!aNewPoint.IsSame(aRefPt, Precision::Confusion(), Precision::PConfusion()))
+          {
+            aNewPoint.SetValue(!IsReversed, aUquad, aVquad);
+            sline->Add(aNewPoint);
+            aPrePointExist = PrePoint_SEAMV;
+            PrePoint = aNewPoint;
+          }
+          else
+          {
+            if(sline->NbPoints() == 1)
+            {
+              //FIRST point of the sline is the pole of the quadric.
+              //Therefore, there is no point in decomposition.
+
+              PrePoint = aRefPt;
+              aPrePointExist = PrePoint_SEAMV;
+            }
+          }
+        }
+        else if(aPrePointExist == PrePoint_POLESEAMU)
+        {//Check if WLine goes through pole
+          
+          aPrePointExist = PrePoint_NONE;
 
           //aPQuad is Pole
           gp_Pnt aPQuad;
@@ -2656,9 +2941,11 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
                                         Precision::PConfusion());
 
           if(!anExtr.IsDone())
+          {
             break;
+          }
 
-          if(anExtr.SquareDistance() < aTol*aTol)
+          if(anExtr.SquareDistance() < theTolTang*theTolTang)
           { //Pole is an intersection point
             //(lies in the quadric and the parametric surface)
 
@@ -2690,6 +2977,12 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
               gp_Vec aVecDu, aVecDv;
               thePSurf->D1(aU0, aV0, aPtemp, aVecDu, aVecDv);
 
+              //Transforms parametric surface in coordinate-system of the quadric
+              gp_Trsf aTr;
+              aTr.SetTransformation((theQuad.TypeQuadric() == GeomAbs_Sphere) ?
+                                      theQuad.Sphere().Position() :
+                                      theQuad.Cone().Position());
+
               //Derivatives of transformed thePSurf
               aVecDu.Transform(aTr);
               aVecDv.Transform(aTr);
@@ -2880,22 +3173,112 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
               aNewPoint.SetValue(!IsReversed, aUquad, aVquad);
               
               sline->Add(aNewPoint);
-              PrePointExist = PrePoint_POLE;
+              aPrePointExist = PrePoint_POLE;
               PrePoint = aNewPoint;
             } // if(!aNewPoint.IsSame(aRefPt, Precision::Confusion()))
             else
             {
+              aPrePointExist = PrePoint_NONE;
+
               if(sline->NbPoints() == 1)
               {
                 //FIRST point of the sline is the pole of the quadric.
                 //Therefore, there is no point in decomposition.
 
                 PrePoint = aRefPt;
-                AnU1=U1;
-                PrePointExist = PrePoint_POLE;
+                aPrePointExist = PrePoint_POLE;
               }
             }
           } //if(anExtr.SquareDistance() < aTol*aTol)
+          else
+          {//Pole is not an intersection point
+            aPrePointExist = PrePoint_SEAMU;
+          }
+        }
+
+        if(aPrePointExist == PrePoint_SEAMU)
+        {//WLine goes through seam
+
+          aPrePointExist = PrePoint_NONE;
+
+          FuncPreciseSeam aF(theQSurf, thePSurf, Standard_True);
+          math_Vector aTol(1, 3), aStartPoint(1,3);
+          
+          //Parameters on parametric surface
+          Standard_Real aUp = 0.0, aVp = 0.0;
+          if(IsReversed)
+          {
+            aSSLine->Value(k).ParametersOnS1(aUp, aVp);
+          }
+          else
+          {
+            aSSLine->Value(k).ParametersOnS2(aUp, aVp);
+          }
+
+          aTol(1) = thePSurf->UResolution(theArcTol);
+          aTol(2) = thePSurf->VResolution(theArcTol);
+          aTol(3) = theQSurf->VResolution(theArcTol);
+          aStartPoint(1) = 0.5*(aU0 + aUp);
+          aStartPoint(2) = 0.5*(aV0 + aVp);
+          aStartPoint(3) = 0.5*(aVQuadRef + V1);
+
+          math_FunctionSetRoot aSRF(aF, aTol);
+          aSRF.Perform(aF, aStartPoint);
+
+          if(!aSRF.IsDone())
+          {
+            break;
+          }
+
+          // Now aStartPoint is useless. Therefore, we use it for keeping
+          // new point.
+          aSRF.Root(aStartPoint);
+          
+          //On parametric
+          aU0 = aStartPoint(1);
+          aV0 = aStartPoint(2);
+
+          //On quadric
+          Standard_Real aUquad = 0.0;
+          Standard_Real aVquad = aStartPoint(3); 
+          const gp_Pnt aPQuad(theQSurf->Value(aUquad, aVquad));
+          const gp_Pnt aP0(thePSurf->Value(aU0, aV0));
+
+          {
+            //Adjust found U-paramter to previous point of the Walking-line
+            Standard_Real aDeltaPar = aUQuadRef-aUquad;
+            const Standard_Real anIncr = Sign(aPeriod, aDeltaPar);
+            while((aDeltaPar > aHalfPeriod) || (aDeltaPar < -aHalfPeriod))
+            {
+              aUquad += anIncr;
+              aDeltaPar = aUQuadRef-aUquad;
+            }
+          }
+
+          IntSurf_PntOn2S aNewPoint;
+          if(IsReversed)
+            aNewPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aU0, aV0, aUquad, aVquad);
+          else
+            aNewPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aUquad, aVquad, aU0, aV0);
+
+          if(!aNewPoint.IsSame(aRefPt, Precision::Confusion(), Precision::PConfusion()))
+          {
+            aNewPoint.SetValue(!IsReversed, aUquad, aVquad);
+            sline->Add(aNewPoint);
+            aPrePointExist = PrePoint_SEAMU;
+            PrePoint = aNewPoint;
+          }
+          else
+          {
+            if(sline->NbPoints() == 1)
+            {
+              //FIRST point of the sline is the pole of the quadric.
+              //Therefore, there is no point in decomposition.
+
+              PrePoint = aRefPt;
+              aPrePointExist = PrePoint_SEAMU;
+            }
+          }
         }
 
         ////
@@ -2904,7 +3287,6 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
 
       sline->Add(aSSLine->Value(k));
       PrePoint = aSSLine->Value(k);
-      AnU1=U1;
     } //for(Standard_Integer k = aFindex; k <= aLindex; k++)
 
     //Creation of new line as part of existing theLine.
@@ -2980,7 +3362,7 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
       aTPntL.SetParameters(U1,V1,U2,V2);
       aTPntL.SetParameter(sline->NbPoints());
       wline->AddVertex(aTPntL);
-      wline->SetLastPoint(sline->NbPoints());
+      wline->SetLastPoint(wline->NbVertex());
 
       IntPatch_SequenceOfLine segm;
       Standard_Boolean isSplited = SplitOnSegments(wline,Standard_False,
@@ -3075,12 +3457,15 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
         aRLine->AddVertex(aTPnt);
       }
 
-      aRLine->SetFirstPoint(1);
-      aRLine->SetLastPoint(sline->NbPoints());
+      if(aLPar - aFPar > Precision::PConfusion())
+      {
+        aRLine->SetFirstPoint(1);
+        aRLine->SetLastPoint(aRLine->NbVertex());
 
-      anArc->Trim(aFPar, aLPar, theArcTol);
+        anArc->Trim(aFPar, aLPar, theArcTol);
 
-      theLines.Append(aRLine);
+        theLines.Append(aRLine);
+      }
     }
 
     if(isDecomposited)

src/IntTools/IntTools_FaceFace.cxx

@@ -648,6 +648,28 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
     }
   }
 
+#ifdef INTTOOLS_FACEFACE_DEBUG
+    if(!myListOfPnts.IsEmpty()) {
+      char aBuff[10000];
+      const IntSurf_PntOn2S& aPt = myListOfPnts.First();
+      Standard_Real u1, v1, u2, v2;
+      aPt.Parameters(u1, v1, u2, v2);
+
+      Sprintf(aBuff,"bopcurves <face1 face2> -2d");
+      IntSurf_ListIteratorOfListOfPntOn2S IterLOP1(myListOfPnts);
+      for(;IterLOP1.More(); IterLOP1.Next())
+      {
+        const IntSurf_PntOn2S& aPt = IterLOP1.Value();
+        Standard_Real u1, v1, u2, v2;
+        aPt.Parameters(u1, v1, u2, v2);
+
+        Sprintf(aBuff, "%s -p %+10.20f %+10.20f %+10.20f %+10.20f", aBuff, u1, v1, u2, v2);
+      }
+
+      cout << aBuff << endl;
+    }
+#endif
+
   const Standard_Boolean isGeomInt = isTreatAnalityc(aF1, aF2);
   myIntersector.Perform(myHS1, dom1, myHS2, dom2, TolArc, TolTang, 
                                   myListOfPnts, RestrictLine, isGeomInt);
@@ -888,28 +910,6 @@ void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
     }
     L = aWLine;
 
-#ifdef INTTOOLS_FACEFACE_DEBUG
-    if(!myListOfPnts.IsEmpty()) {
-      char aBuff[10000];
-      const IntSurf_PntOn2S& aPt = myListOfPnts.First();
-      Standard_Real u1, v1, u2, v2;
-      aPt.Parameters(u1, v1, u2, v2);
-
-      Sprintf(aBuff,"bopcurves <face1 face2> -2d");
-      IntSurf_ListIteratorOfListOfPntOn2S IterLOP1(myListOfPnts);
-      for(;IterLOP1.More(); IterLOP1.Next())
-      {
-        const IntSurf_PntOn2S& aPt = IterLOP1.Value();
-        Standard_Real u1, v1, u2, v2;
-        aPt.Parameters(u1, v1, u2, v2);
-
-        Sprintf(aBuff, "%s -p %+10.20f %+10.20f %+10.20f %+10.20f", aBuff, u1, v1, u2, v2);
-      }
-
-      cout << aBuff << endl;
-    }
-#endif
-
     Standard_Integer nbp = aWLine->NbPnts();
     const IntSurf_PntOn2S& p1 = aWLine->Point(1);
     const IntSurf_PntOn2S& p2 = aWLine->Point(nbp);