0026927: Make Intersection algorithm more adaptive

Proof of concept solution.

diff --git a/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx b/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx
index 8673a17984daa69d577ef7b73e9367012902fb6c..3752099d92cb64d08af06df3213910c711e69312 100644 (file)
--- a/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx
+++ b/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx
@@ -57,7 +57,8 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
                                   const Standard_Real theU2l,
                                   const Standard_Real theTol2D,
                                   const Standard_Real thePeriodOfSurf2,
-                                  const Standard_Boolean isTheReverse);
+                                  const Standard_Boolean isTheReverse,
+                                  const Standard_Boolean isUseCurvatureCheck = Standard_False);
 
 //=======================================================================
 //function : MinMax
@@ -1969,6 +1970,125 @@ static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1,
   return Standard_True;
 }
 
+#include <NCollection_LocalArray.hxx>
+#include <PLib.hxx>
+//=======================================================================
+//function : EvalCurv
+//purpose  : Evaluate curvature in dim-dimension point.
+//=======================================================================
+static Standard_Real EvalCurv(const Standard_Real dim, 
+                              const Standard_Real* V1,
+                              const Standard_Real* V2)
+{
+  // Really V1 and V2 are arrays of size dim;
+  // V1 is first derivative, V2 is second derivative
+  // of n-dimension curve
+  // Curvature is curv = |V1^V2|/|V1|^3
+  // V1^V2 is outer product of two vectors:
+  // P(i,j) = V1(i)*V2(j) - V1(j)*V2(i);
+  Standard_Real mp = 0.;
+  Standard_Integer i, j;
+  Standard_Real p;
+  for(i = 1; i < dim; ++i)
+  {
+    for(j = 0; j < i; ++j)
+    {
+      p = V1[i]*V2[j] - V1[j]*V2[i];
+      mp += p*p;
+    }
+  }
+  //mp *= 2.; //P(j,i) = -P(i,j);
+  mp = Sqrt(mp);
+  //
+  Standard_Real q = 0.;
+  for(i = 0; i < dim; ++i)
+  {
+    q += V1[i]*V1[i];
+  }
+  q = Sqrt(q);
+  //
+  Standard_Real curv = mp / (q*q*q);
+
+  return curv;
+}
+
+//=======================================================================
+//function : computeCurvatureArray
+//purpose  : Evaluate curvature in dim-dimension array.
+//=======================================================================
+static void computeCurvatureArray(const NCollection_LocalArray<Standard_Real>& theCoords,
+                                  const NCollection_Array1<Standard_Real>& thePars,
+                                  const Standard_Integer theDim, 
+                                  NCollection_Array1<Standard_Real> &aCurv)
+{
+  Standard_Real Val[21], Par[3], Res[21];
+  Standard_Integer i, j, k, l, m, ic;
+  Standard_Real aMaxCurv = 0.;
+  Standard_Integer dim = theDim;
+  //
+  i = aCurv.Lower();
+  for(j = 0; j < 3; ++j)
+  {
+    k = i+j;
+    ic = (k - aCurv.Lower()) * dim;
+    l = dim*j;
+    for(m = 0; m < dim; ++m)
+    {
+      Val[l + m] = theCoords[ic + m];
+    }
+    Par[j] = thePars(k);
+  }
+  PLib::EvalLagrange(Par[0], 2, 2, dim, *Val, *Par, *Res);
+  //
+  aCurv(i) = EvalCurv(dim, &Res[dim], &Res[2*dim]);
+  //
+  if(aCurv(i) > aMaxCurv)
+  {
+    aMaxCurv = aCurv(i);
+  }
+  //
+  for(i = aCurv.Lower()+1; i < aCurv.Upper(); ++i)
+  {
+    for(j = 0; j < 3; ++j)
+    {
+      k = i+j-1;
+      ic = (k - aCurv.Lower()) * dim;
+      l = dim*j;
+      for(m = 0; m < dim; ++m)
+      {
+        Val[l + m] = theCoords[ic + m];
+      }
+      Par[j] = thePars(k);
+    }
+    PLib::EvalLagrange(Par[1], 2, 2, dim, *Val, *Par, *Res);
+    //
+    aCurv(i) = EvalCurv(dim, &Res[dim], &Res[2*dim]);
+    if(aCurv(i) > aMaxCurv)
+    {
+      aMaxCurv = aCurv(i);
+    }
+  }
+  //
+  i = aCurv.Upper();
+  for(j = 0; j < 3; ++j)
+  {
+    k = i+j-2;
+    ic = (k - aCurv.Lower()) * dim;
+    l = dim*j;
+    for(m = 0; m < dim; ++m)
+    {
+      Val[l + m] = theCoords[ic + m];
+    }
+    Par[j] = thePars(k);
+  }
+  PLib::EvalLagrange(Par[2], 2, 2, dim, *Val, *Par, *Res);
+  //
+  aCurv(i) = EvalCurv(dim, &Res[dim], &Res[2*dim]);
+  if(aCurv(i) > aMaxCurv)
+  {
+    aMaxCurv = aCurv(i);
+  }
+}
 //=======================================================================
 //function : SeekAdditionalPoints
 //purpose  : 
@@ -1985,7 +2105,8 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
                                   const Standard_Real theU2l,
                                   const Standard_Real theTol2D,
                                   const Standard_Real thePeriodOfSurf2,
-                                  const Standard_Boolean isTheReverse)
+                                  const Standard_Boolean isTheReverse,
+                                  const Standard_Boolean isUseCurvatureCheck)
 {
   if(theLine.IsNull())
     return;
@@ -2018,8 +2139,10 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
   Standard_Integer aLastPointIndex = theEndPointOnLine;
   Standard_Real U1prec = 0.0, V1prec = 0.0, U2prec = 0.0, V2prec = 0.0;
 
+  // Build uniform distribution on parameter space.
   Standard_Integer aNbPointsPrev = 0;
-  while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev))
+  while(aNbPoints < theMinNbPoints && // Number of points check.
+        aNbPoints != aNbPointsPrev)   // Changes on iteration check.
   {
     aNbPointsPrev = aNbPoints;
     for(Standard_Integer fp = theStartPointOnLine, lp = 0; fp < aLastPointIndex; fp = lp + 1)
@@ -2066,7 +2189,7 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
       const gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
 
 #ifdef OCCT_DEBUG
-      //cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << endl;
+      cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << endl;
 #endif
 
       IntSurf_PntOn2S anIP;
@@ -2084,12 +2207,136 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
       aNbPoints++;
       aLastPointIndex++;
     }
+  }
 
-    if(aNbPoints >= theMinNbPoints)
+  if (isUseCurvatureCheck)
+  {
+    // Use curvature for computing additional points.
+    Standard_Integer aDim = 7;
+    aNbPoints = theLine->NbPoints();
+    NCollection_LocalArray<Standard_Real> aCoords(theLine->NbPoints()*aDim);
+    NCollection_Array1<Standard_Real> aParam(1, aNbPoints);
+    NCollection_Array1<Standard_Real> aCurv(1, theLine->NbPoints());
+    Standard_Integer i, j = 0;
+    Standard_Real u,v;
+    Standard_Real aFullDist = 0.0;
     {
-      return;
+      const IntSurf_PntOn2S &aPntOn2S    =  theLine->Value(1);
+
+      aCoords[j++] = aPntOn2S.Value().X();
+      aCoords[j++] = aPntOn2S.Value().Y();
+      aCoords[j++] = aPntOn2S.Value().Z();
+
+      aPntOn2S.ParametersOnS1(u, v);
+      aCoords[j++] = u;
+      aCoords[j++] = v;
+
+      aPntOn2S.ParametersOnS2(u, v);
+      aCoords[j++] = u;
+      aCoords[j++] = v;
+    }
+    for(i = 2; i <= theLine->NbPoints(); ++i)
+    {
+      j = (i - 1) * aDim;
+      const IntSurf_PntOn2S &aPntOn2S    =  theLine->Value(i);
+      const IntSurf_PntOn2S &aPntOn2Prev =  theLine->Value(i - 1);
+
+      aCoords[j++] = aPntOn2S.Value().X();
+      aCoords[j++] = aPntOn2S.Value().Y();
+      aCoords[j++] = aPntOn2S.Value().Z();
+
+      aPntOn2S.ParametersOnS1(u, v);
+      aCoords[j++] = u;
+      aCoords[j++] = v;
+
+      aPntOn2S.ParametersOnS2(u, v);
+      aCoords[j++] = u;
+      aCoords[j++] = v;
+
+      Standard_Real aCurrentDist = aPntOn2S.Value().Distance(aPntOn2Prev.Value());
+      aFullDist += aCurrentDist;
+      aParam(i) = aCurrentDist;
+
+    }
+
+    for(i = 1; i <= theLine->NbPoints(); ++i)
+      aParam(i) /= aFullDist;
+
+    computeCurvatureArray(aCoords, aParam, 7, aCurv);
+
+    Standard_Integer anOldMax = theLine->NbPoints();
+    Standard_Integer aShift = 0;
+    Standard_Integer aRealIdx = 0;
+    for(i = 2; i <= anOldMax; ++i)
+    {
+      const Standard_Real aMax = Max(aCurv(i), aCurv(i - 1));
+      const Standard_Real aMin = Min(aCurv(i), aCurv(i - 1));
+      const Standard_Real aDisbalance = aMax / (aMin + 1.0e-16);
+      if (aDisbalance > 5 &&
+          aDisbalance < 1000)
+      {
+        Standard_Real U1f = 0.0, V1f = 0.0; //first point in 1st suraface
+        Standard_Real U1l = 0.0, V1l = 0.0; //last  point in 1st suraface
+
+        Standard_Real U2f = 0.0, V2f = 0.0; //first point in 2nd suraface
+        Standard_Real U2l = 0.0, V2l = 0.0; //last  point in 2nd suraface
+
+        aRealIdx = i + aShift;
+
+        if(isTheReverse)
+        {
+          theLine->Value(aRealIdx - 1).ParametersOnS2(U1f, V1f);
+          theLine->Value(aRealIdx).ParametersOnS2(U1l, V1l);
+
+          theLine->Value(aRealIdx - 1).ParametersOnS1(U2f, V2f);
+          theLine->Value(aRealIdx).ParametersOnS1(U2l, V2l);
+        }
+        else
+        {
+          theLine->Value(aRealIdx - 1).ParametersOnS1(U1f, V1f);
+          theLine->Value(aRealIdx).ParametersOnS1(U1l, V1l);
+
+          theLine->Value(aRealIdx - 1).ParametersOnS2(U2f, V2f);
+          theLine->Value(aRealIdx).ParametersOnS2(U2l, V2l);
+        }
+        /*cout << aRealIdx - 1 << " " << aRealIdx << endl;
+        cout << theLine->Value(aRealIdx - 1).Value().X() << " "
+             << theLine->Value(aRealIdx - 1).Value().Y() << " "
+             << theLine->Value(aRealIdx - 1).Value().Z() << endl;
+        cout << theLine->Value(aRealIdx).Value().X() << " "
+             << theLine->Value(aRealIdx).Value().Y() << " "
+             << theLine->Value(aRealIdx).Value().Z() << endl;*/
+        Standard_Integer aNewPntNumber = Min(RealToInt(aDisbalance * 3), 25);
+        for(j = 1; j <= aNewPntNumber; ++j)
+        {
+          U1prec = U1f + j * (U1l - U1f) / (aNewPntNumber + 1);
+
+          if(!CylCylComputeParameters(U1prec, theWLIndex, theCoeffs, U2prec, V1prec, V2prec))
+            continue;
+
+          InscribePoint(theU2f, theU2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False);
+
+          const gp_Pnt aP1(theQuad1.Value(U1prec, V1prec)), aP2(theQuad2.Value(U2prec, V2prec));
+          const gp_Pnt aPInt(0.5 * (aP1.XYZ() + aP2.XYZ() ) );
+
+          IntSurf_PntOn2S anIP;
+          if(isTheReverse)
+          {
+            anIP.SetValue(aPInt, U2prec, V2prec, U1prec, V1prec);
+          }
+          else
+          {
+            anIP.SetValue(aPInt, U1prec, V1prec, U2prec, V2prec);
+          }
+
+          theLine->InsertBefore(aRealIdx, anIP);
+          aRealIdx++;
+          aShift++;
+        }
+      }
     }
   }
+
 }
 
 //=======================================================================
@@ -3050,7 +3297,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
             SeekAdditionalPoints( theQuad1, theQuad2, aWLine[i]->Curve(), 
                                   anEquationCoeffs, i, aNbPoints, 1,
                                   aWLine[i]->NbPnts(), aUSurf2f, aUSurf2l,
-                                  theTol2D, aPeriod, isTheReverse);
+                                  theTol2D, aPeriod, isTheReverse, Standard_True);
 
             aWLine[i]->ComputeVertexParameters(theTol3D);
             theSlin.Append(aWLine[i]);