0029807: [Regression to 7.0.0] Impossible to cut cone from prism
[occt.git] / src / IntPatch / IntPatch_SpecialPoints.cxx
index b180d9a..7cda071 100644 (file)
@@ -16,6 +16,7 @@
 #include <IntPatch_SpecialPoints.hxx>
 
 #include <Adaptor3d_HSurface.hxx>
+#include <ElCLib.hxx>
 #include <Extrema_ExtPS.hxx>
 #include <Extrema_GenLocateExtPS.hxx>
 #include <Geom_ConicalSurface.hxx>
@@ -35,10 +36,12 @@ class FuncPreciseSeam: public math_FunctionSetWithDerivatives
 public:
   FuncPreciseSeam(const Handle(Adaptor3d_HSurface)& theQSurf, // quadric
                   const Handle(Adaptor3d_HSurface)& thePSurf, // another surface
-                  const Standard_Boolean isTheUSeam): 
+                  const Standard_Boolean isTheUSeam,
+                  const Standard_Real theIsoParameter): 
         myQSurf(theQSurf),
         myPSurf(thePSurf),
-        mySeamCoordInd(isTheUSeam? 1 : 0) // Defines, U- or V-seam is used
+        mySeamCoordInd(isTheUSeam? 1 : 0), // Defines, U- or V-seam is used
+        myIsoParameter(theIsoParameter)
   {
   };
   
@@ -58,7 +61,7 @@ public:
     try
     {
       const Standard_Integer anIndX = theX.Lower(), anIndF = theF.Lower();
-      Standard_Real aUV[] = {0.0, 0.0};
+      Standard_Real aUV[] = {myIsoParameter, myIsoParameter};
       aUV[mySeamCoordInd] = theX(anIndX+2);
       const gp_Pnt aP1(myPSurf->Value(theX(anIndX), theX(anIndX+1)));
       const gp_Pnt aP2(myQSurf->Value(aUV[0], aUV[1]));
@@ -81,7 +84,7 @@ public:
       const Standard_Integer anIndX = theX.Lower(),
                              anIndRD = theD.LowerRow(),
                              anIndCD = theD.LowerCol();
-      Standard_Real aUV[] = {0.0, 0.0};
+      Standard_Real aUV[] = {myIsoParameter, myIsoParameter};
       aUV[mySeamCoordInd] = theX(anIndX+2);
 
       gp_Pnt aPt;
@@ -131,11 +134,31 @@ private:
   const Handle(Adaptor3d_HSurface)& myQSurf;
   const Handle(Adaptor3d_HSurface)& myPSurf;
 
-  // 1 for U-coordinate, 0 - for V one. 
+  //! 1 for U-coordinate, 0 - for V one. 
   const Standard_Integer mySeamCoordInd;
+
+  //! Constant parameter of iso-line
+  const Standard_Real myIsoParameter;
 };
 
 //=======================================================================
+//function : GetTangent
+//purpose  : Computes tangent having the given parameter.
+//            See calling method(s) for detailed information
+//=======================================================================
+static inline void GetTangent(const Standard_Real theConeSemiAngle,
+                              const Standard_Real theParameter,
+                              gp_XYZ& theResult)
+{
+  const Standard_Real aW2 = theParameter*theParameter;
+  const Standard_Real aCosUn = (1.0 - aW2) / (1.0 + aW2);
+  const Standard_Real aSinUn = 2.0*theParameter / (1.0 + aW2);
+
+  const Standard_Real aTanA = Tan(theConeSemiAngle);
+  theResult.SetCoord(aTanA*aCosUn, aTanA*aSinUn, 1.0);
+}
+
+//=======================================================================
 //function : IsPointOnSurface
 //purpose  : Checks if thePt is in theSurf (with given tolerance).
 //            Returns the foot of projection (theProjPt) and its parameters
@@ -281,6 +304,7 @@ Standard_Boolean IntPatch_SpecialPoints::
                                         const Handle(Adaptor3d_HSurface)& thePSurf,
                                         const IntSurf_PntOn2S& theRefPt,
                                         const Standard_Boolean theIsU,
+                                        const Standard_Real theIsoParameter,
                                         const math_Vector& theToler,
                                         const math_Vector& theInitPoint,
                                         const math_Vector& theInfBound,
@@ -292,7 +316,7 @@ Standard_Boolean IntPatch_SpecialPoints::
   IntSurf::SetPeriod(theIsReversed ? thePSurf : theQSurf,
                      theIsReversed ? theQSurf : thePSurf, anArrOfPeriod);
 
-  FuncPreciseSeam aF(theQSurf, thePSurf, theIsU);
+  FuncPreciseSeam aF(theQSurf, thePSurf, theIsU, theIsoParameter);
 
   math_FunctionSetRoot aSRF(aF, theToler);
   aSRF.Perform(aF, theInitPoint, theInfBound, theSupBound);
@@ -324,6 +348,432 @@ Standard_Boolean IntPatch_SpecialPoints::
 }
 
 //=======================================================================
+//function : ProcessSphere
+//purpose  :
+/*
+The intersection point (including the pole)
+must be satisfied to the following system:
+
+    \left\{\begin{matrix}
+    R*\cos (U_{q})*\cos (V_{q})=S_{x}(U_{s},V_{s})
+    R*\sin (U_{q})*\cos (V_{q})=S_{y}(U_{s},V_{s})
+    R*\sin (V_{q})=S_{z}(U_{s},V_{s})
+    \end{matrix}\right,
+where
+  R is the radius of the sphere;
+  @S_{x}@, @S_{y}@ and @S_{z}@ are X, Y and Z-coordinates of thePSurf;
+  @U_{s}@ and @V_{s}@ are parameters on the parametric surface;
+  @U_{q}@ and @V_{q}@ are equal to theUquad and theVquad correspondingly.
+
+Consequently (from first two equations),
+  \left\{\begin{matrix}
+  \cos (U_{q}) = \frac{S_{x}(U_{s},V_{s})}{R*\cos (V_{q})}
+  \sin (U_{q}) = \frac{S_{y}(U_{s},V_{s})}{R*\cos (V_{q})}
+  \end{matrix}\right.
+
+For pole,
+  V_{q}=\pm \pi /2 \Rightarrow \cos (V_{q}) = 0 (denominator is equal to 0).
+
+Therefore, computation U_{q} directly is impossibly.
+
+Let @V_{q}@ tends to @\pm \pi /2@.
+Then (indeterminate form is evaluated in accordance of L'Hospital rule),
+  \cos (U_{q}) = \lim_{V_{q} \to (\pi /2-0)}
+  \frac{S_{x}(U_{s},V_{s})}{R*\cos (V_{q})}=
+  -\lim_{V_{q} \to (\pi /2-0)}
+  \frac{\frac{\partial S_{x}}
+  {\partial U_{s}}*\frac{\mathrm{d} U_{s}}
+  {\mathrm{d} V_{q}}+\frac{\partial S_{x}}
+  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}
+  {\mathrm{d} V_{q}}}{R*\sin (V_{q})} =
+  -\frac{1}{R}*\frac{\mathrm{d} U_{s}}
+  {\mathrm{d} V_{q}}*(\frac{\partial S_{x}}
+  {\partial U_{s}}+\frac{\partial S_{x}}
+  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}
+  {\mathrm{d} U_{s}}) =
+  -\frac{1}{R}*\frac{\mathrm{d} V_{s}}
+  {\mathrm{d} V_{q}}*(\frac{\partial S_{x}}
+  {\partial U_{s}}*\frac{\mathrm{d} U_{s}}
+  {\mathrm{d} V_{s}}+\frac{\partial S_{x}}
+  {\partial V_{s}}).
+
+Analogicaly for @\sin (U_{q})@ (@S_{x}@ is substituted to @S_{y}@).
+
+Let mean, that
+  \cos (U_{q}) \left | _{V_{q} \to (-\pi /2+0)} = \cos (U_{q}) \left | _{V_{q} \to (\pi /2-0)}
+  \sin (U_{q}) \left | _{V_{q} \to (-\pi /2+0)} = \sin (U_{q}) \left | _{V_{q} \to (\pi /2-0)}
+
+From the 3rd equation of the system, we obtain
+  \frac{\mathrm{d} (R*\sin (V_{q}))}{\mathrm{d} V_{q}} =
+  \frac{\mathrm{d} S_{z}(U_{s},V_{s})}{\mathrm{d} V_{q}}
+or
+  R*\cos (V_{q}) = \frac{\partial S_{z}}{\partial U_{s}}*
+  \frac{\mathrm{d} U_{s}} {\mathrm{d} V_{q}}+\frac{\partial S_{z}}
+  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}{\mathrm{d} V_{q}}.
+
+If @V_{q}=\pm \pi /2@, then
+  \frac{\partial S_{z}}{\partial U_{s}}*
+  \frac{\mathrm{d} U_{s}} {\mathrm{d} V_{q}}+\frac{\partial S_{z}}
+  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}{\mathrm{d} V_{q}} = 0.
+
+Consequently, if @\frac{\partial S_{z}}{\partial U_{s}} \neq 0 @ then
+  \frac{\mathrm{d} U_{s}}{\mathrm{d} V_{s}} =
+  -\frac{\frac{\partial S_{z}}{\partial V_{s}}}
+  {\frac{\partial S_{z}}{\partial U_{s}}}.
+
+If @ \frac{\partial S_{z}}{\partial V_{s}} \neq 0 @ then
+  \frac{\mathrm{d} V_{s}}{\mathrm{d} U_{s}} =
+  -\frac{\frac{\partial S_{z}}{\partial U_{s}}}
+  {\frac{\partial S_{z}}{\partial V_{s}}}
+
+Cases, when @ \frac{\partial S_{z}}{\partial U_{s}} =
+\frac{\partial S_{z}}{\partial V_{s}} = 0 @ are not consider here.
+The reason is written below.
+*/
+//=======================================================================
+Standard_Boolean IntPatch_SpecialPoints::ProcessSphere(const IntSurf_PntOn2S& thePtIso,
+                                                       const gp_Vec& theDUofPSurf,
+                                                       const gp_Vec& theDVofPSurf,
+                                                       const Standard_Boolean theIsReversed,
+                                                       const Standard_Real theVquad,
+                                                       Standard_Real& theUquad,
+                                                       Standard_Boolean& theIsIsoChoosen)
+{
+  theIsIsoChoosen = Standard_False;
+
+  //Vector with {@ \cos (U_{q}) @, @ \sin (U_{q}) @} coordinates.
+  //Ask to pay attention to the fact that this vector is always normalized.
+  gp_Vec2d aV1;
+
+  if ((Abs(theDUofPSurf.Z()) < Precision::PConfusion()) &&
+      (Abs(theDVofPSurf.Z()) < Precision::PConfusion()))
+  {
+    //Example of this case is an intersection of a plane with a sphere
+    //when the plane tangents the sphere in some pole (i.e. only one 
+    //intersection point, not line). In this case, U-coordinate of the
+    //sphere is undefined (can be realy anything).
+    //Another reason is that we have tangent zone around the pole
+    //(see bug #26576).
+    //Computation of correct value of theUquad is impossible.
+    //Therefore, (in order to return something) we will consider
+    //the intersection line goes along some isoline in neighborhood
+    //of the pole.
+
+#ifdef INTPATCH_ADDSPECIALPOINTS_DEBUG
+    cout << "Cannot find UV-coordinate for quadric in the pole."
+      " See considered comment above. IntPatch_SpecialPoints.cxx,"
+      " ProcessSphere(...)" << endl;
+#endif
+    Standard_Real aUIso = 0.0, aVIso = 0.0;
+    if (theIsReversed)
+      thePtIso.ParametersOnS2(aUIso, aVIso);
+    else
+      thePtIso.ParametersOnS1(aUIso, aVIso);
+
+    theUquad = aUIso;
+    theIsIsoChoosen = Standard_True;
+  }
+  else
+  {
+    if (Abs(theDUofPSurf.Z()) > Abs(theDVofPSurf.Z()))
+    {
+      const Standard_Real aDusDvs = theDVofPSurf.Z() / theDUofPSurf.Z();
+      aV1.SetCoord(theDUofPSurf.X()*aDusDvs - theDVofPSurf.X(),
+                   theDUofPSurf.Y()*aDusDvs - theDVofPSurf.Y());
+    }
+    else
+    {
+      const Standard_Real aDvsDus = theDUofPSurf.Z() / theDVofPSurf.Z();
+      aV1.SetCoord(theDVofPSurf.X()*aDvsDus - theDUofPSurf.X(),
+                   theDVofPSurf.Y()*aDvsDus - theDUofPSurf.Y());
+    }
+
+    aV1.Normalize();
+
+    if (Abs(aV1.X()) > Abs(aV1.Y()))
+      theUquad = Sign(asin(aV1.Y()), theVquad);
+    else
+      theUquad = Sign(acos(aV1.X()), theVquad);
+  }
+
+  return Standard_True;
+}
+
+//=======================================================================
+//function : ProcessCone
+//purpose  : 
+/*
+The intersection point (including the pole)
+must be satisfied to the following system:
+
+  \left\{\begin {matrix}
+  (V_{q}\sin(a) + R)*\cos(U_{q})) = S_{x}(U_{s}, V_{s})\\
+  (V_{q}\sin(a) + R)*\sin(U_{q})) = S_{y}(U_{s}, V_{s})\\
+  V_{q}\cos(a) = S_{z}(U_{s}, V_{s})
+  \end {matrix}\right,
+where 
+  R is the radius of the cone;
+  a is its semi-angle;
+  @S_{x}@, @S_{y}@ and @S_{z}@ are X, Y and Z-coordinates of thePSurf;
+  @U_{s}@ and @V_{s}@ are parameters on the parametric surface;
+  @U_{q}@ and @V_{q}@ are equal to theUquad and theVquad correspondingly.
+
+Consequently (from first two equations), 
+  \left\{\begin{matrix}
+  \cos(U_{q})=\frac{S_{x}(U_{s},V_{s})}{(V_{q}\sin(a)+R)}\\
+  \sin(U_{q})=\frac{S_{y}(U_{s}, V_{s})}{(V_{q}\sin(a)+R)}
+  \end{matrix}\right.
+
+For pole, the denominator of these two equations is equal to 0.
+Therefore, computation U_{q} directly is impossibly.
+
+Let @V_{q}@ tends to @\frac{-R}{\sin(a)})@.
+Then (indeterminate form is evaluated in accordance of L'Hospital rule),
+
+  \cos (U_{q}) = 
+  \lim_{V_{q} \to (\frac{-R}{\sin(a)})}\frac{S_{x}(U_{s},V_{s})}{(V_{q}\sin(a)+R)}=
+  \frac{1}{\sin(a)}* \lim_{V_{q} \to (\frac{-R}{\sin(a)})}\frac{dU_{s}}{dV_{q}}*
+  (\frac{\partial S_{x}}{\partial U_{s}}+\frac{\partial S_{x}}{\partial V_{s}}*
+  \frac{dV_{s}}{dU_{s}})=
+  \frac{1}{\sin(a)}* \lim_{V_{q} \to (\frac{-R}{\sin(a)})}\frac{dV_{s}}{dV_{q}}*
+  (\frac{\partial S_{x}}{\partial U_{s}}*
+  \frac{dU_{s}}{dV_{s}}+\frac{\partial S_{x}}{\partial V_{s}})
+
+Analogically for @\sin (U_{q})@ (@S_{x}@ is substituted to @S_{y}@).
+
+After differentiating 3rd equation of the system, we will obtain
+  \cos(a)=\frac{dS_{z}}{dV_{q}}=\frac{dU_{s}}{dV_{q}}*
+  (\frac{\partial S_{z}}{\partial U_{s}}+\frac{\partial S_{z}}{\partial V_{s}}*
+  \frac{dV_{s}}{dU_{s}})
+or
+  \frac{dU_{s}}{dV_{q}}=\frac{\cos(a)}{\frac{\partial S_{z}}{\partial U_{s}}+
+  \frac{\partial S_{z}}{\partial V_{s}}*\frac{dV_{s}}{dU_{s}}}
+
+After substituting we will obtain
+  \cos (U_{q}) = 
+  \cot(a)*\frac{\frac{\partial S_{x}}{\partial U_{s}}+\frac{\partial S_{x}}
+  {\partial V_{s}}*\frac{dV_{s}}{dU_{s}}}{\frac{\partial S_{z}}
+  {\partial U_{s}}+\frac{\partial S_{z}}{\partial V_{s}}*\frac{dV_{s}}{dU_{s}}}
+
+  \sin (U_{q}) =
+  \cot(a)*\frac{\frac{\partial S_{y}}{\partial U_{s}}+\frac{\partial S_{y}}
+  {\partial V_{s}}*\frac{dV_{s}}{dU_{s}}}{\frac{\partial S_{z}}
+  {\partial U_{s}}+\frac{\partial S_{z}}{\partial V_{s}}*\frac{dV_{s}}{dU_{s}}}
+
+So, we have obtained vector with coordinates {@ \cos (U_{q}) @, @ \sin (U_{q}) @}.
+Ask to pay attention to the fact that this vector is always normalized.
+And after normalization this vector will have coordinates
+  {\cos (U_{q}), \sin (U_{q})} = {dS_{x}, dS_{y}}.Normalized().
+
+It means that we have to compute a tangent to the intersection curve in
+the cone apex point. After that, just take its X- and Y-coordinates.
+
+However, we have to compute derivative @\frac{dV_{s}}{dU_{s}}@ in order
+to compute this vector. In order to find this derivative, we use the
+information about direction of tangent to the intersection curve.
+This tangent will be directed along the cone generatrix obtained by intersection
+of the cone with a plane tangent to 2nd (intersected) surface.
+*/
+//=======================================================================
+Standard_Boolean IntPatch_SpecialPoints::ProcessCone(const IntSurf_PntOn2S& thePtIso,
+                                                     const gp_Vec& theDUofPSurf,
+                                                     const gp_Vec& theDVofPSurf,
+                                                     const gp_Cone& theCone,
+                                                     const Standard_Boolean theIsReversed,
+                                                     Standard_Real& theUquad,
+                                                     Standard_Boolean& theIsIsoChoosen)
+{
+  theIsIsoChoosen = Standard_False;
+
+  // A plane tangent to 2nd (intersected) surface.
+  // Its normal.
+  const gp_XYZ aTgPlaneZ(theDUofPSurf.Crossed(theDVofPSurf).XYZ());
+  const Standard_Real aSqModTg = aTgPlaneZ.SquareModulus();
+  if (aSqModTg < Precision::SquareConfusion())
+  {
+    theIsIsoChoosen = Standard_True;
+  }
+
+  gp_XYZ aTgILine[2];
+  const Standard_Integer aNbTangent = !theIsIsoChoosen? 
+                          GetTangentToIntLineForCone(theCone.SemiAngle(),
+                                                     aTgPlaneZ.Divided(Sqrt(aSqModTg)),
+                                                     aTgILine) : 0;
+
+  if (aNbTangent == 0)
+  {
+    theIsIsoChoosen = Standard_True;
+  }
+  else
+  {
+    const Standard_Real aPeriod = M_PI + M_PI;
+    Standard_Real aUIso = 0.0, aVIso = 0.0;
+    if (theIsReversed)
+      thePtIso.ParametersOnS2(aUIso, aVIso);
+    else
+      thePtIso.ParametersOnS1(aUIso, aVIso);
+
+    aUIso = ElCLib::InPeriod(aUIso, 0.0, aPeriod);
+
+    // Sought U-parameter in the apex point
+
+    // For 2 possible parameters value,
+    // one will be chosen which is nearer
+    // to aUIso. Following variables will help to chose.
+    Standard_Real aMinDelta = RealLast();
+    for (Standard_Integer anIdx = 0; anIdx < aNbTangent; anIdx++)
+    {
+      // Vector {@\cos(a), \sin(a)@}
+      gp_Vec2d aVecCS(aTgILine[anIdx].X(), aTgILine[anIdx].Y());
+      const Standard_Real aSqMod = aVecCS.SquareMagnitude();
+      if (aSqMod < Precision::SquareConfusion())
+      {
+        theIsIsoChoosen = Standard_True;
+        break;
+      }
+
+      // Normalize
+      aVecCS.Divide(Sqrt(aSqMod));
+
+      // Angle in range [0, PI/2]
+      Standard_Real anUq = (Abs(aVecCS.X()) < Abs(aVecCS.Y())) ? ACos(Abs(aVecCS.X())) : ASin(Abs(aVecCS.Y()));
+
+      // Convert angles to the range [0, 2*PI]
+      if (aVecCS.Y() < 0.0)
+      {
+        if (aVecCS.X() > 0.0)
+        {
+          anUq = -anUq;
+        }
+        else
+        {
+          anUq += M_PI;
+        }
+      }
+      else if (aVecCS.X() < 0.0)
+      {
+        anUq = M_PI - anUq;
+      }
+
+      //Select the parameter the nearest to aUIso
+      anUq = ElCLib::InPeriod(anUq, 0.0, aPeriod);
+      Standard_Real aDelta = Abs(anUq - aUIso);
+      if (aDelta > M_PI)
+        aDelta = aPeriod - aDelta;
+
+      if (aDelta < aMinDelta)
+      {
+        aMinDelta = aDelta;
+        theUquad = anUq;
+      }
+    }
+  }
+
+  if (theIsIsoChoosen)
+  {
+#ifdef INTPATCH_ADDSPECIALPOINTS_DEBUG
+    cout << "Cannot find UV-coordinate for quadric in the pole."
+      " IntPatch_AddSpecialPoints.cxx, ProcessCone(...)" << endl;
+#endif
+    theIsIsoChoosen = Standard_True;
+
+    Standard_Real aUIso = 0.0, aVIso = 0.0;
+    if (theIsReversed)
+      thePtIso.ParametersOnS2(aUIso, aVIso);
+    else
+      thePtIso.ParametersOnS1(aUIso, aVIso);
+
+    theUquad = aUIso;
+    return Standard_True;
+  }
+  else
+  {
+    return Standard_True;
+  }
+
+  //return Standard_False;
+}
+
+//=======================================================================
+//function : GetTangentToIntLineForCone
+//purpose  : The following conditions must be satisfied:
+//1. The cone is represented in its canonical form.
+//2. The plane goes through the cone apex and has the normal vector thePlnNormal.
+//3. Vector thePlnNormal has already been normalized
+/*
+Let us enter the new coordinate system where the origin will be in the cone apex
+and axes are the same as in World-Coordinate-System (WCS).
+There the horizontal plane (which is parallel XY-plane) with the origin
+(0, 0, 1) will intersect the cone by the circle with center (0, 0, 1),
+direction {0, 0, 1} and radius tg(a) where a is the cone semi-angle.
+Its equation will be
+\left\{\begin{matrix}
+x(U_{n}) = \tan(a)*\cos(U_{n}) = \tan(a)*\frac{1-\tan^{2}(U_{n}/2)}{1+\tan^{2}(U_{n}/2)}\\
+y(U_{n}) = \tan(a)*\sin (U_{n}) = \tan(a)*\frac{2*\tan(U_{n}/2)}{1+\tan^{2}(U_{n}/2)}\\
+z(U_{n}) = 1
+\end{matrix}\right.
+
+The given plane has (in this coordinate system) location (0, 0, 0) and
+the same normal thePlnNormal=={nx,ny,nz}. Its equation is:
+nx*x+ny*y+nz*z==0
+
+After substitution circle's equation to the plane's equation
+we will obtain a quadratic equation
+aA*w^2 + 2*aB*w + aC = 0.
+*/
+//=======================================================================
+Standard_Integer IntPatch_SpecialPoints::GetTangentToIntLineForCone(const Standard_Real theConeSemiAngle,
+                                                                    const gp_XYZ& thePlnNormal,
+                                                                    gp_XYZ theResult[2])
+{
+  const Standard_Real aNullTol = Epsilon(1.0);
+  const Standard_Real aTanA = Tan(theConeSemiAngle);
+  const Standard_Real aA = thePlnNormal.Z() / aTanA - thePlnNormal.X();
+  const Standard_Real aB = thePlnNormal.Y();
+  const Standard_Real aC = thePlnNormal.Z() / aTanA + thePlnNormal.X();
+
+  if (Abs(aA) < aNullTol)
+  {
+    if (Abs(aB) > aNullTol)
+    {
+      //The plane goes along the cone generatrix.
+      GetTangent(theConeSemiAngle, -aC / (aB + aB), theResult[0]);
+      return 1;
+    }
+
+    //The cone and the plane have only one common point.
+    //It is the cone apex.
+    return 0;
+  }
+
+  //Discriminant of this equation is equal to 
+  Standard_Real aDiscr = thePlnNormal.Z() / Sin(theConeSemiAngle);
+  aDiscr = 1.0 - aDiscr*aDiscr;
+
+  if (Abs(aDiscr) < aNullTol)
+  {
+    //The plane goes along the cone generatrix.
+    // Attention! Mathematically, this cond. is equivalent to 
+    // above processed one (Abs(aA) < aNullTol && (Abs(aB) > aNullTol)).
+    // However, we separate this branch in order to eliminate numerical
+    // instability.
+
+    GetTangent(theConeSemiAngle, -aB / aA, theResult[0]);
+    return 1;
+  }
+  else if (aDiscr > 0.0)
+  {
+    const Standard_Real aRD = Sqrt(aDiscr);
+    GetTangent(theConeSemiAngle, (-aB+aRD)/aA, theResult[0]);
+    GetTangent(theConeSemiAngle, (-aB-aRD)/aA, theResult[1]);
+    return 2;
+  }
+
+  // We will never come here.
+  return 0;
+}
+
+//=======================================================================
 //function : AddSingularPole
 //purpose  : theQSurf is the surface possibly containing special point, 
 //            thePSurf is another surface to intersect.
@@ -333,22 +783,11 @@ Standard_Boolean IntPatch_SpecialPoints::
                       AddSingularPole(const Handle(Adaptor3d_HSurface)& theQSurf,
                                       const Handle(Adaptor3d_HSurface)& thePSurf,
                                       const IntSurf_PntOn2S& thePtIso,
-                                      const Standard_Real theTol,
                                       IntPatch_Point& theVertex,
-                                      IntSurf_PntOn2S& theAddedPoint,                                      
+                                      IntSurf_PntOn2S& theAddedPoint,
                                       const Standard_Boolean theIsReversed,
                                       const Standard_Boolean theIsReqRefCheck)
 {
-  const Standard_Real aUpPeriod = thePSurf->IsUPeriodic() ? thePSurf->UPeriod() : 0.0;
-  const Standard_Real aUqPeriod = theQSurf->IsUPeriodic() ? theQSurf->UPeriod() : 0.0;
-  const Standard_Real aVpPeriod = thePSurf->IsVPeriodic() ? thePSurf->VPeriod() : 0.0;
-  const Standard_Real aVqPeriod = theQSurf->IsVPeriodic() ? theQSurf->VPeriod() : 0.0;
-
-  const Standard_Real anArrOfPeriod[4] = {theIsReversed? aUpPeriod : aUqPeriod,
-                                          theIsReversed? aVpPeriod : aVqPeriod,
-                                          theIsReversed? aUqPeriod : aUpPeriod,
-                                          theIsReversed? aVqPeriod : aVpPeriod};
-
   //On parametric
   Standard_Real aU0 = 0.0, aV0 = 0.0;
   //aPQuad is Pole
@@ -379,13 +818,13 @@ Standard_Boolean IntPatch_SpecialPoints::
   }
 
   theQSurf->D0(aUquad, aVquad, aPQuad);
-
-  if (theIsReqRefCheck && (aPQuad.SquareDistance(theVertex.Value()) >= theTol*theTol))
+  const Standard_Real aTol = theVertex.Tolerance();
+  if (theIsReqRefCheck && (aPQuad.SquareDistance(theVertex.Value()) >= aTol*aTol))
   {
     return Standard_False;
   }
 
-  if(!IsPointOnSurface(thePSurf, aPQuad, theTol, aP0, aU0, aV0))
+  if (!IsPointOnSurface(thePSurf, aPQuad, aTol, aP0, aU0, aV0))
   {
     return Standard_False;
   }
@@ -398,25 +837,20 @@ Standard_Boolean IntPatch_SpecialPoints::
   else
     theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aUquad, aVquad, aU0, aV0);
 
-  Standard_Boolean isSame = Standard_False;
-
-  if (theAddedPoint.IsSame(theVertex.PntOn2S(), Precision::Confusion()))
-  {
-    isSame = Standard_True;
-  }
+  const Standard_Boolean isSame = theAddedPoint.IsSame(theVertex.PntOn2S(),
+                                                       Precision::Confusion());
 
   //Found pole does not exist in the Walking-line
   //It must be added there (with correct 2D-parameters)
 
-  //2D-parameters of theparametric surface have already been found (aU0, aV0).
+  //2D-parameters of thePSurf surface have already been found (aU0, aV0).
   //Let find 2D-parameters on the quadric.
 
-  //The algorithm depends on the type of the quadric. Here we consider a Sphere only.
-  //Analogical result can be made for another types (e.g. cone, but formulas will
-  //be different) in case of need.
+  //The algorithm depends on the type of the quadric.
+  //Here we consider a Sphere and cone only.
 
-  //First of all, we need in adjusting thePSurf in the coordinate system of the Sphere
-  //(in order to make the equation of the sphere maximal simple). However, as it will be
+  //First of all, we need in adjusting thePSurf in the coordinate system of the Sphere/Cone
+  //(in order to make its equation maximal simple). However, as it will be
   //shown later, thePSurf is used in algorithm in order to get its derivatives.
   //Therefore, for improving performance, transformation of these vectors is enough
   //(there is no point in transformation of full surface).
@@ -439,146 +873,19 @@ Standard_Boolean IntPatch_SpecialPoints::
 
   if(theQSurf->GetType() == GeomAbs_Sphere)
   {
-    //The intersection point (including the pole)
-    //must be satisfied to the following system:
-    
-    //    \left\{\begin{matrix}
-    //    R*\cos (U_{q})*\cos (V_{q})=S_{x}(U_{s},V_{s})
-    //    R*\sin (U_{q})*\cos (V_{q})=S_{y}(U_{s},V_{s})
-    //    R*\sin (V_{q})=S_{z}(U_{s},V_{s})
-    //    \end{matrix}\right,
-    //where 
-    //  R is the radius of the sphere;
-    //  @S_{x}@, @S_{y}@ and @S_{z}@ are X, Y and Z-coordinates of thePSurf;
-    //  @U_{s}@ and @V_{s}@ are equal to aU0 and aV0 corespondingly;
-    //  @U_{q}@ and @V_{q}@ are equal to aUquad and aVquad corespondingly.
-
-    //Consequently (from first two equations), 
-    //  \left\{\begin{matrix}
-    //  \cos (U_{q}) = \frac{S_{x}(U_{s},V_{s})}{R*\cos (V_{q})}
-    //  \sin (U_{q}) = \frac{S_{y}(U_{s},V_{s})}{R*\cos (V_{q})}
-    //  \end{matrix}\right.
-
-    //For pole, 
-    //  V_{q}=\pm \pi /2 \Rightarrow \cos (V_{q}) = 0 (denominator is equal to 0).
-
-    //Therefore, computation U_{q} directly is impossibly.
-    //
-    //Let @V_{q}@ tends to @\pm \pi /2@.
-    //Then (indeterminate form is evaluated in accordance of L'Hospital rule),
-    //  \cos (U_{q}) = \lim_{V_{q} \to (\pi /2-0)} 
-    //  \frac{S_{x}(U_{s},V_{s})}{R*\cos (V_{q})}= 
-    //  -\lim_{V_{q} \to (\pi /2-0)}
-    //  \frac{\frac{\partial S_{x}}
-    //  {\partial U_{s}}*\frac{\mathrm{d} U_{s}} 
-    //  {\mathrm{d} V_{q}}+\frac{\partial S_{x}} 
-    //  {\partial V_{s}}*\frac{\mathrm{d} V_{s}} 
-    //  {\mathrm{d} V_{q}}}{R*\sin (V_{q})} =  
-    //  -\frac{1}{R}*\frac{\mathrm{d} U_{s}}
-    //  {\mathrm{d} V_{q}}*(\frac{\partial S_{x}} 
-    //  {\partial U_{s}}+\frac{\partial S_{x}}
-    //  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}
-    //  {\mathrm{d} U_{s}}) =
-    //  -\frac{1}{R}*\frac{\mathrm{d} V_{s}}
-    //  {\mathrm{d} V_{q}}*(\frac{\partial S_{x}} 
-    //  {\partial U_{s}}*\frac{\mathrm{d} U_{s}}
-    //  {\mathrm{d} V_{s}}+\frac{\partial S_{x}}
-    //  {\partial V_{s}}).
-
-    //Analogicaly for @\sin (U_{q})@ (@S_{x}@ is substituted to @S_{y}@).
-
-    //Let mean, that
-    //  \cos (U_{q}) \left | _{V_{q} \to (-\pi /2+0)} = \cos (U_{q}) \left | _{V_{q} \to (\pi /2-0)}
-    //  \sin (U_{q}) \left | _{V_{q} \to (-\pi /2+0)} = \sin (U_{q}) \left | _{V_{q} \to (\pi /2-0)}
-
-    //From the 3rd equation of the system, we obtain
-    //  \frac{\mathrm{d} (R*\sin (V_{q}))}{\mathrm{d} V_{q}} =
-    //  \frac{\mathrm{d} S_{z}(U_{s},V_{s})}{\mathrm{d} V_{q}}
-    //or
-    //  R*\cos (V_{q}) = \frac{\partial S_{z}}{\partial U_{s}}*
-    //  \frac{\mathrm{d} U_{s}} {\mathrm{d} V_{q}}+\frac{\partial S_{z}}
-    //  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}{\mathrm{d} V_{q}}.
-
-    //If @V_{q}=\pm \pi /2@, then
-    //  \frac{\partial S_{z}}{\partial U_{s}}*
-    //  \frac{\mathrm{d} U_{s}} {\mathrm{d} V_{q}}+\frac{\partial S_{z}}
-    //  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}{\mathrm{d} V_{q}} = 0.
-
-    //Consequently, if @\frac{\partial S_{z}}{\partial U_{s}} \neq 0 @ then
-    //  \frac{\mathrm{d} U_{s}}{\mathrm{d} V_{s}} =
-    //  -\frac{\frac{\partial S_{z}}{\partial V_{s}}}
-    //  {\frac{\partial S_{z}}{\partial U_{s}}}.
-
-    //If @ \frac{\partial S_{z}}{\partial V_{s}} \neq 0 @ then
-    //  \frac{\mathrm{d} V_{s}}{\mathrm{d} U_{s}} =
-    //  -\frac{\frac{\partial S_{z}}{\partial U_{s}}}
-    //  {\frac{\partial S_{z}}{\partial V_{s}}}
-
-    //Cases, when @ \frac{\partial S_{z}}{\partial U_{s}} = 
-    //\frac{\partial S_{z}}{\partial V_{s}} = 0 @ are not consider here.
-    //The reason is written below.
-
-    //Vector with {@ \cos (U_{q}) @, @ \sin (U_{q}) @} coordinates.
-    //Ask to pay attention to the fact that this vector is always normalyzed.
-    gp_Vec2d aV1;
-
-    if( (Abs(aVecDu.Z()) < Precision::PConfusion()) &&
-        (Abs(aVecDv.Z()) < Precision::PConfusion()))
+    if (!ProcessSphere(thePtIso, aVecDu, aVecDv, theIsReversed, 
+                        aVquad, aUquad, isIsoChoosen))
     {
-      //Example of this case is an intersection of a plane with a sphere
-      //when the plane tangents the sphere in some pole (i.e. only one 
-      //intersection point, not line). In this case, U-coordinate of the
-      //sphere is undefined (can be realy anything).
-      //Another reason is that we have tangent zone around the pole
-      //(see bug #26576).
-      //Computation of correct value of aUquad is impossible.
-      //Therefore, (in oreder to return something) we will consider
-      //the intersection line goes along some isoline in neighbourhood
-      //of the pole.
-
-#ifdef INTPATCH_ADDSPECIALPOINTS_DEBUG
-      cout << "Cannot find UV-coordinate for quadric in the pole."
-        " See considered comment above. IntPatch_AddSpecialPoints.cxx,"
-        " AddSingularPole(...)" << endl;
-#endif
-      Standard_Real aUIso = 0.0, aVIso = 0.0;
-      if(theIsReversed)
-        thePtIso.ParametersOnS2(aUIso, aVIso);
-      else
-        thePtIso.ParametersOnS1(aUIso, aVIso);
-
-      aUquad = aUIso;
-      isIsoChoosen = Standard_True;
-    }
-    else
-    {
-      if(Abs(aVecDu.Z()) > Abs(aVecDv.Z()))
-      {
-        const Standard_Real aDusDvs = aVecDv.Z()/aVecDu.Z();
-        aV1.SetCoord( aVecDu.X()*aDusDvs - aVecDv.X(),
-                      aVecDu.Y()*aDusDvs - aVecDv.Y());
-      }
-      else
-      {
-        const Standard_Real aDvsDus = aVecDu.Z()/aVecDv.Z();
-        aV1.SetCoord( aVecDv.X()*aDvsDus - aVecDu.X(),
-                      aVecDv.Y()*aDvsDus - aVecDu.Y());
-      }
-
-      aV1.Normalize();
-
-      if(Abs(aV1.X()) > Abs(aV1.Y()))
-        aUquad = Sign(asin(aV1.Y()), aVquad);
-      else
-        aUquad = Sign(acos(aV1.X()), aVquad);
+      return Standard_False;
     }
   }
   else //if(theQSurf->GetType() == GeomAbs_Cone)
   {
-    // This case is not processed. However,
-    // it can be done using the same algorithm
-    // as for sphere (formulas will be different).
-    return Standard_False;
+    if (!ProcessCone(thePtIso, aVecDu, aVecDv, theQSurf->Cone(),
+                      theIsReversed, aUquad, isIsoChoosen))
+    {
+      return Standard_False;
+    }
   }
 
   if(theIsReversed)
@@ -594,6 +901,16 @@ Standard_Boolean IntPatch_SpecialPoints::
 
   if (!isIsoChoosen)
   {
+    Standard_Real anArrOfPeriod[4];
+    if (theIsReversed)
+    {
+      IntSurf::SetPeriod(thePSurf, theQSurf, anArrOfPeriod);
+    }
+    else
+    {
+      IntSurf::SetPeriod(theQSurf, thePSurf, anArrOfPeriod);
+    }
+
     AdjustPointAndVertex(theVertex.PntOn2S(), anArrOfPeriod, theAddedPoint);
   }
   else
@@ -606,7 +923,32 @@ Standard_Boolean IntPatch_SpecialPoints::
 
 //=======================================================================
 //function : ContinueAfterSpecialPoint
-//purpose  : 
+//purpose  : If the last point of the line is the pole of the quadric then
+//            the Walking-line has been broken in this point.
+//           However, new line must start from this point. Here we must
+//            find 2D-coordinates of "this new" point.
+/*
+The inters. line in the neighborhood of the Apex/Pole(s) can be
+approximated by the intersection result of the Cone/Sphere with
+the plane going through the Apex/Pole and being tangent to the
+2nd intersected surface. This intersection result is well known.
+
+In case of sphere, the inters. result is a circle.
+If we go along this circle and across the Pole then U-parameter of
+the sphere (@U_{q}@) will change to +/-PI.
+
+In case of cone, the inters. result is two intersected lines (which
+can be merged to one in a special case when the plane goes along
+some generatrix of the cone). The direction of these lines
+are computed by GetTangentToIntLineForCone(...) method).
+
+When the real (not lines) inters. curve goes through the cone apex then
+two variants are possible:
+a) The tangent line to the inters. curve will be left. In this case
+U-parameter of the cone (@U_{q}@) will be change to +/-PI.
+b) Another line (as inters. result of cone + plane) will tangent
+to the inters. curve. In this case @U_{q}@ must be recomputed.
+*/
 //=======================================================================
 Standard_Boolean IntPatch_SpecialPoints::
                   ContinueAfterSpecialPoint(const Handle(Adaptor3d_HSurface)& theQSurf,
@@ -620,36 +962,63 @@ Standard_Boolean IntPatch_SpecialPoints::
   if(theSPType == IntPatch_SPntNone)
     return Standard_False;
 
-  //If the last point of the line is the pole of the quadric.
-  //In this case, Walking-line has been broken in this point.
-  //However, new line must start from this point. Here we must
-  //find its 2D-coordinates.
-
-  //For sphere and cone, some intersection point is satisfied to the system
-  //  \cos(U_{q}) = S_{x}(U_{s},V_{s})/F(V_{q}) 
-  //  \sin(U_{q}) = S_{y}(U_{s},V_{s})/F(V_{q}) 
-
-  //where 
-  //  @S_{x}@, @S_{y}@ are X and Y-coordinates of thePSurf;
-  //  @U_{s}@ and @V_{s}@ are UV-parameters on thePSurf;
-  //  @U_{q}@ and @V_{q}@ are UV-parameters on theQSurf;
-  //  @F(V_{q}) @ is some function, which value independs on @U_{q}@
-  //              (form of this function depends on the type of the quadric).
-
-  //When we go through the pole/apex, the function @F(V_{q}) @ changes sign.
-  //Therefore, some cases are possible, when only @\cos(U_{q}) @ or
-  //only @ \sin(U_{q}) @ change sign.
-
-  //Consequently, when the line goes throug the pole, @U_{q}@ can be
-  //changed on @\pi /2 @ (but not less).
-    
   if(theNewPoint.IsSame(theRefPt, Precision::Confusion(), theTol2D))
   {
     return Standard_False;
   }
 
-  //Here, in case of pole/apex adding, we forbid "jumping" between two neighbor
-  //Walking-point with step greater than pi/4
+  if ((theSPType == IntPatch_SPntPole) && (theQSurf->GetType() == GeomAbs_Cone))
+  {
+    //Check if the condition b) is satisfied.
+    //Repeat the same steps as in 
+    //IntPatch_SpecialPoints::AddSingularPole(...) method.
+
+    //On parametric
+    Standard_Real aU0 = 0.0, aV0 = 0.0;
+    //On quadric
+    Standard_Real aUquad = 0.0, aVquad = 0.0;
+
+    if (theIsReversed)
+      theNewPoint.Parameters(aU0, aV0, aUquad, aVquad);
+    else
+      theNewPoint.Parameters(aUquad, aVquad, aU0, aV0);
+
+    gp_Pnt aPtemp;
+    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(theQSurf->Cone().Position());
+
+    //Derivatives of transformed thePSurf
+    aVecDu.Transform(aTr);
+    aVecDv.Transform(aTr);
+
+    Standard_Boolean isIsoChoosen = Standard_False;
+    ProcessCone(theRefPt, aVecDu, aVecDv, theQSurf->Cone(),
+                theIsReversed, aUquad, isIsoChoosen);
+
+    theNewPoint.SetValue(!theIsReversed, aUquad, aVquad);
+  }
+
+  //As it has already been said, in case of going through the Pole/Apex,
+  //U-parameter of the quadric surface will change to +/-PI. This rule has some
+  //exceptions:
+  //1. When 2nd surface has C0-continuity in the point common with the Apex/Pole.
+  //    In this case, the tangent line to the intersection curve after the Apex/Pole
+  //    must be totally recomputed according to the new derivatives of the 2nd surface.
+  //    Currently, it is not implemented but will be able to be done after the
+  //    corresponding demand.
+  //2. The inters. curve has C1 continuity but huge curvature in the point common with 
+  //    the Apex/Pole. Existing inters. algorithm does not allow putting many points
+  //    near to the Apex/Pole in order to cover this "sharp" piece of the inters. curve.
+  //    Therefore, we use adjusting U-parameter of the quadric surface with
+  //    period PI/2 instead of 2PI. It does not have any mathematical idea
+  //    but allows creating WLine with more or less uniform distributed points.
+  //    In other words, we forbid "jumping" between two neighbor Walking-points
+  //    with step greater than PI/4.
+
   const Standard_Real aPeriod = (theSPType == IntPatch_SPntPole)? M_PI_2 : 2.0*M_PI;
 
   const Standard_Real aUpPeriod = thePSurf->IsUPeriodic() ? thePSurf->UPeriod() : 0.0;