0024915: Wrong intersection curves between two cylinders

Existing method of Cylinder-Cylinder intersection computing is based on finding the analytic line (as a function of one argument) and converting one into the walking-line with set of equidistant (along the line parameter) points.

The main advantage of applied method is using adaptively computed step. Necessary step is computed into every point of the obtained walking-line. At that we receive final walking-line directly (without preliminary analytic line) and we determine moments more precisely, when it should be split (see IntPatch_ImpImpIntersection_4.gxx).

The main disadvantages is bad working this method for non-trimmed cylinders (with infinite bounds), because step value is depend on the boundaries values.

More over, new method always returns walking-line, while intersection result can be an analytic curve (lines, circle, ellipse). That is NO good. Therefore, analytic curve is computed by existing method.

In conclusion, in spite of covering almost all more often meeting cases, new method has limited application. Then we should use the existing old method.

Additionally, method MinMax() is added (see Standard_Real.hxx file). It uses into new algorithm.

Some test cases is changed according to their new behavior.

Test case for issue CR24915 is added.

Into GeometryTest_APICommands.cxx only tabulations were chaged.

"Extending" of isolines (see Geom2dHatch_Hatcher.cxx).

Small correction of test case for issue CR24915.

src/IntPatch/IntPatch_ImpImpIntersection.cdl

@@ -54,7 +54,8 @@ is
     Perform (me: in out;
              S1: HSurface from Adaptor3d; D1: TopolTool from Adaptor3d;
              S2: HSurface from Adaptor3d; D2: TopolTool from Adaptor3d;
-             TolArc,TolTang: Real from Standard)
+             TolArc,TolTang: Real from Standard;
+             isTheTrimmed: Boolean from Standard = Standard_False)
 
 	raises ConstructionError from Standard    
 

src/IntPatch/IntPatch_ImpImpIntersection_1.gxx

@@ -14,6 +14,8 @@
 // Alternatively, this file may be used under the terms of Open CASCADE
 // commercial license or contractual agreement.
 
+#include <Bnd_Box2d.hxx>
+
 static Standard_Boolean IntPP (const IntSurf_Quadric&,
 			       const IntSurf_Quadric&,
 			       const Standard_Real,
@@ -77,6 +79,18 @@ static Standard_Boolean IntCyCy(const IntSurf_Quadric&,
 				IntPatch_SequenceOfLine&,
 				IntPatch_SequenceOfPoint&);
 
+static Standard_Boolean IntCyCyTrim(const IntSurf_Quadric& theQuad1,
+                                    const IntSurf_Quadric& theQuad2,
+                                    const Standard_Real theTol3D,
+                                    const Standard_Real theTol2D,
+                                    const Bnd_Box2d& theUVSurf1,
+                                    const Bnd_Box2d& theUVSurf2,
+                                    const Standard_Boolean isTheReverse,
+                                    Standard_Boolean& isTheEmpty,
+                                    IntPatch_SequenceOfLine& theSlin,
+                                    IntPatch_SequenceOfPoint& theSPnt);
+
+
 static Standard_Boolean IntCySp(const IntSurf_Quadric&,
 				const IntSurf_Quadric&,
 				const Standard_Real,

src/IntPatch/IntPatch_ImpImpIntersection_2.gxx

@@ -50,8 +50,10 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
                                           const Handle(Adaptor3d_HSurface)&  S2,
                                           const Handle(Adaptor3d_TopolTool)& D2,
                                           const Standard_Real TolArc,
-                                          const Standard_Real TolTang) {
+                                          const Standard_Real TolTang,
+                                          const Standard_Boolean isTheTrimmed) {
   done = Standard_False;
+  Standard_Boolean isTrimmed = isTheTrimmed;
   spnt.Clear();
   slin.Clear();
 
@@ -142,13 +144,61 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
       break;
     }
     //
-    case 22: { // Cylinder/Cylinder
-      if (!IntCyCy(quad1, quad2, TolTang, empt, SameSurf, multpoint, slin, spnt)) {
-        return;
+    case 22:
+      { // Cylinder/Cylinder
+        Standard_Boolean isDONE = Standard_False;
+        
+        if(!isTrimmed)
+        {
+          isDONE = IntCyCy(quad1, quad2, TolTang, empt,
+                            SameSurf, multpoint, slin, spnt);
+        }
+        else
+        {
+          Bnd_Box2d aBox1, aBox2;
+
+          const Standard_Real aU1f = S1->FirstUParameter();
+          const Standard_Real aU1l = S1->LastUParameter();
+          const Standard_Real aU2f = S2->FirstUParameter();
+          const Standard_Real aU2l = S2->LastUParameter();
+
+          aBox1.Add(gp_Pnt2d(aU1f, S1->FirstVParameter()));
+          aBox1.Add(gp_Pnt2d(aU1l, S1->LastVParameter()));
+          aBox2.Add(gp_Pnt2d(aU2f, S2->FirstVParameter()));
+          aBox2.Add(gp_Pnt2d(aU2l, S2->LastVParameter()));
+
+          const Standard_Real a2DTol = Min( S1->UResolution(TolTang),
+                                              S2->UResolution(TolTang));
+
+          Standard_Boolean isReversed = ((aU2l - aU2f) < (aU1l - aU1f));
+
+          if(isReversed)
+          {
+            isDONE = IntCyCyTrim(quad2, quad1, TolTang, a2DTol, aBox2, aBox1,
+                                                  isReversed, empt, slin, spnt);
+          }
+          else
+          {
+            isDONE = IntCyCyTrim(quad1, quad2, TolTang, a2DTol, aBox1, aBox2,
+                                                  isReversed, empt, slin, spnt);
+          }
+
+          if(!isDONE)
+          {
+            isDONE = IntCyCy(quad1, quad2, TolTang, empt,
+                            SameSurf, multpoint, slin, spnt);
+            isTrimmed = Standard_False;
+          }
+        }
+
+        if (!isDONE)
+        {
+          return;
+        }
+
+        bEmpty = empt;
+        break;
       }
-      bEmpty = empt;
-      break;
-    }
     //
     case 23:
     case 32: { // Cylinder/Cone
@@ -238,136 +288,140 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
     return;
   }
   //
-  if (!SameSurf) {
-    AFunc.SetQuadric(quad2);
-    AFunc.Set(S1);
+
+  if(!isTrimmed)
+  {
+    if (!SameSurf) {
+      AFunc.SetQuadric(quad2);
+      AFunc.Set(S1);
     
-    solrst.Perform(AFunc, D1, TolArc, TolTang);
-    if (!solrst.IsDone()) {
-      return;
+      solrst.Perform(AFunc, D1, TolArc, TolTang);
+      if (!solrst.IsDone()) {
+        return;
+      }
+
+      if (solrst.AllArcSolution() && typs1 == typs2) {
+        all1 = Standard_True;
+      }
+      nbpt = solrst.NbPoints();
+      nbseg= solrst.NbSegments();
+      for (i=1; i<= nbpt; i++) {
+        pnt1.Append(solrst.Point(i));
+      }
+      for (i=1; i<= nbseg; i++) {
+        edg1.Append(solrst.Segment(i));
+      }
+      nosolonS1 = (nbpt == 0) && (nbseg == 0);
+
+      if (nosolonS1 && all1) {  // cas de face sans restrictions
+        all1 = Standard_False;
+      }
+    }//if (!SameSurf) {
+    else {
+      nosolonS1 = Standard_True;
     }
 
-    if (solrst.AllArcSolution() && typs1 == typs2) {
-      all1 = Standard_True;
-    }
-    nbpt = solrst.NbPoints();
-    nbseg= solrst.NbSegments();
-    for (i=1; i<= nbpt; i++) {
-      pnt1.Append(solrst.Point(i));
-    }
-    for (i=1; i<= nbseg; i++) {
-      edg1.Append(solrst.Segment(i));
-    }
-    nosolonS1 = (nbpt == 0) && (nbseg == 0);
+    if (!SameSurf) {
+      AFunc.SetQuadric(quad1);
+      AFunc.Set(S2);
 
-    if (nosolonS1 && all1) {  // cas de face sans restrictions
-      all1 = Standard_False;
-    }
-  }//if (!SameSurf) {
-  else {
-    nosolonS1 = Standard_True;
-  }
-
-  if (!SameSurf) {
-    AFunc.SetQuadric(quad1);
-    AFunc.Set(S2);
-
-    solrst.Perform(AFunc, D2, TolArc, TolTang);
-    if (!solrst.IsDone()) {
-      return;
-    }
+      solrst.Perform(AFunc, D2, TolArc, TolTang);
+      if (!solrst.IsDone()) {
+        return;
+      }
     
-    if (solrst.AllArcSolution() && typs1 == typs2) {
-      all2 = Standard_True;
-    }
-    nbpt = solrst.NbPoints();
-    nbseg= solrst.NbSegments();
-    for (i=1; i<= nbpt; i++) {
-      pnt2.Append(solrst.Point(i));
-    }
+      if (solrst.AllArcSolution() && typs1 == typs2) {
+        all2 = Standard_True;
+      }
+      nbpt = solrst.NbPoints();
+      nbseg= solrst.NbSegments();
+      for (i=1; i<= nbpt; i++) {
+        pnt2.Append(solrst.Point(i));
+      }
     
-    for (i=1; i<= nbseg; i++) {
-      edg2.Append(solrst.Segment(i));
-    }
-    nosolonS2 = (nbpt == 0) && (nbseg == 0);
+      for (i=1; i<= nbseg; i++) {
+        edg2.Append(solrst.Segment(i));
+      }
+      nosolonS2 = (nbpt == 0) && (nbseg == 0);
 
-    if (nosolonS2 && all2) {  // cas de face sans restrictions
-      all2 = Standard_False;
-    }
-  }// if (!SameSurf) {
-  else {
-    nosolonS2 = Standard_True;
-  }
-  //
-  if (SameSurf || (all1 && all2)) {
-    // faces "paralleles" parfaites
-    empt = Standard_False;
-    tgte = Standard_True;
-    slin.Clear();
-    spnt.Clear();
-
-    gp_Pnt Ptreference;
-
-    switch (typs1) {
-    case GeomAbs_Plane:      {
-      Ptreference = (S1->Plane()).Location();
-    }
-      break;
-    case GeomAbs_Cylinder: {
-      Ptreference = ElSLib::Value(0.,0.,S1->Cylinder());
-    }
-      break;
-    case GeomAbs_Sphere:      {
-      Ptreference = ElSLib::Value(M_PI/4.,M_PI/4.,S1->Sphere());
-    }
-      break;
-    case GeomAbs_Cone:      {
-      Ptreference = ElSLib::Value(0.,10.,S1->Cone());
-    }
-      break;
-    case GeomAbs_Torus:      {
-      Ptreference = ElSLib::Value(0.,0.,S1->Torus());
-    }
-      break;
-    default: 
-      break;
+      if (nosolonS2 && all2) {  // cas de face sans restrictions
+        all2 = Standard_False;
+      }
+    }// if (!SameSurf) {
+    else {
+      nosolonS2 = Standard_True;
     }
     //
-    oppo = quad1.Normale(Ptreference).Dot(quad2.Normale(Ptreference)) < 0.0;
-    done = Standard_True;
-    return;
-  }// if (SameSurf || (all1 && all2)) {
+    if (SameSurf || (all1 && all2)) {
+      // faces "paralleles" parfaites
+      empt = Standard_False;
+      tgte = Standard_True;
+      slin.Clear();
+      spnt.Clear();
 
-  if (!nosolonS1 || !nosolonS2) {
-    empt = Standard_False;
-    // C est la qu il faut commencer a bosser...
-    PutPointsOnLine(S1,S2,pnt1, slin, Standard_True, D1, quad1,quad2,
-                    multpoint,TolArc);
+      gp_Pnt Ptreference;
+
+      switch (typs1) {
+      case GeomAbs_Plane:      {
+        Ptreference = (S1->Plane()).Location();
+      }
+        break;
+      case GeomAbs_Cylinder: {
+        Ptreference = ElSLib::Value(0.,0.,S1->Cylinder());
+      }
+        break;
+      case GeomAbs_Sphere:      {
+        Ptreference = ElSLib::Value(M_PI/4.,M_PI/4.,S1->Sphere());
+      }
+        break;
+      case GeomAbs_Cone:      {
+        Ptreference = ElSLib::Value(0.,10.,S1->Cone());
+      }
+        break;
+      case GeomAbs_Torus:      {
+        Ptreference = ElSLib::Value(0.,0.,S1->Torus());
+      }
+        break;
+      default: 
+        break;
+      }
+      //
+      oppo = quad1.Normale(Ptreference).Dot(quad2.Normale(Ptreference)) < 0.0;
+      done = Standard_True;
+      return;
+    }// if (SameSurf || (all1 && all2)) {
+
+    if (!nosolonS1 || !nosolonS2) {
+      empt = Standard_False;
+      // C est la qu il faut commencer a bosser...
+      PutPointsOnLine(S1,S2,pnt1, slin, Standard_True, D1, quad1,quad2,
+                      multpoint,TolArc);
     
-    PutPointsOnLine(S1,S2,pnt2, slin, Standard_False,D2, quad2,quad1,
-                    multpoint,TolArc);
+      PutPointsOnLine(S1,S2,pnt2, slin, Standard_False,D2, quad2,quad1,
+                      multpoint,TolArc);
 
-    if (edg1.Length() != 0) {
-      ProcessSegments(edg1,slin,quad1,quad2,Standard_True,TolArc);
-    }
+      if (edg1.Length() != 0) {
+        ProcessSegments(edg1,slin,quad1,quad2,Standard_True,TolArc);
+      }
 
-    if (edg2.Length() != 0) {
-      ProcessSegments(edg2,slin,quad1,quad2,Standard_False,TolArc);
-    }
+      if (edg2.Length() != 0) {
+        ProcessSegments(edg2,slin,quad1,quad2,Standard_False,TolArc);
+      }
 
-    if (edg1.Length() !=0 || edg2.Length() !=0) {
-      //      ProcessRLine(slin,S1,S2,TolArc);
-      ProcessRLine(slin,quad1,quad2,TolArc);
+      if (edg1.Length() !=0 || edg2.Length() !=0) {
+        //      ProcessRLine(slin,S1,S2,TolArc);
+        ProcessRLine(slin,quad1,quad2,TolArc);
+      }
+    }//if (!nosolonS1 || !nosolonS2) {
+    else {
+      empt = ((slin.Length()==0) && (spnt.Length()==0));
     }
-  }//if (!nosolonS1 || !nosolonS2) {
-  else {
-    empt = ((slin.Length()==0) && (spnt.Length()==0));
   }
-  //
-  Standard_Integer nblin, aNbPnt;
+  
+  Standard_Integer  nblin = slin.Length(),
+                    aNbPnt = spnt.Length();
   //
   //modified by NIZNHY-PKV Tue Sep 06 10:03:35 2011f
-  aNbPnt=spnt.Length();
   if (aNbPnt) {
     IntPatch_SequenceOfPoint aSIP;
     //
@@ -401,7 +455,6 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
   }//  if (aNbPnt) {
   //modified by NIZNHY-PKV Tue Sep 06 10:18:20 2011t
   //
-  nblin = slin.Length();
   for(i=1; i<=nblin; i++) {
     IntPatch_IType thetype = slin.Value(i)->ArcType();
     if(  (thetype ==  IntPatch_Ellipse)

src/IntPatch/IntPatch_Intersection.cdl

@@ -143,6 +143,16 @@ is
                     typs1, typs2: SurfaceType from GeomAbs)
 
       is private;
+      
+    GeomGeomPerfomTrimSurf( me: in out;
+                            S1: HSurface from Adaptor3d; D1: TopolTool from Adaptor3d;
+                            S2: HSurface from Adaptor3d; D2: TopolTool from Adaptor3d;
+                            TolArc,TolTang: Real from Standard;
+                            LOfPnts: in out ListOfPntOn2S from IntSurf;
+                            RestrictLine: Boolean from Standard;
+                            typs1, typs2: SurfaceType from GeomAbs)
+
+      is private;
                     
     GeomParamPerfom(me: in out;
                     S1: HSurface from Adaptor3d; D1: TopolTool from Adaptor3d;

src/IntPatch/IntPatch_Intersection.cxx

@@ -1137,7 +1137,16 @@ void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)&  theS1,
   }
   else if(ts1 == 1)
   {
-    GeomGeomPerfom(theS1, theD1, theS2, theD2, TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2);
+    if(theD1->DomainIsInfinite() || theD2->DomainIsInfinite())
+    {
+      GeomGeomPerfom(theS1, theD1, theS2, theD2, TolArc, 
+                      TolTang, ListOfPnts, RestrictLine, typs1, typs2);
+    }
+    else
+    {
+      GeomGeomPerfomTrimSurf(theS1, theD1, theS2, theD2,
+              TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2);
+    }
   }
 }
 
@@ -1357,16 +1366,17 @@ void IntPatch_Intersection::GeomGeomPerfom(const Handle(Adaptor3d_HSurface)& the
 }
 
 //=======================================================================
-////function : GeomParamPerfom
+//function : GeomParamPerfom
 //purpose  : 
 //=======================================================================
-void IntPatch_Intersection::GeomParamPerfom(const Handle(Adaptor3d_HSurface)&  theS1,
-                                            const Handle(Adaptor3d_TopolTool)& theD1,
-                                            const Handle(Adaptor3d_HSurface)&  theS2,
-                                            const Handle(Adaptor3d_TopolTool)& theD2,
-                                            const Standard_Boolean isNotAnalitical,
-                                            const GeomAbs_SurfaceType typs1,
-                                            const GeomAbs_SurfaceType typs2)
+void IntPatch_Intersection::
+  GeomParamPerfom(const Handle(Adaptor3d_HSurface)&  theS1,
+                  const Handle(Adaptor3d_TopolTool)& theD1,
+                  const Handle(Adaptor3d_HSurface)&  theS2,
+                  const Handle(Adaptor3d_TopolTool)& theD2,
+                  const Standard_Boolean isNotAnalitical,
+                  const GeomAbs_SurfaceType typs1,
+                  const GeomAbs_SurfaceType typs2)
 {
   IntPatch_ImpPrmIntersection interip;
   if (myIsStartPnt)
@@ -1439,6 +1449,54 @@ void IntPatch_Intersection::GeomParamPerfom(const Handle(Adaptor3d_HSurface)&  t
   }
 }
 
+//=======================================================================
+//function : GeomGeomPerfomTrimSurf
+//purpose  : This function returns ready walking-line (which is not need
+//            in convertation) as an intersection line between two
+//            trimmed surfaces.
+//=======================================================================
+void IntPatch_Intersection::
+  GeomGeomPerfomTrimSurf( const Handle(Adaptor3d_HSurface)& theS1,
+                          const Handle(Adaptor3d_TopolTool)& theD1,
+                          const Handle(Adaptor3d_HSurface)& theS2,
+                          const Handle(Adaptor3d_TopolTool)& theD2,
+                          const Standard_Real theTolArc,
+                          const Standard_Real theTolTang,
+                          IntSurf_ListOfPntOn2S& theListOfPnts,
+                          const Standard_Boolean RestrictLine,
+                          const GeomAbs_SurfaceType theTyps1,
+                          const GeomAbs_SurfaceType theTyps2)
+{
+  IntSurf_Quadric Quad1,Quad2;
+
+  if((theTyps1 == GeomAbs_Cylinder) && (theTyps2 == GeomAbs_Cylinder))
+  {
+    IntPatch_ImpImpIntersection anInt;
+    anInt.Perform(theS1, theD1, theS2, theD2, theTolArc, theTolTang, Standard_True);
+
+    done = anInt.IsDone();
+
+    const Standard_Integer aNbLin = anInt.NbLines();
+    const Standard_Integer aNbPts = anInt.NbPnts();
+
+    for(Standard_Integer aLID = 1; aLID <= aNbLin; aLID++)
+    {
+      const Handle(IntPatch_Line)& aLine = anInt.Line(aLID);
+      slin.Append(aLine);
+    }
+
+    for(Standard_Integer aPID = 1; aPID <= aNbPts; aPID++)
+    {
+      const IntPatch_Point& aPoint = anInt.Point(aPID);
+      spnt.Append(aPoint);
+    }
+  }
+  else
+  {
+    GeomGeomPerfom(theS1, theD1, theS2, theD2, theTolArc, theTolTang, theListOfPnts, RestrictLine, theTyps1, theTyps2);
+  }
+}
+
 
 void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
                                     const Handle(Adaptor3d_TopolTool)& D1,