src/ChFi2d/ChFi2d_AnaFilletAlgo.cxx

   1 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   2 //
   3 // This file is part of Open CASCADE Technology software library.
   4 //
   5 // This library is free software; you can redistribute it and/or modify it under
   6 // the terms of the GNU Lesser General Public License version 2.1 as published
   7 // by the Free Software Foundation, with special exception defined in the file
   8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
   9 // distribution for complete text of the license and disclaimer of any warranty.
  10 //
  11 // Alternatively, this file may be used under the terms of Open CASCADE
  12 // commercial license or contractual agreement.
  13
  14 #include <ChFi2d_AnaFilletAlgo.hxx>
  15
  16 #include <gp_Ax3.hxx>
  17
  18 #include <Standard_TypeMismatch.hxx>
  19
  20 #include <BRepBuilderAPI_MakeEdge.hxx>
  21 #include <BRepBuilderAPI_MakeWire.hxx>
  22 #include <BRepBuilderAPI_MakeFace.hxx>
  23
  24 #include <GeomAPI_ExtremaCurveCurve.hxx>
  25 #include <IntAna2d_AnaIntersection.hxx>
  26 #include <ShapeAnalysis_Wire.hxx>
  27 #include <Geom_Circle.hxx>
  28
  29 #include <BRepAdaptor_Curve.hxx>
  30 #include <BRep_Tool.hxx>
  31
  32 #include <TopoDS.hxx>
  33 #include <TopoDS_Iterator.hxx>
  34
  35 #include <ProjLib.hxx>
  36 #include <TopExp.hxx>
  37 #include <ElSLib.hxx>
  38
  39 // Compute the flag: CW || CCW
  40 static Standard_Boolean isCW(const BRepAdaptor_Curve& AC)
  41 {
  42   const Standard_Real f = AC.FirstParameter();
  43   const Standard_Real l = AC.LastParameter();
  44   Handle(Geom_Circle) circle = Handle(Geom_Circle)::DownCast(AC.Curve().Curve());
  45   gp_Pnt start = AC.Value(f);
  46   gp_Pnt end = AC.Value(l);
  47   gp_Pnt center = AC.Circle().Location();
  48   gp_Ax3 plane = AC.Circle().Position();
  49
  50   // Get point on circle at half angle
  51   gp_Pnt m;
  52   circle->D0(0.5 * (f + l), m);
  53
  54   // Compare angles between vectors to middle point and to the end point.
  55   gp_Vec startv(center, start), endv(center, end), middlev(center, m);
  56   double middlea = startv.AngleWithRef(middlev, plane.Direction());
  57   while(middlea < 0.0)
  58     middlea += 2.0 * M_PI;
  59   double enda = startv.AngleWithRef(endv, plane.Direction());
  60   while(enda < 0.0)
  61     enda += 2.0 * M_PI;
  62
  63   Standard_Boolean is_cw = middlea > enda ? Standard_True : Standard_False;
  64   return is_cw;
  65 }
  66
  67 // Equality of points computed through square distance between the points.
  68 static Standard_Boolean IsEqual(const gp_Pnt& p1, const gp_Pnt& p2)
  69 {
  70   return p1.SquareDistance(p2) < Precision::SquareConfusion();
  71 }
  72 static Standard_Boolean IsEqual(const gp_Pnt2d& p1, const gp_Pnt2d& p2)
  73 {
  74   return p1.SquareDistance(p2) < Precision::SquareConfusion();
  75 }
  76
  77 // An empty constructor.
  78 // Use the method Init() to initialize the class.
  79 ChFi2d_AnaFilletAlgo::ChFi2d_AnaFilletAlgo()
  80 : segment1(Standard_False),
  81   x11(0.0),
  82   y11(0.0),
  83   x12(0.0),
  84   y12(0.0),
  85   xc1(0.0),
  86   yc1(0.0),
  87   radius1(0.0),
  88   cw1(Standard_False),
  89   segment2(Standard_False),
  90   x21(0.0),
  91   y21(0.0),
  92   x22(0.0),
  93   y22(0.0),
  94   xc2(0.0),
  95   yc2(0.0),
  96   radius2(0.0),
  97   cw2(Standard_False)
  98 {
  99 }
 100
 101 // An constructor.
 102 // It expects two edges having a common point of type:
 103 // - segment
 104 // - arc of circle.
 105 ChFi2d_AnaFilletAlgo::ChFi2d_AnaFilletAlgo(const TopoDS_Wire& theWire,
 106                                            const gp_Pln&      thePlane)
 107 : plane(thePlane),
 108   segment1(Standard_False),
 109   x11(0.0),
 110   y11(0.0),
 111   x12(0.0),
 112   y12(0.0),
 113   xc1(0.0),
 114   yc1(0.0),
 115   radius1(0.0),
 116   cw1(Standard_False),
 117   segment2(Standard_False),
 118   x21(0.0),
 119   y21(0.0),
 120   x22(0.0),
 121   y22(0.0),
 122   xc2(0.0),
 123   yc2(0.0),
 124   radius2(0.0),
 125   cw2(Standard_False)
 126 {
 127   Init(theWire, thePlane);
 128 }
 129
 130 // A constructor.
 131 // It expects two edges having a common point of type:
 132 // - segment
 133 // - arc of circle.
 134 ChFi2d_AnaFilletAlgo::ChFi2d_AnaFilletAlgo(const TopoDS_Edge& theEdge1,
 135                                            const TopoDS_Edge& theEdge2,
 136                                            const gp_Pln& thePlane)
 137 : plane(thePlane),
 138   segment1(Standard_False),
 139   x11(0.0),
 140   y11(0.0),
 141   x12(0.0),
 142   y12(0.0),
 143   xc1(0.0),
 144   yc1(0.0),
 145   radius1(0.0),
 146   cw1(Standard_False),
 147   segment2(Standard_False),
 148   x21(0.0),
 149   y21(0.0),
 150   x22(0.0),
 151   y22(0.0),
 152   xc2(0.0),
 153   yc2(0.0),
 154   radius2(0.0),
 155   cw2(Standard_False)
 156 {
 157   // Make a wire consisting of two edges.
 158   Init(theEdge1, theEdge2, thePlane);
 159 }
 160
 161 // Initializes the class by a wire consisting of two edges.
 162 void ChFi2d_AnaFilletAlgo::Init(const TopoDS_Wire& theWire, const gp_Pln& thePlane)
 163 {
 164   plane = thePlane;
 165   TopoDS_Iterator itr(theWire);
 166   for (; itr.More(); itr.Next())
 167   {
 168     if (e1.IsNull())
 169       e1 = TopoDS::Edge(itr.Value());
 170     else if (e2.IsNull())
 171       e2 = TopoDS::Edge(itr.Value());
 172   }
 173   if (e1.IsNull() || e2.IsNull())
 174     throw Standard_TypeMismatch("The algorithm expects a wire consisting of two linear or circular edges.");
 175
 176   // Left neighbour.
 177   BRepAdaptor_Curve AC1(e1);
 178   if (AC1.GetType() != GeomAbs_Line && AC1.GetType() != GeomAbs_Circle)
 179     throw Standard_TypeMismatch("A segment or an arc of circle is expected.");
 180
 181   TopoDS_Vertex v1, v2;
 182   TopExp::Vertices(e1, v1, v2, Standard_True);
 183   if (v1.IsNull() || v2.IsNull())
 184     throw Standard_Failure("An infinite edge.");
 185
 186   gp_Pnt P1 = BRep_Tool::Pnt(v1);
 187   gp_Pnt P2 = BRep_Tool::Pnt(v2);
 188   gp_Pnt2d p1 = ProjLib::Project(thePlane, P1);
 189   gp_Pnt2d p2 = ProjLib::Project(thePlane, P2);
 190   p1.Coord(x11, y11);
 191   p2.Coord(x12, y12);
 192
 193   segment1 = true;
 194   if (AC1.GetType() == GeomAbs_Circle)
 195   {
 196     segment1 = false;
 197     gp_Circ c = AC1.Circle();
 198
 199     gp_Pnt2d loc = ProjLib::Project(thePlane, c.Location());
 200     loc.Coord(xc1, yc1);
 201
 202     radius1 = c.Radius();
 203     cw1 = isCW(AC1);
 204   }
 205
 206   // Right neighbour.
 207   BRepAdaptor_Curve AC2(e2);
 208   if (AC2.GetType() != GeomAbs_Line && AC2.GetType() != GeomAbs_Circle)
 209     throw Standard_TypeMismatch("A segment or an arc of circle is expected.");
 210
 211   TopExp::Vertices(e2, v1, v2, Standard_True);
 212   if (v1.IsNull() || v2.IsNull())
 213     throw Standard_Failure("An infinite edge.");
 214
 215   P1 = BRep_Tool::Pnt(v1);
 216   P2 = BRep_Tool::Pnt(v2);
 217   p1 = ProjLib::Project(thePlane, P1);
 218   p2 = ProjLib::Project(thePlane, P2);
 219   p1.Coord(x21, y21);
 220   p2.Coord(x22, y22);
 221
 222   segment2 = true;
 223   if (AC2.GetType() == GeomAbs_Circle)
 224   {
 225     segment2 = false;
 226     gp_Circ c = AC2.Circle();
 227
 228     gp_Pnt2d loc = ProjLib::Project(thePlane, c.Location());
 229     loc.Coord(xc2, yc2);
 230
 231     radius2 = c.Radius();
 232     cw2 = isCW(AC2);
 233   }
 234 }
 235
 236 // Initializes the class by two edges.
 237 void ChFi2d_AnaFilletAlgo::Init(const TopoDS_Edge& theEdge1, const TopoDS_Edge& theEdge2,
 238                                 const gp_Pln& thePlane)
 239 {
 240   // Make a wire consisting of two edges.
 241
 242   // Get common point.
 243   TopoDS_Vertex v11, v12, v21, v22;
 244   TopExp::Vertices(theEdge1, v11, v12, Standard_True);
 245   TopExp::Vertices(theEdge2, v21, v22, Standard_True);
 246   if (v11.IsNull() || v12.IsNull() || v21.IsNull() || v22.IsNull())
 247     throw Standard_Failure("An infinite edge.");
 248
 249   gp_Pnt p11 = BRep_Tool::Pnt(v11);
 250   gp_Pnt p12 = BRep_Tool::Pnt(v12);
 251   gp_Pnt p21 = BRep_Tool::Pnt(v21);
 252   gp_Pnt p22 = BRep_Tool::Pnt(v22);
 253
 254   gp_Pnt pcommon;
 255   if (IsEqual(p11, p21) || IsEqual(p11, p22))
 256   {
 257     pcommon = p11;
 258   }
 259   else if (IsEqual(p12, p21) || IsEqual(p12, p22))
 260   {
 261     pcommon = p12;
 262   }
 263   else
 264     throw Standard_Failure("The edges have no common point.");
 265
 266   // Reverse the edges in case of need (to construct a wire).
 267   Standard_Boolean is1stReversed(Standard_False), is2ndReversed(Standard_False);
 268   if (IsEqual(pcommon, p11))
 269     is1stReversed = Standard_True;
 270   else if (IsEqual(pcommon, p22))
 271     is2ndReversed = Standard_True;
 272
 273   // Make a wire.
 274   BRepBuilderAPI_MakeWire mkWire;
 275   if (is1stReversed)
 276     mkWire.Add(TopoDS::Edge(theEdge1.Reversed()));
 277   else
 278     mkWire.Add(theEdge1);
 279   if (is2ndReversed)
 280     mkWire.Add(TopoDS::Edge(theEdge2.Reversed()));
 281   else
 282     mkWire.Add(theEdge2);
 283   if (!mkWire.IsDone())
 284     throw Standard_Failure("Can't make a wire.");
 285
 286   const TopoDS_Wire& W = mkWire.Wire();
 287   Init(W, thePlane);
 288 }
 289
 290 // Calculates a fillet.
 291 Standard_Boolean ChFi2d_AnaFilletAlgo::Perform(const Standard_Real radius)
 292 {
 293   Standard_Boolean bRet(false);
 294   if (e1.IsNull() || e2.IsNull() ||
 295       radius < Precision::Confusion())
 296   {
 297     return bRet;
 298   }
 299
 300   // Fillet definition.
 301   Standard_Real xc = 0.0, yc = 0.0;
 302   Standard_Real start = 0.0, end = 0.0;             // parameters on neighbours
 303   Standard_Real xstart = DBL_MAX, ystart = DBL_MAX; // point on left neighbour
 304   Standard_Real xend = DBL_MAX, yend = DBL_MAX;     // point on right neighbour
 305   Standard_Boolean cw = Standard_False;
 306
 307   // Analytical algorithm works for non-intersecting arcs only.
 308   // Check arcs on self-intersection.
 309   Standard_Boolean isCut(Standard_False);
 310   if (!segment1 || !segment2)
 311   {
 312     BRepBuilderAPI_MakeWire mkWire(e1, e2);
 313     if (mkWire.IsDone())
 314     {
 315       const TopoDS_Wire& W = mkWire.Wire();
 316       BRepBuilderAPI_MakeFace mkFace(plane);
 317       if (mkFace.IsDone())
 318       {
 319         const TopoDS_Face& F = mkFace.Face();
 320         ShapeAnalysis_Wire analyzer(W, F, Precision::Confusion());
 321         if (analyzer.CheckSelfIntersection() == Standard_True)
 322         {
 323           // Cut the edges at the point of intersection.
 324           isCut = Standard_True;
 325           if (!Cut(plane, e1, e2))
 326           {
 327             return Standard_False;
 328           }
 329         }
 330       }
 331     }
 332   }// a case of segment - segment
 333
 334   // Choose the case.
 335   BRepAdaptor_Curve AC1(e1), AC2(e2);
 336   if (segment1 && segment2)
 337   {
 338     bRet = SegmentFilletSegment(radius, xc, yc, cw, start, end);
 339   }
 340   else if (segment1 && !segment2)
 341   {
 342     bRet = SegmentFilletArc(radius, xc, yc, cw, start, end, xend, yend);
 343   }
 344   else if (!segment1 && segment2)
 345   {
 346     bRet = ArcFilletSegment(radius, xc, yc, cw, start, end, xstart, ystart);
 347   }
 348   else if (!segment1 && !segment2)
 349   {
 350     bRet = ArcFilletArc(radius, xc, yc, cw, start, end);
 351   }
 352
 353   if (!bRet)
 354     return Standard_False;
 355
 356   // Invert the fillet for left-handed plane.
 357   if (plane.Position().Direct() == Standard_False)
 358     cw = !cw;
 359
 360   // Construct a fillet.
 361   // Make circle.
 362   gp_Pnt center = ElSLib::Value(xc, yc, plane);
 363   const gp_Dir& normal = plane.Position().Direction();
 364   gp_Circ circ(gp_Ax2(center, cw ? -normal : normal), radius);
 365
 366   // Fillet may only shrink a neighbour edge, it can't prolongate it.
 367   const Standard_Real delta1 = AC1.LastParameter() - AC1.FirstParameter();
 368   const Standard_Real delta2 = AC2.LastParameter() - AC2.FirstParameter();
 369   if (!isCut && (start > delta1 || end > delta2))
 370   {
 371     // Check a case when a neighbour edge almost disappears:
 372     // try to reduce the fillet radius for a little (1.e-5 mm).
 373     const Standard_Real little = 100.0 * Precision::Confusion();
 374     const Standard_Real d1 = fabs(start - delta1);
 375     const Standard_Real d2 = fabs(end   - delta2);
 376     if (d1 < little || d2 < little)
 377     {
 378       if (segment1 && segment2)
 379       {
 380         bRet = SegmentFilletSegment(radius - little, xc, yc, cw, start, end);
 381       }
 382       else if (segment1 && !segment2)
 383       {
 384         bRet = SegmentFilletArc(radius - little, xc, yc, cw, start, end, xend, yend);
 385       }
 386       else if (!segment1 && segment2)
 387       {
 388         bRet = ArcFilletSegment(radius - little, xc, yc, cw, start, end, xstart, ystart);
 389       }
 390       else if (!segment1 && !segment2)
 391       {
 392         bRet = ArcFilletArc(radius - little, xc, yc, cw, start, end);
 393       }
 394       if (bRet)
 395       {
 396         // Invert the fillet for left-handed planes.
 397         if (plane.Position().Direct() == Standard_False)
 398           cw = !cw;
 399
 400         // Make the circle again.
 401         center = ElSLib::Value(xc, yc, plane);
 402         circ.SetLocation(center);
 403         circ.SetRadius(radius - little);
 404       }
 405       else
 406       {
 407         return Standard_False;
 408       }
 409     }
 410     else
 411     {
 412       return Standard_False;
 413     }
 414   }
 415   if (bRet)
 416   {
 417     // start: (xstart, ystart) - pstart.
 418     gp_Pnt pstart;
 419     if (xstart != DBL_MAX)
 420     {
 421       pstart = ElSLib::Value(xstart, ystart, plane);
 422     }
 423     else
 424     {
 425       if (e1.Orientation() == TopAbs_FORWARD)
 426         pstart = AC1.Value(AC1.LastParameter() - start);
 427       else
 428         pstart = AC1.Value(AC1.FirstParameter() + start);
 429     }
 430     // end: (xend, yend) -> pend.
 431     gp_Pnt pend;
 432     if (xend != DBL_MAX)
 433     {
 434       pend = ElSLib::Value(xend, yend, plane);
 435     }
 436     else
 437     {
 438       if (e2.Orientation() == TopAbs_FORWARD)
 439         pend = AC2.Value(AC2.FirstParameter() + end);
 440       else
 441         pend = AC2.Value(AC2.LastParameter() - end);
 442     }
 443
 444     // Make arc.
 445     BRepBuilderAPI_MakeEdge mkEdge(circ, pstart, pend);
 446     bRet = mkEdge.IsDone();
 447     if (bRet)
 448     {
 449       fillet = mkEdge.Edge();
 450
 451       // Limit the neighbours.
 452       // Left neighbour.
 453       gp_Pnt p1, p2;
 454       shrinke1.Nullify();
 455       if (e1.Orientation() == TopAbs_FORWARD)
 456       {
 457         p1 = AC1.Value(AC1.FirstParameter());
 458         p2 = pstart;
 459       }
 460       else
 461       {
 462         p1 = pstart;
 463         p2 = AC1.Value(AC1.LastParameter());
 464       }
 465       if (segment1)
 466       {
 467         BRepBuilderAPI_MakeEdge mkSegment1;
 468         mkSegment1.Init(AC1.Curve().Curve(), p1, p2);
 469         if (mkSegment1.IsDone())
 470           shrinke1 = mkSegment1.Edge();
 471       }
 472       else
 473       {
 474         BRepBuilderAPI_MakeEdge mkCirc1;
 475         mkCirc1.Init(AC1.Curve().Curve(), p1, p2);
 476         if (mkCirc1.IsDone())
 477           shrinke1 = mkCirc1.Edge();
 478       }
 479
 480       // Right neighbour.
 481       shrinke2.Nullify();
 482       if (e1.Orientation() == TopAbs_FORWARD)
 483       {
 484         p1 = pend;
 485         p2 = AC2.Value(AC2.LastParameter());
 486       }
 487       else
 488       {
 489         p1 = AC2.Value(AC2.FirstParameter());
 490         p2 = pend;
 491       }
 492       if (segment2)
 493       {
 494         BRepBuilderAPI_MakeEdge mkSegment2;
 495         mkSegment2.Init(AC2.Curve().Curve(), p1, p2);
 496         if (mkSegment2.IsDone())
 497           shrinke2 = mkSegment2.Edge();
 498       }
 499       else
 500       {
 501         BRepBuilderAPI_MakeEdge mkCirc2;
 502         mkCirc2.Init(AC2.Curve().Curve(), p1, p2);
 503         if (mkCirc2.IsDone())
 504           shrinke2 = mkCirc2.Edge();
 505       }
 506
 507       bRet = !shrinke1.IsNull() && !shrinke2.IsNull();
 508     }// fillet edge is done
 509   }// shrinking is good
 510
 511   return bRet;
 512 }
 513
 514 // Retrieves a result (fillet and shrinked neighbours).
 515 const TopoDS_Edge& ChFi2d_AnaFilletAlgo::Result(TopoDS_Edge& theE1, TopoDS_Edge& theE2)
 516 {
 517   theE1 = shrinke1;
 518   theE2 = shrinke2;
 519   return fillet;
 520 }
 521
 522 // WW5 method to compute fillet.
 523 // It returns a constructed fillet definition:
 524 //     center point (xc, yc)
 525 //     point on the 1st segment (xstart, ystart)
 526 //     point on the 2nd segment (xend, yend)
 527 //     is the arc of fillet clockwise (cw = true) or counterclockwise (cw = false).
 528 Standard_Boolean ChFi2d_AnaFilletAlgo::SegmentFilletSegment(const Standard_Real radius,
 529                                                             Standard_Real& xc, Standard_Real& yc,
 530                                                             Standard_Boolean& cw,
 531                                                             Standard_Real& start, Standard_Real& end)
 532 {
 533   // Make normalized vectors at p12.
 534   gp_Pnt2d p11(x11, y11);
 535   gp_Pnt2d p12(x12, y12);
 536   gp_Pnt2d p22(x22, y22);
 537
 538   // Check length of segments.
 539   if (IsEqual(p12, p11) || IsEqual(p12, p22))
 540   {
 541     return Standard_False;
 542   }
 543
 544   // Make vectors.
 545   gp_Vec2d v1(p12, p11);
 546   gp_Vec2d v2(p12, p22);
 547   v1.Normalize();
 548   v2.Normalize();
 549
 550   // Make bisectrissa.
 551   gp_Vec2d bisec = 0.5 * (v1 + v2);
 552
 553   // Check bisectrissa.
 554   if (bisec.SquareMagnitude() < Precision::SquareConfusion())
 555     return Standard_False;
 556
 557   // Normalize the bisectrissa.
 558   bisec.Normalize();
 559
 560   // Angle at bisectrissa.
 561   Standard_Real beta = v1.Angle(bisec);
 562
 563   // Length along the bisectrissa till the center of fillet.
 564   Standard_Real L = radius / sin(fabs(beta));
 565
 566   // Center point of fillet.
 567   gp_Pnt2d pc = p12.Translated(L * bisec);
 568   pc.Coord(xc, yc);
 569
 570   // Shrinking length along segments.
 571   start = sqrt(L * L - radius * radius);
 572   end = start;
 573
 574   // Orientation of fillet.
 575   cw = beta > 0.0;
 576   return Standard_True;
 577 }
 578
 579 // A function constructs a fillet between a segment and an arc.
 580 Standard_Boolean ChFi2d_AnaFilletAlgo::SegmentFilletArc(const Standard_Real radius,
 581                                                         Standard_Real& xc, Standard_Real& yc,
 582                                                         Standard_Boolean& cw,
 583                                                         Standard_Real& start, Standard_Real& end,
 584                                                         Standard_Real& xend, Standard_Real& yend)
 585 {
 586   // Make a line parallel to the segment at the side of center point of fillet.
 587   // This side may be defined through making a bisectrissa for vectors at p12 (or p21).
 588
 589   // Make 2D points.
 590   gp_Pnt2d p12(x12, y12);
 591   gp_Pnt2d p11(x11, y11);
 592   gp_Pnt2d pc2(xc2, yc2);
 593
 594   // Check length of segment.
 595   if (p11.SquareDistance(p12) < gp::Resolution())
 596     return Standard_False;
 597
 598   // Make 2D vectors.
 599   gp_Vec2d v1(p12, p11);
 600   gp_Vec2d v2(p12, pc2);
 601
 602   // Rotate the arc vector to become tangential at p21.
 603   if (cw2)
 604     v2.Rotate(+M_PI_2);
 605   else
 606     v2.Rotate(-M_PI_2);
 607
 608   // If vectors coincide (segment and arc are tangent),
 609   // the algorithm doesn't work...
 610   Standard_Real angle = v1.Angle(v2);
 611   if (fabs(angle) < Precision::Angular())
 612     return Standard_False;
 613
 614   // Make a bissectrisa of vectors at p12.
 615   v2.Normalize();
 616   v1.Normalize();
 617   gp_Vec2d bisec = 0.5 * (v1 + v2);
 618
 619   // If segment and arc look in opposite direction,
 620   // no fillet is possible.
 621   if (bisec.SquareMagnitude() < gp::Resolution())
 622     return Standard_False;
 623
 624   // Define an appropriate point to choose center of fillet.
 625   bisec.Normalize();
 626   gp_Pnt2d nearp = p12.Translated(radius * bisec);
 627   gp_Lin2d nearl(p12, bisec);
 628
 629   // Make a line parallel to segment and
 630   // passing near the "near" point.
 631   gp_Vec2d d1(v1);
 632   gp_Lin2d line(p11, -d1);
 633   d1.Rotate(M_PI_2);
 634   line.Translate(radius * d1);
 635   if (line.Distance(nearp) > radius)
 636     line.Translate(-2.0 * radius * d1);
 637
 638   // Make a circle of radius of the arc +/- fillet radius.
 639   gp_Ax2d axes(pc2, gp::DX2d());
 640   gp_Circ2d circ(axes, radius2 + radius);
 641   if (radius2 > radius && circ.Distance(nearp) > radius)
 642     circ.SetRadius(radius2 - radius);
 643
 644   // Calculate intersection of the line and the circle.
 645   IntAna2d_AnaIntersection intersector(line, circ);
 646   if (!intersector.IsDone() || !intersector.NbPoints())
 647     return Standard_False;
 648
 649   // Find center point of fillet.
 650   Standard_Integer i;
 651   Standard_Real minDist = DBL_MAX;
 652   for (i = 1; i <= intersector.NbPoints(); ++i)
 653   {
 654     const IntAna2d_IntPoint& intp = intersector.Point(i);
 655     const gp_Pnt2d& p = intp.Value();
 656
 657     Standard_Real d = nearl.Distance(p);
 658     if (d < minDist)
 659     {
 660       minDist = d;
 661       p.Coord(xc, yc);
 662     }
 663   }
 664
 665   // Shrink of segment.
 666   gp_Pnt2d pc(xc, yc);
 667   Standard_Real L2 = pc.SquareDistance(p12);
 668   const Standard_Real Rf2 = radius * radius;
 669   start = sqrt(L2 - Rf2);
 670
 671   // Shrink of arc.
 672   gp_Vec2d pcc(pc2, pc);
 673   end = fabs(gp_Vec2d(pc2, p12).Angle(pcc));
 674
 675   // Duplicate the information on shrink the arc:
 676   // calculate a point on the arc coinciding with the end of fillet.
 677   line.SetLocation(pc2);
 678   line.SetDirection(pcc);
 679   circ.SetLocation(pc2);
 680   circ.SetRadius(radius2);
 681   intersector.Perform(line, circ);
 682   if (!intersector.IsDone() || !intersector.NbPoints())
 683     return Standard_False;
 684
 685   xend = DBL_MAX;
 686   yend = DBL_MAX;
 687   for (i = 1; i <= intersector.NbPoints(); ++i)
 688   {
 689     const IntAna2d_IntPoint& intp = intersector.Point(i);
 690     const gp_Pnt2d& p = intp.Value();
 691
 692     const Standard_Real d2 = p.SquareDistance(pc);
 693     if (fabs(d2 - Rf2) < Precision::Confusion())
 694     {
 695       p.Coord(xend, yend);
 696       break;
 697     }
 698   }
 699
 700   // Orientation of the fillet.
 701   angle = v1.Angle(v2);
 702   cw = angle > 0.0;
 703   return Standard_True;
 704 }
 705
 706 // A function constructs a fillet between an arc and a segment.
 707 Standard_Boolean ChFi2d_AnaFilletAlgo::ArcFilletSegment(const Standard_Real radius,
 708                                                         Standard_Real& xc, Standard_Real& yc,
 709                                                         Standard_Boolean& cw,
 710                                                         Standard_Real& start, Standard_Real& end,
 711                                                         Standard_Real& xstart, Standard_Real& ystart)
 712 {
 713   // Make a line parallel to the segment at the side of center point of fillet.
 714   // This side may be defined through making a bisectrissa for vectors at p12 (or p21).
 715
 716   // Make 2D points.
 717   gp_Pnt2d p12(x12, y12);
 718   gp_Pnt2d p22(x22, y22);
 719   gp_Pnt2d pc1(xc1, yc1);
 720
 721   // Check length of segment.
 722   if (p12.SquareDistance(p22) < gp::Resolution())
 723     return Standard_False;
 724
 725   // Make 2D vectors.
 726   gp_Vec2d v1(p12, pc1);
 727   gp_Vec2d v2(p12, p22);
 728
 729   // Rotate the arc vector to become tangential at p21.
 730   if (cw1)
 731     v1.Rotate(-M_PI_2);
 732   else
 733     v1.Rotate(+M_PI_2);
 734
 735   // If vectors coincide (segment and arc are tangent),
 736   // the algorithm doesn't work...
 737   Standard_Real angle = v1.Angle(v2);
 738   if (fabs(angle) < Precision::Angular())
 739     return Standard_False;
 740
 741   // Make a bisectrissa of vectors at p12.
 742   v1.Normalize();
 743   v2.Normalize();
 744   gp_Vec2d bisec = 0.5 * (v1 + v2);
 745
 746   // If segment and arc look in opposite direction,
 747   // no fillet is possible.
 748   if (bisec.SquareMagnitude() < gp::Resolution())
 749     return Standard_False;
 750
 751   // Define an appropriate point to choose center of fillet.
 752   bisec.Normalize();
 753   gp_Pnt2d nearPoint = p12.Translated(radius * bisec);
 754   gp_Lin2d nearLine(p12, bisec);
 755
 756   // Make a line parallel to segment and
 757   // passing near the "near" point.
 758   gp_Vec2d aD2Vec(v2);
 759   gp_Lin2d line(p22, -aD2Vec);
 760   aD2Vec.Rotate(M_PI_2);
 761   line.Translate(radius * aD2Vec);
 762   if (line.Distance(nearPoint) > radius)
 763     line.Translate(-2.0 * radius * aD2Vec);
 764
 765   // Make a circle of radius of the arc +/- fillet radius.
 766   gp_Ax2d axes(pc1, gp::DX2d());
 767   gp_Circ2d circ(axes, radius1 + radius);
 768   if (radius1 > radius && circ.Distance(nearPoint) > radius)
 769     circ.SetRadius(radius1 - radius);
 770
 771   // Calculate intersection of the line and the big circle.
 772   IntAna2d_AnaIntersection intersector(line, circ);
 773   if (!intersector.IsDone() || !intersector.NbPoints())
 774     return Standard_False;
 775
 776   // Find center point of fillet.
 777   Standard_Integer i;
 778   Standard_Real minDist = DBL_MAX;
 779   for (i = 1; i <= intersector.NbPoints(); ++i)
 780   {
 781     const IntAna2d_IntPoint& intp = intersector.Point(i);
 782     const gp_Pnt2d& p = intp.Value();
 783
 784     Standard_Real d = nearLine.Distance(p);
 785     if (d < minDist)
 786     {
 787       minDist = d;
 788       p.Coord(xc, yc);
 789     }
 790   }
 791
 792   // Shrink of segment.
 793   gp_Pnt2d pc(xc, yc);
 794   Standard_Real L2 = pc.SquareDistance(p12);
 795   const Standard_Real Rf2 = radius * radius;
 796   end = sqrt(L2 - Rf2);
 797
 798   // Shrink of arc.
 799   gp_Vec2d pcc(pc1, pc);
 800   start = fabs(gp_Vec2d(pc1, p12).Angle(pcc));
 801
 802   // Duplicate the information on shrink the arc:
 803   // calculate a point on the arc coinciding with the start of fillet.
 804   line.SetLocation(pc1);
 805   line.SetDirection(pcc);
 806   circ.SetLocation(pc1);
 807   circ.SetRadius(radius1);
 808   intersector.Perform(line, circ);
 809   if (!intersector.IsDone() || !intersector.NbPoints())
 810     return Standard_False;
 811
 812   xstart = DBL_MAX;
 813   ystart = DBL_MAX;
 814   for (i = 1; i <= intersector.NbPoints(); ++i)
 815   {
 816     const IntAna2d_IntPoint& intp = intersector.Point(i);
 817     const gp_Pnt2d& p = intp.Value();
 818
 819     const Standard_Real d2 = p.SquareDistance(pc);
 820     if (fabs(d2 - Rf2) < Precision::SquareConfusion())
 821     {
 822       p.Coord(xstart, ystart);
 823       break;
 824     }
 825   }
 826
 827   // Orientation of the fillet.
 828   angle = v2.Angle(v1);
 829   cw = angle < 0.0;
 830   return Standard_True;
 831 }
 832
 833 // WW5 method to compute fillet: arc - arc.
 834 // It returns a constructed fillet definition:
 835 //     center point (xc, yc)
 836 //     shrinking parameter of the 1st circle (start)
 837 //     shrinking parameter of the 2nd circle (end)
 838 //     if the arc of fillet clockwise (cw = true) or counterclockwise (cw = false).
 839 Standard_Boolean ChFi2d_AnaFilletAlgo::ArcFilletArc(const Standard_Real radius,
 840                                                     Standard_Real& xc, Standard_Real& yc,
 841                                                     Standard_Boolean& cw,
 842                                                     Standard_Real& start, Standard_Real& end)
 843 {
 844   // Make points.
 845   const gp_Pnt2d pc1(xc1, yc1);
 846   const gp_Pnt2d pc2(xc2, yc2);
 847   const gp_Pnt2d p12(x12, y12);
 848
 849   // Make vectors at p12.
 850   gp_Vec2d v1(pc1, p12);
 851   gp_Vec2d v2(pc2, p12);
 852
 853   // Rotate the vectors so that they are tangent to circles at p12.
 854   if (cw1)
 855     v1.Rotate(+M_PI_2);
 856   else
 857     v1.Rotate(-M_PI_2);
 858   if (cw2)
 859     v2.Rotate(-M_PI_2);
 860   else
 861     v2.Rotate(+M_PI_2);
 862
 863   // Make a "check" point for choosing an offset circle.
 864   v1.Normalize();
 865   v2.Normalize();
 866   gp_Vec2d bisec = 0.5 * (v1 + v2);
 867   if (bisec.SquareMagnitude() < gp::Resolution())
 868     return Standard_False;
 869
 870   const gp_Pnt2d checkp = p12.Translated(radius * bisec);
 871   const gp_Lin2d checkl(p12, bisec);
 872
 873   // Make two circles of radius r1 +/- r and r2 +/- r
 874   // with center point equal to pc1 and pc2.
 875   // Arc 1.
 876   gp_Ax2d axes(pc1, gp::DX2d());
 877   gp_Circ2d c1(axes, radius1 + radius);
 878   if (radius1 > radius && c1.Distance(checkp) > radius)
 879     c1.SetRadius(radius1 - radius);
 880   // Arc 2.
 881   axes.SetLocation(pc2);
 882   gp_Circ2d c2(axes, radius2 + radius);
 883   if (radius2 > radius && c2.Distance(checkp) > radius)
 884     c2.SetRadius(radius2 - radius);
 885
 886   // Calculate an intersection point of these two circles
 887   // and choose the one closer to the "check" point.
 888   IntAna2d_AnaIntersection intersector(c1, c2);
 889   if (!intersector.IsDone() || !intersector.NbPoints())
 890     return Standard_False;
 891
 892   // Find center point of fillet.
 893   gp_Pnt2d pc;
 894   Standard_Real minDist = DBL_MAX;
 895   for (int i = 1; i <= intersector.NbPoints(); ++i)
 896   {
 897     const IntAna2d_IntPoint& intp = intersector.Point(i);
 898     const gp_Pnt2d& p = intp.Value();
 899
 900     Standard_Real d = checkp.SquareDistance(p);
 901     if (d < minDist)
 902     {
 903       minDist = d;
 904       pc = p;
 905     }
 906   }
 907   pc.Coord(xc, yc);
 908
 909   // Orientation of fillet.
 910   Standard_Real angle = v1.Angle(v2);
 911   if (fabs(angle) < Precision::Angular())
 912   {
 913     angle = gp_Vec2d(pc, pc1).Angle(gp_Vec2d(pc, pc2));
 914     cw = angle < 0.0;
 915   }
 916   else
 917   {
 918     cw = angle > 0.0;
 919   }
 920
 921   // Shrinking of circles.
 922   start = fabs(gp_Vec2d(pc1, p12).Angle(gp_Vec2d(pc1, pc)));
 923   end = fabs(gp_Vec2d(pc2, p12).Angle(gp_Vec2d(pc2, pc)));
 924   return Standard_True;
 925 }
 926
 927 // Cuts intersecting edges of a contour.
 928 Standard_Boolean ChFi2d_AnaFilletAlgo::Cut(const gp_Pln& thePlane, TopoDS_Edge& theE1, TopoDS_Edge& theE2)
 929 {
 930   gp_Pnt p;
 931   Standard_Boolean found(Standard_False);
 932   Standard_Real param1 = 0.0, param2 = 0.0;
 933   Standard_Real f1, l1, f2, l2;
 934   Handle(Geom_Curve) c1 = BRep_Tool::Curve(theE1, f1, l1);
 935   Handle(Geom_Curve) c2 = BRep_Tool::Curve(theE2, f2, l2);
 936   GeomAPI_ExtremaCurveCurve extrema(c1, c2, f1, l1, f2, l2);
 937   if (extrema.NbExtrema())
 938   {
 939     Standard_Integer i, nb = extrema.NbExtrema();
 940     for (i = 1; i <= nb; ++i)
 941     {
 942       const Standard_Real d = extrema.Distance(i);
 943       if (d < Precision::Confusion())
 944       {
 945         extrema.Parameters(i, param1, param2);
 946         if (fabs(l1 - param1) > Precision::Confusion() &&
 947             fabs(f2 - param2) > Precision::Confusion())
 948         {
 949           found = Standard_True;
 950           extrema.Points(i, p, p);
 951           break;
 952         }
 953       }
 954     }
 955   }
 956
 957   if (found)
 958   {
 959     BRepBuilderAPI_MakeEdge mkEdge1(c1, f1, param1);
 960     if (mkEdge1.IsDone())
 961     {
 962       theE1 = mkEdge1.Edge();
 963
 964       BRepBuilderAPI_MakeEdge mkEdge2(c2, param2, l2);
 965       if (mkEdge2.IsDone())
 966       {
 967         theE2 = mkEdge2.Edge();
 968
 969         gp_Pnt2d p2d = ProjLib::Project(thePlane, p);
 970         p2d.Coord(x12, y12);
 971         x21 = x12;
 972         y21 = y12;
 973         return Standard_True;
 974       }
 975     }
 976   }
 977   return Standard_False;
 978 }