X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FINTERP_KERNEL%2FGeometric2D%2FInterpKernelGeo2DEdgeArcCircle.cxx;h=f73b9e154295082f7cbe64569e2786a5132d962e;hb=9727e779d56acece98be02cdccd0f99cc5ef0fa2;hp=fc3fce7edfb48c74ad2d967435259ac4d7c261f0;hpb=62227e928392a71a014b2799d946401caf574dca;p=tools%2Fmedcoupling.git diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx index fc3fce7ed..f73b9e154 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2016 CEA/DEN, EDF R&D +// Copyright (C) 2007-2019 CEA/DEN, EDF R&D // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -26,6 +26,7 @@ #include #include +#include using namespace INTERP_KERNEL; @@ -131,11 +132,11 @@ bool ArcCArcCIntersector::internalAreColinears(const EdgeArcCircle& a1, const Ed a2.getCenter(centerL); radiusL=a2.getRadius(); a1.getCenter(centerB); radiusB=a1.getRadius(); } - // dividing from the begining by radiusB^2 to keep precision + // dividing from the beginning by radiusB^2 to keep precision 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 +153,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,10 +168,10 @@ 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) +void ArcCArcCIntersector::areOverlappedOrOnlyColinears(bool& obviousNoIntersection, bool& areOverlapped) { _dist=Node::distanceBtw2Pt(getE1().getCenter(),getE2().getCenter()); double radius1=getE1().getRadius(); double radius2=getE2().getRadius(); @@ -192,22 +193,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 +226,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 +248,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())); } @@ -313,13 +321,11 @@ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristi } return ret;*/ -ArcCSegIntersector::ArcCSegIntersector(const EdgeArcCircle& e1, const EdgeLin& e2, bool reverse):CrossTypeEdgeIntersector(e1,e2,reverse) -{ -} - -void ArcCSegIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped) +ArcCSegIntersector::ArcCSegIntersector(const EdgeArcCircle& e1, const EdgeLin& e2, bool reverse): + CrossTypeEdgeIntersector(e1,e2,reverse), + _deltaRoot_div_dr(0.), + _i1S2E(false),_i1E2E(false) { - areOverlapped=false;//No overlapping by construction const double *center=getE1().getCenter(); _dx=(*(_e2.getEndNode()))[0]-(*(_e2.getStartNode()))[0]; _dy=(*(_e2.getEndNode()))[1]-(*(_e2.getStartNode()))[1]; @@ -327,11 +333,41 @@ void ArcCSegIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, _cross= ((*(_e2.getStartNode()))[0]-center[0])*((*(_e2.getEndNode()))[1]-center[1])- ((*(_e2.getStartNode()))[1]-center[1])*((*(_e2.getEndNode()))[0]-center[0]); - _determinant=getE1().getRadius()*getE1().getRadius()/_drSq-_cross*_cross/(_drSq*_drSq); - if(_determinant>-2*QuadraticPlanarPrecision::getPrecision())//QuadraticPlanarPrecision::getPrecision()*QuadraticPlanarPrecision::getPrecision()*_drSq*_drSq/(2.*_dx*_dx)) +} + +/** + See http://mathworld.wolfram.com/Circle-LineIntersection.html + _cross is 'D', the computation is done with the translation to put back the circle at the origin +*/ +void ArcCSegIntersector::areOverlappedOrOnlyColinears(bool& obviousNoIntersection, bool& areOverlapped) +{ + areOverlapped=false;//No overlapping by construction + + // Similar optimisation than SegSegIntersector::areOverlappedOrOnlyColinears() + bool dnu1, dnu2; + identifyEarlyIntersection(dnu1, dnu2, _i1S2E, _i1E2E); + + const double R = getE1().getRadius(); + + // We need to compute d = R*R-_cross*_cross/_drSq + // In terms of numerical precision, this can trigger 'catastrophic cancellation' and is hence better expressed as: + double _dr = sqrt(_drSq); + double diff = (R-_cross/_dr), add=(R+_cross/_dr); + // Ah ah: we will be taking a square root later. If we want the user to be able to use an epsilon finer than 1.0e-8, then we need + // to prevent ourselves going below machine precision (typ. 1.0e-16 for double). + const double eps_machine = std::numeric_limits::epsilon(); + diff = fabs(diff/R) < eps_machine ? 0.0 : diff; + add = fabs(add/R) < eps_machine ? 0.0 : add; + double d = add*diff; + // Compute deltaRoot_div_dr := sqrt(delta)/dr, where delta has the meaning of Wolfram. + // Then 2*deltaRoot_div_dr is the distance between the two intersection points of the line with the circle. This is what we compare to eps. + // We compute it in such a way that it can be used in boolean tests too (a very negative value means we're far apart from intersection) + _deltaRoot_div_dr = Node::sign(d)*sqrt(fabs(d)); + + if( 2*_deltaRoot_div_dr > -QuadraticPlanarPrecision::getPrecision()) obviousNoIntersection=false; else - obviousNoIntersection=true; + obviousNoIntersection=true; } /*! @@ -351,29 +387,71 @@ std::list< IntersectElement > ArcCSegIntersector::getIntersectionsCharacteristic { std::list< IntersectElement > ret; const double *center=getE1().getCenter(); - if(!(fabs(_determinant)<(2.*QuadraticPlanarPrecision::getPrecision())))//QuadraticPlanarPrecision::getPrecision()*QuadraticPlanarPrecision::getPrecision()*_drSq*_drSq/(2.*_dx*_dx)) - { - double determinant=EdgeArcCircle::SafeSqrt(_determinant); + if(!(2*fabs(_deltaRoot_div_dr) < QuadraticPlanarPrecision::getPrecision())) // see comments in areOverlappedOrOnlyColinears() + { // Two intersection nodes + // -> if a common node found, there is a chance that this is the only one (i.e. second intersection point is outside e1 and e2) + if(_earlyInter) + { + // Check tangent vector of the arc circle at the common node with the linear segment. + // There we can tell if the arc of circle is 'moving away' from the seg, or if it might intersect it twice + const Node &n(*_earlyInter->getNodeOnly()); + const double * center(getE1().getCenter()); + + double tang[2] = {-(n[1]-center[1]), n[0]-center[0]}; // (-y, x) is the tangent vector in the trigo direction with (x,y) = (center->node) + bool invSeg = _i1S2E || _i1E2E; + double linEdge[2] = {invSeg ? (-_dx) : _dx, invSeg ? (-_dy) : _dy}; + if(tang[1]*linEdge[0]-tang[0]*linEdge[1] < 0) + { + ret.push_back(*_earlyInter); + return ret; + } + } + + double determinant=fabs(_deltaRoot_div_dr)/sqrt(_drSq); double x1=(_cross*_dy/_drSq+Node::sign(_dy)*_dx*determinant)+center[0]; double y1=(-_cross*_dx/_drSq+fabs(_dy)*determinant)+center[1]; Node *intersect1=new Node(x1,y1); intersect1->declareOn(); - bool i1_1S=_e1.getStartNode()->isEqual(*intersect1); - bool i1_1E=_e1.getEndNode()->isEqual(*intersect1); - bool i1_2S=_e2.getStartNode()->isEqual(*intersect1); - bool i1_2E=_e2.getEndNode()->isEqual(*intersect1); - ret.push_back(IntersectElement(getE1().getCharactValue(*intersect1),getE2().getCharactValue(*intersect1),i1_1S,i1_1E,i1_2S,i1_2E,intersect1,_e1,_e2,keepOrder())); - // double x2=(_cross*_dy/_drSq-Node::sign(_dy)*_dx*determinant)+center[0]; double y2=(-_cross*_dx/_drSq-fabs(_dy)*determinant)+center[1]; Node *intersect2=new Node(x2,y2); intersect2->declareOn(); - bool i2_1S=_e1.getStartNode()->isEqual(*intersect2); - bool i2_1E=_e1.getEndNode()->isEqual(*intersect2); - bool i2_2S=_e2.getStartNode()->isEqual(*intersect2); - bool i2_2E=_e2.getEndNode()->isEqual(*intersect2); - ret.push_back(IntersectElement(getE1().getCharactValue(*intersect2),getE2().getCharactValue(*intersect2),i2_1S,i2_1E,i2_2S,i2_2E,intersect2,_e1,_e2,keepOrder())); + + bool isN1(false), isN2(false); + if (_earlyInter) + { + // Which node do we actually already found? Assume this is the closest ... + const Node &iN = *(_earlyInter->getNodeOnly()); + const Node &n1(*intersect1), &n2(*intersect2); + double d1 = std::max(fabs(iN[0]-n1[0]), fabs(iN[1]-n1[1])); + double d2 = std::max(fabs(iN[0]-n2[0]), fabs(iN[1]-n2[1])); + isN1 = d1 < d2; isN2 = !isN1; + if (isN1) intersect1->decrRef(); + if (isN2) intersect2->decrRef(); + ret.push_back(*_earlyInter); + } + if (!isN1) + { + bool i1_1S=_e1.getStartNode()->isEqual(*intersect1); + bool i1_1E=_e1.getEndNode()->isEqual(*intersect1); + bool i1_2S=_e2.getStartNode()->isEqual(*intersect1); + bool i1_2E=_e2.getEndNode()->isEqual(*intersect1); + ret.push_back(IntersectElement(getE1().getCharactValue(*intersect1),getE2().getCharactValue(*intersect1),i1_1S,i1_1E,i1_2S,i1_2E,intersect1,_e1,_e2,keepOrder())); + } + if(!isN2) + { + bool i2_1S=_e1.getStartNode()->isEqual(*intersect2); + bool i2_1E=_e1.getEndNode()->isEqual(*intersect2); + bool i2_2S=_e2.getStartNode()->isEqual(*intersect2); + bool i2_2E=_e2.getEndNode()->isEqual(*intersect2); + ret.push_back(IntersectElement(getE1().getCharactValue(*intersect2),getE2().getCharactValue(*intersect2),i2_1S,i2_1E,i2_2S,i2_2E,intersect2,_e1,_e2,keepOrder())); + } } else//tangent intersection { + if (_earlyInter) + { + ret.push_back(*_earlyInter); + return ret; + } double x=(_cross*_dy)/_drSq+center[0]; double y=(-_cross*_dx)/_drSq+center[1]; Node *intersect3=new Node(x,y); intersect3->declareOnTangent(); @@ -466,7 +544,7 @@ void EdgeArcCircle::unApplySimilarity(double xBary, double yBary, double dimChar /*! * 'eps' is expected to be > 0. * 'conn' is of size 3. conn[0] is start id, conn[1] is end id and conn[2] is middle id. - * 'offset' is typically the number of nodes already existing in global 2D curve mesh. Additionnal coords 'addCoo' ids will be put after the already existing. + * 'offset' is typically the number of nodes already existing in global 2D curve mesh. Additional coords 'addCoo' ids will be put after the already existing. */ void EdgeArcCircle::tesselate(const int *conn, int offset, double eps, std::vector& newConn, std::vector& addCoo) const { @@ -522,10 +600,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) { @@ -814,28 +889,3 @@ void EdgeArcCircle::updateBounds() if(IsIn2Pi(_angle0,_angle,M_PI)) _bounds[0]=_center[0]-_radius; } - -void EdgeArcCircle::fillGlobalInfoAbs(bool direction, const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, - std::vector& edgesThis, std::vector& addCoo, std::map mapAddCoo) const -{ - int tmp[2]; - _start->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp); - _end->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp+1); - if(direction) - { - edgesThis.push_back(tmp[0]); - edgesThis.push_back(tmp[1]); - } - else - { - edgesThis.push_back(tmp[1]); - edgesThis.push_back(tmp[0]); - } -} - -void EdgeArcCircle::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 -{ - _start->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther); - _end->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther); -}