src/ChFi2d/ChFi2d_FilletAlgo.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_FilletAlgo.hxx>
  15
  16 #include <GeomProjLib.hxx>
  17 #include <BRep_Tool.hxx>
  18 #include <Precision.hxx>
  19 #include <ElSLib.hxx>
  20 #include <ElCLib.hxx>
  21
  22 #include <Geom2dAPI_ProjectPointOnCurve.hxx>
  23 #include <GeomAPI_ProjectPointOnCurve.hxx>
  24 #include <Geom2dAPI_InterCurveCurve.hxx>
  25
  26 #include <TopoDS.hxx>
  27 #include <TopoDS_Iterator.hxx>
  28 #include <TColStd_ListIteratorOfListOfReal.hxx>
  29
  30 #include <gp_Circ.hxx>
  31 #include <Geom_Circle.hxx>
  32 #include <Geom2d_Line.hxx>
  33
  34 #include <BRepBuilderAPI_MakeEdge.hxx>
  35 #include <BRepAdaptor_Curve.hxx>
  36
  37 ChFi2d_FilletAlgo::ChFi2d_FilletAlgo()
  38 : myStart1(0.0),
  39   myEnd1  (0.0),
  40   myStart2(0.0),
  41   myEnd2  (0.0),
  42   myRadius(0.0),
  43   myStartSide    (Standard_False),
  44   myEdgesExchnged(Standard_False),
  45   myDegreeOfRecursion(0)
  46 {
  47 }
  48
  49 ChFi2d_FilletAlgo::ChFi2d_FilletAlgo(const TopoDS_Wire& theWire, const gp_Pln& thePlane)
  50 : myStart1(0.0),
  51   myEnd1  (0.0),
  52   myStart2(0.0),
  53   myEnd2  (0.0),
  54   myRadius(0.0),
  55   myStartSide    (Standard_False),
  56   myEdgesExchnged(Standard_False),
  57   myDegreeOfRecursion(0)
  58 {
  59   Init(theWire, thePlane);
  60 }
  61
  62 ChFi2d_FilletAlgo::ChFi2d_FilletAlgo(const TopoDS_Edge& theEdge1,
  63                                      const TopoDS_Edge& theEdge2,
  64                                      const gp_Pln& thePlane)
  65 : myEdge1(theEdge1),
  66   myEdge2(theEdge2),
  67   myStart1(0.0),
  68   myEnd1  (0.0),
  69   myStart2(0.0),
  70   myEnd2  (0.0),
  71   myRadius(0.0),
  72   myStartSide    (Standard_False),
  73   myEdgesExchnged(Standard_False),
  74   myDegreeOfRecursion(0)
  75 {
  76   Init(theEdge1, theEdge2, thePlane);
  77 }
  78
  79 void ChFi2d_FilletAlgo::Init(const TopoDS_Wire& theWire, const gp_Pln& thePlane)
  80 {
  81   TopoDS_Edge theEdge1, theEdge2;
  82   TopoDS_Iterator itr(theWire);
  83   for (; itr.More(); itr.Next())
  84   {
  85     if (theEdge1.IsNull())
  86       theEdge1 = TopoDS::Edge(itr.Value());
  87     else if (theEdge2.IsNull())
  88       theEdge2 = TopoDS::Edge(itr.Value());
  89     else
  90       break;
  91   }
  92   if (theEdge1.IsNull() || theEdge2.IsNull())
  93     throw Standard_ConstructionError("The fillet algorithms expects a wire consisting of two edges.");
  94   Init(theEdge1, theEdge2, thePlane);
  95 }
  96
  97 void ChFi2d_FilletAlgo::Init(const TopoDS_Edge& theEdge1,
  98                              const TopoDS_Edge& theEdge2,
  99                              const gp_Pln& thePlane)
 100 {
 101   myPlane = new Geom_Plane(thePlane);
 102
 103   myEdgesExchnged = Standard_False;
 104
 105   BRepAdaptor_Curve aBAC1(theEdge1);
 106   BRepAdaptor_Curve aBAC2(theEdge2);
 107   if (aBAC1.GetType() < aBAC2.GetType())
 108   { // first curve must be more complicated
 109     myEdge1 = theEdge2;
 110     myEdge2 = theEdge1;
 111     myEdgesExchnged = Standard_True;
 112   }
 113   else
 114   {
 115     myEdge1 = theEdge1;
 116     myEdge2 = theEdge2;
 117   }
 118
 119   Handle(Geom_Curve) aCurve1 = BRep_Tool::Curve(myEdge1, myStart1, myEnd1);
 120   Handle(Geom_Curve) aCurve2 = BRep_Tool::Curve(myEdge2, myStart2, myEnd2);
 121
 122   myCurve1 = GeomProjLib::Curve2d(aCurve1, myStart1, myEnd1, myPlane);
 123   myCurve2 = GeomProjLib::Curve2d(aCurve2, myStart2, myEnd2, myPlane);
 124
 125   while (myCurve1->IsPeriodic() && myStart1 >= myEnd1)
 126     myEnd1 += myCurve1->Period();
 127   while (myCurve2->IsPeriodic() && myStart2 >= myEnd2)
 128     myEnd2 += myCurve2->Period();
 129
 130   if (aBAC1.GetType() == aBAC2.GetType())
 131   {
 132     if (myEnd2 - myStart2 < myEnd1 - myStart1)
 133     { // first curve must be parametrically shorter
 134       TopoDS_Edge anEdge = myEdge1;
 135       myEdge1 = myEdge2;
 136       myEdge2 = anEdge;
 137       Handle(Geom2d_Curve) aCurve = myCurve1;
 138       myCurve1 = myCurve2;
 139       myCurve2 = aCurve;
 140       Standard_Real a = myStart1;
 141       myStart1 = myStart2;
 142       myStart2 = a;
 143       a = myEnd1;
 144       myEnd1 = myEnd2;
 145       myEnd2 = a;
 146       myEdgesExchnged = Standard_True;
 147     }
 148   }
 149 }
 150
 151 //! This function returns true if linear segment from start point of the
 152 //! fillet arc to the end point is intersected by the first or second
 153 //! curve: in this case fillet is invalid.
 154 static Standard_Boolean IsRadiusIntersected(const Handle(Geom2d_Curve)& theCurve, const Standard_Real theCurveMin, const double theCurveMax,
 155                                             const gp_Pnt2d theStart, const gp_Pnt2d theEnd, const Standard_Boolean theStartConnected)
 156 {
 157   //Check the given start and end if they are identical. If yes
 158   //return false
 159   if (theStart.SquareDistance(theEnd) < Precision::SquareConfusion())
 160   {
 161     return Standard_False;
 162   }
 163   Handle(Geom2d_Line) line = new Geom2d_Line(theStart, gp_Dir2d(gp_Vec2d(theStart, theEnd)));
 164   Geom2dAPI_InterCurveCurve anInter(theCurve, line, Precision::Confusion());
 165   Standard_Integer a;
 166   gp_Pnt2d aPoint;
 167   for(a = anInter.NbPoints(); a > 0; a--)
 168   {
 169     aPoint = anInter.Point(a);
 170     Geom2dAPI_ProjectPointOnCurve aProjInt(aPoint, theCurve, theCurveMin, theCurveMax);
 171     if (aProjInt.NbPoints() < 1 || aProjInt.LowerDistanceParameter() > Precision::Confusion())
 172       continue; // point is not on edge
 173
 174     if (aPoint.Distance(theStart) < Precision::Confusion())
 175     {
 176       if (!theStartConnected)
 177         return Standard_True;
 178     }
 179     if (aPoint.Distance(theEnd) < Precision::Confusion())
 180       return Standard_True;
 181     if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), Precision::Angular()))
 182       return Standard_True;
 183   }
 184   Handle(Geom2d_Curve) aCurve = theCurve;
 185   for(a = anInter.NbSegments(); a > 0; a--)
 186   {
 187     //anInter.Segment(a, aCurve); //not implemented (bug in OCC)
 188     aPoint = aCurve->Value(aCurve->FirstParameter());
 189
 190     Geom2dAPI_ProjectPointOnCurve aProjInt(aPoint, theCurve, theCurveMin, theCurveMax);
 191     if (aProjInt.NbPoints() && aProjInt.LowerDistanceParameter() < Precision::Confusion())
 192     { // point is on edge
 193       if (aPoint.Distance(theStart) < Precision::Confusion())
 194         if (!theStartConnected)
 195           return Standard_True;
 196       if (aPoint.Distance(theEnd) < Precision::Confusion())
 197         return Standard_True;
 198       if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), Precision::Angular()))
 199         return Standard_True;
 200     }
 201     aPoint = aCurve->Value(aCurve->LastParameter());
 202
 203     aProjInt.Init(aPoint, theCurve, theCurveMin, theCurveMax);
 204     if (aProjInt.NbPoints() && aProjInt.LowerDistanceParameter() < Precision::Confusion())
 205     { // point is on edge
 206       if (aPoint.Distance(theStart) < Precision::Confusion())
 207         if (!theStartConnected)
 208           return Standard_True;
 209       if (aPoint.Distance(theEnd) < Precision::Confusion())
 210         return Standard_True;
 211       if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), Precision::Angular()))
 212         return Standard_True;
 213     }
 214   }
 215   return Standard_False;
 216 }
 217
 218 void ChFi2d_FilletAlgo::FillPoint(FilletPoint* thePoint, const Standard_Real theLimit)
 219 {
 220
 221   // on the intersection point
 222   Standard_Boolean aValid = Standard_False;
 223   Standard_Real aStep = Precision::Confusion();
 224   gp_Pnt2d aCenter, aPoint; // center of fillet and point on curve1
 225   Standard_Real aParam = thePoint->getParam();
 226   if (theLimit < aParam) aStep = -aStep;
 227   for(aValid = Standard_False; !aValid; aParam += aStep)
 228   {
 229     if ((aParam - aStep - theLimit) * (aParam - theLimit) <= 0)
 230       break; // limit was exceeded
 231     aStep *= 2;
 232     gp_Vec2d aVec;
 233     myCurve1->D1(aParam, aPoint, aVec);
 234     if (aVec.SquareMagnitude() < Precision::Confusion())
 235       continue;
 236
 237     gp_Vec2d aPerp(((myStartSide)?-1:1) * aVec.Y(), ((myStartSide)?1:-1) * aVec.X());
 238     aPerp.Normalize();
 239     aPerp.Multiply(myRadius);
 240     aCenter = aPoint.Translated(aPerp);
 241
 242
 243     Geom2dAPI_ProjectPointOnCurve aProjInt(aPoint, myCurve2, myStart2, myEnd2);
 244     if (aProjInt.NbPoints() == 0 || aPoint.Distance(aProjInt.NearestPoint()) > Precision::Confusion())
 245     {
 246       aValid = Standard_True;
 247       break;
 248     }
 249   }
 250   if (aValid)
 251   {
 252     thePoint->setParam(aParam);
 253     thePoint->setCenter(aCenter);
 254     aValid = !IsRadiusIntersected(myCurve2, myStart2, myEnd2, aPoint, aCenter, Standard_True);
 255   }
 256
 257   Geom2dAPI_ProjectPointOnCurve aProj(aCenter, myCurve2);
 258   int a, aNB = aProj.NbPoints();
 259   for(a = aNB; a > 0; a--)
 260   {
 261     if (aPoint.SquareDistance(aProj.Point(a)) < Precision::Confusion())
 262       continue;
 263
 264     Standard_Boolean aValid2 = aValid;
 265     if (aValid2)
 266       aValid2 = !IsRadiusIntersected(myCurve1, myStart1, myEnd1, aCenter, aProj.Point(a), Standard_False);
 267
 268     // checking the right parameter
 269     Standard_Real aParamProj = aProj.Parameter(a);
 270     while(myCurve2->IsPeriodic() && aParamProj < myStart2)
 271       aParamProj += myCurve2->Period();
 272
 273     const Standard_Real d = aProj.Distance(a);
 274     thePoint->appendValue(d * d - myRadius * myRadius, (aParamProj >= myStart2 && aParamProj <= myEnd2 && aValid2));
 275     if (Abs(d - myRadius) < Precision::Confusion())
 276       thePoint->setParam2(aParamProj);
 277   }
 278 }
 279
 280 void ChFi2d_FilletAlgo::FillDiff(FilletPoint* thePoint, Standard_Real theDiffStep, Standard_Boolean theFront)
 281 {
 282   Standard_Real aDelta = theFront?(theDiffStep):(-theDiffStep);
 283   FilletPoint* aDiff = new FilletPoint(thePoint->getParam() + aDelta);
 284   FillPoint(aDiff, aDelta * 999.);
 285   if (!thePoint->calculateDiff(aDiff))
 286   {
 287     aDiff->setParam(thePoint->getParam() - aDelta);
 288     FillPoint(aDiff,  - aDelta * 999);
 289     thePoint->calculateDiff(aDiff);
 290   }
 291   delete aDiff;
 292 }
 293
 294 // returns true, if at least one result was found
 295 Standard_Boolean ChFi2d_FilletAlgo::Perform(const Standard_Real theRadius)
 296 {
 297   myDegreeOfRecursion = 0;
 298   myResultParams.Clear();
 299   myResultOrientation.Clear();
 300
 301   Standard_Real aNBSteps;
 302   Geom2dAdaptor_Curve aGAC(myCurve1);
 303   switch (aGAC.GetType())
 304   {
 305   case GeomAbs_Line:
 306     aNBSteps = 2;
 307     break;
 308   case GeomAbs_Circle:
 309     aNBSteps = 4;
 310     break;
 311   case GeomAbs_Ellipse:
 312     aNBSteps = 5;
 313     break;
 314   case GeomAbs_BSplineCurve:
 315     aNBSteps = 2 + aGAC.Degree() * aGAC.NbPoles();
 316     break;
 317   default: // unknown: maximum
 318     aNBSteps = 100;
 319   }
 320   //cout<<"aNBSteps = "<<aNBSteps<<endl;
 321
 322   myRadius = theRadius;
 323   Standard_Real aParam, aStep, aDStep;
 324   aStep = (myEnd1 - myStart1) / aNBSteps;
 325   aDStep = 1.e-4 * aStep;
 326
 327   int aCycle;
 328   for(aCycle = 2, myStartSide = Standard_False; aCycle; myStartSide = !myStartSide, aCycle--)
 329   {
 330     FilletPoint *aLeft = NULL, *aRight;
 331
 332     for(aParam = myStart1 + aStep; aParam < myEnd1 || Abs(myEnd1 - aParam) < Precision::Confusion(); aParam += aStep)
 333     {
 334       if (!aLeft)
 335       {
 336         aLeft = new FilletPoint(aParam - aStep);
 337         FillPoint(aLeft, aParam);
 338         FillDiff(aLeft, aDStep, Standard_True);
 339       }
 340
 341       aRight = new FilletPoint(aParam);
 342       FillPoint(aRight, aParam - aStep);
 343       FillDiff(aRight, aDStep, Standard_False);
 344
 345       aLeft->FilterPoints(aRight);
 346       PerformNewton(aLeft, aRight);
 347
 348       delete aLeft;
 349       aLeft = aRight;
 350     }//for
 351     delete aLeft;
 352   }//for
 353
 354   return !myResultParams.IsEmpty();
 355 }
 356
 357 Standard_Boolean ChFi2d_FilletAlgo::ProcessPoint(FilletPoint* theLeft, FilletPoint* theRight, Standard_Real theParameter)
 358 {
 359   if (theParameter >= theLeft->getParam() && theParameter < theRight->getParam())
 360   {
 361     Standard_Real aDX = (theRight->getParam() - theLeft->getParam());
 362     if (theParameter - theLeft->getParam() < aDX/100.0)
 363     {
 364       theParameter = theLeft->getParam() + aDX/100.0;
 365     }
 366     if (theRight->getParam() - theParameter < aDX/100.0)
 367     {
 368       theParameter = theRight->getParam() - aDX/100.0;
 369     }
 370
 371     // Protection on infinite loops.
 372     myDegreeOfRecursion++;
 373     Standard_Real diffx = 0.001 * aDX;
 374     if (myDegreeOfRecursion > 100)
 375     {
 376       diffx *= 10.0;
 377       if (myDegreeOfRecursion > 1000)
 378       {
 379         diffx *= 10.0;
 380         if (myDegreeOfRecursion > 3000)
 381         {
 382           return Standard_True;
 383         }
 384       }
 385     }
 386
 387     FilletPoint* aPoint1 = theLeft->Copy();
 388     FilletPoint* aPoint2 = new FilletPoint(theParameter);
 389     FillPoint(aPoint2, aPoint1->getParam());
 390     FillDiff(aPoint2, diffx, Standard_True);
 391
 392     aPoint1->FilterPoints(aPoint2);
 393     PerformNewton(aPoint1, aPoint2);
 394     aPoint2->FilterPoints(theRight);
 395     PerformNewton(aPoint2, theRight);
 396
 397     delete aPoint1;
 398     delete aPoint2;
 399     return Standard_True;
 400   }
 401
 402   return Standard_False;
 403 }
 404
 405 void ChFi2d_FilletAlgo::PerformNewton(FilletPoint* theLeft, FilletPoint* theRight)
 406 {
 407   int a;
 408   // check the left: if this is solution store it and remove it from the list of researching points of theLeft
 409   a = theLeft->hasSolution(myRadius);
 410   if (a)
 411   {
 412     if (theLeft->isValid(a))
 413     {
 414       myResultParams.Append(theLeft->getParam());
 415       myResultOrientation.Append(myStartSide);
 416     }
 417     return;
 418   }
 419
 420   Standard_Real aDX = theRight->getParam() - theLeft->getParam();
 421   if (aDX < 1.e-6 * Precision::Confusion())
 422   {
 423     a = theRight->hasSolution(myRadius);
 424     if (a && theRight->isValid(a))
 425     {
 426       myResultParams.Append(theRight->getParam());
 427       myResultOrientation.Append(myStartSide);
 428     }
 429     return;
 430   }
 431   for(a = 1; a <= theLeft->getNBValues(); a++)
 432   {
 433     int aNear = theLeft->getNear(a);
 434
 435     Standard_Real aA = (theRight->getDiff(aNear) - theLeft->getDiff(a)) / aDX;
 436     Standard_Real aB = theLeft->getDiff(a) - aA * theLeft->getParam();
 437     Standard_Real aC = theLeft->getValue(a) - theLeft->getDiff(a) * theLeft->getParam() + aA * theLeft->getParam() * theLeft->getParam() / 2.0;
 438     Standard_Real aDet = aB * aB - 2.0 * aA * aC;
 439
 440     if (Abs(aA) < Precision::Confusion())
 441     { // linear case
 442       //cout<<"###"<<endl;
 443       if (Abs(aB) > 10e-20)
 444       {
 445         Standard_Real aX0 = - aC / aB; // use extremum
 446         if (aX0 > theLeft->getParam() && aX0 < theRight->getParam())
 447           ProcessPoint(theLeft, theRight, aX0);
 448       }
 449       else
 450       {
 451         ProcessPoint(theLeft, theRight, theLeft->getParam() + aDX / 2.0); // linear division otherwise
 452       }
 453     }
 454     else
 455     {
 456       if (Abs(aB) > Abs(aDet * 1000000.))
 457       { // possible floating point operations accurancy errors
 458         //cout<<"*";
 459         ProcessPoint(theLeft, theRight, theLeft->getParam() + aDX / 2.0); // linear division otherwise
 460       }
 461       else
 462       {
 463         if (aDet > 0)
 464         { // two solutions
 465           aDet = sqrt(aDet);
 466           Standard_Boolean aRes = ProcessPoint(theLeft, theRight, (- aB + aDet) / aA);
 467           if (!aRes)
 468             aRes = ProcessPoint(theLeft, theRight, (- aB - aDet) / aA);
 469           if (!aRes)
 470             ProcessPoint(theLeft, theRight, theLeft->getParam() + aDX / 2.0); // linear division otherwise
 471         }
 472         else
 473         {
 474           //cout<<"%%%"<<endl;
 475           Standard_Real aX0 = - aB / aA; // use extremum
 476           if (aX0 > theLeft->getParam() && aX0 < theRight->getParam())
 477             ProcessPoint(theLeft, theRight, aX0);
 478           else
 479             ProcessPoint(theLeft, theRight, theLeft->getParam() + aDX / 2.0); // linear division otherwise
 480         }
 481       }
 482     }
 483   }//for
 484 }
 485
 486 // returns number of possible solutions.
 487 int ChFi2d_FilletAlgo::NbResults(const gp_Pnt& thePoint)
 488 {
 489   Standard_Real aX, aY;
 490   gp_Pnt2d aTargetPoint2d;
 491   ElSLib::PlaneParameters(myPlane->Pln().Position(), thePoint, aX, aY);
 492   aTargetPoint2d.SetCoord(aX, aY);
 493
 494   //iterate through all possible solutions.
 495   int i = 1, nb = 0;
 496   TColStd_ListIteratorOfListOfReal anIter(myResultParams);
 497   for(; anIter.More(); anIter.Next(), i++)
 498   {
 499     myStartSide = (myResultOrientation.Value(i)) ? Standard_True : Standard_False;
 500     FilletPoint *aPoint = new FilletPoint(anIter.Value());
 501     FillPoint(aPoint, anIter.Value() + 1.);
 502     if (aPoint->hasSolution(myRadius))
 503       nb++;
 504     delete aPoint;
 505   }//for
 506
 507   return nb;
 508 }
 509
 510 // returns result (fillet edge, modified edge1, modified edge2), neares to the given point <thePoint>
 511 TopoDS_Edge ChFi2d_FilletAlgo::Result(const gp_Pnt& thePoint, TopoDS_Edge& theEdge1, TopoDS_Edge& theEdge2, const int iSolution)
 512 {
 513   TopoDS_Edge aResult;
 514   gp_Pnt2d aTargetPoint2d;
 515   Standard_Real aX, aY;
 516   ElSLib::PlaneParameters(myPlane->Pln().Position(), thePoint, aX, aY);
 517   aTargetPoint2d.SetCoord(aX, aY);
 518
 519   // choose the nearest circle
 520   Standard_Real aDistance = 0.0, aP;
 521   FilletPoint *aNearest;
 522   int a, iSol = 1;
 523   TColStd_ListIteratorOfListOfReal anIter(myResultParams);
 524   for(aNearest = NULL, a = 1; anIter.More(); anIter.Next(), a++)
 525   {
 526     myStartSide = (myResultOrientation.Value(a))?Standard_True:Standard_False;
 527     FilletPoint *aPoint = new FilletPoint(anIter.Value());
 528     FillPoint(aPoint, anIter.Value() + 1.);
 529     if (!aPoint->hasSolution(myRadius))
 530     {
 531       delete aPoint;
 532       continue;
 533     }
 534     aP = DBL_MAX;
 535     if (iSolution == -1)
 536     {
 537       aP = Abs(aPoint->getCenter().Distance(aTargetPoint2d) - myRadius);
 538     }
 539     else if (iSolution == iSol)
 540     {
 541       aP = 0.0;
 542     }
 543     if (!aNearest || aP < aDistance)
 544     {
 545       aNearest = aPoint;
 546       aDistance = aP;
 547     }
 548     else
 549     {
 550       delete aPoint;
 551     }
 552     if (iSolution == iSol)
 553       break;
 554     iSol++;
 555   }//for
 556
 557   if (!aNearest)
 558     return aResult;
 559
 560   // create circle edge
 561   gp_Pnt aCenter = ElSLib::PlaneValue(aNearest->getCenter().X(), aNearest->getCenter().Y(), myPlane->Pln().Position());
 562   Handle(Geom_Circle) aCircle = new Geom_Circle(gp_Ax2(aCenter, myPlane->Pln().Axis().Direction()), myRadius);
 563   gp_Pnt2d aPoint2d1, aPoint2d2;
 564   myCurve1->D0(aNearest->getParam(), aPoint2d1);
 565   myCurve2->D0(aNearest->getParam2(), aPoint2d2);
 566   gp_Pnt aPoint1 = ElSLib::PlaneValue(aPoint2d1.X(), aPoint2d1.Y(), myPlane->Pln().Position());
 567   gp_Pnt aPoint2 = ElSLib::PlaneValue(aPoint2d2.X(), aPoint2d2.Y(), myPlane->Pln().Position());
 568
 569   GeomAPI_ProjectPointOnCurve aProj(thePoint, aCircle);
 570   Standard_Real aTargetParam = aProj.LowerDistanceParameter();
 571   gp_Pnt aPointOnCircle = aProj.NearestPoint();
 572
 573   // There is a bug in Open CASCADE in calculation of nearest point to a circle near the parameter 0.0
 574   // Therefore I check this extrema point manually:
 575   gp_Pnt p0 = ElCLib::Value(0.0, aCircle->Circ());
 576   if (p0.Distance(thePoint) < aPointOnCircle.Distance(thePoint))
 577   {
 578     aTargetParam = 0.0;
 579     aPointOnCircle = p0;
 580   }
 581
 582   aProj.Perform(aPoint1);
 583   Standard_Real aParam1 = aProj.LowerDistanceParameter();
 584   aProj.Perform(aPoint2);
 585   Standard_Real aParam2 = aProj.LowerDistanceParameter();
 586   Standard_Boolean aIsOut = ((aParam1 < aTargetParam && aParam2 < aTargetParam) || (aParam1 > aTargetParam && aParam2 > aTargetParam));
 587   if (aParam1 > aParam2)
 588     aIsOut = !aIsOut;
 589   BRepBuilderAPI_MakeEdge aBuilder(aCircle->Circ(), aIsOut ? aParam2 : aParam1, aIsOut? aParam1 : aParam2);
 590   aResult = aBuilder.Edge();
 591
 592   // divide edges
 593   Standard_Real aStart, anEnd;
 594   Handle(Geom_Curve) aCurve = BRep_Tool::Curve(myEdge1, aStart, anEnd);
 595   gp_Vec aDir;
 596   aCurve->D1(aNearest->getParam(), aPoint1, aDir);
 597
 598   gp_Vec aCircleDir;
 599   aCircle->D1(aParam1, aPoint1, aCircleDir);
 600
 601   if ((aCircleDir.Angle(aDir) > M_PI / 2.0) ^ aIsOut)
 602     aStart = aNearest->getParam();
 603   else
 604     anEnd = aNearest->getParam();
 605
 606   //Check the case when start and end are identical. This happens
 607   //when the edge decreases to size 0. Old ww5 allows such
 608   //cases. So we are again bug compatible
 609   if (fabs(aStart - anEnd) < Precision::Confusion())
 610     anEnd = aStart + Precision::Confusion();
 611   //Divide edge
 612   BRepBuilderAPI_MakeEdge aDivider1(aCurve, aStart, anEnd);
 613   if (myEdgesExchnged)
 614     theEdge2 = aDivider1.Edge();
 615   else
 616     theEdge1 = aDivider1.Edge();
 617
 618   aCurve = BRep_Tool::Curve(myEdge2, aStart, anEnd);
 619   aCurve->D1(aNearest->getParam2(), aPoint2, aDir);
 620
 621   aCircle->D1(aParam2, aPoint2, aCircleDir);
 622
 623   if ((aCircleDir.Angle(aDir) > M_PI / 2.0) ^ (!aIsOut))
 624     aStart = aNearest->getParam2();
 625   else
 626     anEnd = aNearest->getParam2();
 627
 628   //Check the case when start and end are identical. This happens
 629   //when the edge decreases to size 0. Old ww5 allows such
 630   //cases. So we are again bug compatible
 631   if (fabs(aStart - anEnd) < Precision::Confusion())
 632     anEnd = aStart + Precision::Confusion();
 633   BRepBuilderAPI_MakeEdge aDivider2(aCurve, aStart, anEnd);
 634   if (myEdgesExchnged)
 635     theEdge1 = aDivider2.Edge();
 636   else
 637     theEdge2 = aDivider2.Edge();
 638
 639   delete aNearest;
 640   return aResult;
 641 }
 642
 643 FilletPoint::FilletPoint(const Standard_Real theParam)
 644 : myParam (theParam),
 645   myParam2(0.0)
 646 {
 647 }
 648
 649 void FilletPoint::appendValue(Standard_Real theValue, Standard_Boolean theValid)
 650 {
 651   Standard_Integer a;
 652   for(a = 1; a <= myV.Length(); a++)
 653   {
 654     if (theValue < myV.Value(a))
 655     {
 656       myV.InsertBefore(a, theValue);
 657       myValid.InsertBefore(a, theValid);
 658       return;
 659     }
 660   }
 661   myV.Append(theValue);
 662   myValid.Append(theValid);
 663 }
 664
 665 Standard_Boolean FilletPoint::calculateDiff(FilletPoint* thePoint)
 666 {
 667   Standard_Integer a;
 668   Standard_Boolean aDiffsSet = (myD.Length() != 0);
 669   Standard_Real aDX = thePoint->getParam() - myParam, aDY = 0.0;
 670   if (thePoint->myV.Length() == myV.Length())
 671   { // absolutely the same points
 672     for(a = 1; a <= myV.Length(); a++)
 673     {
 674       aDY = thePoint->myV.Value(a) - myV.Value(a);
 675       if (aDiffsSet)
 676         myD.SetValue(a, aDY / aDX);
 677       else
 678         myD.Append(aDY / aDX);
 679     }
 680     return Standard_True;
 681   }
 682   // between the diffeerent points searching for nearest analogs
 683   Standard_Integer b;
 684   for(a = 1; a <= myV.Length(); a++)
 685   {
 686     for(b = 1; b <= thePoint->myV.Length(); b++)
 687     {
 688       if (b == 1 || Abs(thePoint->myV.Value(b) - myV.Value(a)) < Abs(aDY))
 689         aDY = thePoint->myV.Value(b) - myV.Value(a);
 690     }
 691     if (aDiffsSet)
 692     {
 693       if (Abs(aDY / aDX) < Abs(myD.Value(a)))
 694         myD.SetValue(a, aDY / aDX);
 695     }
 696     else
 697     {
 698       myD.Append(aDY / aDX);
 699     }
 700   }//for
 701
 702   return Standard_False;
 703 }
 704
 705 void FilletPoint::FilterPoints(FilletPoint* thePoint)
 706 {
 707   Standard_Integer a, b;
 708   TColStd_SequenceOfReal aDiffs;
 709   Standard_Real aY, aY2, aDX = thePoint->getParam() - myParam;
 710   for(a = 1; a <= myV.Length(); a++)
 711   {
 712     // searching for near point from thePoint
 713     Standard_Integer aNear = 0;
 714     Standard_Real aDiff = aDX * 10000.;
 715     aY = myV.Value(a) + myD.Value(a) * aDX;
 716     for(b = 1; b <= thePoint->myV.Length(); b++)
 717     {
 718       // calculate hypothesis value of the Y2 with the constant first and second derivative
 719       aY2 = aY + aDX * (thePoint->myD.Value(b) - myD.Value(a)) / 2.0;
 720       if (aNear == 0 || Abs(aY2 - thePoint->myV.Value(b)) < Abs(aDiff))
 721       {
 722         aNear = b;
 723         aDiff = aY2 - thePoint->myV.Value(b);
 724       }
 725     }//for b...
 726
 727     if (aNear)
 728     {
 729       if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
 730       {// the same sign at the same sides of the interval
 731         if (myV.Value(a) * myD.Value(a) > 0)
 732         {
 733           if (Abs(myD.Value(a)) > Precision::Confusion())
 734             aNear = 0;
 735         }
 736         else
 737         {
 738           if (Abs(myV.Value(a)) > Abs(thePoint->myV.Value(aNear)))
 739             if (thePoint->myV.Value(aNear) * thePoint->myD.Value(aNear) < 0 && Abs(thePoint->myD.Value(aNear)) > Precision::Confusion())
 740             {
 741               aNear = 0;
 742             }
 743         }
 744       }
 745     }//if aNear
 746
 747     if (aNear)
 748     {
 749       if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
 750       {
 751         if ((myV.Value(a) + myD.Value(a) * aDX) * myV.Value(a) > Precision::Confusion() &&
 752             (thePoint->myV.Value(aNear) + thePoint->myD.Value(aNear) * aDX) * thePoint->myV.Value(aNear) > Precision::Confusion())
 753         {
 754           aNear = 0;
 755         }
 756       }
 757     }//if aNear
 758
 759     if (aNear)
 760     {
 761       if (Abs(aDiff / aDX) > 1.e+7)
 762       {
 763         aNear = 0;
 764       }
 765     }
 766
 767     if (aNear == 0)
 768     {   // there is no near: remove it from the list
 769       myV.Remove(a);
 770       myD.Remove(a);
 771       myValid.Remove(a);
 772       a--;
 773     }
 774     else
 775     {
 776       Standard_Boolean aFound = Standard_False;
 777       for(b = 1; b <= myNear.Length(); b++)
 778       {
 779         if (myNear.Value(b) == aNear)
 780         {
 781           if (Abs(aDiffs.Value(b)) < Abs(aDiff))
 782           { // return this 'near'
 783             aFound = Standard_True;
 784             myV.Remove(a);
 785             myD.Remove(a);
 786             myValid.Remove(a);
 787             a--;
 788             break;
 789           }
 790           else
 791           { // remove the old 'near'
 792             myV.Remove(b);
 793             myD.Remove(b);
 794             myValid.Remove(b);
 795             myNear.Remove(b);
 796             aDiffs.Remove(b);
 797             a--;
 798             break;
 799           }
 800         }
 801       }//for b...
 802       if (!aFound)
 803       {
 804         myNear.Append(aNear);
 805         aDiffs.Append(aDiff);
 806       }
 807     }//else
 808   }//for a...
 809 }
 810
 811 FilletPoint* FilletPoint::Copy()
 812 {
 813   FilletPoint* aCopy = new FilletPoint(myParam);
 814   Standard_Integer a;
 815   for(a = 1; a <= myV.Length(); a++)
 816   {
 817     aCopy->myV.Append(myV.Value(a));
 818     aCopy->myD.Append(myD.Value(a));
 819     aCopy->myValid.Append(myValid.Value(a));
 820   }
 821   return aCopy;
 822 }
 823
 824 int FilletPoint::hasSolution(const Standard_Real theRadius)
 825 {
 826   Standard_Integer a;
 827   for(a = 1; a <= myV.Length(); a++)
 828   {
 829     if (Abs(sqrt(Abs(Abs(myV.Value(a)) + theRadius * theRadius)) - theRadius) < Precision::Confusion())
 830       return a;
 831   }
 832   return 0;
 833 }
 834
 835 void FilletPoint::remove(int theIndex)
 836 {
 837   myV.Remove(theIndex);
 838   myD.Remove(theIndex);
 839   myValid.Remove(theIndex);
 840   myNear.Remove(theIndex);
 841 }