src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx

   1 // Copyright (C) 2007-2021  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
  22 #include "InterpKernelGeo2DEdgeLin.hxx"
  23 #include "InterpKernelException.hxx"
  24 #include "InterpKernelGeo2DNode.hxx"
  25 #include "NormalizedUnstructuredMesh.hxx"
  26
  27 #include <sstream>
  28 #include <algorithm>
  29 #include <limits>
  30
  31 using namespace INTERP_KERNEL;
  32
  33 ArcCArcCIntersector::ArcCArcCIntersector(const EdgeArcCircle& e1, const EdgeArcCircle& e2):SameTypeEdgeIntersector(e1,e2),_dist(0.)
  34 {
  35 }
  36
  37 bool ArcCArcCIntersector::haveTheySameDirection() const
  38 {
  39   return (getE1().getAngle()>0. &&  getE2().getAngle()>0.) || (getE1().getAngle()<0. &&  getE2().getAngle()<0.);
  40 }
  41
  42 bool ArcCArcCIntersector::areColinears() const
  43 {
  44   double radiusL,radiusB;
  45   double centerL[2],centerB[2];
  46   double tmp,cst;
  47   return internalAreColinears(getE1(),getE2(),tmp,cst,radiusL,centerL,radiusB,centerB);
  48 }
  49
  50 /*!
  51  * Precondition 'start' and 'end' are on the same curve than this.
  52  */
  53 void ArcCArcCIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
  54 {
  55   bool obvious1,obvious2;
  56   obviousCaseForCurvAbscisse(start,whereStart,commonNode,obvious1);
  57   obviousCaseForCurvAbscisse(end,whereEnd,commonNode,obvious2);
  58   if(obvious1 && obvious2)
  59     return ;
  60   double angleInRadStart=getAngle(start);
  61   double angleInRadEnd=getAngle(end);
  62   if(obvious1 || obvious2)
  63     {
  64       if(obvious1)
  65         {
  66           if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
  67             whereEnd=INSIDE;
  68           else
  69             whereEnd=OUT_AFTER;
  70           return ;
  71         }
  72       else
  73         {
  74           if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart))
  75             whereStart=INSIDE;
  76           else
  77             whereStart=OUT_BEFORE;
  78           return ;
  79         }
  80     }
  81   if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart))
  82     {
  83       whereStart=INSIDE;
  84       if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
  85         whereEnd=INSIDE;
  86       else
  87         whereEnd=OUT_AFTER;
  88     }
  89   else
  90     {//we are out in start.
  91       if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
  92         {
  93           whereStart=OUT_BEFORE;
  94           whereEnd=INSIDE;
  95         }
  96       else
  97         {
  98           if(EdgeArcCircle::IsIn2Pi(getE2().getAngle0(),getE2().getAngle(),getE1().getAngle0()))
  99             {//_e2 contains stictly _e1
 100               whereStart=OUT_BEFORE;
 101               whereEnd=OUT_AFTER;
 102             }
 103           else
 104             {//_e2 is outside from _e1
 105               whereStart=OUT_BEFORE;
 106               whereEnd=OUT_BEFORE;
 107             }
 108         }
 109     }
 110 }
 111
 112 /*!
 113  * Return angle between ]-Pi;Pi[
 114  */
 115 double ArcCArcCIntersector::getAngle(Node *node) const
 116 {
 117   return EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(((*node)[0]-getE1().getCenter()[0])/getE1().getRadius(),((*node)[1]-getE1().getCenter()[1])/getE1().getRadius());
 118 }
 119
 120 bool ArcCArcCIntersector::internalAreColinears(const EdgeArcCircle& a1, const EdgeArcCircle& a2, double& distBetweenCenters, double& cst,
 121                                                double& radiusL, double centerL[2], double& radiusB, double centerB[2])
 122 {
 123   double lgth1=fabs(a1.getAngle()*a1.getRadius());
 124   double lgth2=fabs(a2.getAngle()*a2.getRadius());
 125   if(lgth1<lgth2)
 126     {//a1 is the little one ('L') and a2 the big one ('B')
 127       a1.getCenter(centerL); radiusL=a1.getRadius();
 128       a2.getCenter(centerB); radiusB=a2.getRadius();
 129     }
 130   else
 131     {
 132       a2.getCenter(centerL); radiusL=a2.getRadius();
 133       a1.getCenter(centerB); radiusB=a1.getRadius();
 134     }
 135   // dividing from the beginning by radiusB^2 to keep precision
 136   distBetweenCenters=Node::distanceBtw2PtSq(centerL,centerB);
 137   cst=distBetweenCenters/(radiusB*radiusB);
 138   cst+=radiusL*radiusL/(radiusB*radiusB);
 139   return Node::areDoubleEqualsWPRight(cst,1.,2.);
 140 }
 141
 142 bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeArcCircle& a2)
 143 {
 144   double radiusL,radiusB;
 145   double centerL[2],centerB[2];
 146   double tmp(0.),cst(0.);
 147   if(!internalAreColinears(a1,a2,tmp,cst,radiusL,centerL,radiusB,centerB))
 148     return false;
 149   //
 150   double angle0L,angleL;
 151   Bounds *merge=a1.getBounds().nearlyAmIIntersectingWith(a2.getBounds());
 152   merge->getInterceptedArc(centerL,radiusL,angle0L,angleL);
 153   delete merge;
 154   //
 155   tmp=sqrt(tmp);
 156   if(Node::areDoubleEqualsWPLeft(tmp,0.,10*std::max(radiusL,radiusB)))
 157     return Node::areDoubleEquals(radiusL,radiusB);
 158   double phi=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect((centerL[0]-centerB[0])/tmp,(centerL[1]-centerB[1])/tmp);
 159   double cst2=2*radiusL*tmp/(radiusB*radiusB);
 160   double cmpContainer[4];
 161   int sizeOfCmpContainer=2;
 162   cmpContainer[0]=cst+cst2*cos(phi-angle0L);
 163   cmpContainer[1]=cst+cst2*cos(phi-angle0L+angleL);
 164   double a=EdgeArcCircle::NormalizeAngle(phi-angle0L);
 165   if(EdgeArcCircle::IsIn2Pi(angle0L,angleL,a))
 166     cmpContainer[sizeOfCmpContainer++]=cst+cst2;
 167   a=EdgeArcCircle::NormalizeAngle(phi-angle0L+M_PI);
 168   if(EdgeArcCircle::IsIn2Pi(angle0L,angleL,a))
 169     cmpContainer[sizeOfCmpContainer++]=cst-cst2;
 170   a=*std::max_element(cmpContainer,cmpContainer+sizeOfCmpContainer);
 171   return Node::areDoubleEqualsWPRight(a,1.,2.);
 172 }
 173
 174 void ArcCArcCIntersector::areOverlappedOrOnlyColinears(bool& obviousNoIntersection, bool& areOverlapped)
 175 {
 176   _dist=Node::distanceBtw2Pt(getE1().getCenter(),getE2().getCenter());
 177   double radius1=getE1().getRadius(); double radius2=getE2().getRadius();
 178   if(_dist>radius1+radius2+QuadraticPlanarPrecision::getPrecision() || _dist+std::min(radius1,radius2)+QuadraticPlanarPrecision::getPrecision()<std::max(radius1,radius2))
 179     {
 180       obviousNoIntersection=true;
 181       areOverlapped=false;
 182       return ;
 183     }
 184   if(areArcsOverlapped(getE1(),getE2()))//(Node::areDoubleEquals(_dist,0.) && Node::areDoubleEquals(radius1,radius2))
 185     {
 186       obviousNoIntersection=false;
 187       areOverlapped=true;
 188     }
 189   else
 190     {
 191       obviousNoIntersection=false;
 192       areOverlapped=false;
 193     }
 194 }
 195
 196 /**
 197  Heart of the algorithm for arc/arc intersection.
 198  See http://mathworld.wolfram.com/Circle-CircleIntersection.html
 199  The computation is done in the coordinate system where Ox is the line between the 2 circle centers.
 200 */
 201 std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristicVal() const
 202 {
 203   std::list< IntersectElement > ret;
 204   const double *center1=getE1().getCenter();
 205   const double *center2=getE2().getCenter();
 206   double radius1=getE1().getRadius(); double radius2=getE2().getRadius();
 207   double d1_1=(_dist*_dist-radius2*radius2+radius1*radius1)/(2.*_dist);  // computation of 'x' on wolfram
 208   double u[2];//u is normalized vector from center1 to center2.
 209   u[0]=(center2[0]-center1[0])/_dist; u[1]=(center2[1]-center1[1])/_dist;
 210   double d1_1y=EdgeArcCircle::SafeSqrt(radius1*radius1-d1_1*d1_1);      // computation of 'y' on wolfram
 211   double angleE1=EdgeArcCircle::NormalizeAngle(getE1().getAngle0()+getE1().getAngle());
 212   double angleE2=EdgeArcCircle::NormalizeAngle(getE2().getAngle0()+getE2().getAngle());
 213   if(!Node::areDoubleEquals(d1_1y,0))
 214     {
 215       //2 intersections
 216       double v1[2],v2[2];
 217       // coming back to our coordinate system:
 218       v1[0]=u[0]*d1_1-u[1]*d1_1y; v1[1]=u[1]*d1_1+u[0]*d1_1y;
 219       v2[0]=u[0]*d1_1+u[1]*d1_1y; v2[1]=u[1]*d1_1-u[0]*d1_1y;
 220       Node *node1=new Node(center1[0]+v1[0],center1[1]+v1[1]); node1->declareOn();
 221       Node *node2=new Node(center1[0]+v2[0],center1[1]+v2[1]); node2->declareOn();
 222       double angle1_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1);
 223       double angle2_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v2[0]/radius1,v2[1]/radius1);
 224       double v3[2],v4[2];
 225       v3[0]=center1[0]-center2[0]+v1[0]; v3[1]=center1[1]-center2[1]+v1[1];
 226       v4[0]=center1[0]-center2[0]+v2[0]; v4[1]=center1[1]-center2[1]+v2[1];
 227       double angle1_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v3[0]/radius2,v3[1]/radius2);
 228       double angle2_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v4[0]/radius2,v4[1]/radius2);
 229       // Check whether intersection points are exactly ON the other arc or not
 230       //   -> the curvilinear distance (=radius*angle) must below eps
 231       bool e1_1S=Node::areDoubleEqualsWPLeft(angle1_1,getE1().getAngle0(),radius1);
 232       bool e1_1E=Node::areDoubleEqualsWPLeft(angle1_1,angleE1,radius1);
 233       bool e1_2S=Node::areDoubleEqualsWPLeft(angle1_2,getE2().getAngle0(),radius1);
 234       bool e1_2E=Node::areDoubleEqualsWPLeft(angle1_2,angleE2,radius1);
 235       //
 236       bool e2_1S=Node::areDoubleEqualsWPLeft(angle2_1,getE1().getAngle0(),radius2);
 237       bool e2_1E=Node::areDoubleEqualsWPLeft(angle2_1,angleE1,radius2);
 238       bool e2_2S=Node::areDoubleEqualsWPLeft(angle2_2,getE2().getAngle0(),radius2);
 239       bool e2_2E=Node::areDoubleEqualsWPLeft(angle2_2,angleE2,radius2);
 240       ret.push_back(IntersectElement(angle1_1,angle1_2,e1_1S,e1_1E,e1_2S,e1_2E,node1,_e1,_e2,keepOrder()));
 241       ret.push_back(IntersectElement(angle2_1,angle2_2,e2_1S,e2_1E,e2_2S,e2_2E,node2,_e1,_e2,keepOrder()));
 242     }
 243   else
 244     {
 245       //tangent intersection
 246       double v1[2],v2[2];
 247       v1[0]=d1_1*u[0]; v1[1]=d1_1*u[1];
 248       v2[0]=center1[0]-center2[0]+v1[0]; v2[1]=center1[1]-center2[1]+v1[1];
 249       double angle0_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1);
 250       double angle0_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v2[0]/radius2,v2[1]/radius2);
 251       bool e0_1S=Node::areDoubleEqualsWPLeft(angle0_1,getE1().getAngle0(),radius1);
 252       bool e0_1E=Node::areDoubleEqualsWPLeft(angle0_1,angleE1,radius1);
 253       bool e0_2S=Node::areDoubleEqualsWPLeft(angle0_2,getE2().getAngle0(),radius2);
 254       bool e0_2E=Node::areDoubleEqualsWPLeft(angle0_2,angleE2,radius2);
 255       Node *node=new Node(center1[0]+d1_1*u[0],center1[1]+d1_1*u[1]); node->declareOnTangent();
 256       ret.push_back(IntersectElement(angle0_1,angle0_2,e0_1S,e0_1E,e0_2S,e0_2E,node,_e1,_e2,keepOrder()));
 257     }
 258   return ret;
 259 }
 260 /*double angle0_2;
 261   double signDeltaAngle2;
 262   double d1_2;
 263   if(u[1]<0.)
 264   angle0_1=-angle0_1;
 265   if(d1_1>=0.)
 266   {
 267   if(_dist>radius1)
 268   {
 269   angle0_2=angle0_1+M_PI;
 270   signDeltaAngle2=-1.;
 271   }
 272   else
 273   {
 274   angle0_2=angle0_1;
 275   signDeltaAngle2=1.;
 276   }
 277   }
 278   else
 279   {
 280   angle0_1+=M_PI;
 281   angle0_2=angle0_1;
 282   signDeltaAngle2=1.;
 283   }
 284   angle0_1=NormalizeAngle(angle0_1);
 285   angle0_2=NormalizeAngle(angle0_2);
 286   double angleE1=NormalizeAngle(getE1().getAngle0()+getE1().getAngle());
 287   double angleE2=NormalizeAngle(getE2().getAngle0()+getE2().getAngle());
 288   if(!(Node::areDoubleEquals(d1_1,radius1) || Node::areDoubleEquals(d1_1,-radius1)) )
 289   {
 290   //2 intersections
 291   double deltaAngle1=EdgeArcCircle::SafeAcos(fabs(d1_1)/radius1); //owns to 0;Pi/2 by construction
 292   double deltaAngle2=EdgeArcCircle::SafeAcos(fabs(d1_2)/radius2); //owns to 0;Pi/2 by construction
 293   double angle1_1=NormalizeAngle(angle0_1+deltaAngle1);// Intersection 1 seen for _e1
 294   double angle2_1=NormalizeAngle(angle0_1-deltaAngle1);// Intersection 2 seen for _e1
 295   double angle1_2=NormalizeAngle(angle0_2+signDeltaAngle2*deltaAngle2);// Intersection 1 seen for _e2
 296   double angle2_2=NormalizeAngle(angle0_2-signDeltaAngle2*deltaAngle2);// Intersection 2 seen for _e2
 297   //
 298   bool e1_1S=Node::areDoubleEqualsWP(angle1_1,getE1().getAngle0(),radius1);
 299   bool e1_1E=Node::areDoubleEqualsWP(angle1_1,angleE1,radius1);
 300   bool e1_2S=Node::areDoubleEqualsWP(angle1_2,getE2().getAngle0(),radius1);
 301   bool e1_2E=Node::areDoubleEqualsWP(angle1_2,angleE2,radius1);
 302   //
 303   bool e2_1S=Node::areDoubleEqualsWP(angle2_1,getE1().getAngle0(),radius2);
 304   bool e2_1E=Node::areDoubleEqualsWP(angle2_1,angleE1,radius2);
 305   bool e2_2S=Node::areDoubleEqualsWP(angle2_2,getE2().getAngle0(),radius2);
 306   bool e2_2E=Node::areDoubleEqualsWP(angle2_2,angleE2,radius2);
 307   Node *node1=new Node(center1[0]+radius1*cos(angle1_1),center1[0]+radius1*sin(angle1_1)); node1->declareOn();
 308   Node *node2=new Node(center1[0]+radius1*cos(angle2_1),center1[0]+radius1*sin(angle2_1)); node2->declareOn();
 309   ret.push_back(IntersectElement(angle1_1,angle1_2,e1_1S,e1_1E,e1_2S,e1_2E,node1,_e1,_e2,keepOrder()));
 310   ret.push_back(IntersectElement(angle2_1,angle2_2,e2_1S,e2_1E,e2_2S,e2_2E,node2,_e1,_e2,keepOrder()));
 311   }
 312   else
 313   //tangent intersection
 314   {
 315   bool e0_1S=Node::areDoubleEqualsWP(angle0_1,getE1().getAngle0(),radius1);
 316   bool e0_1E=Node::areDoubleEqualsWP(angle0_1,angleE1,radius1);
 317   bool e0_2S=Node::areDoubleEqualsWP(angle0_2,getE2().getAngle0(),radius2);
 318   bool e0_2E=Node::areDoubleEqualsWP(angle0_2,angleE2,radius2);
 319   Node *node=new Node(center1[0]+radius1*cos(angle0_1),center1[0]+radius1*sin(angle0_1)); node->declareOnTangent();
 320   ret.push_back(IntersectElement(angle0_1,angle0_2,e0_1S,e0_1E,e0_2S,e0_2E,node,_e1,_e2,keepOrder()));
 321   }
 322   return ret;*/
 323
 324 ArcCSegIntersector::ArcCSegIntersector(const EdgeArcCircle& e1, const EdgeLin& e2, bool reverse):
 325     CrossTypeEdgeIntersector(e1,e2,reverse),
 326     _deltaRoot_div_dr(0.),
 327     _i1S2E(false),_i1E2E(false)
 328 {
 329   const double *center=getE1().getCenter();
 330   _dx=(*(_e2.getEndNode()))[0]-(*(_e2.getStartNode()))[0];
 331   _dy=(*(_e2.getEndNode()))[1]-(*(_e2.getStartNode()))[1];
 332   _drSq=_dx*_dx+_dy*_dy;
 333   _cross=
 334       ((*(_e2.getStartNode()))[0]-center[0])*((*(_e2.getEndNode()))[1]-center[1])-
 335       ((*(_e2.getStartNode()))[1]-center[1])*((*(_e2.getEndNode()))[0]-center[0]);
 336 }
 337
 338 /**
 339   See http://mathworld.wolfram.com/Circle-LineIntersection.html
 340   _cross is 'D', the computation is done with the translation to put back the circle at the origin
 341 */
 342 void ArcCSegIntersector::areOverlappedOrOnlyColinears(bool& obviousNoIntersection, bool& areOverlapped)
 343 {
 344   areOverlapped=false;//No overlapping by construction
 345
 346   // Similar optimisation than SegSegIntersector::areOverlappedOrOnlyColinears()
 347   bool dnu1, dnu2;
 348   identifyEarlyIntersection(dnu1, dnu2, _i1S2E, _i1E2E);
 349
 350   const double R = getE1().getRadius();
 351
 352   // We need to compute d = R*R-_cross*_cross/_drSq
 353   // In terms of numerical precision, this can trigger 'catastrophic cancellation' and is hence better expressed as:
 354   double _dr = sqrt(_drSq);
 355   double diff = (R-_cross/_dr), add=(R+_cross/_dr);
 356   // 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
 357   // to prevent ourselves going below machine precision (typ. 1.0e-16 for double).
 358   const double eps_machine = std::numeric_limits<double>::epsilon();
 359   diff = fabs(diff/R) < eps_machine ? 0.0 : diff;
 360   add  = fabs(add/R)  < eps_machine ? 0.0 : add;
 361   double d = add*diff;
 362   // Compute deltaRoot_div_dr := sqrt(delta)/dr, where delta has the meaning of Wolfram.
 363   // 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.
 364   // 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)
 365   _deltaRoot_div_dr = Node::sign(d)*sqrt(fabs(d));
 366
 367   if( 2*_deltaRoot_div_dr > -QuadraticPlanarPrecision::getPrecision())
 368     obviousNoIntersection=false;
 369   else
 370     obviousNoIntersection=true;
 371 }
 372
 373 /*!
 374  * By construction, no chance that an arc of circle and line to be colinear.
 375  */
 376 bool ArcCSegIntersector::areColinears() const
 377 {
 378   return false;
 379 }
 380
 381 void ArcCSegIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
 382 {
 383   throw Exception("Internal error. Should never been called : no overlapping possible between arc of circle and a segment.");
 384 }
 385
 386 std::list< IntersectElement > ArcCSegIntersector::getIntersectionsCharacteristicVal() const
 387 {
 388   std::list< IntersectElement > ret;
 389   const double *center=getE1().getCenter();
 390   if(!(2*fabs(_deltaRoot_div_dr) < QuadraticPlanarPrecision::getPrecision())) // see comments in areOverlappedOrOnlyColinears()
 391     {   // Two intersection nodes
 392       // -> 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)
 393       if(_earlyInter)
 394         {
 395           // Check tangent vector of the arc circle at the common node with the linear segment.
 396           // There we can tell if the arc of circle is 'moving away' from the seg, or if it might intersect it twice
 397           const Node &n(*_earlyInter->getNodeOnly());
 398
 399           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)
 400           bool invSeg = _i1S2E || _i1E2E;
 401           double linEdge[2] = {invSeg ? (-_dx) : _dx, invSeg ? (-_dy) : _dy};
 402           if(tang[1]*linEdge[0]-tang[0]*linEdge[1] < 0)
 403             {
 404               ret.push_back(*_earlyInter);
 405               return ret;
 406             }
 407         }
 408
 409       double determinant=fabs(_deltaRoot_div_dr)/sqrt(_drSq);
 410       double x1=(_cross*_dy/_drSq+Node::sign(_dy)*_dx*determinant)+center[0];
 411       double y1=(-_cross*_dx/_drSq+fabs(_dy)*determinant)+center[1];
 412       Node *intersect1=new Node(x1,y1); intersect1->declareOn();
 413       double x2=(_cross*_dy/_drSq-Node::sign(_dy)*_dx*determinant)+center[0];
 414       double y2=(-_cross*_dx/_drSq-fabs(_dy)*determinant)+center[1];
 415       Node *intersect2=new Node(x2,y2); intersect2->declareOn();
 416
 417       bool isN1(false), isN2(false);
 418       if (_earlyInter)
 419         {
 420           // Which node do we actually already found? Assume this is the closest ...
 421           const Node &iN = *(_earlyInter->getNodeOnly());
 422           const Node &n1(*intersect1), &n2(*intersect2);
 423           double d1 = std::max(fabs(iN[0]-n1[0]), fabs(iN[1]-n1[1]));
 424           double d2 = std::max(fabs(iN[0]-n2[0]), fabs(iN[1]-n2[1]));
 425           isN1 = d1 < d2; isN2 = !isN1;
 426           if (isN1) intersect1->decrRef();
 427           if (isN2) intersect2->decrRef();
 428           ret.push_back(*_earlyInter);
 429         }
 430       if (!isN1)
 431         {
 432           bool i1_1S=_e1.getStartNode()->isEqual(*intersect1);
 433           bool i1_1E=_e1.getEndNode()->isEqual(*intersect1);
 434           bool i1_2S=_e2.getStartNode()->isEqual(*intersect1);
 435           bool i1_2E=_e2.getEndNode()->isEqual(*intersect1);
 436           ret.push_back(IntersectElement(getE1().getCharactValue(*intersect1),getE2().getCharactValue(*intersect1),i1_1S,i1_1E,i1_2S,i1_2E,intersect1,_e1,_e2,keepOrder()));
 437         }
 438       if(!isN2)
 439         {
 440           bool i2_1S=_e1.getStartNode()->isEqual(*intersect2);
 441           bool i2_1E=_e1.getEndNode()->isEqual(*intersect2);
 442           bool i2_2S=_e2.getStartNode()->isEqual(*intersect2);
 443           bool i2_2E=_e2.getEndNode()->isEqual(*intersect2);
 444           ret.push_back(IntersectElement(getE1().getCharactValue(*intersect2),getE2().getCharactValue(*intersect2),i2_1S,i2_1E,i2_2S,i2_2E,intersect2,_e1,_e2,keepOrder()));
 445         }
 446     }
 447   else//tangent intersection
 448     {
 449       if (_earlyInter)
 450         {
 451           ret.push_back(*_earlyInter);
 452           return ret;
 453         }
 454       double x=(_cross*_dy)/_drSq+center[0];
 455       double y=(-_cross*_dx)/_drSq+center[1];
 456       Node *intersect3=new Node(x,y); intersect3->declareOnTangent();
 457       bool i_1S=_e1.getStartNode()->isEqual(*intersect3);
 458       bool i_1E=_e1.getEndNode()->isEqual(*intersect3);
 459       bool i_2S=_e2.getStartNode()->isEqual(*intersect3);
 460       bool i_2E=_e2.getEndNode()->isEqual(*intersect3);
 461       ret.push_back(IntersectElement(_e1.getCharactValue(*intersect3),_e2.getCharactValue(*intersect3),i_1S,i_1E,i_2S,i_2E,intersect3,_e1,_e2,keepOrder()));
 462     }
 463   return ret;
 464 }
 465
 466 EdgeArcCircle::EdgeArcCircle(std::istream& lineInXfig)
 467 {
 468   const unsigned NB_OF_SKIP_FIELDS=15;
 469   std::string tmpS;
 470   for(unsigned i=0;i<NB_OF_SKIP_FIELDS;i++)
 471     lineInXfig >> tmpS;
 472   _start=new Node(lineInXfig);
 473   Node *middle=new Node(lineInXfig);
 474   _end=new Node(lineInXfig);
 475   GetArcOfCirclePassingThru(*_start,*middle,*_end,_center,_radius,_angle,_angle0);
 476   middle->decrRef();
 477   updateBounds();
 478 }
 479
 480 EdgeArcCircle::EdgeArcCircle(Node *start, Node *middle, Node *end, bool direction):Edge(start,end, direction)
 481 {
 482   GetArcOfCirclePassingThru(*_start,*middle,*_end,_center,_radius,_angle,_angle0);
 483   updateBounds();
 484 }
 485
 486 EdgeArcCircle::EdgeArcCircle(double sX, double sY, double mX, double mY, double eX, double eY):Edge(sX,sY,eX,eY)
 487 {
 488   double middle[2]; middle[0]=mX; middle[1]=mY;
 489   GetArcOfCirclePassingThru(*_start,middle,*_end,_center,_radius,_angle,_angle0);
 490   updateBounds();
 491 }
 492
 493 /*!
 494  * @param angle0 in ]-Pi;Pi[
 495  * @param deltaAngle in ]-2.*Pi;2.*Pi[
 496  */
 497 EdgeArcCircle::EdgeArcCircle(Node *start, Node *end, const double *center, double radius, double angle0, double deltaAngle, bool direction):Edge(start,end,direction),_angle0(angle0),_angle(deltaAngle),
 498     _radius(radius)
 499 {
 500   _center[0]=center[0];
 501   _center[1]=center[1];
 502   updateBounds();
 503 }
 504
 505 void EdgeArcCircle::changeMiddle(Node *newMiddle)
 506 {
 507   GetArcOfCirclePassingThru(*_start,*newMiddle,*_end,_center,_radius,_angle,_angle0);
 508   updateBounds();
 509 }
 510
 511 Edge *EdgeArcCircle::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) const
 512 {
 513   double sx=((*start)[0]-_center[0])/_radius;
 514   double sy=((*start)[1]-_center[1])/_radius;
 515   double ex=((*end)[0]-_center[0])/_radius;
 516   double ey=((*end)[1]-_center[1])/_radius;
 517   double angle0=GetAbsoluteAngleOfNormalizedVect(direction?sx:ex,direction?sy:ey);
 518   double deltaAngle=GetAbsoluteAngleOfNormalizedVect(sx*ex+sy*ey,sx*ey-sy*ex);
 519   if(deltaAngle>0. && _angle<0.)
 520     deltaAngle-=2.*M_PI;
 521   else if(deltaAngle<0. && _angle>0.)
 522     deltaAngle+=2.*M_PI;
 523   deltaAngle=direction?deltaAngle:-deltaAngle;
 524   return new EdgeArcCircle(start,end,_center,_radius,angle0,deltaAngle,direction);
 525 }
 526
 527 void EdgeArcCircle::applySimilarity(double xBary, double yBary, double dimChar)
 528 {
 529   Edge::applySimilarity(xBary,yBary,dimChar);
 530   _radius/=dimChar;
 531   _center[0]=(_center[0]-xBary)/dimChar;
 532   _center[1]=(_center[1]-yBary)/dimChar;
 533 }
 534
 535 void EdgeArcCircle::unApplySimilarity(double xBary, double yBary, double dimChar)
 536 {
 537   Edge::unApplySimilarity(xBary,yBary,dimChar);
 538   _radius*=dimChar;
 539   _center[0]=_center[0]*dimChar+xBary;
 540   _center[1]=_center[1]*dimChar+yBary;
 541 }
 542
 543 /*!
 544  * 'eps' is expected to be > 0.
 545  * 'conn' is of size 3. conn[0] is start id, conn[1] is end id and conn[2] is middle id.
 546  * '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.
 547  */
 548 void EdgeArcCircle::tesselate(const mcIdType *conn, mcIdType offset, double eps, std::vector<mcIdType>& newConn, std::vector<double>& addCoo) const
 549 {
 550   newConn.push_back(INTERP_KERNEL::NORM_POLYL);
 551   int nbOfSubDiv=(int)(fabs(_angle)/eps);
 552   if(nbOfSubDiv<=2)
 553     {
 554       newConn.push_back(conn[0]); newConn.push_back(conn[2]); newConn.push_back(conn[1]);
 555       return ;
 556     }
 557   double signOfAngle=_angle>0.?1.:-1.;
 558   mcIdType offset2=offset+(ToIdType(addCoo.size()))/2;
 559   newConn.push_back(conn[0]);
 560   for(int i=1;i<nbOfSubDiv;i++,offset2++)
 561     {
 562       double angle=_angle0+i*eps*signOfAngle;
 563       newConn.push_back(offset2);
 564       addCoo.push_back(_center[0]+_radius*cos(angle)); addCoo.push_back(_center[1]+_radius*sin(angle));
 565     }
 566   newConn.push_back(conn[1]);
 567 }
 568
 569 EdgeArcCircle *EdgeArcCircle::BuildFromNodes(Node *start, Node *middle, Node *end)
 570 {
 571   EdgeLin *e1,*e2;
 572   e1=new EdgeLin(start,middle);
 573   e2=new EdgeLin(middle,end);
 574   SegSegIntersector inters(*e1,*e2);
 575   bool colinearity=inters.areColinears();
 576   delete e1; delete e2;
 577   if(colinearity)
 578     {
 579       start->decrRef(); middle->decrRef(); end->decrRef();
 580       return 0;
 581     }
 582   else
 583     {
 584       EdgeArcCircle *ret=new EdgeArcCircle(start,middle,end);
 585       start->decrRef(); middle->decrRef(); end->decrRef();
 586       return ret;
 587     }
 588 }
 589
 590 /*!
 591  * Given an \b NON normalized vector 'vect', returns its norm 'normVect' and its
 592  * angle in ]-Pi,Pi] relative to Ox axe.
 593  */
 594 double EdgeArcCircle::GetAbsoluteAngle(const double *vect, double& normVect)
 595 {
 596   normVect=Node::norm(vect);
 597   return GetAbsoluteAngleOfNormalizedVect(vect[0]/normVect,vect[1]/normVect);
 598 }
 599
 600 /*!
 601  * Given a \b normalized vector defined by (ux,uy) returns its angle in ]-Pi;Pi].
 602  * Actually in the current implementation, the vector does not even need to be normalized ...
 603  */
 604 double EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(double ux, double uy)
 605 {
 606   return atan2(uy, ux);
 607 }
 608
 609 void EdgeArcCircle::GetArcOfCirclePassingThru(const double *start, const double *middle, const double *end,
 610                                               double *center, double& radius, double& angleInRad, double& angleInRad0)
 611 {
 612   double delta=(middle[0]-start[0])*(end[1]-middle[1])-(end[0]-middle[0])*(middle[1]-start[1]);
 613   double b1=(middle[1]*middle[1]+middle[0]*middle[0]-start[0]*start[0]-start[1]*start[1])/2;
 614   double b2=(end[1]*end[1]+end[0]*end[0]-middle[0]*middle[0]-middle[1]*middle[1])/2;
 615   center[0]=((end[1]-middle[1])*b1+(start[1]-middle[1])*b2)/delta;
 616   center[1]=((middle[0]-end[0])*b1+(middle[0]-start[0])*b2)/delta;
 617   radius=SafeSqrt((start[0]-center[0])*(start[0]-center[0])+(start[1]-center[1])*(start[1]-center[1]));
 618   angleInRad0=GetAbsoluteAngleOfNormalizedVect((start[0]-center[0])/radius,(start[1]-center[1])/radius);
 619   double angleInRadM=GetAbsoluteAngleOfNormalizedVect((middle[0]-center[0])/radius,(middle[1]-center[1])/radius);
 620   angleInRad=GetAbsoluteAngleOfNormalizedVect(((start[0]-center[0])*(end[0]-center[0])+(start[1]-center[1])*(end[1]-center[1]))/(radius*radius),
 621       ((start[0]-center[0])*(end[1]-center[1])-(start[1]-center[1])*(end[0]-center[0]))/(radius*radius));
 622   if(IsAngleNotIn(angleInRad0,angleInRad,angleInRadM))
 623     angleInRad=angleInRad<0?2*M_PI+angleInRad:angleInRad-2*M_PI;
 624 }
 625
 626 void EdgeArcCircle::dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const
 627 {
 628   stream << "5 1 0 1 ";
 629   fillXfigStreamForLoc(stream);
 630   stream << " 7 50 -1 -1 0.000 0 ";
 631   if( (direction && (-_angle)>=0) || (!direction && (-_angle)<0))
 632     stream << '0';//'0'
 633   else
 634     stream << '1';//'1'
 635   stream << " 1 0 ";
 636   stream << box.fitXForXFigD(_center[0],resolution) << " " << box.fitYForXFigD(_center[1],resolution) << " ";
 637   direction?_start->dumpInXfigFile(stream,resolution,box):_end->dumpInXfigFile(stream,resolution,box);
 638   Node *middle=buildRepresentantOfMySelf();
 639   middle->dumpInXfigFile(stream,resolution,box);
 640   middle->decrRef();
 641   direction?_end->dumpInXfigFile(stream,resolution,box):_start->dumpInXfigFile(stream,resolution,box);
 642   stream << std::endl << "1 1 2.00 120.00 180.00" << std::endl;
 643 }
 644
 645 void EdgeArcCircle::update(Node *m)
 646 {
 647   GetArcOfCirclePassingThru(*_start,*m,*_end,_center,_radius,_angle,_angle0);
 648   updateBounds();
 649 }
 650
 651 /*!
 652  * This methods computes :
 653  * \f[
 654  * \int_{Current Edge} -ydx
 655  * \f]
 656  */
 657 double EdgeArcCircle::getAreaOfZone() const
 658 {
 659   return -_radius*_radius*(sin(_angle)-_angle)/2.+((*_start)[0]-(*_end)[0])*((*_start)[1]+(*_end)[1])/2.;
 660 }
 661
 662 double EdgeArcCircle::getCurveLength() const
 663 {
 664   return fabs(_angle*_radius);
 665 }
 666
 667 void EdgeArcCircle::getBarycenter(double *bary) const
 668 {
 669   bary[0]=_center[0]+_radius*cos(_angle0+_angle/2.);
 670   bary[1]=_center[1]+_radius*sin(_angle0+_angle/2.);
 671 }
 672
 673 /*!
 674  * \f[
 675  * bary[0]=\int_{Current Edge} -yxdx
 676  * \f]
 677  * \f[
 678  * bary[1]=\int_{Current Edge} -\frac{y^{2}}{2}dx
 679  * \f]
 680  * To compute these 2 expressions in this class we have :
 681  * \f[
 682  * x=x_{0}+Radius \cdot cos(\theta)
 683  * \f]
 684  * \f[
 685  * y=y_{0}+Radius \cdot sin(\theta)
 686  * \f]
 687  * \f[
 688  * dx=-Radius \cdot sin(\theta) \cdot d\theta
 689  * \f]
 690  */
 691 void EdgeArcCircle::getBarycenterOfZone(double *bary) const
 692 {
 693   double x0=_center[0];
 694   double y0=_center[1];
 695   double angle1=_angle0+_angle;
 696   double tmp1=sin(angle1);
 697   double tmp0=sin(_angle0);
 698   double tmp2=_radius*_radius*_radius;
 699   double tmp3=cos(angle1);
 700   double tmp4=cos(_angle0);
 701   bary[0]=_radius*x0*y0*(tmp4-tmp3)+_radius*_radius*(y0*(cos(2*_angle0)-cos(2*angle1))/4.+
 702       x0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.))
 703       +tmp2*(tmp1*tmp1*tmp1-tmp0*tmp0*tmp0)/3.;
 704   bary[1]=y0*y0*_radius*(tmp4-tmp3)/2.+_radius*_radius*y0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.)
 705         +tmp2*(tmp4-tmp3+(tmp3*tmp3*tmp3-tmp4*tmp4*tmp4)/3.)/2.;
 706 }
 707
 708 /**
 709  * Compute the "middle" of two points on the arc of circle.
 710  * The order (p1,p2) or (p2,p1) doesn't matter. p1 and p2 have to be localized on the edge defined by this.
 711  * \param[out] mid the point located half-way between p1 and p2 on the arc defined by this.
 712  * \sa getMiddleOfPointsOriented() a generalisation working also when p1 and p2 are not on the arc.
 713  */
 714 void EdgeArcCircle::getMiddleOfPoints(const double *p1, const double *p2, double *mid) const
 715 {
 716   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);
 717   double angle1(GetAbsoluteAngleOfNormalizedVect(dx1,dy1)),angle2(GetAbsoluteAngleOfNormalizedVect(dx2,dy2));
 718   //
 719   double myDelta1(angle1-_angle0),myDelta2(angle2-_angle0);
 720   if(_angle>0.)
 721     { myDelta1=myDelta1>-QuadraticPlanarPrecision::getPrecision()?myDelta1:myDelta1+2.*M_PI; myDelta2=myDelta2>-QuadraticPlanarPrecision::getPrecision()?myDelta2:myDelta2+2.*M_PI; }
 722   else
 723     { myDelta1=myDelta1<QuadraticPlanarPrecision::getPrecision()?myDelta1:myDelta1-2.*M_PI; myDelta2=myDelta2<QuadraticPlanarPrecision::getPrecision()?myDelta2:myDelta2-2.*M_PI; }
 724   ////
 725   mid[0]=_center[0]+_radius*cos(_angle0+(myDelta1+myDelta2)/2.);
 726   mid[1]=_center[1]+_radius*sin(_angle0+(myDelta1+myDelta2)/2.);
 727 }
 728
 729 /**
 730  * Compute the "middle" of two points on the arc of circle.
 731  * Walk on the circle from p1 to p2 using the rotation direction indicated by this->_angle (i.e. by the orientation of the arc).
 732  * This function is sensitive to the ordering of p1 and p2.
 733  * \param[out] mid the point located half-way between p1 and p2
 734  * \sa getMiddleOfPoints() to be used when the order of p1 and p2 is not relevant.
 735  */
 736 void EdgeArcCircle::getMiddleOfPointsOriented(const double *p1, const double *p2, double *mid) const
 737 {
 738   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);
 739   double angle1(GetAbsoluteAngleOfNormalizedVect(dx1,dy1)),angle2(GetAbsoluteAngleOfNormalizedVect(dx2,dy2));
 740
 741   if (angle1 <= 0.0)
 742     angle1 += 2.*M_PI;
 743   if (angle2 <= 0.0)
 744     angle2 += 2.*M_PI;
 745
 746   double avg;
 747   if((_angle>0. && angle1 <= angle2) || (_angle<=0. && angle1 >= angle2))
 748     avg = (angle1+angle2)/2.;
 749   else
 750     avg = (angle1+angle2)/2. - M_PI;
 751
 752   mid[0]=_center[0]+_radius*cos(avg);
 753   mid[1]=_center[1]+_radius*sin(avg);
 754 }
 755
 756
 757 /*!
 758  * Characteristic value used is angle in ]_Pi;Pi[ from axe 0x.
 759  */
 760 bool EdgeArcCircle::isIn(double characterVal) const
 761 {
 762   return IsIn2Pi(_angle0,_angle,characterVal);
 763 }
 764
 765 Node *EdgeArcCircle::buildRepresentantOfMySelf() const
 766 {
 767   return new Node(_center[0]+_radius*cos(_angle0+_angle/2.),_center[1]+_radius*sin(_angle0+_angle/2.));
 768 }
 769
 770 /*!
 771  * Characteristic value used is angle in ]_Pi;Pi[ from axe 0x.
 772  * 'val1' and 'val2' have been detected previously as owning to this.
 773  */
 774 bool EdgeArcCircle::isLower(double val1, double val2) const
 775 {
 776   double myDelta1=val1-_angle0;
 777   double myDelta2=val2-_angle0;
 778   if(_angle>0.)
 779     {
 780       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.
 781       myDelta2=myDelta2>-(_radius*QuadraticPlanarPrecision::getPrecision())?myDelta2:myDelta2+2.*M_PI;
 782       return myDelta1<myDelta2;
 783     }
 784   else
 785     {
 786       myDelta1=myDelta1<(_radius*QuadraticPlanarPrecision::getPrecision())?myDelta1:myDelta1-2.*M_PI;
 787       myDelta2=myDelta2<(_radius*QuadraticPlanarPrecision::getPrecision())?myDelta2:myDelta2-2.*M_PI;
 788       return myDelta2<myDelta1;
 789     }
 790 }
 791
 792 /*!
 793  * For Arc circle the caract value is angle with Ox between -Pi and Pi.
 794  */
 795 double EdgeArcCircle::getCharactValue(const Node& node) const
 796 {
 797   double dx=(node[0]-_center[0])/_radius;
 798   double dy=(node[1]-_center[1])/_radius;
 799   return GetAbsoluteAngleOfNormalizedVect(dx,dy);
 800 }
 801
 802 double EdgeArcCircle::getCharactValueBtw0And1(const Node& node) const
 803 {
 804   double dx=(node[0]-_center[0])/_radius;
 805   double dy=(node[1]-_center[1])/_radius;
 806   double angle=GetAbsoluteAngleOfNormalizedVect(dx,dy);
 807   //
 808   double myDelta=angle-_angle0;
 809   if(_angle>0.)
 810     myDelta=myDelta>=0.?myDelta:myDelta+2.*M_PI;
 811   else
 812     myDelta=myDelta<=0.?myDelta:myDelta-2.*M_PI;
 813   return myDelta/_angle;
 814 }
 815
 816 double EdgeArcCircle::getDistanceToPoint(const double *pt) const
 817 {
 818   double angle=Node::computeAngle(_center,pt);
 819   if(IsIn2Pi(_angle0,_angle,angle))
 820     return fabs(Node::distanceBtw2Pt(_center,pt)-_radius);
 821   else
 822     {
 823       double dist1=Node::distanceBtw2Pt(*_start,pt);
 824       double dist2=Node::distanceBtw2Pt(*_end,pt);
 825       return std::min(dist1,dist2);
 826     }
 827 }
 828
 829 bool EdgeArcCircle::isNodeLyingOn(const double *coordOfNode) const
 830 {
 831   double dist=Node::distanceBtw2Pt(_center,coordOfNode);
 832   if(Node::areDoubleEquals(dist,_radius))
 833     {
 834       double angle=Node::computeAngle(_center,coordOfNode);
 835       return IsIn2Pi(_angle0,_angle,angle);
 836     }
 837   else
 838     return false;
 839 }
 840
 841 /*!
 842  * Idem IsAngleNotIn except that here 'start' in ]-Pi;Pi[ and delta in ]-2*Pi;2Pi[.
 843  * @param angleIn in ]-Pi;Pi[.
 844  */
 845 bool EdgeArcCircle::IsIn2Pi(double start, double delta, double angleIn)
 846 {
 847   double myDelta=angleIn-start;
 848   if(delta>0.)
 849     {
 850       myDelta=myDelta>=0.?myDelta:myDelta+2.*M_PI;
 851       return myDelta>0. && myDelta<delta;
 852     }
 853   else
 854     {
 855       myDelta=myDelta<=0.?myDelta:myDelta-2.*M_PI;
 856       return myDelta<0. && myDelta>delta;
 857     }
 858 }
 859
 860 /*!
 861  * Given the arc 'a' defined by 'start' angle and a 'delta' [-Pi;Pi] states for the angle 'angleIn' [-Pi;Pi] if it owns or not 'a'.
 862  */
 863 bool EdgeArcCircle::IsAngleNotIn(double start, double delta, double angleIn)
 864 {
 865   double tmp=start;
 866   if(tmp<0.)
 867     tmp+=2*M_PI;
 868   double tmp2=angleIn;
 869   if(tmp2<0.)
 870     tmp2+=2*M_PI;
 871   if(tmp+delta>=2.*M_PI)
 872     return (tmp2<tmp) && (tmp2>tmp+delta-2*M_PI);
 873   else if(tmp+delta>=0.)
 874     return (tmp2<std::min(tmp,tmp+delta) || tmp2>std::max(tmp,tmp+delta));
 875   else
 876     return (tmp2>tmp) && (tmp2<(tmp+delta+2.*M_PI));
 877 }
 878
 879 void EdgeArcCircle::updateBounds()
 880 {
 881   _bounds.setValues(std::min((*_start)[0],(*_end)[0]),std::max((*_start)[0],(*_end)[0]),std::min((*_start)[1],(*_end)[1]),std::max((*_start)[1],(*_end)[1]));
 882   if(IsIn2Pi(_angle0,_angle,M_PI/2))
 883     _bounds[3]=_center[1]+_radius;
 884   if(IsIn2Pi(_angle0,_angle,-M_PI/2))
 885     _bounds[2]=_center[1]-_radius;
 886   if(IsIn2Pi(_angle0,_angle,0.))
 887     _bounds[1]=_center[0]+_radius;
 888   if(IsIn2Pi(_angle0,_angle,M_PI))
 889     _bounds[0]=_center[0]-_radius;
 890 }