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