src/GEOMImpl/GEOMImpl_Fillet1d.cxx

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