From c64f1fb7cab0e3e22c97825751fbdc6352c51660 Mon Sep 17 00:00:00 2001 From: aml Date: Wed, 25 Nov 2015 21:38:20 +0300 Subject: [PATCH] 0026927: Make Intersection algorithm more adaptive Proof of concept solution. --- .../IntPatch_ImpImpIntersection_4.gxx | 261 +++++++++++++++++- 1 file changed, 254 insertions(+), 7 deletions(-) diff --git a/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx b/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx index 8673a17984..3752099d92 100644 --- 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 +#include +//======================================================================= +//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& theCoords, + const NCollection_Array1& thePars, + const Standard_Integer theDim, + NCollection_Array1 &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 aCoords(theLine->NbPoints()*aDim); + NCollection_Array1 aParam(1, aNbPoints); + NCollection_Array1 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]); -- 2.39.5