GEOMImpl_Fillet1d.cxx

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