1 // Copyright (C) 2007-2008 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 #include "EdgeArcCircle.hxx"
20 #include "EdgeLin.hxx"
21 #include "InterpKernelException.hxx"
28 using namespace INTERP_KERNEL;
30 ArcCArcCIntersector::ArcCArcCIntersector(const EdgeArcCircle& e1, const EdgeArcCircle& e2):SameTypeEdgeIntersector(e1,e2),_dist(0.)
34 bool ArcCArcCIntersector::haveTheySameDirection() const
36 return (getE1().getAngle()>0. && getE2().getAngle()>0.) || (getE1().getAngle()<0. && getE2().getAngle()<0.);
40 * Precondition 'start' and 'end' are on the same curve than this.
42 void ArcCArcCIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
44 bool obvious1,obvious2;
45 obviousCaseForCurvAbscisse(start,whereStart,commonNode,obvious1);
46 obviousCaseForCurvAbscisse(end,whereEnd,commonNode,obvious2);
47 if(obvious1 && obvious2)
49 double angleInRadStart=getAngle(start);
50 double angleInRadEnd=getAngle(end);
51 if(obvious1 || obvious2)
55 if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
63 if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart))
66 whereStart=OUT_BEFORE;
70 if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart))
73 if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
79 {//we are out in start.
80 if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
82 whereStart=OUT_BEFORE;
87 if(EdgeArcCircle::isIn2Pi(getE2().getAngle0(),getE2().getAngle(),getE1().getAngle0()))
88 {//_e2 contains stictly _e1
89 whereStart=OUT_BEFORE;
93 {//_e2 is outside from _e1
94 whereStart=OUT_BEFORE;
102 * Return angle between ]-Pi;Pi[
104 double ArcCArcCIntersector::getAngle(Node *node) const
106 return EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(((*node)[0]-getE1().getCenter()[0])/getE1().getRadius(),((*node)[1]-getE1().getCenter()[1])/getE1().getRadius());
109 bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeArcCircle& a2)
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());
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();
122 a2.getCenter(centerL); radiusL=a2.getRadius(); angle0L=a2.getAngle0(); angleL=a2.getAngle();
123 a1.getCenter(centerB); radiusB=a1.getRadius();
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.))
132 Bounds *merge=a1.getBounds().nearlyAmIIntersectingWith(a2.getBounds());
133 merge->getInterceptedArc(centerL,radiusL,angle0L,angleL);
137 if(Node::areDoubleEqualsWP(tmp,0.,1/(10*std::max(radiusL,radiusB))))
139 if(Node::areDoubleEquals(radiusL,radiusB))
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.);
160 void ArcCArcCIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped)
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))
166 obviousNoIntersection=true;
170 if(areArcsOverlapped(getE1(),getE2()))//(Node::areDoubleEquals(_dist,0.) && Node::areDoubleEquals(radius1,radius2))
172 obviousNoIntersection=false;
177 obviousNoIntersection=false;
182 std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristicVal() const
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))
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);
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);
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);
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()));
224 //tangent intersection
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()));
240 double signDeltaAngle2;
248 angle0_2=angle0_1+M_PI;
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)) )
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
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);
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()));
292 //tangent intersection
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()));
303 ArcCSegIntersector::ArcCSegIntersector(const EdgeArcCircle& e1, const EdgeLin& e2, bool reverse):CrossTypeEdgeIntersector(e1,e2,reverse)
307 void ArcCSegIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped)
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;
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;
321 obviousNoIntersection=true;
324 void ArcCSegIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
326 throw Exception("Internal error. Should never been called : no overlapping possible between arc of circle and a segment.");
329 std::list< IntersectElement > ArcCSegIntersector::getIntersectionsCharacteristicVal() const
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))
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()));
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()));
354 else//tangent intersection
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()));
368 EdgeArcCircle::EdgeArcCircle(std::istream& lineInXfig)
370 const unsigned NB_OF_SKIP_FIELDS=15;
372 for(unsigned i=0;i<NB_OF_SKIP_FIELDS;i++)
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);
382 EdgeArcCircle::EdgeArcCircle(Node *start, Node *middle, Node *end, bool direction):Edge(start,end, direction)
384 getArcOfCirclePassingThru(*_start,*middle,*_end,_center,_radius,_angle,_angle0);
388 EdgeArcCircle::EdgeArcCircle(double sX, double sY, double mX, double mY, double eX, double eY):Edge(sX,sY,eX,eY)
390 double middle[2]; middle[0]=mX; middle[1]=mY;
391 getArcOfCirclePassingThru(*_start,middle,*_end,_center,_radius,_angle,_angle0);
396 * @param angle0 in ]-Pi;Pi[
397 * @param deltaAngle in ]-2.*Pi;2.*Pi[
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)
402 _center[0]=center[0];
403 _center[1]=center[1];
407 void EdgeArcCircle::changeMiddle(Node *newMiddle)
409 getArcOfCirclePassingThru(*_start,*newMiddle,*_end,_center,_radius,_angle,_angle0);
413 Edge *EdgeArcCircle::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) const
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.)
423 else if(deltaAngle<0. && _angle>0.)
425 deltaAngle=direction?deltaAngle:-deltaAngle;
426 return new EdgeArcCircle(start,end,_center,_radius,angle0,deltaAngle,direction);
429 void EdgeArcCircle::applySimilarity(double xBary, double yBary, double dimChar)
431 Edge::applySimilarity(xBary,yBary,dimChar);
433 _center[0]=(_center[0]-xBary)/dimChar;
434 _center[1]=(_center[1]-yBary)/dimChar;
438 * Given an \b NON normalized vector 'vect', returns its norm 'normVect' and its
439 * angle in ]-Pi,Pi] relative to Ox axe.
441 double EdgeArcCircle::getAbsoluteAngle(const double *vect, double& normVect)
443 normVect=Node::norm(vect);
444 return getAbsoluteAngleOfNormalizedVect(vect[0]/normVect,vect[1]/normVect);
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.)
454 double EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(double ux, double uy)
456 //When arc is lower than 0.707 Using Asin
459 double ret=safeAcos(ux);
467 double ret=safeAsin(uy);
477 void EdgeArcCircle::getArcOfCirclePassingThru(const double *start, const double *middle, const double *end,
478 double *center, double& radius, double& angleInRad, double& angleInRad0)
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;
494 void EdgeArcCircle::dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const
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))
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);
509 direction?_end->dumpInXfigFile(stream,resolution,box):_start->dumpInXfigFile(stream,resolution,box);
513 void EdgeArcCircle::update(Node *m)
515 getArcOfCirclePassingThru(*_start,*m,*_end,_center,_radius,_angle,_angle0);
520 * This methods computes :
522 * \int_{Current Edge} -ydx
525 double EdgeArcCircle::getAreaOfZone() const
527 return -_radius*_radius*(sin(_angle)-_angle)/2.+((*_start)[0]-(*_end)[0])*((*_start)[1]+(*_end)[1])/2.;
530 double EdgeArcCircle::getCurveLength() const
532 return fabs(_angle*_radius);
535 void EdgeArcCircle::getBarycenter(double *bary) const
537 bary[0]=_center[0]+_radius*cos(_angle0+_angle/2.);
538 bary[1]=_center[1]+_radius*sin(_angle0+_angle/2.);
543 * bary[0]=\int_{Current Edge} -yxdx
546 * bary[1]=\int_{Current Edge} -\frac{y^{2}}{2}dx
548 * To compute these 2 expressions in this class we have :
550 * x=x_{0}+Radius \cdot cos(\theta)
553 * y=y_{0}+Radius \cdot sin(\theta)
556 * dx=-Radius \cdot sin(\theta) \cdot d\theta
559 void EdgeArcCircle::getBarycenterOfZone(double *bary) const
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.;
577 * Characteristic value used is angle in ]_Pi;Pi[ from axe 0x.
579 bool EdgeArcCircle::isIn(double characterVal) const
581 return isIn2Pi(_angle0,_angle,characterVal);
584 Node *EdgeArcCircle::buildRepresentantOfMySelf() const
586 return new Node(_center[0]+_radius*cos(_angle0+_angle/2.),_center[1]+_radius*sin(_angle0+_angle/2.));
590 * Characteristic value used is angle in ]_Pi;Pi[ from axe 0x.
591 * 'val1' and 'val2' have been detected previously as owning to this.
593 bool EdgeArcCircle::isLower(double val1, double val2) const
595 double myDelta1=val1-_angle0;
596 double myDelta2=val2-_angle0;
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;
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;
612 * For Arc circle the caract value is angle with Ox between -Pi and Pi.
614 double EdgeArcCircle::getCharactValue(const Node& node) const
616 double dx=(node[0]-_center[0])/_radius;
617 double dy=(node[1]-_center[1])/_radius;
618 return getAbsoluteAngleOfNormalizedVect(dx,dy);
621 double EdgeArcCircle::getDistanceToPoint(const double *pt) const
623 double angle=Node::computeAngle(_center,pt);
624 if(isIn2Pi(_angle0,_angle,angle))
625 return fabs(Node::distanceBtw2Pt(_center,pt)-_radius);
628 double dist1=Node::distanceBtw2Pt(*_start,pt);
629 double dist2=Node::distanceBtw2Pt(*_end,pt);
630 return std::min(dist1,dist2);
634 bool EdgeArcCircle::isNodeLyingOn(const double *coordOfNode) const
636 double dist=Node::distanceBtw2Pt(_center,coordOfNode);
637 if(Node::areDoubleEquals(dist,_radius))
639 double angle=Node::computeAngle(_center,coordOfNode);
640 return isIn2Pi(_angle0,_angle,angle);
647 * Idem isAngleNotIn except that here 'start' in ]-Pi;Pi[ and delta in ]-2*Pi;2Pi[.
648 * @param angleIn in ]-Pi;Pi[.
650 bool EdgeArcCircle::isIn2Pi(double start, double delta, double angleIn)
652 double myDelta=angleIn-start;
655 myDelta=myDelta>=0.?myDelta:myDelta+2.*M_PI;
656 return myDelta>0. && myDelta<delta;
660 myDelta=myDelta<=0.?myDelta:myDelta-2.*M_PI;
661 return myDelta<0. && myDelta>delta;
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'.
668 bool EdgeArcCircle::isAngleNotIn(double start, double delta, double angleIn)
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));
681 return (tmp2>tmp) && (tmp2<(tmp+delta+2.*M_PI));
684 void EdgeArcCircle::updateBounds()
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;