src/GEOMImpl/GEOMImpl_Fillet1d.cxx

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