From 893e88f2a8635d9081ff95c39d9fab23d8dff2a8 Mon Sep 17 00:00:00 2001 From: abn Date: Mon, 9 Apr 2018 17:12:30 +0200 Subject: [PATCH] [Backport from 8.5] Bug fix: points located on the arc where wrongly detected by arc/arc intersector. --- .../InterpKernelGeo2DEdgeArcCircle.cxx | 48 ++++++++++--------- .../InterpKernelGeo2DEdgeArcCircle.hxx | 8 ++-- .../Geometric2D/InterpKernelGeo2DNode.hxx | 4 +- .../QuadraticPlanarInterpTest.hxx | 3 ++ .../QuadraticPlanarInterpTest4.cxx | 33 +++++++++++++ 5 files changed, 70 insertions(+), 26 deletions(-) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx index fc3fce7ed..7f1788eac 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx @@ -135,7 +135,7 @@ bool ArcCArcCIntersector::internalAreColinears(const EdgeArcCircle& a1, const Ed distBetweenCenters=Node::distanceBtw2PtSq(centerL,centerB); cst=distBetweenCenters/(radiusB*radiusB); cst+=radiusL*radiusL/(radiusB*radiusB); - return Node::areDoubleEqualsWP(cst,1.,2.); + return Node::areDoubleEqualsWPRight(cst,1.,2.); } bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeArcCircle& a2) @@ -152,7 +152,7 @@ bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeA delete merge; // tmp=sqrt(tmp); - if(Node::areDoubleEqualsWP(tmp,0.,1/(10*std::max(radiusL,radiusB)))) + if(Node::areDoubleEqualsWPLeft(tmp,0.,10*std::max(radiusL,radiusB))) return Node::areDoubleEquals(radiusL,radiusB); double phi=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect((centerL[0]-centerB[0])/tmp,(centerL[1]-centerB[1])/tmp); double cst2=2*radiusL*tmp/(radiusB*radiusB); @@ -167,7 +167,7 @@ bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeA if(EdgeArcCircle::IsIn2Pi(angle0L,angleL,a)) cmpContainer[sizeOfCmpContainer++]=cst-cst2; a=*std::max_element(cmpContainer,cmpContainer+sizeOfCmpContainer); - return Node::areDoubleEqualsWP(a,1.,2.); + return Node::areDoubleEqualsWPRight(a,1.,2.); } void ArcCArcCIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped) @@ -192,22 +192,28 @@ void ArcCArcCIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind } } +/** + Heart of the algorithm for arc/arc intersection. + See http://mathworld.wolfram.com/Circle-CircleIntersection.html + The computation is done in the coordinate system where Ox is the line between the 2 circle centers. +*/ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristicVal() const { std::list< IntersectElement > ret; const double *center1=getE1().getCenter(); const double *center2=getE2().getCenter(); double radius1=getE1().getRadius(); double radius2=getE2().getRadius(); - double d1_1=(_dist*_dist-radius2*radius2+radius1*radius1)/(2.*_dist); + double d1_1=(_dist*_dist-radius2*radius2+radius1*radius1)/(2.*_dist); // computation of 'x' on wolfram double u[2];//u is normalized vector from center1 to center2. u[0]=(center2[0]-center1[0])/_dist; u[1]=(center2[1]-center1[1])/_dist; - double d1_1y=EdgeArcCircle::SafeSqrt(radius1*radius1-d1_1*d1_1); + double d1_1y=EdgeArcCircle::SafeSqrt(radius1*radius1-d1_1*d1_1); // computation of 'y' on wolfram double angleE1=EdgeArcCircle::NormalizeAngle(getE1().getAngle0()+getE1().getAngle()); double angleE2=EdgeArcCircle::NormalizeAngle(getE2().getAngle0()+getE2().getAngle()); if(!Node::areDoubleEquals(d1_1y,0)) { //2 intersections double v1[2],v2[2]; + // coming back to our coordinate system: v1[0]=u[0]*d1_1-u[1]*d1_1y; v1[1]=u[1]*d1_1+u[0]*d1_1y; v2[0]=u[0]*d1_1+u[1]*d1_1y; v2[1]=u[1]*d1_1-u[0]*d1_1y; Node *node1=new Node(center1[0]+v1[0],center1[1]+v1[1]); node1->declareOn(); @@ -219,16 +225,17 @@ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristi v4[0]=center1[0]-center2[0]+v2[0]; v4[1]=center1[1]-center2[1]+v2[1]; double angle1_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v3[0]/radius2,v3[1]/radius2); double angle2_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v4[0]/radius2,v4[1]/radius2); + // Check whether intersection points are exactly ON the other arc or not + // -> the curvilinear distance (=radius*angle) must below eps + bool e1_1S=Node::areDoubleEqualsWPLeft(angle1_1,getE1().getAngle0(),radius1); + bool e1_1E=Node::areDoubleEqualsWPLeft(angle1_1,angleE1,radius1); + bool e1_2S=Node::areDoubleEqualsWPLeft(angle1_2,getE2().getAngle0(),radius1); + bool e1_2E=Node::areDoubleEqualsWPLeft(angle1_2,angleE2,radius1); // - bool e1_1S=Node::areDoubleEqualsWP(angle1_1,getE1().getAngle0(),radius1); - bool e1_1E=Node::areDoubleEqualsWP(angle1_1,angleE1,radius1); - bool e1_2S=Node::areDoubleEqualsWP(angle1_2,getE2().getAngle0(),radius1); - bool e1_2E=Node::areDoubleEqualsWP(angle1_2,angleE2,radius1); - // - bool e2_1S=Node::areDoubleEqualsWP(angle2_1,getE1().getAngle0(),radius2); - bool e2_1E=Node::areDoubleEqualsWP(angle2_1,angleE1,radius2); - bool e2_2S=Node::areDoubleEqualsWP(angle2_2,getE2().getAngle0(),radius2); - bool e2_2E=Node::areDoubleEqualsWP(angle2_2,angleE2,radius2); + bool e2_1S=Node::areDoubleEqualsWPLeft(angle2_1,getE1().getAngle0(),radius2); + bool e2_1E=Node::areDoubleEqualsWPLeft(angle2_1,angleE1,radius2); + bool e2_2S=Node::areDoubleEqualsWPLeft(angle2_2,getE2().getAngle0(),radius2); + bool e2_2E=Node::areDoubleEqualsWPLeft(angle2_2,angleE2,radius2); ret.push_back(IntersectElement(angle1_1,angle1_2,e1_1S,e1_1E,e1_2S,e1_2E,node1,_e1,_e2,keepOrder())); ret.push_back(IntersectElement(angle2_1,angle2_2,e2_1S,e2_1E,e2_2S,e2_2E,node2,_e1,_e2,keepOrder())); } @@ -240,10 +247,10 @@ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristi v2[0]=center1[0]-center2[0]+v1[0]; v2[1]=center1[1]-center2[1]+v1[1]; double angle0_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1); double angle0_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v2[0]/radius2,v2[1]/radius2); - bool e0_1S=Node::areDoubleEqualsWP(angle0_1,getE1().getAngle0(),radius1); - bool e0_1E=Node::areDoubleEqualsWP(angle0_1,angleE1,radius1); - bool e0_2S=Node::areDoubleEqualsWP(angle0_2,getE2().getAngle0(),radius2); - bool e0_2E=Node::areDoubleEqualsWP(angle0_2,angleE2,radius2); + bool e0_1S=Node::areDoubleEqualsWPLeft(angle0_1,getE1().getAngle0(),radius1); + bool e0_1E=Node::areDoubleEqualsWPLeft(angle0_1,angleE1,radius1); + bool e0_2S=Node::areDoubleEqualsWPLeft(angle0_2,getE2().getAngle0(),radius2); + bool e0_2E=Node::areDoubleEqualsWPLeft(angle0_2,angleE2,radius2); Node *node=new Node(center1[0]+d1_1*u[0],center1[1]+d1_1*u[1]); node->declareOnTangent(); ret.push_back(IntersectElement(angle0_1,angle0_2,e0_1S,e0_1E,e0_2S,e0_2E,node,_e1,_e2,keepOrder())); } @@ -522,10 +529,7 @@ double EdgeArcCircle::GetAbsoluteAngle(const double *vect, double& normVect) /*! * Given a \b normalized vector defined by (ux,uy) returns its angle in ]-Pi;Pi]. - * So before using this method ux*ux+uy*uy should as much as possible close to 1. - * This methods is quite time consuming in order to keep as much as possible precision. - * It is NOT ALWAYS possible to do that only in one call of acos. Sometimes call to asin is necessary - * due to imperfection of acos near 0. and Pi (cos x ~ 1-x*x/2.) + * Actually in the current implementation, the vector does not even need to be normalized ... */ double EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(double ux, double uy) { diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx index 31f5fafc1..19559f9bc 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx @@ -126,11 +126,13 @@ namespace INTERP_KERNEL void fillGlobalInfoAbs2(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, std::vector& edgesOther, std::vector& addCoo, std::map& mapAddCoo) const; protected: - //!Value between -2Pi and 2Pi - double _angle; - //!Value between -Pi and Pi + //! Absolute angle where the arc starts. Value between -Pi and Pi double _angle0; + //! Angular span of the arc. Value between -2Pi and 2Pi + double _angle; + //! Radius of the arc of circle double _radius; + //! Center of the arc of circle double _center[2]; }; } diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.hxx index 970b0f328..cbe1e0471 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.hxx @@ -90,7 +90,9 @@ namespace INTERP_KERNEL static double norm(const double *vect) { return sqrt(vect[0]*vect[0]+vect[1]*vect[1]); } static bool areDoubleEquals(double a, double b) { return fabs(a-b) < QuadraticPlanarPrecision::getPrecision(); } //! idem areDoubleEquals except that precision of comparison is modified. - static bool areDoubleEqualsWP(double a, double b, double k) { return fabs(a-b) < k*QuadraticPlanarPrecision::getPrecision(); } + static bool areDoubleEqualsWPLeft(double a, double b, double k) { return k*fabs(a-b) < QuadraticPlanarPrecision::getPrecision(); } + //! idem areDoubleEquals except that precision of comparison is modified. + static bool areDoubleEqualsWPRight(double a, double b, double k) { return fabs(a-b) < k*QuadraticPlanarPrecision::getPrecision(); } static double distanceBtw2Pt(const double *a, const double *b) { return sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])); } static double distanceBtw2PtSq(const double *a, const double *b) { return (a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1]); } // diff --git a/src/INTERP_KERNELTest/QuadraticPlanarInterpTest.hxx b/src/INTERP_KERNELTest/QuadraticPlanarInterpTest.hxx index b5cb15293..87b0ba653 100644 --- a/src/INTERP_KERNELTest/QuadraticPlanarInterpTest.hxx +++ b/src/INTERP_KERNELTest/QuadraticPlanarInterpTest.hxx @@ -110,6 +110,8 @@ namespace INTERP_TEST CPPUNIT_TEST( checkIsInOrOut ); CPPUNIT_TEST( checkGetMiddleOfPoints ); CPPUNIT_TEST( checkGetMiddleOfPointsOriented ); + CPPUNIT_TEST( checkArcArcIntersection1 ); + CPPUNIT_TEST_SUITE_END(); public: void setUp(); @@ -202,6 +204,7 @@ namespace INTERP_TEST void checkIsInOrOut(); void checkGetMiddleOfPoints(); void checkGetMiddleOfPointsOriented(); + void checkArcArcIntersection1(); private: INTERP_KERNEL::QuadraticPolygon *buildQuadraticPolygonCoarseInfo(const double *coords, const int *conn, int lgth); diff --git a/src/INTERP_KERNELTest/QuadraticPlanarInterpTest4.cxx b/src/INTERP_KERNELTest/QuadraticPlanarInterpTest4.cxx index 019eeb84e..67f0a1104 100644 --- a/src/INTERP_KERNELTest/QuadraticPlanarInterpTest4.cxx +++ b/src/INTERP_KERNELTest/QuadraticPlanarInterpTest4.cxx @@ -1682,4 +1682,37 @@ void QuadraticPlanarInterpTest::checkMakePartitionAbs1() pol1.dumpInXfigFileWithOther(pol2,"tony.fig"); } +/** + * Arc/arc intersection was buggy: a point was detected OFF when used in a linear segment, but + * detected ON when used in an (almost flat) arc of circle. + */ +void QuadraticPlanarInterpTest::checkArcArcIntersection1() +{ + double eps=1.0e-8; + INTERP_KERNEL::QuadraticPlanarPrecision::setPrecision(eps); + INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision::setArcDetectionPrecision(eps); + + Node *n0=new Node(6.37533,38.8928); Node *n3=new Node(6.29194,39.2789); + Node *n1=new Node(6.13158,38.8308); Node *n4=new Node(6.31919,38.7607); + Node *n2=new Node(6.25346,38.8618); Node *n5=new Node(6.38778,39.0241); + + Node *n6=new Node(6.2534549999999998, 38.861800000000002); // to have a linear edge e1 + + //EdgeArcCircle *e1=new EdgeArcCircle(n0, n2, n6, true); // to have a linear edge e1 + EdgeArcCircle *e1=new EdgeArcCircle(n0, n2, n1, true); + EdgeArcCircle *e2=new EdgeArcCircle(n3, n5, n4, true); + + MergePoints merge; + QuadraticPolygon c1,c2; + e1->intersectWith(e2,merge,c1,c2); + CPPUNIT_ASSERT_EQUAL(2,c1.size()); + CPPUNIT_ASSERT_EQUAL(2,c2.size()); + //clean-up + n0->decrRef(); n1->decrRef(); n2->decrRef(); n3->decrRef(); n4->decrRef(); n5->decrRef(); n6->decrRef(); + e1->decrRef(); e2->decrRef(); } + + +} + + -- 2.39.2