0024772: Modeling Algorithms - Intersection of cylinder and cone produces incorrect results

GeomInt/GeomInt_IntSS.cxx - handle flat cone case
IntPatch/IntPatch_Intersection.hxx, IntPatch/IntPatch_Intersection.cxx - method for preparing surfaces is added
bugs/lowalgos/intss/bug24772 - test case is changed according new behavior

src/Adaptor3d/Adaptor3d_TopolTool.cxx

@@ -30,34 +30,6 @@ IMPLEMENT_STANDARD_RTTIEXT(Adaptor3d_TopolTool,Standard_Transient)
 
 #define myInfinite Precision::Infinite()
 
-static void GetConeApexParam(const gp_Cone& C, Standard_Real& U, Standard_Real& V)
-{
-  const gp_Ax3& Pos = C.Position();
-  Standard_Real Radius = C.RefRadius();
-  Standard_Real SAngle = C.SemiAngle();
-  const gp_Pnt& P = C.Apex();
- 
-  gp_Trsf T;
-  T.SetTransformation (Pos);
-  gp_Pnt Ploc = P.Transformed (T);
-
-  if(Ploc.X() ==0.0  &&  Ploc.Y()==0.0 ) {
-    U = 0.0;
-  }
-  else if ( -Radius > Ploc.Z()* Tan(SAngle) ) {
-    // the point is at the `wrong` side of the apex
-    U = atan2(-Ploc.Y(), -Ploc.X());
-  }
-  else {
-    U = atan2(Ploc.Y(),Ploc.X());
-  }
-  if      (U < -1.e-16)  U += (M_PI+M_PI);
-  else if (U < 0)        U = 0;
-
-  V =  sin(SAngle) * ( Ploc.X() * cos(U) + Ploc.Y() * sin(U) - Radius)
-    + cos(SAngle) * Ploc.Z();
-}
-
 
 Adaptor3d_TopolTool::Adaptor3d_TopolTool ()
 : myNbSamplesU(-1),
@@ -1376,3 +1348,39 @@ Standard_Boolean Adaptor3d_TopolTool::IsUniformSampling() const
     return Standard_False;
   return Standard_True;
 }
+
+//=======================================================================
+//function : GetConeApexParam
+//purpose  : Computes the cone's apex parameters
+//=======================================================================
+void Adaptor3d_TopolTool::GetConeApexParam (const gp_Cone& theC, Standard_Real& theU, Standard_Real& theV)
+{
+  const gp_Ax3& Pos = theC.Position();
+  Standard_Real Radius = theC.RefRadius();
+  Standard_Real SAngle = theC.SemiAngle();
+  const gp_Pnt& P = theC.Apex();
+
+  gp_Trsf T;
+  T.SetTransformation(Pos);
+  gp_Pnt Ploc = P.Transformed(T);
+
+  if (Ploc.X() == 0.0 && Ploc.Y() == 0.0)
+  {
+    theU = 0.0;
+  }
+  else if (-Radius > Ploc.Z() * Tan(SAngle))
+  {
+    // the point is at the `wrong` side of the apex
+    theU = atan2(-Ploc.Y(), -Ploc.X());
+  }
+  else
+  {
+    theU = atan2(Ploc.Y(), Ploc.X());
+  }
+
+  if (theU < -1.e-16)  theU += (M_PI + M_PI);
+  else if (theU < 0)   theU = 0;
+
+  theV = sin(SAngle) * (Ploc.X() * cos(theU) + Ploc.Y() * sin(theU) - Radius)
+         + cos(SAngle) * Ploc.Z();
+}

src/Adaptor3d/Adaptor3d_TopolTool.hxx

@@ -146,6 +146,12 @@ public:
   //! Returns true if provide uniform sampling of points.
   Standard_EXPORT virtual Standard_Boolean IsUniformSampling() const;
 
+  //! Computes the cone's apex parameters.
+  //! @param[in] theC conical surface
+  //! @param[in] theU U parameter of cone's apex
+  //! @param[in] theV V parameter of cone's apex
+  Standard_EXPORT static void GetConeApexParam (const gp_Cone& theC, Standard_Real& theU, Standard_Real& theV);
+
   DEFINE_STANDARD_RTTIEXT(Adaptor3d_TopolTool,Standard_Transient)
 
 protected:

src/GeomInt/GeomInt_IntSS.cxx

@@ -77,38 +77,78 @@ void GeomInt_IntSS::Perform(const Handle(Geom_Surface)& S1,
   myTolReached2d = myTolReached3d = 0.0;
   myNbrestr = 0;
   sline.Clear();
-  Handle(Adaptor3d_TopolTool) dom1 = new Adaptor3d_TopolTool(myHS1);
-  Handle(Adaptor3d_TopolTool) dom2 = new Adaptor3d_TopolTool(myHS2);
-  myLConstruct.Load(dom1,dom2,myHS1,myHS2);
   
   Standard_Real TolArc = Tol;
   Standard_Real TolTang = Tol;
-  Standard_Real UVMaxStep = IntPatch_Intersection::DefineUVMaxStep(myHS1, dom1, myHS2, dom2);
   Standard_Real Deflection = 0.1;
   if (myHS1->GetType() == GeomAbs_BSplineSurface && myHS2->GetType() == GeomAbs_BSplineSurface)
   {
     Deflection /= 10.;
   }
 
+  Handle(Adaptor3d_TopolTool) dom1 = new Adaptor3d_TopolTool (myHS1);
+  Handle(Adaptor3d_TopolTool) dom2 = new Adaptor3d_TopolTool (myHS2);
 
-  myIntersector.SetTolerances(TolArc,TolTang,UVMaxStep,Deflection);
+  NCollection_Vector< Handle(Adaptor3d_Surface)> aVecHS1;
+  NCollection_Vector< Handle(Adaptor3d_Surface)> aVecHS2;
 
-  if(myHS1 == myHS2) {
-    myIntersector.Perform(myHS1,dom1,TolArc,TolTang);
+  if (myHS1 == myHS2)
+  {
+    aVecHS1.Append (myHS1);
+    aVecHS2.Append (myHS2);
   }
-  else if (!useStart) {
-    myIntersector.Perform(myHS1,dom1,myHS2,dom2,TolArc,TolTang);
-  }
-  else {
-    myIntersector.Perform(myHS1,dom1,myHS2,dom2,U1,V1,U2,V2,TolArc,TolTang);
+  else
+  {
+    myIntersector.PrepareSurfaces (myHS1, dom1, myHS2, dom2, Tol, aVecHS1, aVecHS2);
   }
 
-  // ============================================================
-  if (myIntersector.IsDone()) {
-    const Standard_Integer nblin = myIntersector.NbLines();
-    for (Standard_Integer i=1; i<= nblin; i++)
+  for (Standard_Integer aNumOfHS1 = 0; aNumOfHS1 < aVecHS1.Length(); aNumOfHS1++)
+  {
+    const Handle(Adaptor3d_Surface)& aHS1 = aVecHS1.Value (aNumOfHS1);
+
+    for (Standard_Integer aNumOfHS2 = 0; aNumOfHS2 < aVecHS2.Length(); aNumOfHS2++)
     {
-      MakeCurve(i,dom1,dom2,Tol,Approx,ApproxS1,ApproxS2);
+      const Handle(Adaptor3d_Surface)& aHS2 = aVecHS2.Value (aNumOfHS2);
+
+      Handle(Adaptor3d_TopolTool) aDom1 = new Adaptor3d_TopolTool (aHS1);
+      Handle(Adaptor3d_TopolTool) aDom2 = new Adaptor3d_TopolTool (aHS2);
+
+      myLConstruct.Load (aDom1 ,aDom2,
+                         Handle(GeomAdaptor_Surface)::DownCast (aHS1),
+                         Handle(GeomAdaptor_Surface)::DownCast (aHS2));
+
+      Standard_Real UVMaxStep = IntPatch_Intersection::DefineUVMaxStep (aHS1, aDom1, aHS2, aDom2);
+
+      myIntersector.SetTolerances (TolArc, TolTang, UVMaxStep, Deflection);
+
+      if (aHS1 == aHS2)
+      {
+        myIntersector.Perform (aHS1, aDom1, TolArc, TolTang);
+      }
+      else if (!useStart)
+      {
+        myIntersector.Perform (aHS1, aDom1, aHS2, aDom2, TolArc, TolTang);
+      }
+      else
+      {
+        TopAbs_State aState1 = aDom1->Classify (gp_Pnt2d (U1, V1), Tol);
+        TopAbs_State aState2 = aDom2->Classify (gp_Pnt2d (U2, V2), Tol);
+
+        if ((aState1 == TopAbs_IN || aState1 == TopAbs_ON) &&
+            (aState2 == TopAbs_IN || aState2 == TopAbs_ON))
+        {
+          myIntersector.Perform (aHS1, aDom1, aHS2, aDom2, U1, V1, U2, V2, TolArc, TolTang);
+        }
+      }
+
+      // ============================================================
+      if (myIntersector.IsDone()) {
+        const Standard_Integer nblin = myIntersector.NbLines();
+        for (Standard_Integer i = 1; i <= nblin; i++)
+        {
+          MakeCurve (i, aDom1, aDom2, Tol, Approx, ApproxS1, ApproxS2);
+        }
+      }
     }
   }
 }

src/IntPatch/IntPatch_Intersection.cxx

@@ -1923,3 +1923,72 @@ Standard_Real IntPatch_Intersection::DefineUVMaxStep(
   return anUVMaxStep;
 }
 
+//=======================================================================
+//function : splitCone
+//purpose  : Splits cone by the apex
+//=======================================================================
+static void splitCone(
+  const Handle(Adaptor3d_Surface)& theS,
+  const Handle(Adaptor3d_TopolTool)& theD,
+  const Standard_Real theTol,
+  NCollection_Vector< Handle(Adaptor3d_Surface)>& theVecHS)
+{
+  if (theS->GetType() != GeomAbs_Cone)
+  {
+    throw Standard_NoSuchObject("IntPatch_Intersection : Surface is not Cone");
+  }
+
+  gp_Cone aCone = theS->Cone();
+
+  Standard_Real aU0, aV0;
+  Adaptor3d_TopolTool::GetConeApexParam (aCone, aU0, aV0);
+
+  TopAbs_State aState = theD->Classify (gp_Pnt2d (aU0, aV0), theTol);
+
+  if (aState == TopAbs_IN || aState == TopAbs_ON)
+  {
+    const Handle(Adaptor3d_Surface) aHSDn = theS->VTrim (theS->FirstVParameter(), aV0, Precision::PConfusion());
+    const Handle(Adaptor3d_Surface) aHSUp = theS->VTrim (aV0, theS->LastVParameter(), Precision::PConfusion());
+
+    theVecHS.Append (aHSDn);
+    theVecHS.Append (aHSUp);
+  }
+  else
+  {
+    theVecHS.Append (theS);
+  }
+}
+
+//=======================================================================
+//function : PrepareSurfaces
+//purpose  : Prepares surfaces for intersection
+//=======================================================================
+void IntPatch_Intersection::PrepareSurfaces(
+  const Handle(Adaptor3d_Surface)& theS1,
+  const Handle(Adaptor3d_TopolTool)& theD1,
+  const Handle(Adaptor3d_Surface)& theS2,
+  const Handle(Adaptor3d_TopolTool)& theD2,
+  const Standard_Real theTol,
+  NCollection_Vector< Handle(Adaptor3d_Surface)>& theVecHS1,
+  NCollection_Vector< Handle(Adaptor3d_Surface)>& theVecHS2)
+{
+  if ((theS1->GetType() == GeomAbs_Cone) && (Abs (M_PI / 2. - Abs (theS1->Cone().SemiAngle())) < theTol))
+  {
+    splitCone (theS1, theD1, theTol, theVecHS1);
+  }
+  else
+  {
+    theVecHS1.Append (theS1);
+  }
+
+  if ((theS2->GetType() == GeomAbs_Cone) && (Abs (M_PI / 2. - Abs (theS2->Cone().SemiAngle())) < theTol))
+  {
+    splitCone (theS2, theD2, theTol, theVecHS2);
+  }
+  else
+  {
+    theVecHS2.Append (theS2);
+  }
+}
+
+

src/IntPatch/IntPatch_Intersection.hxx

@@ -22,6 +22,7 @@
 #include <IntPatch_SequenceOfLine.hxx>
 #include <IntSurf_ListOfPntOn2S.hxx>
 #include <GeomAbs_SurfaceType.hxx>
+#include <NCollection_Vector.hxx>
 
 class Adaptor3d_TopolTool;
 
@@ -149,6 +150,15 @@ public:
     const Handle(Adaptor3d_Surface)&  theS2,
     const Handle(Adaptor3d_TopolTool)& theD2);
 
+  //! Prepares surfaces for intersection
+  Standard_EXPORT static void PrepareSurfaces(
+      const Handle(Adaptor3d_Surface)& theS1,
+      const Handle(Adaptor3d_TopolTool)& theD1,
+      const Handle(Adaptor3d_Surface)& theS2,
+      const Handle(Adaptor3d_TopolTool)& theD2,
+      const Standard_Real Tol,
+      NCollection_Vector< Handle(Adaptor3d_Surface)>& theSeqHS1,
+      NCollection_Vector< Handle(Adaptor3d_Surface)>& theSeqHS2);
 
 protected:
 

tests/lowalgos/intss/bug24772

@@ -8,38 +8,61 @@ puts ""
 
 restore [locate_data_file bug24772_s1.draw] s1
 restore [locate_data_file bug24772_s2.draw] s2
+restore [locate_data_file bug24772_s3.draw] s3
 
 smallview +X+Y
 fit
 zoom 6
 
-set bug_info [intersect r s1 s2]
+# case of unmatched axes
+set bug_info_1 [intersect r1 s1 s2]
 
-if {[llength $bug_info] != 4} {
-#  puts "ERROR: OCC24722 is reproduced."
+# case of matched axes
+set bug_info_2 [intersect r2 s2 s3]
+
+if {[llength $bug_info_1] != 3} {
+  puts "ERROR: OCC24722 is reproduced."
 } else {
-  # snapshot r_1
+  # snapshot r1_1
   clear
   display s1
   display s2
-  display r_1
-  xwd $imagedir/${casename}_r_1.png
-  # snapshot r_2
+  display r1_1
+  xwd $imagedir/${casename}_r1_1.png
+  # snapshot r1_2
   clear
   display s1
   display s2
-  display r_2
-  xwd $imagedir/${casename}_r_2.png
-  # snapshot r_3
+  display r1_2
+  xwd $imagedir/${casename}_r1_2.png
+  # snapshot r1_3
   clear
   display s1
   display s2
-  display r_3
-  xwd $imagedir/${casename}_r_3.png
-  # snapshot r_4
-  clear
-  display s1
-  display s2
-  display r_4
-  xwd $imagedir/${casename}_r_4.png
+  display r1_3
+  xwd $imagedir/${casename}_r1_3.png
+}
+
+
+if {[llength $bug_info_2] != 3} {
+  puts "ERROR: OCC24722 is reproduced."
+} else {
+  # snapshot r2_1
+  clear
+  display s2
+  display s3
+  display r2_1
+  xwd $imagedir/${casename}_r2_1.png
+  # snapshot r2_2
+  clear
+  display s2
+  display s3
+  display r2_2
+  xwd $imagedir/${casename}_r2_2.png
+  # snapshot r2_3
+  clear
+  display s2
+  display s3
+  display r2_3
+  xwd $imagedir/${casename}_r2_3.png
 }