src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.cxx

   1 //  Copyright (C) 2007-2008  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.
   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 #include "EdgeArcCircle.hxx"
  20 #include "EdgeLin.hxx"
  21 #include "InterpKernelException.hxx"
  22 #include "Node.hxx"
  23
  24 #include <sstream>
  25 #include <algorithm>
  26
  27 using namespace std;
  28 using namespace INTERP_KERNEL;
  29
  30 ArcCArcCIntersector::ArcCArcCIntersector(const EdgeArcCircle& e1, const EdgeArcCircle& e2):SameTypeEdgeIntersector(e1,e2),_dist(0.)
  31 {
  32 }
  33
  34 bool ArcCArcCIntersector::haveTheySameDirection() const
  35 {
  36   return (getE1().getAngle()>0. &&  getE2().getAngle()>0.) || (getE1().getAngle()<0. &&  getE2().getAngle()<0.);
  37 }
  38
  39 /*!
  40  * Precondition 'start' and 'end' are on the same curve than this.
  41  */
  42 void ArcCArcCIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
  43 {
  44   bool obvious1,obvious2;
  45   obviousCaseForCurvAbscisse(start,whereStart,commonNode,obvious1);
  46   obviousCaseForCurvAbscisse(end,whereEnd,commonNode,obvious2);
  47   if(obvious1 && obvious2)
  48     return ;
  49   double angleInRadStart=getAngle(start);
  50   double angleInRadEnd=getAngle(end);
  51   if(obvious1 || obvious2)
  52     {
  53       if(obvious1)
  54         {
  55           if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
  56             whereEnd=INSIDE;
  57           else
  58             whereEnd=OUT_AFTER;
  59           return ;
  60         }
  61       else
  62         {
  63           if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart))
  64             whereStart=INSIDE;
  65           else
  66             whereStart=OUT_BEFORE;
  67           return ;
  68         }
  69     }
  70   if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart))
  71     {
  72       whereStart=INSIDE;
  73       if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
  74         whereEnd=INSIDE;
  75       else
  76         whereEnd=OUT_AFTER;
  77     }
  78   else
  79     {//we are out in start.
  80       if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
  81         {
  82           whereStart=OUT_BEFORE;
  83           whereEnd=INSIDE;
  84         }
  85       else
  86         {
  87           if(EdgeArcCircle::isIn2Pi(getE2().getAngle0(),getE2().getAngle(),getE1().getAngle0()))
  88             {//_e2 contains stictly _e1
  89               whereStart=OUT_BEFORE;
  90               whereEnd=OUT_AFTER;
  91             }
  92           else
  93             {//_e2 is outside from _e1
  94               whereStart=OUT_BEFORE;
  95               whereEnd=OUT_BEFORE;
  96             }
  97         }
  98     }
  99 }
 100
 101 /*!
 102  * Return angle between ]-Pi;Pi[
 103  */
 104 double ArcCArcCIntersector::getAngle(Node *node) const
 105 {
 106   return EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(((*node)[0]-getE1().getCenter()[0])/getE1().getRadius(),((*node)[1]-getE1().getCenter()[1])/getE1().getRadius());
 107 }
 108
 109 bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeArcCircle& a2)
 110 {
 111   double centerL[2],radiusL,angle0L,angleL;
 112   double centerB[2],radiusB;
 113   double lgth1=fabs(a1.getAngle()*a1.getRadius());
 114   double lgth2=fabs(a2.getAngle()*a2.getRadius());
 115   if(lgth1<lgth2)
 116     {//a1 is the little one ('L') and a2 the big one ('B')
 117       a1.getCenter(centerL); radiusL=a1.getRadius(); angle0L=a1.getAngle0(); angleL=a1.getAngle();
 118       a2.getCenter(centerB); radiusB=a2.getRadius();
 119     }
 120   else
 121     {
 122       a2.getCenter(centerL); radiusL=a2.getRadius(); angle0L=a2.getAngle0(); angleL=a2.getAngle();
 123       a1.getCenter(centerB); radiusB=a1.getRadius();
 124     }
 125   // dividing from the begining by radiusB^2 to keep precision
 126   double tmp=Node::distanceBtw2PtSq(centerL,centerB);
 127   double cst=tmp/(radiusB*radiusB);
 128   cst+=radiusL*radiusL/(radiusB*radiusB);
 129   if(!Node::areDoubleEqualsWP(cst,1.,2.))
 130     return false;
 131   //
 132   Bounds *merge=a1.getBounds().nearlyAmIIntersectingWith(a2.getBounds());
 133   merge->getInterceptedArc(centerL,radiusL,angle0L,angleL);
 134   delete merge;
 135   //
 136   tmp=sqrt(tmp);
 137   if(Node::areDoubleEqualsWP(tmp,0.,1/(10*std::max(radiusL,radiusB))))
 138     {
 139       if(Node::areDoubleEquals(radiusL,radiusB))
 140         return true;
 141       else
 142         return false;
 143     }
 144   double phi=EdgeArcCircle::getAbsoluteAngleOfNormalizedVect((centerL[0]-centerB[0])/tmp,(centerL[1]-centerB[1])/tmp);
 145   double cst2=2*radiusL*tmp/(radiusB*radiusB);
 146   double cmpContainer[4];
 147   int sizeOfCmpContainer=2;
 148   cmpContainer[0]=cst+cst2*cos(phi-angle0L);
 149   cmpContainer[1]=cst+cst2*cos(phi-angle0L+angleL);
 150   double a=EdgeArcCircle::normalizeAngle(phi-angle0L);
 151   if(EdgeArcCircle::isIn2Pi(angle0L,angleL,a))
 152     cmpContainer[sizeOfCmpContainer++]=cst+cst2;
 153   a=EdgeArcCircle::normalizeAngle(phi-angle0L+M_PI);
 154   if(EdgeArcCircle::isIn2Pi(angle0L,angleL,a))
 155     cmpContainer[sizeOfCmpContainer++]=cst-cst2;
 156   a=*max_element(cmpContainer,cmpContainer+sizeOfCmpContainer);
 157   return Node::areDoubleEqualsWP(a,1.,2.);
 158 }
 159
 160 void ArcCArcCIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped)
 161 {
 162   _dist=Node::distanceBtw2Pt(getE1().getCenter(),getE2().getCenter());
 163   double radius1=getE1().getRadius(); double radius2=getE2().getRadius();
 164   if(_dist>radius1+radius2+QUADRATIC_PLANAR::_precision || _dist+std::min(radius1,radius2)+QUADRATIC_PLANAR::_precision<std::max(radius1,radius2))
 165     {
 166       obviousNoIntersection=true;
 167       areOverlapped=false;
 168       return ;
 169     }
 170   if(areArcsOverlapped(getE1(),getE2()))//(Node::areDoubleEquals(_dist,0.) && Node::areDoubleEquals(radius1,radius2))
 171     {
 172       obviousNoIntersection=false;
 173       areOverlapped=true;
 174     }
 175   else
 176     {
 177       obviousNoIntersection=false;
 178       areOverlapped=false;
 179     }
 180 }
 181
 182 std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristicVal() const
 183 {
 184   std::list< IntersectElement > ret;
 185   const double *center1=getE1().getCenter();
 186   const double *center2=getE2().getCenter();
 187   double radius1=getE1().getRadius(); double radius2=getE2().getRadius();
 188   double d1_1=(_dist*_dist-radius2*radius2+radius1*radius1)/(2.*_dist);
 189   double u[2];//u is normalized vector from center1 to center2.
 190   u[0]=(center2[0]-center1[0])/_dist; u[1]=(center2[1]-center1[1])/_dist;
 191   double d1_1y=EdgeArcCircle::safeSqrt(radius1*radius1-d1_1*d1_1);
 192   double angleE1=EdgeArcCircle::normalizeAngle(getE1().getAngle0()+getE1().getAngle());
 193   double angleE2=EdgeArcCircle::normalizeAngle(getE2().getAngle0()+getE2().getAngle());
 194   if(!Node::areDoubleEquals(d1_1y,0))
 195     {
 196       //2 intersections
 197       double v1[2],v2[2];
 198       v1[0]=u[0]*d1_1-u[1]*d1_1y; v1[1]=u[1]*d1_1+u[0]*d1_1y;
 199       v2[0]=u[0]*d1_1+u[1]*d1_1y; v2[1]=u[1]*d1_1-u[0]*d1_1y;
 200       Node *node1=new Node(center1[0]+v1[0],center1[1]+v1[1]); node1->declareOn();
 201       Node *node2=new Node(center1[0]+v2[0],center1[1]+v2[1]); node2->declareOn();
 202       double angle1_1=EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1);
 203       double angle2_1=EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(v2[0]/radius1,v2[1]/radius1);
 204       double v3[2],v4[2];
 205       v3[0]=center1[0]-center2[0]+v1[0]; v3[1]=center1[1]-center2[1]+v1[1];
 206       v4[0]=center1[0]-center2[0]+v2[0]; v4[1]=center1[1]-center2[1]+v2[1];
 207       double angle1_2=EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(v3[0]/radius2,v3[1]/radius2);
 208       double angle2_2=EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(v4[0]/radius2,v4[1]/radius2);
 209       //
 210       bool e1_1S=Node::areDoubleEqualsWP(angle1_1,getE1().getAngle0(),radius1);
 211       bool e1_1E=Node::areDoubleEqualsWP(angle1_1,angleE1,radius1);
 212       bool e1_2S=Node::areDoubleEqualsWP(angle1_2,getE2().getAngle0(),radius1);
 213       bool e1_2E=Node::areDoubleEqualsWP(angle1_2,angleE2,radius1);
 214       //
 215       bool e2_1S=Node::areDoubleEqualsWP(angle2_1,getE1().getAngle0(),radius2);
 216       bool e2_1E=Node::areDoubleEqualsWP(angle2_1,angleE1,radius2);
 217       bool e2_2S=Node::areDoubleEqualsWP(angle2_2,getE2().getAngle0(),radius2);
 218       bool e2_2E=Node::areDoubleEqualsWP(angle2_2,angleE2,radius2);
 219       ret.push_back(IntersectElement(angle1_1,angle1_2,e1_1S,e1_1E,e1_2S,e1_2E,node1,_e1,_e2,keepOrder()));
 220       ret.push_back(IntersectElement(angle2_1,angle2_2,e2_1S,e2_1E,e2_2S,e2_2E,node2,_e1,_e2,keepOrder()));
 221     }
 222   else
 223     {
 224       //tangent intersection
 225       double v1[2],v2[2];
 226       v1[0]=d1_1*u[0]; v1[1]=d1_1*u[1];
 227       v2[0]=center1[0]-center2[0]+v1[0]; v2[1]=center1[1]-center2[1]+v1[1];
 228       double angle0_1=EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1);
 229       double angle0_2=EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(v2[0]/radius2,v2[1]/radius2);
 230       bool e0_1S=Node::areDoubleEqualsWP(angle0_1,getE1().getAngle0(),radius1);
 231       bool e0_1E=Node::areDoubleEqualsWP(angle0_1,angleE1,radius1);
 232       bool e0_2S=Node::areDoubleEqualsWP(angle0_2,getE2().getAngle0(),radius2);
 233       bool e0_2E=Node::areDoubleEqualsWP(angle0_2,angleE2,radius2);
 234       Node *node=new Node(center1[0]+d1_1*u[0],center1[1]+d1_1*u[1]); node->declareOnTangent();
 235       ret.push_back(IntersectElement(angle0_1,angle0_2,e0_1S,e0_1E,e0_2S,e0_2E,node,_e1,_e2,keepOrder()));
 236     }
 237   return ret;
 238 }
 239 /*double angle0_2;
 240   double signDeltaAngle2;
 241   double d1_2;
 242   if(u[1]<0.)
 243   angle0_1=-angle0_1;
 244   if(d1_1>=0.)
 245   {
 246   if(_dist>radius1)
 247   {
 248   angle0_2=angle0_1+M_PI;
 249   signDeltaAngle2=-1.;
 250   }
 251   else
 252   {
 253   angle0_2=angle0_1;
 254   signDeltaAngle2=1.;
 255   }
 256   }
 257   else
 258   {
 259   angle0_1+=M_PI;
 260   angle0_2=angle0_1;
 261   signDeltaAngle2=1.;
 262   }
 263   angle0_1=normalizeAngle(angle0_1);
 264   angle0_2=normalizeAngle(angle0_2);
 265   double angleE1=normalizeAngle(getE1().getAngle0()+getE1().getAngle());
 266   double angleE2=normalizeAngle(getE2().getAngle0()+getE2().getAngle());
 267   if(!(Node::areDoubleEquals(d1_1,radius1) || Node::areDoubleEquals(d1_1,-radius1)) )
 268   {
 269   //2 intersections
 270   double deltaAngle1=EdgeArcCircle::safeAcos(fabs(d1_1)/radius1); //owns to 0;Pi/2 by construction
 271   double deltaAngle2=EdgeArcCircle::safeAcos(fabs(d1_2)/radius2); //owns to 0;Pi/2 by construction
 272   double angle1_1=normalizeAngle(angle0_1+deltaAngle1);// Intersection 1 seen for _e1
 273   double angle2_1=normalizeAngle(angle0_1-deltaAngle1);// Intersection 2 seen for _e1
 274   double angle1_2=normalizeAngle(angle0_2+signDeltaAngle2*deltaAngle2);// Intersection 1 seen for _e2
 275   double angle2_2=normalizeAngle(angle0_2-signDeltaAngle2*deltaAngle2);// Intersection 2 seen for _e2
 276   //
 277   bool e1_1S=Node::areDoubleEqualsWP(angle1_1,getE1().getAngle0(),radius1);
 278   bool e1_1E=Node::areDoubleEqualsWP(angle1_1,angleE1,radius1);
 279   bool e1_2S=Node::areDoubleEqualsWP(angle1_2,getE2().getAngle0(),radius1);
 280   bool e1_2E=Node::areDoubleEqualsWP(angle1_2,angleE2,radius1);
 281   //
 282   bool e2_1S=Node::areDoubleEqualsWP(angle2_1,getE1().getAngle0(),radius2);
 283   bool e2_1E=Node::areDoubleEqualsWP(angle2_1,angleE1,radius2);
 284   bool e2_2S=Node::areDoubleEqualsWP(angle2_2,getE2().getAngle0(),radius2);
 285   bool e2_2E=Node::areDoubleEqualsWP(angle2_2,angleE2,radius2);
 286   Node *node1=new Node(center1[0]+radius1*cos(angle1_1),center1[0]+radius1*sin(angle1_1)); node1->declareOn();
 287   Node *node2=new Node(center1[0]+radius1*cos(angle2_1),center1[0]+radius1*sin(angle2_1)); node2->declareOn();
 288   ret.push_back(IntersectElement(angle1_1,angle1_2,e1_1S,e1_1E,e1_2S,e1_2E,node1,_e1,_e2,keepOrder()));
 289   ret.push_back(IntersectElement(angle2_1,angle2_2,e2_1S,e2_1E,e2_2S,e2_2E,node2,_e1,_e2,keepOrder()));
 290   }
 291   else
 292   //tangent intersection
 293   {
 294   bool e0_1S=Node::areDoubleEqualsWP(angle0_1,getE1().getAngle0(),radius1);
 295   bool e0_1E=Node::areDoubleEqualsWP(angle0_1,angleE1,radius1);
 296   bool e0_2S=Node::areDoubleEqualsWP(angle0_2,getE2().getAngle0(),radius2);
 297   bool e0_2E=Node::areDoubleEqualsWP(angle0_2,angleE2,radius2);
 298   Node *node=new Node(center1[0]+radius1*cos(angle0_1),center1[0]+radius1*sin(angle0_1)); node->declareOnTangent();
 299   ret.push_back(IntersectElement(angle0_1,angle0_2,e0_1S,e0_1E,e0_2S,e0_2E,node,_e1,_e2,keepOrder()));
 300   }
 301   return ret;*/
 302
 303 ArcCSegIntersector::ArcCSegIntersector(const EdgeArcCircle& e1, const EdgeLin& e2, bool reverse):CrossTypeEdgeIntersector(e1,e2,reverse)
 304 {
 305 }
 306
 307 void ArcCSegIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped)
 308 {
 309   areOverlapped=false;//No overlapping by contruction
 310   const double *center=getE1().getCenter();
 311   _dx=(*(_e2.getEndNode()))[0]-(*(_e2.getStartNode()))[0];
 312   _dy=(*(_e2.getEndNode()))[1]-(*(_e2.getStartNode()))[1];
 313   _drSq=_dx*_dx+_dy*_dy;
 314   _cross=
 315     ((*(_e2.getStartNode()))[0]-center[0])*((*(_e2.getEndNode()))[1]-center[1])-
 316     ((*(_e2.getStartNode()))[1]-center[1])*((*(_e2.getEndNode()))[0]-center[0]);
 317   _determinant=getE1().getRadius()*getE1().getRadius()/_drSq-_cross*_cross/(_drSq*_drSq);
 318   if(_determinant>-2*QUADRATIC_PLANAR::_precision)//QUADRATIC_PLANAR::_precision*QUADRATIC_PLANAR::_precision*_drSq*_drSq/(2.*_dx*_dx))
 319     obviousNoIntersection=false;
 320   else
 321     obviousNoIntersection=true;
 322 }
 323
 324 void ArcCSegIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
 325 {
 326   throw Exception("Internal error. Should never been called : no overlapping possible between arc of circle and a segment.");
 327 }
 328
 329 std::list< IntersectElement > ArcCSegIntersector::getIntersectionsCharacteristicVal() const
 330 {
 331   std::list< IntersectElement > ret;
 332   const double *center=getE1().getCenter();
 333   if(!(fabs(_determinant)<(2.*QUADRATIC_PLANAR::_precision)))//QUADRATIC_PLANAR::_precision*QUADRATIC_PLANAR::_precision*_drSq*_drSq/(2.*_dx*_dx))
 334     {
 335       double determinant=EdgeArcCircle::safeSqrt(_determinant);
 336       double x1=(_cross*_dy/_drSq+Node::sign(_dy)*_dx*determinant)+center[0];
 337       double y1=(-_cross*_dx/_drSq+fabs(_dy)*determinant)+center[1];
 338       Node *intersect1=new Node(x1,y1); intersect1->declareOn();
 339       bool i1_1S=_e1.getStartNode()->isEqual(*intersect1);
 340       bool i1_1E=_e1.getEndNode()->isEqual(*intersect1);
 341       bool i1_2S=_e2.getStartNode()->isEqual(*intersect1);
 342       bool i1_2E=_e2.getEndNode()->isEqual(*intersect1);
 343       ret.push_back(IntersectElement(getE1().getCharactValue(*intersect1),getE2().getCharactValue(*intersect1),i1_1S,i1_1E,i1_2S,i1_2E,intersect1,_e1,_e2,keepOrder()));
 344       //
 345       double x2=(_cross*_dy/_drSq-Node::sign(_dy)*_dx*determinant)+center[0];
 346       double y2=(-_cross*_dx/_drSq-fabs(_dy)*determinant)+center[1];
 347       Node *intersect2=new Node(x2,y2); intersect2->declareOn();
 348       bool i2_1S=_e1.getStartNode()->isEqual(*intersect2);
 349       bool i2_1E=_e1.getEndNode()->isEqual(*intersect2);
 350       bool i2_2S=_e2.getStartNode()->isEqual(*intersect2);
 351       bool i2_2E=_e2.getEndNode()->isEqual(*intersect2);
 352       ret.push_back(IntersectElement(getE1().getCharactValue(*intersect2),getE2().getCharactValue(*intersect2),i2_1S,i2_1E,i2_2S,i2_2E,intersect2,_e1,_e2,keepOrder()));
 353     }
 354   else//tangent intersection
 355     {
 356       double x=(_cross*_dy)/_drSq+center[0];
 357       double y=(-_cross*_dx)/_drSq+center[1];
 358       Node *intersect=new Node(x,y); intersect->declareOnTangent();
 359       bool i_1S=_e1.getStartNode()->isEqual(*intersect);
 360       bool i_1E=_e1.getEndNode()->isEqual(*intersect);
 361       bool i_2S=_e2.getStartNode()->isEqual(*intersect);
 362       bool i_2E=_e2.getEndNode()->isEqual(*intersect);
 363       ret.push_back(IntersectElement(_e1.getCharactValue(*intersect),_e2.getCharactValue(*intersect),i_1S,i_1E,i_2S,i_2E,intersect,_e1,_e2,keepOrder()));
 364     }
 365   return ret;
 366 }
 367
 368 EdgeArcCircle::EdgeArcCircle(std::istream& lineInXfig)
 369 {
 370   const unsigned NB_OF_SKIP_FIELDS=15;
 371   std::string tmpS;
 372   for(unsigned i=0;i<NB_OF_SKIP_FIELDS;i++)
 373     lineInXfig >> tmpS;
 374   _start=new Node(lineInXfig);
 375   Node *middle=new Node(lineInXfig);
 376   _end=new Node(lineInXfig);
 377   getArcOfCirclePassingThru(*_start,*middle,*_end,_center,_radius,_angle,_angle0);
 378   middle->decrRef();
 379   updateBounds();
 380 }
 381
 382 EdgeArcCircle::EdgeArcCircle(Node *start, Node *middle, Node *end, bool direction):Edge(start,end, direction)
 383 {
 384   getArcOfCirclePassingThru(*_start,*middle,*_end,_center,_radius,_angle,_angle0);
 385   updateBounds();
 386 }
 387
 388 EdgeArcCircle::EdgeArcCircle(double sX, double sY, double mX, double mY, double eX, double eY):Edge(sX,sY,eX,eY)
 389 {
 390   double middle[2]; middle[0]=mX; middle[1]=mY;
 391   getArcOfCirclePassingThru(*_start,middle,*_end,_center,_radius,_angle,_angle0);
 392   updateBounds();
 393 }
 394
 395 /*!
 396  * @param angle0 in ]-Pi;Pi[
 397  * @param deltaAngle in ]-2.*Pi;2.*Pi[
 398  */
 399 EdgeArcCircle::EdgeArcCircle(Node *start, Node *end, const double *center, double radius, double angle0, double deltaAngle, bool direction):Edge(start,end,direction),_angle(deltaAngle),
 400                                                                                                                                             _angle0(angle0),_radius(radius)
 401 {
 402   _center[0]=center[0];
 403   _center[1]=center[1];
 404   updateBounds();
 405 }
 406
 407 void EdgeArcCircle::changeMiddle(Node *newMiddle)
 408 {
 409   getArcOfCirclePassingThru(*_start,*newMiddle,*_end,_center,_radius,_angle,_angle0);
 410   updateBounds();
 411 }
 412
 413 Edge *EdgeArcCircle::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) const
 414 {
 415   double sx=((*start)[0]-_center[0])/_radius;
 416   double sy=((*start)[1]-_center[1])/_radius;
 417   double ex=((*end)[0]-_center[0])/_radius;
 418   double ey=((*end)[1]-_center[1])/_radius;
 419   double angle0=getAbsoluteAngleOfNormalizedVect(direction?sx:ex,direction?sy:ey);
 420   double deltaAngle=getAbsoluteAngleOfNormalizedVect(sx*ex+sy*ey,sx*ey-sy*ex);
 421   if(deltaAngle>0. && _angle<0.)
 422     deltaAngle-=2.*M_PI;
 423   else if(deltaAngle<0. && _angle>0.)
 424     deltaAngle+=2.*M_PI;
 425   deltaAngle=direction?deltaAngle:-deltaAngle;
 426   return new EdgeArcCircle(start,end,_center,_radius,angle0,deltaAngle,direction);
 427 }
 428
 429 void EdgeArcCircle::applySimilarity(double xBary, double yBary, double dimChar)
 430 {
 431   Edge::applySimilarity(xBary,yBary,dimChar);
 432   _radius/=dimChar;
 433   _center[0]=(_center[0]-xBary)/dimChar;
 434   _center[1]=(_center[1]-yBary)/dimChar;
 435 }
 436
 437 /*!
 438  * Given an \b NON normalized vector 'vect', returns its norm 'normVect' and its
 439  * angle in ]-Pi,Pi] relative to Ox axe.
 440  */
 441 double EdgeArcCircle::getAbsoluteAngle(const double *vect, double& normVect)
 442 {
 443   normVect=Node::norm(vect);
 444   return getAbsoluteAngleOfNormalizedVect(vect[0]/normVect,vect[1]/normVect);
 445 }
 446
 447 /*!
 448  * Given a \b normalized vector defined by (ux,uy) returns its angle in ]-Pi;Pi].
 449  * So before using this method ux*ux+uy*uy should as much as possible close to 1.
 450  * This methods is quite time consuming in order to keep as much as possible precision.
 451  * It is NOT ALWAYS possible to do that only in one call of acos. Sometimes call to asin is necessary
 452  * due to imperfection of acos near 0. and Pi (cos x ~ 1-x*x/2.)
 453  */
 454 double EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(double ux, double uy)
 455 {
 456   //When arc is lower than 0.707 Using Asin
 457   if(fabs(ux)<0.707)
 458     {
 459       double ret=safeAcos(ux);
 460       if(uy>0.)
 461         return ret;
 462       ret=-ret;
 463       return ret;
 464     }
 465   else
 466     {
 467       double ret=safeAsin(uy);
 468       if(ux>0.)
 469         return ret;
 470       if(ret>0.)
 471         return M_PI-ret;
 472       else
 473         return -M_PI-ret;
 474     }
 475 }
 476
 477 void EdgeArcCircle::getArcOfCirclePassingThru(const double *start, const double *middle, const double *end,
 478                                               double *center, double& radius, double& angleInRad, double& angleInRad0)
 479 {
 480   double delta=(middle[0]-start[0])*(end[1]-middle[1])-(end[0]-middle[0])*(middle[1]-start[1]);
 481   double b1=(middle[1]*middle[1]+middle[0]*middle[0]-start[0]*start[0]-start[1]*start[1])/2;
 482   double b2=(end[1]*end[1]+end[0]*end[0]-middle[0]*middle[0]-middle[1]*middle[1])/2;
 483   center[0]=((end[1]-middle[1])*b1+(start[1]-middle[1])*b2)/delta;
 484   center[1]=((middle[0]-end[0])*b1+(middle[0]-start[0])*b2)/delta;
 485   radius=safeSqrt((start[0]-center[0])*(start[0]-center[0])+(start[1]-center[1])*(start[1]-center[1]));
 486   angleInRad0=getAbsoluteAngleOfNormalizedVect((start[0]-center[0])/radius,(start[1]-center[1])/radius);
 487   double angleInRadM=getAbsoluteAngleOfNormalizedVect((middle[0]-center[0])/radius,(middle[1]-center[1])/radius);
 488   angleInRad=getAbsoluteAngleOfNormalizedVect(((start[0]-center[0])*(end[0]-center[0])+(start[1]-center[1])*(end[1]-center[1]))/(radius*radius),
 489                                               ((start[0]-center[0])*(end[1]-center[1])-(start[1]-center[1])*(end[0]-center[0]))/(radius*radius));
 490   if(isAngleNotIn(angleInRad0,angleInRad,angleInRadM))
 491     angleInRad=angleInRad<0?2*M_PI+angleInRad:angleInRad-2*M_PI;
 492 }
 493
 494 void EdgeArcCircle::dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const
 495 {
 496   stream << "5 1 0 1 ";
 497   fillXfigStreamForLoc(stream);
 498   stream << " 7 50 -1 -1 0.000 0 ";
 499   if( (direction && _angle>=0) || (!direction && _angle<0))
 500     stream << '0';//'0'
 501   else
 502     stream << '1';//'1'
 503   stream << " 0 0 ";
 504   stream << box.fitXForXFigD(_center[0],resolution) << " " << box.fitYForXFigD(_center[1],resolution) << " ";
 505   direction?_start->dumpInXfigFile(stream,resolution,box):_end->dumpInXfigFile(stream,resolution,box);
 506   Node *middle=buildRepresentantOfMySelf();
 507   middle->dumpInXfigFile(stream,resolution,box);
 508   middle->decrRef();
 509   direction?_end->dumpInXfigFile(stream,resolution,box):_start->dumpInXfigFile(stream,resolution,box);
 510   stream << endl;
 511 }
 512
 513 void EdgeArcCircle::update(Node *m)
 514 {
 515   getArcOfCirclePassingThru(*_start,*m,*_end,_center,_radius,_angle,_angle0);
 516   updateBounds();
 517 }
 518
 519 /*!
 520  * This methods computes :
 521  * \f[
 522  * \int_{Current Edge} -ydx
 523  * \f]
 524  */
 525 double EdgeArcCircle::getAreaOfZone() const
 526 {
 527   return -_radius*_radius*(sin(_angle)-_angle)/2.+((*_start)[0]-(*_end)[0])*((*_start)[1]+(*_end)[1])/2.;
 528 }
 529
 530 double EdgeArcCircle::getCurveLength() const
 531 {
 532   return fabs(_angle*_radius);
 533 }
 534
 535 void EdgeArcCircle::getBarycenter(double *bary) const
 536 {
 537   bary[0]=_center[0]+_radius*cos(_angle0+_angle/2.);
 538   bary[1]=_center[1]+_radius*sin(_angle0+_angle/2.);
 539 }
 540
 541 /*!
 542  * \f[
 543  * bary[0]=\int_{Current Edge} -yxdx
 544  * \f]
 545  * \f[
 546  * bary[1]=\int_{Current Edge} -\frac{y^{2}}{2}dx
 547  * \f]
 548  * To compute these 2 expressions in this class we have :
 549  * \f[
 550  * x=x_{0}+Radius \cdot cos(\theta)
 551  * \f]
 552  * \f[
 553  * y=y_{0}+Radius \cdot sin(\theta)
 554  * \f]
 555  * \f[
 556  * dx=-Radius \cdot sin(\theta) \cdot d\theta
 557  * \f]
 558  */
 559 void EdgeArcCircle::getBarycenterOfZone(double *bary) const
 560 {
 561   double x0=_center[0];
 562   double y0=_center[1];
 563   double angle1=_angle0+_angle;
 564   double tmp1=sin(angle1);
 565   double tmp0=sin(_angle0);
 566   double tmp2=_radius*_radius*_radius;
 567   double tmp3=cos(angle1);
 568   double tmp4=cos(_angle0);
 569   bary[0]=_radius*x0*y0*(tmp4-tmp3)+_radius*_radius*(y0*(cos(2*_angle0)-cos(2*angle1))/4.+
 570                                                      x0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.))
 571     +tmp2*(tmp1*tmp1*tmp1-tmp0*tmp0*tmp0)/3.;
 572   bary[1]=y0*y0*_radius*(tmp4-tmp3)/2.+_radius*_radius*y0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.)
 573     +tmp2*(tmp4-tmp3+(tmp3*tmp3*tmp3-tmp4*tmp4*tmp4)/3.)/2.;
 574 }
 575
 576 /*!
 577  * Characteristic value used is angle in ]_Pi;Pi[ from axe 0x.
 578  */
 579 bool EdgeArcCircle::isIn(double characterVal) const
 580 {
 581   return isIn2Pi(_angle0,_angle,characterVal);
 582 }
 583
 584 Node *EdgeArcCircle::buildRepresentantOfMySelf() const
 585 {
 586   return new Node(_center[0]+_radius*cos(_angle0+_angle/2.),_center[1]+_radius*sin(_angle0+_angle/2.));
 587 }
 588
 589 /*!
 590  * Characteristic value used is angle in ]_Pi;Pi[ from axe 0x.
 591  * 'val1' and 'val2' have been detected previously as owning to this.
 592  */
 593 bool EdgeArcCircle::isLower(double val1, double val2) const
 594 {
 595   double myDelta1=val1-_angle0;
 596   double myDelta2=val2-_angle0;
 597   if(_angle>0.)
 598     {
 599       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.
 600       myDelta2=myDelta2>-(_radius*QUADRATIC_PLANAR::_precision)?myDelta2:myDelta2+2.*M_PI;
 601       return myDelta1<myDelta2;
 602     }
 603   else
 604     {
 605       myDelta1=myDelta1<(_radius*QUADRATIC_PLANAR::_precision)?myDelta1:myDelta1-2.*M_PI;
 606       myDelta2=myDelta2<(_radius*QUADRATIC_PLANAR::_precision)?myDelta2:myDelta2-2.*M_PI;
 607       return myDelta2<myDelta1;
 608     }
 609 }
 610
 611 /*!
 612  * For Arc circle the caract value is angle with Ox between -Pi and Pi.
 613  */
 614 double EdgeArcCircle::getCharactValue(const Node& node) const
 615 {
 616   double dx=(node[0]-_center[0])/_radius;
 617   double dy=(node[1]-_center[1])/_radius;
 618   return getAbsoluteAngleOfNormalizedVect(dx,dy);
 619 }
 620
 621 double EdgeArcCircle::getDistanceToPoint(const double *pt) const
 622 {
 623   double angle=Node::computeAngle(_center,pt);
 624   if(isIn2Pi(_angle0,_angle,angle))
 625     return fabs(Node::distanceBtw2Pt(_center,pt)-_radius);
 626   else
 627     {
 628       double dist1=Node::distanceBtw2Pt(*_start,pt);
 629       double dist2=Node::distanceBtw2Pt(*_end,pt);
 630       return std::min(dist1,dist2);
 631     }
 632 }
 633
 634 bool EdgeArcCircle::isNodeLyingOn(const double *coordOfNode) const
 635 {
 636   double dist=Node::distanceBtw2Pt(_center,coordOfNode);
 637   if(Node::areDoubleEquals(dist,_radius))
 638     {
 639       double angle=Node::computeAngle(_center,coordOfNode);
 640       return isIn2Pi(_angle0,_angle,angle);
 641     }
 642   else
 643     return false;
 644 }
 645
 646 /*!
 647  * Idem isAngleNotIn except that here 'start' in ]-Pi;Pi[ and delta in ]-2*Pi;2Pi[.
 648  * @param angleIn in ]-Pi;Pi[.
 649  */
 650 bool EdgeArcCircle::isIn2Pi(double start, double delta, double angleIn)
 651 {
 652   double myDelta=angleIn-start;
 653   if(delta>0.)
 654     {
 655       myDelta=myDelta>=0.?myDelta:myDelta+2.*M_PI;
 656       return myDelta>0. && myDelta<delta;
 657     }
 658   else
 659     {
 660       myDelta=myDelta<=0.?myDelta:myDelta-2.*M_PI;
 661       return myDelta<0. && myDelta>delta;
 662     }
 663 }
 664
 665 /*!
 666  * 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'.
 667  */
 668 bool EdgeArcCircle::isAngleNotIn(double start, double delta, double angleIn)
 669 {
 670   double tmp=start;
 671   if(tmp<0.)
 672     tmp+=2*M_PI;
 673   double tmp2=angleIn;
 674   if(tmp2<0.)
 675     tmp2+=2*M_PI;
 676   if(tmp+delta>=2.*M_PI)
 677     return (tmp2<tmp) && (tmp2>tmp+delta-2*M_PI);
 678   else if(tmp+delta>=0.)
 679     return (tmp2<std::min(tmp,tmp+delta) || tmp2>std::max(tmp,tmp+delta));
 680   else
 681     return (tmp2>tmp) && (tmp2<(tmp+delta+2.*M_PI));
 682 }
 683
 684 void EdgeArcCircle::updateBounds()
 685 {
 686   _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]));
 687   if(isIn2Pi(_angle0,_angle,M_PI/2))
 688     _bounds[3]=_center[1]+_radius;
 689   if(isIn2Pi(_angle0,_angle,-M_PI/2))
 690     _bounds[2]=_center[1]-_radius;
 691   if(isIn2Pi(_angle0,_angle,0.))
 692     _bounds[1]=_center[0]+_radius;
 693   if(isIn2Pi(_angle0,_angle,M_PI))
 694     _bounds[0]=_center[0]-_radius;
 695 }