src/GEOMImpl/GEOMImpl_Fillet1d.cxx

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