From: nbv Date: Thu, 21 Jan 2016 06:40:46 +0000 (+0300) Subject: 0023914: Intersection algorithm produced too many intersection points X-Git-Url: http://git.dev.opencascade.org/gitweb/?a=commitdiff_plain;h=refs%2Fheads%2FCR23914_4;p=occt-copy.git 0023914: Intersection algorithm produced too many intersection points Method IntWalk_PWalking::CalculateStepData(...) has been implemented in order to compute the step in parametric space. Relation between the step of 3D intersection curve and local curvature is described in A Practical Algorithm for Surface/Surface Intersection. Ling Huang Xinxiong Zhu （School of Mechanical Engineering & Automation, Beijing University of Aeronautics & Astronautics, Beijing, 100083） --- diff --git a/src/IntPatch/IntPatch_PrmPrmIntersection.cxx b/src/IntPatch/IntPatch_PrmPrmIntersection.cxx index c68e8d5a86..34cacdce7c 100644 --- a/src/IntPatch/IntPatch_PrmPrmIntersection.cxx +++ b/src/IntPatch/IntPatch_PrmPrmIntersection.cxx @@ -2341,13 +2341,13 @@ void IntPatch_PrmPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur //Try to extend the intersection line to boundary, if it is possibly Standard_Boolean hasBeenAdded = PW.PutToBoundary(Surf1, Surf2); - const Standard_Integer aMinNbPoints = 40; - if(iPWNbPoints < aMinNbPoints) - { - hasBeenAdded = - PW.SeekAdditionalPoints(Surf1, Surf2, aMinNbPoints) || hasBeenAdded; - iPWNbPoints = PW.NbPoints(); - } + //const Standard_Integer aMinNbPoints = 40; + //if(iPWNbPoints < aMinNbPoints) + //{ + // hasBeenAdded = + // PW.SeekAdditionalPoints(Surf1, Surf2, aMinNbPoints) || hasBeenAdded; + // iPWNbPoints = PW.NbPoints(); + //} RejectLine = Standard_False; Point3dDebut = PW.Value(1).Value(); @@ -2581,11 +2581,11 @@ void IntPatch_PrmPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur { if(PW.NbPoints()>2) { - const Standard_Integer aMinNbPoints = 40; - if(PW.NbPoints() < aMinNbPoints) - { - PW.SeekAdditionalPoints(Surf1, Surf2, aMinNbPoints); - } + //const Standard_Integer aMinNbPoints = 40; + //if(PW.NbPoints() < aMinNbPoints) + //{ + // PW.SeekAdditionalPoints(Surf1, Surf2, aMinNbPoints); + //} //----------------------------------------------- //-- Verification a posteriori : diff --git a/src/IntWalk/IntWalk_PWalking.cxx b/src/IntWalk/IntWalk_PWalking.cxx index abbca82257..d1e280128e 100644 --- a/src/IntWalk/IntWalk_PWalking.cxx +++ b/src/IntWalk/IntWalk_PWalking.cxx @@ -819,6 +819,11 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep, Standard_Integer LevelOfIterWithoutAppend = -1; // + const Standard_Real aStepMax[4] = { 0.1*(u1max - u1min), + 0.1*(v1max - v1min), + 0.1*(u2max - u2min), + 0.1*(v2max - v2min)}; + const Standard_Real aTol[4] = { Epsilon(u1max - u1min), Epsilon(v1max - v1min), Epsilon(u2max - u2min), @@ -854,15 +859,37 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep, previousPoint.Parameters(Param(1),Param(2),Param(3),Param(4)); // //--ofv.begin - Standard_Real aIncKey, aEps, dP1, dP2, dP3, dP4; - // - dP1 = sensCheminement * pasuv[0] * previousd1.X() /f; - dP2 = sensCheminement * pasuv[1] * previousd1.Y() /f; - dP3 = sensCheminement * pasuv[2] * previousd2.X() /f; - dP4 = sensCheminement * pasuv[3] * previousd2.Y() /f; - // - aIncKey=5.*(Standard_Real)IncKey; - aEps=1.e-7; + const Standard_Real aEps = Precision::Confusion(); + Standard_Real dP1, dP2, dP3, dP4; + + if((Status == IntWalk_OK) && CalculateStepData(Param, aStepMax)) + { +#ifdef INTWALK_PWALKING_DEBUG + printf("+++++++ The steps have been recomputed by ADAPTIVE algorith. New value:\n" + " dU1 = %10.20f;\n dV1 = %10.20f;\n dU2 = %10.20f;\n dV2 = %10.20f. -------\n", pasuv[0], pasuv[1], pasuv[2], pasuv[3]); +#endif + + dP1 = dP2 = dP3 = dP4 = 1.0; + } + else + {// +#ifdef INTWALK_PWALKING_DEBUG + printf("+++++++ The steps have been kept as PREVIOUS:\n" + " dU1 = %10.20f;\n dV1 = %10.20f;\n dU2 = %10.20f;\n dV2 = %10.20f. -------\n", pasuv[0], pasuv[1], pasuv[2], pasuv[3]); +#endif + + dP1 = previousd1.X() /f; + dP2 = previousd1.Y() /f; + dP3 = previousd2.X() /f; + dP4 = previousd2.Y() /f; + } + + dP1 *= (sensCheminement * pasuv[0]); + dP2 *= (sensCheminement * pasuv[1]); + dP3 *= (sensCheminement * pasuv[2]); + dP4 *= (sensCheminement * pasuv[3]); + + const Standard_Real aIncKey = 5.*(Standard_Real)IncKey; if(ChoixIso == IntImp_UIsoparametricOnCaro1 && Abs(dP1) < aEps) { dP1 *= aIncKey; @@ -882,12 +909,48 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep, { dP4 *= aIncKey; } - //--ofv.end + // Param(1) += dP1; Param(2) += dP2; Param(3) += dP3; Param(4) += dP4; + + // + //if (Param(1) < UFirst1) + //{ + // Param(1) = UFirst1; + //} + //if (ULast1 < Param(1)) + //{ + // Param(1) = ULast1; + //} + //if (Param(2) < VFirst1) + //{ + // Param(2) = VFirst1; + //} + //if (VLast1 < Param(2)) + //{ + // Param(2) = ULast1; + //} + //// + //if (Param(3) < UFirst2) + //{ + // Param(3) = UFirst2; + //} + //if (ULast2 < Param(3)) + //{ + // Param(3) = ULast2; + //} + //if (Param(4) < VFirst2) + //{ + // Param(4) = VFirst2; + //} + //if (VLast2 < Param(4)) + //{ + // Param(4) = ULast2; + //} + //========================== SvParam[0]=Param(1); SvParam[1]=Param(2); @@ -1298,6 +1361,7 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep, if(RejectIndex >= RejectIndexMAX) { + Arrive = Standard_True; break; } @@ -1385,6 +1449,7 @@ void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep, if(RejectIndex >= RejectIndexMAX) { + Arrive = Standard_True; break; } @@ -3250,3 +3315,194 @@ TestArret(const Standard_Boolean DejaReparti, } } +//======================================================================= +//function : CalculateStepData +//purpose : Computes new step (in UV-coordinates) which depends on +// the local curvature of the intersection line. +// If the step has been computed successfuly the method +// returns TRUE. +//======================================================================= +Standard_Boolean IntWalk_PWalking::CalculateStepData(const TColStd_Array1OfReal& theParams, + const Standard_Real theArrMaxStep[]) +{ + //return Standard_False; + + const Standard_Real aCoefStepMin = 100.0; + //sqrt(z*(2-z))*(1-z)/(1-2*z*(2-z)), when z = 0.01 maximal relative error while + //approximation the curve by circle + const Standard_Real aCoef = 0.125; + //As same as in approximation (see ApproxInt_ImpPrmSvSurfaces::Compute(...)) + const Standard_Real aNullValue = Precision::Approximation()* + Precision::Approximation(); + + const Handle(Adaptor3d_HSurface) &aSurf1 = myIntersectionOn2S.Function().AuxillarSurface1(), + &aSurf2 = myIntersectionOn2S.Function().AuxillarSurface2(); + + if((aSurf1->UContinuity() < GeomAbs_C2) || (aSurf2->VContinuity() < GeomAbs_C2)) + { + return Standard_False; + } + + if((aSurf2->UContinuity() < GeomAbs_C2) || (aSurf2->VContinuity() < GeomAbs_C2)) + { + return Standard_False; + } + + //if(aSurf1->IsUPeriodic() || aSurf1->IsVPeriodic()) + // return Standard_False; + + //if(aSurf2->IsUPeriodic() || aSurf2->IsVPeriodic()) + // return Standard_False; + + gp_Pnt aPt; + gp_Vec aDU1, aDV1, aDUU1, aDUV1, aDVV1; + gp_Vec aDU2, aDV2, aDUU2, aDUV2, aDVV2; + + aSurf1->D2(theParams(1), theParams(2), aPt, aDU1, aDV1, aDUU1, aDVV1, aDUV1); + aSurf2->D2(theParams(3), theParams(4), aPt, aDU2, aDV2, aDUU2, aDVV2, aDUV2); + +#ifdef INTWALK_PWALKING_DEBUG + int aTestID = 0; + Standard_Real anExpectedSqRad = -1.0; + switch(aTestID) + { + case 1: + //Intersection between two spherical surfaces: O1(0.0, 0.0, 0.0), R1 = 3 + //and O2(5.0, 0.0, 0.0), R2 = 5.0. + //Considered point has coordinates: (0.9, 0.0, 0.3*sqrt(91.0)). + + aDU1.SetCoord(0.00000000000000000, 0.90000000000000002, 0.00000000000000000); + aDV1.SetCoord(-2.8618176042508372, 0.00000000000000000, 0.90000000000000002); + aDUU1.SetCoord(-0.90000000000000002, 0.00000000000000000, 0.00000000000000000); + aDUV1.SetCoord(0.00000000000000000, -2.8618176042508372, 0.00000000000000000); + aDVV1.SetCoord(-0.90000000000000002, 0.00000000000000000, -2.8618176042508372); + aDU2.SetCoord(0.00000000000000000, -4.0999999999999996, 0.00000000000000000); + aDV2.SetCoord(-2.8618176042508372, 0.00000000000000000, -4.0999999999999996); + aDUU2.SetCoord(4.0999999999999996, 0.00000000000000000, 0.00000000000000000); + aDUV2.SetCoord(0.00000000000000000, -2.8618176042508372, 0.00000000000000000); + aDVV2.SetCoord(4.0999999999999996, 0.00000000000000000, -2.8618176042508372); + anExpectedSqRad = 819.0/100.0; + break; + case 2: + //Intersection between spherical surfaces: O1(0.0, 0.0, 0.0), R1 = 10 + //and the plane 3*x+4*y+z=26. + //Considered point has coordinates: (-1.68, 5.76, 8.0). + + aDU1.SetCoord(-5.76, -1.68, 0.0); + aDV1.SetCoord(2.24, -7.68, 6.0); + aDUU1.SetCoord(1.68, -5.76, 0.0); + aDUV1.SetCoord(7.68, 2.24, 0.0); + aDVV1.SetCoord(1.68, -5.76, -8.0); + aDU2.SetCoord(1.0, 0.0, -3.0); + aDV2.SetCoord(0.0, 1.0, -4.0); + aDUU2.SetCoord(0.0, 0.0, 0.0); + aDUV2.SetCoord(0.0, 0.0, 0.0); + aDVV2.SetCoord(0.0, 0.0, 0.0); + anExpectedSqRad = 74.0; + break; + default: + aTestID = 0; + break; + } +#endif + + const gp_Vec aN1(aDU1.Crossed(aDV1)), aN2(aDU2.Crossed(aDV2)); + //Tangent vactor to the intersection curve + const gp_Vec aCTan(aN1.Crossed(aN2)); + const Standard_Real aSqMagnFDer = aCTan.SquareMagnitude(); + + if(aSqMagnFDer < aNullValue) + return Standard_False; + + Standard_Real aDuS1 = 0.0, aDvS1 = 0.0, aDuS2 = 0.0, aDvS2 = 1.0; + + { + //This algorithm is described in NonSingularProcessing() function + //in ApproxInt_ImpPrmSvSurfaces.gxx file + Standard_Real aSqNMagn = aN1.SquareMagnitude(); + gp_Vec aTgU(aCTan.Crossed(aDU1)), aTgV(aCTan.Crossed(aDV1)); + Standard_Real aDeltaU = aTgV.SquareMagnitude()/aSqNMagn; + Standard_Real aDeltaV = aTgU.SquareMagnitude()/aSqNMagn; + + aDuS1 = Sign(sqrt(aDeltaU), aTgV.Dot(aN1)); + aDvS1 = -Sign(sqrt(aDeltaV), aTgU.Dot(aN1)); + + aSqNMagn = aN2.SquareMagnitude(); + aTgU.SetXYZ(aCTan.Crossed(aDU2).XYZ()); + aTgV.SetXYZ(aCTan.Crossed(aDV2).XYZ()); + aDeltaU = aTgV.SquareMagnitude()/aSqNMagn; + aDeltaV = aTgU.SquareMagnitude()/aSqNMagn; + + aDuS2 = Sign(sqrt(aDeltaU), aTgV.Dot(aN2)); + aDvS2 = -Sign(sqrt(aDeltaV), aTgU.Dot(aN2)); + } + + //According to "Marching along surface/surface intersection curves with an adaptive step length" + //by Tz.E.Stoyagov (http://www.sciencedirect.com/science/article/pii/016783969290046R) + //we obtain the system: + // {A*a+B*b=F1 + // {B*a+C*b=F2 + //where a and b should be found. + //After that, 2nd intersection curve derivative can be computed as + // r''(t)=a*N1+b*N2. + + const Standard_Real aA = aN1.Dot(aN1), aB = aN1.Dot(aN2), aC = aN2.Dot(aN2); + const Standard_Real aDetSyst = aB*aB - aA*aC; + + if(Abs(aDetSyst) < aNullValue) + {//Indetermined system solution + return Standard_False; + } + + const Standard_Real aF1 = aDuS1*aDuS1*aDUU1.Dot(aN1) + 2.0*aDuS1*aDvS1*aDUV1.Dot(aN1) + aDvS1*aDvS1*aDVV1.Dot(aN1); + const Standard_Real aF2 = aDuS2*aDuS2*aDUU2.Dot(aN2) + 2.0*aDuS2*aDvS2*aDUV2.Dot(aN2) + aDvS2*aDvS2*aDVV2.Dot(aN2); + + //Principal normal to the intersection curve + const gp_Vec aCNorm((aF1*aC-aF2*aB)/aDetSyst*aN1 + (aA*aF2-aF1*aB)/aDetSyst*aN2); + const Standard_Real aSqMagnSDer = aCNorm.CrossSquareMagnitude(aCTan); + + if(aSqMagnSDer < aNullValue) + {//Intersection curve has null curvature in observed point + return Standard_False; + } + +#ifdef INTWALK_PWALKING_DEBUG + if(aTestID) + { + //square of curvature radius + const Standard_Real aFactSqRad = aSqMagnFDer*aSqMagnFDer*aSqMagnFDer/aSqMagnSDer; + + if(Abs(aFactSqRad - anExpectedSqRad) < Precision::Confusion()) + { + printf("OK: Curvature radius is equal to expected (%5.10g)", anExpectedSqRad); + } + else + { + printf("Error: Curvature radius is not equal to expected: %5.10g != %5.10g", aFactSqRad, anExpectedSqRad); + } + + return Standard_False; + } +#endif + + + //Current step-system ([aDuS1 aDvS1 aDuS2 aDvS2]) provides the length of the + //intersection curve tangent to be equal to Lc=sqrt(aSqMagnFDer). + //We need in tangent length to be equal to Ln=aCoef*sqrt((aSqMagnFDer^3)/aSqMagnSDer). + //For that, new step-system step-system should be (Ln/Lc)*[aDuS1 aDvS1 aDuS2 aDvS2]. + + //(Ln/Lc) + const Standard_Real aRatio = aCoef*aSqMagnFDer/sqrt(aSqMagnSDer); + + pasuv[0] = Min(Abs(aDuS1*aRatio), theArrMaxStep[0]); + pasuv[1] = Min(Abs(aDvS1*aRatio), theArrMaxStep[1]); + pasuv[2] = Min(Abs(aDuS2*aRatio), theArrMaxStep[2]); + pasuv[3] = Min(Abs(aDvS2*aRatio), theArrMaxStep[3]); + + pasuv[0] = Max(pasuv[0], aCoefStepMin*ResoU1); + pasuv[1] = Max(pasuv[1], aCoefStepMin*ResoV1); + pasuv[2] = Max(pasuv[2], aCoefStepMin*ResoU2); + pasuv[3] = Max(pasuv[3], aCoefStepMin*ResoV2); + + return Standard_True; +} \ No newline at end of file diff --git a/src/IntWalk/IntWalk_PWalking.hxx b/src/IntWalk/IntWalk_PWalking.hxx index 35b4aea5e0..91b0f3180f 100644 --- a/src/IntWalk/IntWalk_PWalking.hxx +++ b/src/IntWalk/IntWalk_PWalking.hxx @@ -142,6 +142,7 @@ public: protected: + Standard_EXPORT Standard_Boolean CalculateStepData(const TColStd_Array1OfReal& theParams, const Standard_Real theArrMaxStep[]); diff --git a/src/IntWalk/IntWalk_PWalking.lxx b/src/IntWalk/IntWalk_PWalking.lxx index e1614170c9..3866bd3051 100644 --- a/src/IntWalk/IntWalk_PWalking.lxx +++ b/src/IntWalk/IntWalk_PWalking.lxx @@ -62,22 +62,14 @@ inline const gp_Dir& IntWalk_PWalking::TangentAtLine return tgdir; } -#define REGLAGE 0 - inline void IntWalk_PWalking::AddAPoint(Handle(IntSurf_LineOn2S)& theLine, const IntSurf_PntOn2S& POn2S) { -#if REGLAGE - Standard_Integer n=theLine->NbPoints(); - if(n) { - gp_Vec V(POn2S.Value(),theLine->Value(n).Value()); +#ifdef INTWALK_PWALKING_DEBUG + const gp_Pnt& aP = POn2S.Value(); Standard_Real u1,v1,u2,v2; - Standard_Real U1,V1,U2,V2; POn2S.Parameters(u1,v1,u2,v2); - theLine->Value(n).Parameters(U1,V1,U2,V2); - printf("\n%3d: (%10.5g)(%+12.5g %+12.5g %+12.5g) (%+12.5g %+12.5g) (%+12.5g %+12.5g)",n, - V.Magnitude(),V.X(),V.Y(),V.Z(),U1-u1,V1-v1,U2-u2,V2-v2); - fflush(stdout); - } + printf("Added point: (%+10.20f %+10.20f %+10.20f) (%+10.20f %+10.20f) " + "(%+10.20f %+10.20f)\n", aP.X(), aP.Y(), aP.Z(), u1, v1, u2, v2); #endif theLine->Add(POn2S); }