src/GEOMImpl/GEOMImpl_Fillet1d.cxx

   1 // Copyright (C) 2007-2023  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 //  File   : GEOMImpl_Fillet1d.cxx
  21 //  Module : GEOMImpl
  22
  23 #include "GEOMImpl_Fillet1d.hxx"
  24
  25 #include <BRep_Tool.hxx>
  26 #include <BRepAdaptor_Curve.hxx>
  27 #include <BRepBuilderAPI_MakeEdge.hxx>
  28 #include <ElCLib.hxx>
  29 #include <ElSLib.hxx>
  30 #include <gp_Circ.hxx>
  31 #include <Geom2d_Line.hxx>
  32 #include <Geom2dAPI_ProjectPointOnCurve.hxx>
  33 #include <Geom2dAPI_InterCurveCurve.hxx>
  34 #include <GeomAPI_ProjectPointOnCurve.hxx>
  35 #include <GeomProjLib.hxx>
  36 #include <Geom_Circle.hxx>
  37 #include <Precision.hxx>
  38 #include <TColStd_ListIteratorOfListOfReal.hxx>
  39 #include <IntRes2d_IntersectionSegment.hxx>
  40 #include <TopExp.hxx>
  41
  42 #include <Standard_NotImplemented.hxx>
  43
  44
  45 /**
  46  * This function returns Standard_True if it is possible to divide edge, i.e.
  47  * if one parameter either start or end one is inside the edge. This function
  48  * is used in the method GEOMImpl_Fillet1d::Result.
  49  *
  50  * \param theEdge the edge
  51  * \param theStart the start parameter
  52  * \param theEnd the end parameter
  53  * \return Standard_True if it is possible to split edge;
  54  *         Standard_False otherwise.
  55  */
  56 static Standard_Boolean IsDivideEdge(const TopoDS_Edge   &theEdge,
  57                                      const Standard_Real  theStart,
  58                                      const Standard_Real  theEnd)
  59 {
  60   Standard_Real      aFirst;
  61   Standard_Real      aLast;
  62   Handle(Geom_Curve) aCurve    = BRep_Tool::Curve(theEdge, aFirst, aLast);
  63   gp_Pnt             aPStart   = aCurve->Value(theStart);
  64   gp_Pnt             aPEnd     = aCurve->Value(theEnd);
  65   TopoDS_Vertex      aVFirst   = TopExp::FirstVertex(theEdge);
  66   TopoDS_Vertex      aVLast    = TopExp::LastVertex(theEdge);
  67   Standard_Real      aTolFirst = BRep_Tool::Tolerance(aVFirst);
  68   Standard_Real      aTolLast  = BRep_Tool::Tolerance(aVLast);
  69   Standard_Real      aTolConf  = Precision::Confusion();
  70   gp_Pnt             aPFirst   = BRep_Tool::Pnt(aVFirst);
  71   gp_Pnt             aPLast    = BRep_Tool::Pnt(aVLast);
  72   Standard_Real      aDistSF   = aPStart.Distance(aPFirst);
  73   Standard_Real      aDistSL   = aPStart.Distance(aPLast);
  74   Standard_Real      aDistEF   = aPEnd.Distance(aPFirst);
  75   Standard_Real      aDistEL   = aPEnd.Distance(aPLast);
  76   Standard_Boolean   isSplit   = Standard_True;
  77
  78   if (aDistSF <= aTolFirst + aTolConf ||
  79       aDistSL <= aTolLast  + aTolConf) {
  80     if (aDistEF <= aTolFirst + aTolConf ||
  81         aDistEL <= aTolLast  + aTolConf) {
  82
  83       isSplit = Standard_False;
  84       // in this case the original edge is thrown, and distance (gap) from new arc end
  85       // to a vertex of original wire can reach (aVertexTolerance + Precision::Confusion()).
  86       // Resulting wire is fixed (Mantis issue 0023411) in GEOMImpl_Fillet1dDriver::MakeFillet()
  87     }
  88   }
  89
  90   return isSplit;
  91 }
  92
  93 /**
  94  * class GEOMImpl_Fillet1d
  95  */
  96
  97 //=======================================================================
  98 //function : Constructor
  99 //purpose  :
 100 //=======================================================================
 101 GEOMImpl_Fillet1d::GEOMImpl_Fillet1d(const TopoDS_Edge& theEdge1,
 102                                      const TopoDS_Edge& theEdge2,
 103                                      const gp_Pln& thePlane)
 104 : myEdgesExchnged( Standard_False )
 105 {
 106   myPlane = new Geom_Plane(thePlane);
 107
 108   BRepAdaptor_Curve aBAC1(theEdge1);
 109   BRepAdaptor_Curve aBAC2(theEdge2);
 110   if (aBAC1.GetType() < aBAC2.GetType())
 111   { // first curve must be more complicated
 112     myEdge1 = theEdge2;
 113     myEdge2 = theEdge1;
 114     myEdgesExchnged = Standard_True;
 115   }
 116   else
 117   {
 118     myEdge1 = theEdge1;
 119     myEdge2 = theEdge2;
 120   }
 121
 122   Handle(Geom_Curve) aCurve1 = BRep_Tool::Curve(myEdge1, myStart1, myEnd1);
 123   Handle(Geom_Curve) aCurve2 = BRep_Tool::Curve(myEdge2, myStart2, myEnd2);
 124
 125   myCurve1 = GeomProjLib::Curve2d(aCurve1, myStart1, myEnd1, myPlane);
 126   myCurve2 = GeomProjLib::Curve2d(aCurve2, myStart2, myEnd2, myPlane);
 127
 128   while (myCurve1->IsPeriodic() && myStart1 >= myEnd1)
 129     myEnd1 += myCurve1->Period();
 130   while (myCurve2->IsPeriodic() && myStart2 >= myEnd2)
 131     myEnd2 += myCurve2->Period();
 132
 133   if (aBAC1.GetType() == aBAC2.GetType())
 134   {
 135     if (myEnd2 - myStart2 < myEnd1 - myStart1)
 136     { // first curve must be parametrically shorter
 137       TopoDS_Edge anEdge = myEdge1;
 138       myEdge1 = myEdge2;
 139       myEdge2 = anEdge;
 140       Handle(Geom2d_Curve) aCurve = myCurve1;
 141       myCurve1 = myCurve2;
 142       myCurve2 = aCurve;
 143       Standard_Real a = myStart1;
 144       myStart1 = myStart2;
 145       myStart2 = a;
 146       a = myEnd1;
 147       myEnd1 = myEnd2;
 148       myEnd2 = a;
 149       myEdgesExchnged = Standard_True;
 150     }
 151   }
 152 }
 153
 154 //=======================================================================
 155 //function : isRadiusIntersected
 156 //purpose  : local function
 157 //=======================================================================
 158 static Standard_Boolean isRadiusIntersected(const Handle(Geom2d_Curve)& theCurve,
 159                                             const gp_Pnt2d theStart,
 160                                             const gp_Pnt2d theEnd,
 161                                             const Standard_Boolean theStartConnected)
 162 {
 163   const Standard_Real aTol = Precision::Confusion();
 164   const Standard_Real anAngTol = Precision::Angular();
 165   Handle(Geom2d_Line) aRadiusLine = new Geom2d_Line (theStart, gp_Dir2d(gp_Vec2d(theStart, theEnd)));
 166   Geom2dAPI_InterCurveCurve anInter (theCurve, aRadiusLine, aTol);
 167   Standard_Integer a;
 168   gp_Pnt2d aPoint;
 169   for(a = anInter.NbPoints(); a > 0; a--)
 170   {
 171     aPoint = anInter.Point(a);
 172     if ( aPoint.Distance(theStart) < aTol && !theStartConnected )
 173       return Standard_True;
 174     if (aPoint.Distance(theEnd) < aTol * 200)
 175       return Standard_True;
 176     if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
 177       return Standard_True;
 178   }
 179   Handle(Geom2d_Curve) aCurve;
 180   for(a = anInter.NbSegments(); a > 0; a--)
 181   {
 182     // Porting to DEV version of OCCT 10.02.2017 BEGIN
 183     Standard_NotImplemented::Raise("The treatment of tangential intersection is not implemented");
 184     // Porting to DEV version of OCCT 10.02.2017 END
 185   }
 186   return Standard_False;
 187 }
 188
 189
 190 //=======================================================================
 191 //function : fillPoint
 192 //purpose  :
 193 //=======================================================================
 194 void GEOMImpl_Fillet1d::fillPoint(GEOMImpl_Fillet1dPoint* thePoint)
 195 {
 196   gp_Pnt2d aPoint;
 197   gp_Vec2d aVec;
 198   const Standard_Real aTol = Precision::Confusion();
 199   myCurve1->D1(thePoint->GetParam(), aPoint, aVec);
 200   if (aVec.SquareMagnitude() < aTol)
 201     return;
 202
 203   gp_Vec2d aPerp(((myStartSide)?-1:1) * aVec.Y(), ((myStartSide)?1:-1) * aVec.X());
 204   aPerp.Normalize();
 205   aPerp.Multiply(myRadius);
 206   gp_Pnt2d aCenter = aPoint.Translated(aPerp);
 207   thePoint->SetCenter(aCenter);
 208
 209   // on the intersection point
 210   Standard_Boolean aValid = Standard_True;
 211   Geom2dAPI_ProjectPointOnCurve aProjInt(aPoint, myCurve2);
 212   if (aProjInt.NbPoints() && aPoint.Distance(aProjInt.NearestPoint()) < aTol)
 213     aValid = Standard_False;
 214   else
 215     aValid = !isRadiusIntersected(myCurve2, aPoint, aCenter, Standard_True);
 216
 217   Geom2dAPI_ProjectPointOnCurve aProj(aCenter, myCurve2);
 218   Standard_Integer a, aNB = aProj.NbPoints();
 219   for(a = aNB; a > 0; a--)
 220   {
 221     if (aPoint.Distance(aProj.Point(a)) < aTol)
 222       continue;
 223
 224     Standard_Boolean aValid2 = aValid;
 225     if (aValid2)
 226       aValid2 = !isRadiusIntersected(myCurve1, aCenter, aProj.Point(a), Standard_False);
 227
 228     // checking the right parameter
 229     Standard_Real aParam = aProj.Parameter(a);
 230     while(myCurve2->IsPeriodic() && aParam < myStart2)
 231       aParam += myCurve2->Period();
 232
 233     thePoint->AddValue(aProj.Distance(a) * aProj.Distance(a) - myRadius * myRadius,
 234                        (aParam >= myStart2 && aParam <= myEnd2 && aValid2));
 235     if (fabs(fabs(aProj.Distance(a)) - myRadius) < aTol)
 236       thePoint->SetParam2(aParam);
 237   }
 238 }
 239
 240 //=======================================================================
 241 //function : fillDiff
 242 //purpose  :
 243 //=======================================================================
 244 void GEOMImpl_Fillet1d::fillDiff(GEOMImpl_Fillet1dPoint* thePoint, Standard_Real theDiffStep, Standard_Boolean theFront)
 245 {
 246   GEOMImpl_Fillet1dPoint* aDiff =
 247     new GEOMImpl_Fillet1dPoint(thePoint->GetParam() + (theFront?(theDiffStep):(-theDiffStep)));
 248   fillPoint(aDiff);
 249   if (!thePoint->ComputeDifference(aDiff))
 250   {
 251     aDiff->SetParam(thePoint->GetParam() + (theFront?(-theDiffStep):(theDiffStep)));
 252     fillPoint(aDiff);
 253     thePoint->ComputeDifference(aDiff);
 254   }
 255   delete aDiff;
 256 }
 257
 258 //=======================================================================
 259 //function : Perform
 260 //purpose  :
 261 //=======================================================================
 262 Standard_Boolean GEOMImpl_Fillet1d::Perform(const Standard_Real theRadius)
 263 {
 264   myDegreeOfRecursion = 0;
 265   myResultParams.Clear();
 266   myResultOrientation.Clear();
 267
 268   Standard_Real aNBSteps = 100;
 269   Geom2dAdaptor_Curve aGAC(myCurve1);
 270   switch (aGAC.GetType())
 271   {
 272     case GeomAbs_Line:
 273       aNBSteps = 2;
 274       break;
 275     case GeomAbs_Circle:
 276       aNBSteps = 4;
 277       break;
 278     case GeomAbs_Ellipse:
 279       aNBSteps = 5;
 280       break;
 281     case GeomAbs_BezierCurve:
 282     case GeomAbs_BSplineCurve:
 283       aNBSteps = 2 + aGAC.Degree() * aGAC.NbPoles();
 284       break;
 285     default: // unknown: maximum
 286       aNBSteps = 100;
 287   }
 288
 289   myRadius = theRadius;
 290
 291   // Compute the intervals.
 292   const Standard_Real aTol = Precision::Confusion();
 293   Geom2dAPI_InterCurveCurve anAPIInter(myCurve1, myCurve2, aTol);
 294   const Geom2dInt_GInter &anInter = anAPIInter.Intersector();
 295   Standard_Integer aNb = anInter.NbPoints();
 296   Standard_Integer i;
 297   TColStd_ListOfReal aParams;
 298   TColStd_ListIteratorOfListOfReal anIter;
 299
 300   // Treat intersection points.
 301   for(i = 1; i <= aNb; i++) {
 302     const IntRes2d_IntersectionPoint &aPoint = anInter.Point(i);
 303     Standard_Real                     aParam = aPoint.ParamOnFirst();
 304
 305     // Adjust parameter on periodic curve.
 306     if (myCurve1->IsPeriodic()) {
 307       aParam = ElCLib::InPeriod
 308         (aParam, myStart1, myStart1 + myCurve1->Period());
 309     }
 310
 311     if (aParam > myStart1 + aTol && aParam < myEnd1 - aTol) {
 312       // Add the point in the list in increasing order.
 313       for(anIter.Initialize(aParams); anIter.More(); anIter.Next()) {
 314         if (anIter.Value() > aParam) {
 315           aParams.InsertBefore(aParam, anIter);
 316           break;
 317         }
 318       }
 319
 320       if (!anIter.More()) {
 321         aParams.Append(aParam);
 322       }
 323     }
 324   }
 325
 326   // Treat intersection segments.
 327   aNb = anInter.NbSegments();
 328
 329   for(i = 1; i <= aNb; i++) {
 330     const IntRes2d_IntersectionSegment &aSegment = anInter.Segment(i);
 331
 332     if (aSegment.HasFirstPoint() && aSegment.HasLastPoint()) {
 333       Standard_Real aParam1 = aSegment.FirstPoint().ParamOnFirst();
 334       Standard_Real aParam2 = aSegment.LastPoint().ParamOnFirst();
 335
 336       // Adjust parameters on periodic curve.
 337       if (myCurve1->IsPeriodic()) {
 338         ElCLib::AdjustPeriodic(myStart1, myStart1 + myCurve1->Period(),
 339                                aTol, aParam1, aParam2);
 340       }
 341
 342       if (aParam1 > myStart1 + aTol && aParam1 < myEnd1 - aTol &&
 343           aParam2 > myStart1 + aTol && aParam2 < myEnd1 - aTol) {
 344         // Add the point in the list in increasing order.
 345         const Standard_Real aParam = 0.5*(aParam1 + aParam2);
 346
 347         for(anIter.Initialize(aParams); anIter.More(); anIter.Next()) {
 348           if (anIter.Value() > aParam) {
 349             aParams.InsertBefore(aParam, anIter);
 350             break;
 351           }
 352         }
 353
 354         if (!anIter.More()) {
 355           aParams.Append(aParam);
 356         }
 357       }
 358     }
 359   }
 360
 361   // Add start and end parameters to the list.
 362   aParams.Prepend(myStart1);
 363   aParams.Append(myEnd1);
 364   anIter.Initialize(aParams);
 365
 366   // Perform each interval.
 367   Standard_Real aStart = anIter.Value();
 368
 369   for (anIter.Next(); anIter.More(); anIter.Next()) {
 370     const Standard_Real anEnd = anIter.Value();
 371
 372     // Perform the interval.
 373     performInterval(aStart, anEnd, aNBSteps);
 374     aStart = anEnd;
 375   }
 376
 377   if (myResultParams.Extent())
 378     return Standard_True;
 379
 380   return Standard_False;
 381 }
 382
 383 //=======================================================================
 384 //function : performInterval
 385 //purpose  :
 386 //=======================================================================
 387 void GEOMImpl_Fillet1d::performInterval(const Standard_Real theStart,
 388                                         const Standard_Real theEnd,
 389                                         const Standard_Integer theNBSteps)
 390 {
 391   Standard_Real aParam, aStep, aDStep;
 392   aStep = (theEnd - theStart) / theNBSteps;
 393   aDStep = aStep/1000.;
 394
 395   Standard_Integer aCycle;
 396   for(aCycle = 2, myStartSide = Standard_False; aCycle; myStartSide = !myStartSide, aCycle--)
 397   {
 398     GEOMImpl_Fillet1dPoint *aLeft = NULL, *aRight = NULL;
 399
 400     for(aParam = theStart + aStep; aParam < theEnd || fabs(theEnd - aParam) < Precision::Confusion(); aParam += aStep)
 401     {
 402       if (!aLeft)
 403       {
 404         aLeft = new GEOMImpl_Fillet1dPoint(aParam - aStep);
 405         fillPoint(aLeft);
 406         fillDiff(aLeft, aDStep, Standard_True);
 407       }
 408
 409       aRight = new GEOMImpl_Fillet1dPoint(aParam);
 410       fillPoint(aRight);
 411       fillDiff(aRight, aDStep, Standard_False);
 412
 413       aLeft->FilterPoints(aRight);
 414       performNewton(aLeft, aRight);
 415
 416       delete aLeft;
 417       aLeft = aRight;
 418     }
 419     delete aLeft;
 420   }
 421 }
 422
 423 //=======================================================================
 424 //function : processPoint
 425 //purpose  :
 426 //=======================================================================
 427 Standard_Boolean GEOMImpl_Fillet1d::processPoint(GEOMImpl_Fillet1dPoint* theLeft,
 428                                                  GEOMImpl_Fillet1dPoint* theRight,
 429                                                  Standard_Real           theParameter)
 430 {
 431   if (theParameter >= theLeft->GetParam() && theParameter < theRight->GetParam())
 432   {
 433     Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
 434     if (theParameter - theLeft->GetParam() < aDX / 100.)
 435     {
 436       theParameter = theLeft->GetParam() + aDX / 100.;
 437     }
 438     if (theRight->GetParam() - theParameter < aDX / 100.)
 439     {
 440       theParameter = theRight->GetParam() - aDX / 100.;
 441     }
 442
 443     // Protection on infinite loop.
 444     myDegreeOfRecursion++;
 445     Standard_Real diffx = 0.001 * aDX;
 446     if (myDegreeOfRecursion > 1000)
 447     {
 448         diffx *= 10.0;
 449         if (myDegreeOfRecursion > 10000)
 450         {
 451             diffx *= 10.0;
 452             if (myDegreeOfRecursion > 100000)
 453             {
 454                 return Standard_True;
 455             }
 456         }
 457     }
 458
 459     GEOMImpl_Fillet1dPoint* aPoint1 = theLeft->Copy();
 460     GEOMImpl_Fillet1dPoint* aPoint2 = new GEOMImpl_Fillet1dPoint(theParameter);
 461     fillPoint(aPoint2);
 462     fillDiff(aPoint2, diffx, Standard_True);
 463
 464     aPoint1->FilterPoints(aPoint2);
 465     performNewton(aPoint1, aPoint2);
 466     aPoint2->FilterPoints(theRight);
 467     performNewton(aPoint2, theRight);
 468
 469     delete aPoint1;
 470     delete aPoint2;
 471     return Standard_True;
 472   }
 473
 474   return Standard_False;
 475 }
 476
 477 //=======================================================================
 478 //function : performNewton
 479 //purpose  :
 480 //=======================================================================
 481 void GEOMImpl_Fillet1d::performNewton(GEOMImpl_Fillet1dPoint* theLeft,
 482                                       GEOMImpl_Fillet1dPoint* theRight)
 483 {
 484   Standard_Integer a;
 485   // check the left: if this is solution store it and remove it from the list of researching points of theLeft
 486   a = theLeft->HasSolution(myRadius);
 487   if (a)
 488   {
 489     if (theLeft->IsValid(a))
 490     {
 491       myResultParams.Append(theLeft->GetParam());
 492       myResultOrientation.Append(myStartSide);
 493     }
 494     return;
 495   }
 496
 497   Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
 498   if ( aDX < Precision::Confusion() / 1000000.)
 499   {
 500     a = theRight->HasSolution(myRadius);
 501     if (a)
 502       if (theRight->IsValid(a))
 503       {
 504         myResultParams.Append(theRight->GetParam());
 505         myResultOrientation.Append(myStartSide);
 506       }
 507     return;
 508   }
 509
 510   for(a = 1; a <= theLeft->GetNBValues(); a++)
 511   {
 512     Standard_Integer aNear = theLeft->GetNear(a);
 513
 514     Standard_Real aA = (theRight->GetDiff(aNear) - theLeft->GetDiff(a)) / aDX;
 515     Standard_Real aB = theLeft->GetDiff(a) - aA * theLeft->GetParam();
 516     Standard_Real aC = theLeft->GetValue(a) - theLeft->GetDiff(a) * theLeft->GetParam() +
 517                        aA * theLeft->GetParam() * theLeft->GetParam() / 2.0;
 518     Standard_Real aDet = aB * aB - 2.0 * aA * aC;
 519
 520     if ( fabs(aDet) < gp::Resolution() )
 521       continue;
 522
 523     if (fabs(aA) < Precision::Confusion())
 524     { // linear case
 525       if (fabs(aB) > 10e-20)
 526       {
 527         Standard_Real aX0 = - aC / aB; // use extremum
 528         if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
 529           processPoint(theLeft, theRight, aX0);
 530       }
 531       else
 532       {
 533         processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 534       }
 535     }
 536     else
 537     {
 538       if (fabs(aB) > fabs(aDet * 1000000.))
 539       {  // possible floating point operations accuracy errors
 540         processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 541       }
 542       else
 543       {
 544         if (aDet > 0)
 545         { // two solutions
 546           aDet = sqrt(aDet);
 547           Standard_Boolean aRes = processPoint(theLeft, theRight, (- aB + aDet) / aA);
 548           if (!aRes)
 549             aRes = processPoint(theLeft, theRight, (- aB - aDet) / aA);
 550           if (!aRes)
 551             processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 552         }
 553         else
 554         {
 555           Standard_Real aX0 = - aB / aA; // use extremum
 556           if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
 557             processPoint(theLeft, theRight, aX0);
 558           else
 559             processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 560         }
 561       }
 562     }
 563   }
 564 }
 565
 566 //=======================================================================
 567 //function : Result
 568 //purpose  :
 569 //=======================================================================
 570 TopoDS_Edge GEOMImpl_Fillet1d::Result(const gp_Pnt& thePoint,
 571                                       TopoDS_Edge& theEdge1,
 572                                       TopoDS_Edge& theEdge2)
 573 {
 574   TopoDS_Edge aResult;
 575   gp_Pnt2d aTargetPoint2d;
 576   Standard_Real aX, aY;
 577   ElSLib::PlaneParameters(myPlane->Pln().Position(), thePoint, aX, aY);
 578   aTargetPoint2d.SetCoord(aX, aY);
 579
 580   // choose the nearest circle
 581   Standard_Real aDistance = 0., aP; // todo: aDistance must be explicitly initialized to avoid warning (see below)
 582   GEOMImpl_Fillet1dPoint *aNearest;
 583   Standard_Integer a;
 584   TColStd_ListIteratorOfListOfReal anIter(myResultParams);
 585   for(aNearest = NULL, a = 1; anIter.More(); anIter.Next(), a++)
 586   {
 587     myStartSide = (myResultOrientation.Value(a)) ? Standard_True : Standard_False;
 588     GEOMImpl_Fillet1dPoint *aPoint = new GEOMImpl_Fillet1dPoint(anIter.Value());
 589     fillPoint(aPoint);
 590     if (!aPoint->HasSolution(myRadius))
 591       continue;
 592     aP = fabs(aPoint->GetCenter().Distance(aTargetPoint2d) - myRadius);
 593     if (!aNearest || aP < aDistance) // todo: aDistance must be explicitly initialized to avoid warning (see above)
 594     {
 595       aNearest = aPoint;
 596       aDistance = aP;
 597     }
 598     else
 599     {
 600       delete aPoint;
 601     }
 602    }
 603
 604   if (!aNearest)
 605      return aResult;
 606
 607   // create circle edge
 608   gp_Pnt aCenter = ElSLib::PlaneValue(aNearest->GetCenter().X(),
 609                                       aNearest->GetCenter().Y(),
 610                                       myPlane->Pln().Position());
 611   Handle(Geom_Circle) aCircle =
 612     new Geom_Circle(gp_Ax2(aCenter, myPlane->Pln().Axis().Direction()), myRadius);
 613   gp_Pnt2d aPoint2d1, aPoint2d2;
 614   myCurve1->D0(aNearest->GetParam(), aPoint2d1);
 615   myCurve2->D0(aNearest->GetParam2(), aPoint2d2);
 616   gp_Pnt aPoint1 = ElSLib::PlaneValue(aPoint2d1.X(), aPoint2d1.Y(), myPlane->Pln().Position());
 617   gp_Pnt aPoint2 = ElSLib::PlaneValue(aPoint2d2.X(), aPoint2d2.Y(), myPlane->Pln().Position());
 618
 619   GeomAPI_ProjectPointOnCurve aProj(thePoint, aCircle);
 620   Standard_Real aTarGetParam = aProj.LowerDistanceParameter();
 621   gp_Pnt aPointOnCircle = aProj.NearestPoint();
 622
 623   // Check extrema point manually, because there is a bug in Open CASCADE
 624   //  in calculation of nearest point to a circle near the parameter 0.0
 625   gp_Pnt p0 = ElCLib::Value(0.0, aCircle->Circ());
 626   if (p0.Distance(thePoint) < aPointOnCircle.Distance(thePoint))
 627   {
 628      aTarGetParam = 0.0;
 629      aPointOnCircle = p0;
 630   }
 631
 632   aProj.Perform(aPoint1);
 633   Standard_Real aParam1 = aProj.LowerDistanceParameter();
 634     aProj.Perform(aPoint2);
 635   Standard_Real aParam2 = aProj.LowerDistanceParameter();
 636   Standard_Boolean aIsOut = ((aParam1 < aTarGetParam && aParam2 < aTarGetParam) ||
 637                              (aParam1 > aTarGetParam && aParam2 > aTarGetParam));
 638   if (aParam1 > aParam2)
 639     aIsOut = !aIsOut;
 640   BRepBuilderAPI_MakeEdge aBuilder(aCircle->Circ(),
 641                                    aIsOut ? aParam2 : aParam1,
 642                                    aIsOut? aParam1 : aParam2);
 643   aResult = aBuilder.Edge();
 644
 645   // divide edges
 646   Standard_Real aStart, anEnd;
 647   Handle(Geom_Curve) aCurve = BRep_Tool::Curve(myEdge1, aStart, anEnd);
 648   gp_Vec aDir;
 649   aCurve->D1(aNearest->GetParam(), aPoint1, aDir);
 650
 651   gp_Vec aCircleDir;
 652   aCircle->D1(aParam1, aPoint1, aCircleDir);
 653   if ((aCircleDir.Angle(aDir) > M_PI / 2.0) ^ aIsOut)
 654     aStart = aNearest->GetParam();
 655   else
 656     anEnd = aNearest->GetParam();
 657
 658   if (IsDivideEdge(myEdge1, aStart, anEnd))
 659   {
 660       //Divide edge
 661       BRepBuilderAPI_MakeEdge aDivider1(aCurve, aStart, anEnd);
 662       if (myEdgesExchnged)
 663         theEdge2 = aDivider1.Edge();
 664       else
 665         theEdge1 = aDivider1.Edge();
 666   }
 667
 668   aCurve = BRep_Tool::Curve(myEdge2, aStart, anEnd);
 669   aCurve->D1(aNearest->GetParam2(), aPoint2, aDir);
 670
 671   aCircle->D1(aParam2, aPoint2, aCircleDir);
 672   if ((aCircleDir.Angle(aDir) > M_PI / 2.0) ^ (!aIsOut))
 673     aStart = aNearest->GetParam2();
 674   else
 675     anEnd = aNearest->GetParam2();
 676
 677   if (IsDivideEdge(myEdge2, aStart, anEnd))
 678   {
 679       BRepBuilderAPI_MakeEdge aDivider2(aCurve, aStart, anEnd);
 680       if (myEdgesExchnged)
 681         theEdge1 = aDivider2.Edge();
 682       else
 683         theEdge2 = aDivider2.Edge();
 684   }
 685
 686   delete aNearest;
 687   return aResult;
 688 }
 689
 690 //=======================================================================
 691 //function : AddValue
 692 //purpose  :
 693 //=======================================================================
 694 void GEOMImpl_Fillet1dPoint::AddValue(Standard_Real theValue, Standard_Boolean theValid)
 695 {
 696   Standard_Integer a;
 697   for(a = 1; a <= myV.Length(); a++)
 698   {
 699     if (theValue < myV.Value(a))
 700     {
 701       myV.InsertBefore(a, theValue);
 702       myValid.InsertBefore(a, (Standard_Integer)theValid);
 703       return;
 704     }
 705   }
 706   myV.Append(theValue);
 707   myValid.Append((Standard_Integer)theValid);
 708 }
 709
 710 //=======================================================================
 711 //function : ComputeDifference
 712 //purpose  :
 713 //=======================================================================
 714 Standard_Boolean GEOMImpl_Fillet1dPoint::ComputeDifference(GEOMImpl_Fillet1dPoint* thePoint)
 715 {
 716   Standard_Integer a;
 717   Standard_Boolean aDiffsSet = (myD.Length() != 0);
 718   Standard_Real aDX = thePoint->GetParam() - myParam, aDY = 0.; // todo: aDY must be explicitly initialized to avoid warning (see below)
 719   if (thePoint->myV.Length() == myV.Length())
 720   { // absolutely the same points
 721     for(a = 1; a <= myV.Length(); a++)
 722     {
 723       aDY = thePoint->myV.Value(a) - myV.Value(a);
 724       if ( aDiffsSet )
 725         myD.SetValue(a, fabs(aDX) > gp::Resolution() ? (aDY/aDX) : 0);
 726       else
 727         myD.Append( fabs(aDX) > gp::Resolution() ? (aDY/aDX) : 0);
 728     }
 729     return Standard_True;
 730   }
 731   // between the diffeerent points searching for nearest analogs
 732   Standard_Integer b;
 733   for(a = 1; a <= myV.Length(); a++)
 734   {
 735     for(b = 1; b <= thePoint->myV.Length(); b++)
 736     {
 737       if (b == 1 || fabs(thePoint->myV.Value(b) - myV.Value(a)) < fabs(aDY))
 738         aDY = thePoint->myV.Value(b) - myV.Value(a);
 739     }
 740     if (aDiffsSet)
 741     {
 742       if ( fabs(aDX) > gp::Resolution() && fabs(aDY / aDX) < fabs(myD.Value(a)))
 743         myD.SetValue(a, aDY / aDX);
 744       else
 745         myD.SetValue(a, 0);
 746     }
 747     else
 748     {
 749       myD.Append( fabs(aDX) > gp::Resolution() ? aDY/aDX : 0); // todo: aDY must be explicitly initialized to avoid warning (see above)
 750     }
 751   }
 752
 753   return Standard_False;
 754 }
 755
 756 //=======================================================================
 757 //function : FilterPoints
 758 //purpose  :
 759 //=======================================================================
 760 void GEOMImpl_Fillet1dPoint::FilterPoints(GEOMImpl_Fillet1dPoint* thePoint)
 761 {
 762   Standard_Integer a, b;
 763   TColStd_SequenceOfReal aDiffs;
 764   Standard_Real aY, aY2, aDX = thePoint->GetParam() - myParam;
 765   for(a = 1; a <= myV.Length(); a++)
 766   {
 767     // searching for near point from thePoint
 768     Standard_Integer aNear = 0;
 769     Standard_Real aDiff = aDX * 10000.;
 770     aY = myV.Value(a) + myD.Value(a) * aDX;
 771     for(b = 1; b <= thePoint->myV.Length(); b++)
 772     {
 773       // calculate hypothesis value of the Y2 with the constant first and second derivative
 774       aY2 = aY + aDX * (thePoint->myD.Value(b) - myD.Value(a)) / 2.0;
 775       if (aNear == 0 || fabs(aY2 - thePoint->myV.Value(b)) < fabs(aDiff))
 776       {
 777         aNear = b;
 778         aDiff = aY2 - thePoint->myV.Value(b);
 779       }
 780     }//for b...
 781
 782     if (aNear)
 783     {
 784       if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
 785       {// the same sign at the same sides of the interval
 786         if (myV.Value(a) * myD.Value(a) > 0)
 787         {
 788           if (fabs(myD.Value(a)) > Precision::Confusion())
 789             aNear = 0;
 790         }
 791         else
 792         {
 793           if (fabs(myV.Value(a)) > fabs(thePoint->myV.Value(aNear)))
 794             if (thePoint->myV.Value(aNear) * thePoint->myD.Value(aNear) < 0 &&
 795                 fabs(thePoint->myD.Value(aNear)) > Precision::Confusion())
 796             {
 797               aNear = 0;
 798             }
 799         }
 800       }
 801     }
 802
 803     if (aNear)
 804     {
 805       if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
 806       {
 807         if ((myV.Value(a) + myD.Value(a) * aDX) * myV.Value(a) > Precision::Confusion() &&
 808         (thePoint->myV.Value(aNear) + thePoint->myD.Value(aNear) * aDX) * thePoint->myV.Value(aNear) > Precision::Confusion())
 809         {
 810           aNear = 0;
 811         }
 812       }
 813     }
 814
 815     if (aNear)
 816     {
 817       if (  fabs(aDX) < gp::Resolution() || fabs(aDiff / aDX) > 1.e+7)
 818       {
 819         aNear = 0;
 820       }
 821     }
 822
 823     if (aNear == 0)
 824     {  // there is no near: remove it from the list
 825       myV.Remove(a);
 826       myD.Remove(a);
 827       myValid.Remove(a);
 828       a--;
 829     }
 830     else
 831     {
 832       Standard_Boolean aFound = Standard_False;
 833       for(b = 1; b <= myNear.Length(); b++)
 834       {
 835         if (myNear.Value(b) == aNear)
 836         {
 837           if (fabs(aDiffs.Value(b)) < fabs(aDiff))
 838           { // return this 'near'
 839             aFound = Standard_True;
 840             myV.Remove(a);
 841             myD.Remove(a);
 842             myValid.Remove(a);
 843             a--;
 844             break;
 845           }
 846           else
 847           { // remove the old 'near'
 848             myV.Remove(b);
 849             myD.Remove(b);
 850             myValid.Remove(b);
 851             myNear.Remove(b);
 852             aDiffs.Remove(b);
 853             a--;
 854             break;
 855           }
 856         }
 857       }//for b...
 858       if (!aFound)
 859       {
 860         myNear.Append(aNear);
 861         aDiffs.Append(aDiff);
 862       }
 863     }
 864   }//for a...
 865 }
 866
 867 //=======================================================================
 868 //function : Copy
 869 //purpose  :
 870 //=======================================================================
 871 GEOMImpl_Fillet1dPoint* GEOMImpl_Fillet1dPoint::Copy()
 872 {
 873   GEOMImpl_Fillet1dPoint* aCopy = new GEOMImpl_Fillet1dPoint(myParam);
 874   Standard_Integer a;
 875   for(a = 1; a <= myV.Length(); a++)
 876   {
 877     aCopy->myV.Append(myV.Value(a));
 878     aCopy->myD.Append(myD.Value(a));
 879     aCopy->myValid.Append(myValid.Value(a));
 880   }
 881   return aCopy;
 882 }
 883
 884 //=======================================================================
 885 //function : HasSolution
 886 //purpose  :
 887 //=======================================================================
 888 Standard_Integer GEOMImpl_Fillet1dPoint::HasSolution(const Standard_Real theRadius)
 889 {
 890   Standard_Integer a;
 891   for(a = 1; a <= myV.Length(); a++)
 892   {
 893     if (fabs(sqrt(fabs(fabs(myV.Value(a)) + theRadius * theRadius)) - theRadius) < Precision::Confusion() / 10.)
 894       return a;
 895   }
 896   return 0;
 897 }
 898
 899 //=======================================================================
 900 //function : RemoveSolution
 901 //purpose  :
 902 //=======================================================================
 903 void GEOMImpl_Fillet1dPoint::RemoveSolution(Standard_Integer theIndex)
 904 {
 905   myV.Remove(theIndex);
 906   myD.Remove(theIndex);
 907   myValid.Remove(theIndex);
 908   myNear.Remove(theIndex);
 909 }