X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FINTERP_KERNEL%2FGeometric2D%2FInterpKernelGeo2DEdgeArcCircle.cxx;h=f73b9e154295082f7cbe64569e2786a5132d962e;hb=9727e779d56acece98be02cdccd0f99cc5ef0fa2;hp=ab736a68ba3bb5bd4d6a3d314da19f75536068e9;hpb=1123dccd6613b2e8abba35182759d5c4a11ecc8d;p=tools%2Fmedcoupling.git diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx index ab736a68b..f73b9e154 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2014 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; @@ -38,6 +39,14 @@ bool ArcCArcCIntersector::haveTheySameDirection() const return (getE1().getAngle()>0. && getE2().getAngle()>0.) || (getE1().getAngle()<0. && getE2().getAngle()<0.); } +bool ArcCArcCIntersector::areColinears() const +{ + double radiusL,radiusB; + double centerL[2],centerB[2]; + double tmp,cst; + return internalAreColinears(getE1(),getE2(),tmp,cst,radiusL,centerL,radiusB,centerB); +} + /*! * Precondition 'start' and 'end' are on the same curve than this. */ @@ -108,41 +117,44 @@ double ArcCArcCIntersector::getAngle(Node *node) const return EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(((*node)[0]-getE1().getCenter()[0])/getE1().getRadius(),((*node)[1]-getE1().getCenter()[1])/getE1().getRadius()); } -bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeArcCircle& a2) +bool ArcCArcCIntersector::internalAreColinears(const EdgeArcCircle& a1, const EdgeArcCircle& a2, double& distBetweenCenters, double& cst, + double& radiusL, double centerL[2], double& radiusB, double centerB[2]) { - double centerL[2],radiusL,angle0L,angleL; - double centerB[2],radiusB; double lgth1=fabs(a1.getAngle()*a1.getRadius()); double lgth2=fabs(a2.getAngle()*a2.getRadius()); if(lgth1getInterceptedArc(centerL,radiusL,angle0L,angleL); delete merge; // tmp=sqrt(tmp); - if(Node::areDoubleEqualsWP(tmp,0.,1/(10*std::max(radiusL,radiusB)))) - { - if(Node::areDoubleEquals(radiusL,radiusB)) - return true; - else - return false; - } + 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); double cmpContainer[4]; @@ -156,14 +168,14 @@ 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(); - if(_dist>radius1+radius2+QUADRATIC_PLANAR::_precision || _dist+std::min(radius1,radius2)+QUADRATIC_PLANAR::_precisionradius1+radius2+QuadraticPlanarPrecision::getPrecision() || _dist+std::min(radius1,radius2)+QuadraticPlanarPrecision::getPrecision() 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(); @@ -208,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())); } @@ -229,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())); } @@ -302,25 +321,61 @@ 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]; _drSq=_dx*_dx+_dy*_dy; _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*QUADRATIC_PLANAR::_precision)//QUADRATIC_PLANAR::_precision*QUADRATIC_PLANAR::_precision*_drSq*_drSq/(2.*_dx*_dx)) + ((*(_e2.getStartNode()))[0]-center[0])*((*(_e2.getEndNode()))[1]-center[1])- + ((*(_e2.getStartNode()))[1]-center[1])*((*(_e2.getEndNode()))[0]-center[0]); +} + +/** + 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; +} + +/*! + * By construction, no chance that an arc of circle and line to be colinear. + */ +bool ArcCSegIntersector::areColinears() const +{ + return false; } void ArcCSegIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const @@ -332,29 +387,71 @@ std::list< IntersectElement > ArcCSegIntersector::getIntersectionsCharacteristic { std::list< IntersectElement > ret; const double *center=getE1().getCenter(); - if(!(fabs(_determinant)<(2.*QUADRATIC_PLANAR::_precision)))//QUADRATIC_PLANAR::_precision*QUADRATIC_PLANAR::_precision*_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(); @@ -399,7 +496,7 @@ EdgeArcCircle::EdgeArcCircle(double sX, double sY, double mX, double mY, double * @param deltaAngle in ]-2.*Pi;2.*Pi[ */ EdgeArcCircle::EdgeArcCircle(Node *start, Node *end, const double *center, double radius, double angle0, double deltaAngle, bool direction):Edge(start,end,direction),_angle(deltaAngle), - _angle0(angle0),_radius(radius) + _angle0(angle0),_radius(radius) { _center[0]=center[0]; _center[1]=center[1]; @@ -447,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 { @@ -503,32 +600,11 @@ 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) { - //When arc is lower than 0.707 Using Asin - if(fabs(ux)<0.707) - { - double ret=SafeAcos(ux); - if(uy>0.) - return ret; - ret=-ret; - return ret; - } - else - { - double ret=SafeAsin(uy); - if(ux>0.) - return ret; - if(ret>0.) - return M_PI-ret; - else - return -M_PI-ret; - } + return atan2(uy, ux); } void EdgeArcCircle::GetArcOfCirclePassingThru(const double *start, const double *middle, const double *end, @@ -543,7 +619,7 @@ void EdgeArcCircle::GetArcOfCirclePassingThru(const double *start, const double angleInRad0=GetAbsoluteAngleOfNormalizedVect((start[0]-center[0])/radius,(start[1]-center[1])/radius); double angleInRadM=GetAbsoluteAngleOfNormalizedVect((middle[0]-center[0])/radius,(middle[1]-center[1])/radius); angleInRad=GetAbsoluteAngleOfNormalizedVect(((start[0]-center[0])*(end[0]-center[0])+(start[1]-center[1])*(end[1]-center[1]))/(radius*radius), - ((start[0]-center[0])*(end[1]-center[1])-(start[1]-center[1])*(end[0]-center[0]))/(radius*radius)); + ((start[0]-center[0])*(end[1]-center[1])-(start[1]-center[1])*(end[0]-center[0]))/(radius*radius)); if(IsAngleNotIn(angleInRad0,angleInRad,angleInRadM)) angleInRad=angleInRad<0?2*M_PI+angleInRad:angleInRad-2*M_PI; } @@ -624,12 +700,61 @@ void EdgeArcCircle::getBarycenterOfZone(double *bary) const double tmp3=cos(angle1); double tmp4=cos(_angle0); bary[0]=_radius*x0*y0*(tmp4-tmp3)+_radius*_radius*(y0*(cos(2*_angle0)-cos(2*angle1))/4.+ - x0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.)) - +tmp2*(tmp1*tmp1*tmp1-tmp0*tmp0*tmp0)/3.; + x0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.)) + +tmp2*(tmp1*tmp1*tmp1-tmp0*tmp0*tmp0)/3.; bary[1]=y0*y0*_radius*(tmp4-tmp3)/2.+_radius*_radius*y0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.) - +tmp2*(tmp4-tmp3+(tmp3*tmp3*tmp3-tmp4*tmp4*tmp4)/3.)/2.; + +tmp2*(tmp4-tmp3+(tmp3*tmp3*tmp3-tmp4*tmp4*tmp4)/3.)/2.; +} + +/** + * Compute the "middle" of two points on the arc of circle. + * The order (p1,p2) or (p2,p1) doesn't matter. p1 and p2 have to be localized on the edge defined by this. + * \param[out] mid the point located half-way between p1 and p2 on the arc defined by this. + * \sa getMiddleOfPointsOriented() a generalisation working also when p1 and p2 are not on the arc. + */ +void EdgeArcCircle::getMiddleOfPoints(const double *p1, const double *p2, double *mid) const +{ + double dx1((p1[0]-_center[0])/_radius),dy1((p1[1]-_center[1])/_radius),dx2((p2[0]-_center[0])/_radius),dy2((p2[1]-_center[1])/_radius); + double angle1(GetAbsoluteAngleOfNormalizedVect(dx1,dy1)),angle2(GetAbsoluteAngleOfNormalizedVect(dx2,dy2)); + // + double myDelta1(angle1-_angle0),myDelta2(angle2-_angle0); + if(_angle>0.) + { myDelta1=myDelta1>-QuadraticPlanarPrecision::getPrecision()?myDelta1:myDelta1+2.*M_PI; myDelta2=myDelta2>-QuadraticPlanarPrecision::getPrecision()?myDelta2:myDelta2+2.*M_PI; } + else + { myDelta1=myDelta1_angle (i.e. by the orientation of the arc). + * This function is sensitive to the ordering of p1 and p2. + * \param[out] mid the point located half-way between p1 and p2 + * \sa getMiddleOfPoints() to be used when the order of p1 and p2 is not relevant. + */ +void EdgeArcCircle::getMiddleOfPointsOriented(const double *p1, const double *p2, double *mid) const +{ + double dx1((p1[0]-_center[0])/_radius),dy1((p1[1]-_center[1])/_radius),dx2((p2[0]-_center[0])/_radius),dy2((p2[1]-_center[1])/_radius); + double angle1(GetAbsoluteAngleOfNormalizedVect(dx1,dy1)),angle2(GetAbsoluteAngleOfNormalizedVect(dx2,dy2)); + + if (angle1 <= 0.0) + angle1 += 2.*M_PI; + if (angle2 <= 0.0) + angle2 += 2.*M_PI; + + double avg; + if((_angle>0. && angle1 <= angle2) || (_angle<=0. && angle1 >= angle2)) + avg = (angle1+angle2)/2.; + else + avg = (angle1+angle2)/2. - M_PI; + + mid[0]=_center[0]+_radius*cos(avg); + mid[1]=_center[1]+_radius*sin(avg); } + /*! * Characteristic value used is angle in ]_Pi;Pi[ from axe 0x. */ @@ -653,14 +778,14 @@ bool EdgeArcCircle::isLower(double val1, double val2) const double myDelta2=val2-_angle0; if(_angle>0.) { - myDelta1=myDelta1>-(_radius*QUADRATIC_PLANAR::_precision)?myDelta1:myDelta1+2.*M_PI;//in some cases val1 or val2 are so close to angle0 that myDelta is close to 0. but negative. - myDelta2=myDelta2>-(_radius*QUADRATIC_PLANAR::_precision)?myDelta2:myDelta2+2.*M_PI; + myDelta1=myDelta1>-(_radius*QuadraticPlanarPrecision::getPrecision())?myDelta1:myDelta1+2.*M_PI;//in some cases val1 or val2 are so close to angle0 that myDelta is close to 0. but negative. + myDelta2=myDelta2>-(_radius*QuadraticPlanarPrecision::getPrecision())?myDelta2:myDelta2+2.*M_PI; return myDelta1& 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); -}