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