src/GEOMImpl/GEOMImpl_Fillet1d.cxx

   1 //  Copyright (C) 2007-2010  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
  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
  40 /**
  41  * class GEOMImpl_Fillet1d
  42  */
  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   myResultParams.Clear();
 227   myResultOrientation.Clear();
 228
 229   Standard_Real aNBSteps = 100;
 230   Geom2dAdaptor_Curve aGAC(myCurve1);
 231   switch (aGAC.GetType())
 232   {
 233     case GeomAbs_Line:
 234       aNBSteps = 2;
 235       break;
 236     case GeomAbs_Circle:
 237       aNBSteps = 4;
 238       break;
 239     case GeomAbs_Ellipse:
 240       aNBSteps = 5;
 241       break;
 242     case GeomAbs_BezierCurve:
 243     case GeomAbs_BSplineCurve:
 244       aNBSteps = 2 + aGAC.Degree() * aGAC.NbPoles();
 245       break;
 246     default: // unknown: maximum
 247       aNBSteps = 100;
 248   }
 249
 250   myRadius = theRadius;
 251   Standard_Real aParam, aStep, aDStep;
 252   aStep = (myEnd1 - myStart1) / aNBSteps;
 253   aDStep = aStep/1000.;
 254
 255   Standard_Integer aCycle;
 256   for(aCycle = 2, myStartSide = Standard_False; aCycle; myStartSide = !myStartSide, aCycle--)
 257   {
 258     GEOMImpl_Fillet1dPoint *aLeft = NULL, *aRight = NULL;
 259
 260     for(aParam = myStart1 + aStep; aParam < myEnd1 || fabs(myEnd1 - aParam) < Precision::Confusion(); aParam += aStep)
 261     {
 262       if (!aLeft)
 263       {
 264         aLeft = new GEOMImpl_Fillet1dPoint(aParam - aStep);
 265         fillPoint(aLeft);
 266         fillDiff(aLeft, aDStep, Standard_True);
 267       }
 268
 269       aRight = new GEOMImpl_Fillet1dPoint(aParam);
 270       fillPoint(aRight);
 271       fillDiff(aRight, aDStep, Standard_False);
 272
 273       aLeft->FilterPoints(aRight);
 274       performNewton(aLeft, aRight);
 275
 276       delete aLeft;
 277       aLeft = aRight;
 278     }
 279     delete aLeft;
 280   }
 281
 282   if (myResultParams.Extent())
 283     return Standard_True;
 284
 285   return Standard_False;
 286 }
 287
 288 //=======================================================================
 289 //function : processPoint
 290 //purpose  :
 291 //=======================================================================
 292 Standard_Boolean GEOMImpl_Fillet1d::processPoint(GEOMImpl_Fillet1dPoint* theLeft,
 293                                                  GEOMImpl_Fillet1dPoint* theRight,
 294                                                  Standard_Real           theParameter)
 295 {
 296   if (theParameter >= theLeft->GetParam() && theParameter < theRight->GetParam())
 297   {
 298     Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
 299     if (theParameter - theLeft->GetParam() < aDX / 100.)
 300     {
 301       theParameter = theLeft->GetParam() + aDX / 100.;
 302     }
 303     if (theRight->GetParam() - theParameter < aDX / 100.)
 304     {
 305       theParameter = theRight->GetParam() - aDX / 100.;
 306     }
 307
 308     GEOMImpl_Fillet1dPoint* aPoint1 = theLeft->Copy();
 309     GEOMImpl_Fillet1dPoint* aPoint2 = new GEOMImpl_Fillet1dPoint(theParameter);
 310     fillPoint(aPoint2);
 311     fillDiff(aPoint2, aDX / 1000., Standard_True);
 312
 313     aPoint1->FilterPoints(aPoint2);
 314     performNewton(aPoint1, aPoint2);
 315     aPoint2->FilterPoints(theRight);
 316     performNewton(aPoint2, theRight);
 317
 318     delete aPoint1;
 319     delete aPoint2;
 320     return Standard_True;
 321   }
 322
 323   return Standard_False;
 324 }
 325
 326 //=======================================================================
 327 //function : performNewton
 328 //purpose  :
 329 //=======================================================================
 330 void GEOMImpl_Fillet1d::performNewton(GEOMImpl_Fillet1dPoint* theLeft,
 331                                       GEOMImpl_Fillet1dPoint* theRight)
 332 {
 333   Standard_Integer a;
 334   // check the left: if this is solution store it and remove it from the list of researching points of theLeft
 335   a = theLeft->HasSolution(myRadius);
 336   if (a)
 337   {
 338     if (theLeft->IsValid(a))
 339     {
 340       myResultParams.Append(theLeft->GetParam());
 341       myResultOrientation.Append(myStartSide);
 342     }
 343     return;
 344   }
 345
 346   Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
 347   if ( aDX < Precision::Confusion() / 1000000.)
 348   {
 349     a = theRight->HasSolution(myRadius);
 350     if (a)
 351       if (theRight->IsValid(a))
 352       {
 353         myResultParams.Append(theRight->GetParam());
 354         myResultOrientation.Append(myStartSide);
 355       }
 356     return;
 357   }
 358
 359   for(a = 1; a <= theLeft->GetNBValues(); a++)
 360   {
 361     Standard_Integer aNear = theLeft->GetNear(a);
 362
 363     Standard_Real aA = (theRight->GetDiff(aNear) - theLeft->GetDiff(a)) / aDX;
 364     Standard_Real aB = theLeft->GetDiff(a) - aA * theLeft->GetParam();
 365     Standard_Real aC = theLeft->GetValue(a) - theLeft->GetDiff(a) * theLeft->GetParam() +
 366                        aA * theLeft->GetParam() * theLeft->GetParam() / 2.0;
 367     Standard_Real aDet = aB * aB - 2.0 * aA * aC;
 368
 369     if ( fabs(aDet) < gp::Resolution() )
 370       continue;
 371
 372     if (fabs(aA) < Precision::Confusion())
 373     { // linear case
 374       if (fabs(aB) > 10e-20)
 375       {
 376         Standard_Real aX0 = - aC / aB; // use extremum
 377         if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
 378           processPoint(theLeft, theRight, aX0);
 379       }
 380       else
 381       {
 382         processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 383       }
 384     }
 385     else
 386     {
 387       if (fabs(aB) > fabs(aDet * 1000000.))
 388       {  // possible floating point operations accurancy errors
 389         processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 390       }
 391       else
 392       {
 393         if (aDet > 0)
 394         { // two solutions
 395           aDet = sqrt(aDet);
 396           Standard_Boolean aRes = processPoint(theLeft, theRight, (- aB + aDet) / aA);
 397           if (!aRes)
 398             aRes = processPoint(theLeft, theRight, (- aB - aDet) / aA);
 399           if (!aRes)
 400             processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 401         }
 402         else
 403         {
 404           Standard_Real aX0 = - aB / aA; // use extremum
 405           if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
 406             processPoint(theLeft, theRight, aX0);
 407           else
 408             processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 409         }
 410       }
 411     }
 412   }
 413 }
 414
 415 //=======================================================================
 416 //function : Result
 417 //purpose  :
 418 //=======================================================================
 419 TopoDS_Edge GEOMImpl_Fillet1d::Result(const gp_Pnt& thePoint,
 420                                       TopoDS_Edge& theEdge1,
 421                                       TopoDS_Edge& theEdge2)
 422 {
 423   TopoDS_Edge aResult;
 424   gp_Pnt2d aTargetPoint2d;
 425   Standard_Real aX, aY;
 426   ElSLib::PlaneParameters(myPlane->Pln().Position(), thePoint, aX, aY);
 427   aTargetPoint2d.SetCoord(aX, aY);
 428
 429   // choose the nearest circle
 430   Standard_Real aDistance, aP;
 431   GEOMImpl_Fillet1dPoint *aNearest;
 432   Standard_Integer a;
 433   TColStd_ListIteratorOfListOfReal anIter(myResultParams);
 434   for(aNearest = NULL, a = 1; anIter.More(); anIter.Next(), a++)
 435   {
 436     myStartSide = (myResultOrientation.Value(a)) ? Standard_True : Standard_False;
 437     GEOMImpl_Fillet1dPoint *aPoint = new GEOMImpl_Fillet1dPoint(anIter.Value());
 438     fillPoint(aPoint);
 439     if (!aPoint->HasSolution(myRadius))
 440       continue;
 441     aP = fabs(aPoint->GetCenter().Distance(aTargetPoint2d) - myRadius);
 442     if (!aNearest || aP < aDistance)
 443     {
 444       aNearest = aPoint;
 445       aDistance = aP;
 446     }
 447     else
 448     {
 449       delete aPoint;
 450     }
 451    }
 452
 453   if (!aNearest)
 454      return aResult;
 455
 456   // create circle edge
 457   gp_Pnt aCenter = ElSLib::PlaneValue(aNearest->GetCenter().X(),
 458                                       aNearest->GetCenter().Y(),
 459                                       myPlane->Pln().Position());
 460   Handle(Geom_Circle) aCircle =
 461     new Geom_Circle(gp_Ax2(aCenter, myPlane->Pln().Axis().Direction()), myRadius);
 462   gp_Pnt2d aPoint2d1, aPoint2d2;
 463   myCurve1->D0(aNearest->GetParam(), aPoint2d1);
 464   myCurve2->D0(aNearest->GetParam2(), aPoint2d2);
 465   gp_Pnt aPoint1 = ElSLib::PlaneValue(aPoint2d1.X(), aPoint2d1.Y(), myPlane->Pln().Position());
 466   gp_Pnt aPoint2 = ElSLib::PlaneValue(aPoint2d2.X(), aPoint2d2.Y(), myPlane->Pln().Position());
 467
 468   GeomAPI_ProjectPointOnCurve aProj(thePoint, aCircle);
 469   Standard_Real aTarGetParam = aProj.LowerDistanceParameter();
 470   gp_Pnt aPointOnCircle = aProj.NearestPoint();
 471
 472   // Check extrema point manually, because there is a bug in Open CASCADE
 473   //  in calculation of nearest point to a circle near the parameter 0.0
 474   gp_Pnt p0 = ElCLib::Value(0.0, aCircle->Circ());
 475   if (p0.Distance(thePoint) < aPointOnCircle.Distance(thePoint))
 476   {
 477      aTarGetParam = 0.0;
 478      aPointOnCircle = p0;
 479   }
 480
 481   aProj.Perform(aPoint1);
 482   Standard_Real aParam1 = aProj.LowerDistanceParameter();
 483     aProj.Perform(aPoint2);
 484   Standard_Real aParam2 = aProj.LowerDistanceParameter();
 485   Standard_Boolean aIsOut = ((aParam1 < aTarGetParam && aParam2 < aTarGetParam) ||
 486                              (aParam1 > aTarGetParam && aParam2 > aTarGetParam));
 487   if (aParam1 > aParam2)
 488     aIsOut = !aIsOut;
 489   BRepBuilderAPI_MakeEdge aBuilder(aCircle->Circ(),
 490                                    aIsOut ? aParam2 : aParam1,
 491                                    aIsOut? aParam1 : aParam2);
 492   aResult = aBuilder.Edge();
 493
 494   // divide edges
 495   Standard_Real aStart, anEnd;
 496   Handle(Geom_Curve) aCurve = BRep_Tool::Curve(myEdge1, aStart, anEnd);
 497   gp_Vec aDir;
 498   aCurve->D1(aNearest->GetParam(), aPoint1, aDir);
 499
 500   gp_Vec aCircleDir;
 501   aCircle->D1(aParam1, aPoint1, aCircleDir);
 502   if ((aCircleDir.Angle(aDir) > PI / 2.0) ^ aIsOut)
 503     aStart = aNearest->GetParam();
 504   else
 505     anEnd = aNearest->GetParam();
 506
 507   //Check the case when start and end are identical. This happens
 508   //when the edge decreases to size 0. Old ww5 allows such
 509   //cases. So we are again bug compatible
 510   if (fabs(aStart - anEnd) < Precision::PConfusion()/*gp::Resolution()*/)
 511     anEnd = 0.01;
 512   //Divide edge
 513   BRepBuilderAPI_MakeEdge aDivider1(aCurve, aStart, anEnd);
 514   if (myEdgesExchnged)
 515     theEdge2 = aDivider1.Edge();
 516   else
 517     theEdge1 = aDivider1.Edge();
 518
 519
 520   aCurve = BRep_Tool::Curve(myEdge2, aStart, anEnd);
 521   aCurve->D1(aNearest->GetParam2(), aPoint2, aDir);
 522
 523   aCircle->D1(aParam2, aPoint2, aCircleDir);
 524   if ((aCircleDir.Angle(aDir) > PI / 2.0) ^ (!aIsOut))
 525     aStart = aNearest->GetParam2();
 526   else
 527     anEnd = aNearest->GetParam2();
 528
 529   //Check the case when start and end are identical. This happens
 530   //when the edge decreases to size 0. Old ww5 allows such
 531   //cases. So we are again bug compatible
 532   if (fabs(aStart - anEnd) < Precision::PConfusion()/*gp::Resolution()*/)
 533     anEnd = 0.01;
 534   BRepBuilderAPI_MakeEdge aDivider2(aCurve, aStart, anEnd);
 535   if (myEdgesExchnged)
 536     theEdge1 = aDivider2.Edge();
 537   else
 538     theEdge2 = aDivider2.Edge();
 539
 540   delete aNearest;
 541   return aResult;
 542 }
 543
 544 //=======================================================================
 545 //function : AddValue
 546 //purpose  :
 547 //=======================================================================
 548 void GEOMImpl_Fillet1dPoint::AddValue(Standard_Real theValue, Standard_Boolean theValid)
 549 {
 550   Standard_Integer a;
 551   for(a = 1; a <= myV.Length(); a++)
 552   {
 553     if (theValue < myV.Value(a))
 554     {
 555       myV.InsertBefore(a, theValue);
 556       myValid.InsertBefore(a, (Standard_Integer)theValid);
 557       return;
 558     }
 559   }
 560   myV.Append(theValue);
 561   myValid.Append((Standard_Integer)theValid);
 562 }
 563
 564 //=======================================================================
 565 //function : ComputeDifference
 566 //purpose  :
 567 //=======================================================================
 568 Standard_Boolean GEOMImpl_Fillet1dPoint::ComputeDifference(GEOMImpl_Fillet1dPoint* thePoint)
 569 {
 570   Standard_Integer a;
 571   Standard_Boolean aDiffsSet = (myD.Length() != 0);
 572   Standard_Real aDX = thePoint->GetParam() - myParam, aDY;
 573   if (thePoint->myV.Length() == myV.Length())
 574   { // absolutely the same points
 575     for(a = 1; a <= myV.Length(); a++)
 576     {
 577       aDY = thePoint->myV.Value(a) - myV.Value(a);
 578       if ( aDiffsSet )
 579         myD.SetValue(a, fabs(aDX) > gp::Resolution() ? (aDY/aDX) : 0);
 580       else
 581         myD.Append( fabs(aDX) > gp::Resolution() ? (aDY/aDX) : 0);
 582     }
 583     return Standard_True;
 584   }
 585   // between the diffeerent points searching for nearest analogs
 586   Standard_Integer b;
 587   for(a = 1; a <= myV.Length(); a++)
 588   {
 589     for(b = 1; b <= thePoint->myV.Length(); b++)
 590     {
 591       if (b == 1 || fabs(thePoint->myV.Value(b) - myV.Value(a)) < fabs(aDY))
 592         aDY = thePoint->myV.Value(b) - myV.Value(a);
 593     }
 594     if (aDiffsSet)
 595     {
 596       if ( fabs(aDX) > gp::Resolution() && fabs(aDY / aDX) < fabs(myD.Value(a)))
 597         myD.SetValue(a, aDY / aDX);
 598       else
 599         myD.SetValue(a, 0);
 600     }
 601     else
 602     {
 603       myD.Append( fabs(aDX) > gp::Resolution() ? aDY/aDX : 0);
 604     }
 605   }
 606
 607   return Standard_False;
 608 }
 609
 610 //=======================================================================
 611 //function : FilterPoints
 612 //purpose  :
 613 //=======================================================================
 614 void GEOMImpl_Fillet1dPoint::FilterPoints(GEOMImpl_Fillet1dPoint* thePoint)
 615 {
 616   Standard_Integer a, b;
 617   TColStd_SequenceOfReal aDiffs;
 618   Standard_Real aY, aY2, aDX = thePoint->GetParam() - myParam;
 619   for(a = 1; a <= myV.Length(); a++)
 620   {
 621     // searching for near point from thePoint
 622     Standard_Integer aNear = 0;
 623     Standard_Real aDiff = aDX * 10000.;
 624     aY = myV.Value(a) + myD.Value(a) * aDX;
 625     for(b = 1; b <= thePoint->myV.Length(); b++)
 626     {
 627       // calculate hypothesis value of the Y2 with the constant first and second derivative
 628       aY2 = aY + aDX * (thePoint->myD.Value(b) - myD.Value(a)) / 2.0;
 629       if (aNear == 0 || fabs(aY2 - thePoint->myV.Value(b)) < fabs(aDiff))
 630       {
 631         aNear = b;
 632         aDiff = aY2 - thePoint->myV.Value(b);
 633       }
 634     }//for b...
 635
 636     if (aNear)
 637     {
 638       if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
 639       {// the same sign at the same sides of the interval
 640         if (myV.Value(a) * myD.Value(a) > 0)
 641         {
 642           if (fabs(myD.Value(a)) > Precision::Confusion())
 643             aNear = 0;
 644         }
 645         else
 646         {
 647           if (fabs(myV.Value(a)) > fabs(thePoint->myV.Value(aNear)))
 648             if (thePoint->myV.Value(aNear) * thePoint->myD.Value(aNear) < 0 &&
 649                 fabs(thePoint->myD.Value(aNear)) > Precision::Confusion())
 650             {
 651               aNear = 0;
 652             }
 653         }
 654       }
 655     }
 656
 657     if (aNear)
 658     {
 659       if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
 660       {
 661         if ((myV.Value(a) + myD.Value(a) * aDX) * myV.Value(a) > Precision::Confusion() &&
 662         (thePoint->myV.Value(aNear) + thePoint->myD.Value(aNear) * aDX) * thePoint->myV.Value(aNear) > Precision::Confusion())
 663         {
 664           aNear = 0;
 665         }
 666       }
 667     }
 668
 669     if (aNear)
 670     {
 671       if (  fabs(aDX) < gp::Resolution() || fabs(aDiff / aDX) > 1.e+7)
 672       {
 673         aNear = 0;
 674       }
 675     }
 676
 677     if (aNear == 0)
 678     {  // there is no near: remove it from the list
 679       myV.Remove(a);
 680       myD.Remove(a);
 681       myValid.Remove(a);
 682       a--;
 683     }
 684     else
 685     {
 686       Standard_Boolean aFound = Standard_False;
 687       for(b = 1; b <= myNear.Length(); b++)
 688       {
 689         if (myNear.Value(b) == aNear)
 690         {
 691           if (fabs(aDiffs.Value(b)) < fabs(aDiff))
 692           { // return this 'near'
 693             aFound = Standard_True;
 694             myV.Remove(a);
 695             myD.Remove(a);
 696             myValid.Remove(a);
 697             a--;
 698             break;
 699           }
 700           else
 701           { // remove the old 'near'
 702             myV.Remove(b);
 703             myD.Remove(b);
 704             myValid.Remove(b);
 705             myNear.Remove(b);
 706             aDiffs.Remove(b);
 707             a--;
 708             break;
 709           }
 710         }
 711       }//for b...
 712       if (!aFound)
 713       {
 714         myNear.Append(aNear);
 715         aDiffs.Append(aDiff);
 716       }
 717     }
 718   }//for a...
 719 }
 720
 721 //=======================================================================
 722 //function : Copy
 723 //purpose  :
 724 //=======================================================================
 725 GEOMImpl_Fillet1dPoint* GEOMImpl_Fillet1dPoint::Copy()
 726 {
 727   GEOMImpl_Fillet1dPoint* aCopy = new GEOMImpl_Fillet1dPoint(myParam);
 728   Standard_Integer a;
 729   for(a = 1; a <= myV.Length(); a++)
 730   {
 731     aCopy->myV.Append(myV.Value(a));
 732     aCopy->myD.Append(myD.Value(a));
 733     aCopy->myValid.Append(myValid.Value(a));
 734   }
 735   return aCopy;
 736 }
 737
 738 //=======================================================================
 739 //function : HasSolution
 740 //purpose  :
 741 //=======================================================================
 742 Standard_Integer GEOMImpl_Fillet1dPoint::HasSolution(const Standard_Real theRadius)
 743 {
 744   Standard_Integer a;
 745   for(a = 1; a <= myV.Length(); a++)
 746   {
 747     if (fabs(sqrt(fabs(fabs(myV.Value(a)) + theRadius * theRadius)) - theRadius) < Precision::Confusion() / 10.)
 748       return a;
 749   }
 750   return 0;
 751 }
 752
 753 //=======================================================================
 754 //function : RemoveSolution
 755 //purpose  :
 756 //=======================================================================
 757 void GEOMImpl_Fillet1dPoint::RemoveSolution(Standard_Integer theIndex)
 758 {
 759   myV.Remove(theIndex);
 760   myD.Remove(theIndex);
 761   myValid.Remove(theIndex);
 762   myNear.Remove(theIndex);
 763 }