src/GEOMImpl/GEOMImpl_Fillet1d.cxx

   1 // Copyright (C) 2005  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   2 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   3 //
   4 // This library is free software; you can redistribute it and/or
   5 // modify it under the terms of the GNU Lesser General Public
   6 // License as published by the Free Software Foundation; either
   7 // version 2.1 of the License.
   8 //
   9 // This library is distributed in the hope that it will be useful
  10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 // Lesser General Public License for more details.
  13 //
  14 // You should have received a copy of the GNU Lesser General Public
  15 // License along with this library; if not, write to the Free Software
  16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 //
  18 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  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)
 121       if (!theStartConnected)
 122         return Standard_True;
 123     if (aPoint.Distance(theEnd) < aTol)
 124       return Standard_True;
 125     if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
 126       return Standard_True;
 127   }
 128   Handle(Geom2d_Curve) aCurve;
 129   for(a = anInter.NbSegments(); a > 0; a--)
 130   {
 131     anInter.Segment(a, aCurve);
 132     aPoint = aCurve->Value(aCurve->FirstParameter());
 133     if (aPoint.Distance(theStart) < aTol)
 134       if (!theStartConnected)
 135         return Standard_True;
 136     if (aPoint.Distance(theEnd) < aTol)
 137       return Standard_True;
 138     if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
 139       return Standard_True;
 140     aPoint = aCurve->Value(aCurve->LastParameter());
 141     if (aPoint.Distance(theStart) < aTol)
 142       if (!theStartConnected)
 143         return Standard_True;
 144     if (aPoint.Distance(theEnd) < aTol)
 145       return Standard_True;
 146     if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
 147       return Standard_True;
 148   }
 149   return Standard_False;
 150 }
 151
 152
 153 //=======================================================================
 154 //function : fillPoint
 155 //purpose  :
 156 //=======================================================================
 157 void GEOMImpl_Fillet1d::fillPoint(GEOMImpl_Fillet1dPoint* thePoint)
 158 {
 159   gp_Pnt2d aPoint;
 160   gp_Vec2d aVec;
 161   const Standard_Real aTol = Precision::Confusion();
 162   myCurve1->D1(thePoint->GetParam(), aPoint, aVec);
 163   if (aVec.SquareMagnitude() < aTol)
 164     return;
 165
 166   gp_Vec2d aPerp(((myStartSide)?-1:1) * aVec.Y(), ((myStartSide)?1:-1) * aVec.X());
 167   aPerp.Normalize();
 168   aPerp.Multiply(myRadius);
 169   gp_Pnt2d aCenter = aPoint.Translated(aPerp);
 170   thePoint->SetCenter(aCenter);
 171
 172   // on the intersection point
 173   Standard_Boolean aValid = Standard_True;
 174   Geom2dAPI_ProjectPointOnCurve aProjInt(aPoint, myCurve2);
 175   if (aProjInt.NbPoints() && aPoint.Distance(aProjInt.NearestPoint()) < aTol)
 176     aValid = Standard_False;
 177   else
 178     aValid = !isRadiusIntersected(myCurve2, aPoint, aCenter, Standard_True);
 179
 180   Geom2dAPI_ProjectPointOnCurve aProj(aCenter, myCurve2);
 181   Standard_Integer a, aNB = aProj.NbPoints();
 182   for(a = aNB; a > 0; a--)
 183   {
 184     if (aPoint.Distance(aProj.Point(a)) < aTol)
 185       continue;
 186
 187     Standard_Boolean aValid2 = aValid;
 188     if (aValid2)
 189       aValid2 = !isRadiusIntersected(myCurve1, aCenter, aProj.Point(a), Standard_False);
 190
 191     // checking the right parameter
 192     Standard_Real aParam = aProj.Parameter(a);
 193     while(myCurve2->IsPeriodic() && aParam < myStart2)
 194       aParam += myCurve2->Period();
 195
 196     thePoint->AddValue(aProj.Distance(a) * aProj.Distance(a) - myRadius * myRadius,
 197                        (aParam >= myStart2 && aParam <= myEnd2 && aValid2));
 198     if (fabs(fabs(aProj.Distance(a)) - myRadius) < aTol)
 199       thePoint->SetParam2(aParam);
 200   }
 201 }
 202
 203 //=======================================================================
 204 //function : fillDiff
 205 //purpose  :
 206 //=======================================================================
 207 void GEOMImpl_Fillet1d::fillDiff(GEOMImpl_Fillet1dPoint* thePoint, Standard_Real theDiffStep, Standard_Boolean theFront)
 208 {
 209   GEOMImpl_Fillet1dPoint* aDiff =
 210     new GEOMImpl_Fillet1dPoint(thePoint->GetParam() + (theFront?(theDiffStep):(-theDiffStep)));
 211   fillPoint(aDiff);
 212   if (!thePoint->ComputeDifference(aDiff))
 213   {
 214     aDiff->SetParam(thePoint->GetParam() + (theFront?(-theDiffStep):(theDiffStep)));
 215     fillPoint(aDiff);
 216     thePoint->ComputeDifference(aDiff);
 217   }
 218   delete aDiff;
 219 }
 220
 221 //=======================================================================
 222 //function : Perform
 223 //purpose  :
 224 //=======================================================================
 225 Standard_Boolean GEOMImpl_Fillet1d::Perform(const Standard_Real theRadius)
 226 {
 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   Standard_Real aParam, aStep, aDStep;
 253   aStep = (myEnd1 - myStart1) / aNBSteps;
 254   aDStep = aStep/1000.;
 255
 256   Standard_Integer aCycle;
 257   for(aCycle = 2, myStartSide = Standard_False; aCycle; myStartSide = !myStartSide, aCycle--)
 258   {
 259     GEOMImpl_Fillet1dPoint *aLeft = NULL, *aRight = NULL;
 260
 261     for(aParam = myStart1 + aStep; aParam < myEnd1 || fabs(myEnd1 - aParam) < Precision::Confusion(); aParam += aStep)
 262     {
 263       if (!aLeft)
 264       {
 265         aLeft = new GEOMImpl_Fillet1dPoint(aParam - aStep);
 266         fillPoint(aLeft);
 267         fillDiff(aLeft, aDStep, Standard_True);
 268       }
 269
 270       aRight = new GEOMImpl_Fillet1dPoint(aParam);
 271       fillPoint(aRight);
 272       fillDiff(aRight, aDStep, Standard_False);
 273
 274       aLeft->FilterPoints(aRight);
 275       performNewton(aLeft, aRight);
 276
 277       delete aLeft;
 278       aLeft = aRight;
 279     }
 280     delete aLeft;
 281   }
 282
 283   if (myResultParams.Extent())
 284     return Standard_True;
 285
 286   return Standard_False;
 287 }
 288
 289 //=======================================================================
 290 //function : processPoint
 291 //purpose  :
 292 //=======================================================================
 293 Standard_Boolean GEOMImpl_Fillet1d::processPoint(GEOMImpl_Fillet1dPoint* theLeft,
 294                                                  GEOMImpl_Fillet1dPoint* theRight,
 295                                                  Standard_Real           theParameter)
 296 {
 297   if (theParameter >= theLeft->GetParam() && theParameter < theRight->GetParam())
 298   {
 299     Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
 300     if (theParameter - theLeft->GetParam() < aDX / 100.)
 301     {
 302       theParameter = theLeft->GetParam() + aDX / 100.;
 303     }
 304     if (theRight->GetParam() - theParameter < aDX / 100.)
 305     {
 306       theParameter = theRight->GetParam() - aDX / 100.;
 307     }
 308
 309     GEOMImpl_Fillet1dPoint* aPoint1 = theLeft->Copy();
 310     GEOMImpl_Fillet1dPoint* aPoint2 = new GEOMImpl_Fillet1dPoint(theParameter);
 311     fillPoint(aPoint2);
 312     fillDiff(aPoint2, aDX / 1000., Standard_True);
 313
 314     aPoint1->FilterPoints(aPoint2);
 315     performNewton(aPoint1, aPoint2);
 316     aPoint2->FilterPoints(theRight);
 317     performNewton(aPoint2, theRight);
 318
 319     delete aPoint1;
 320     delete aPoint2;
 321     return Standard_True;
 322   }
 323
 324   return Standard_False;
 325 }
 326
 327 //=======================================================================
 328 //function : performNewton
 329 //purpose  :
 330 //=======================================================================
 331 void GEOMImpl_Fillet1d::performNewton(GEOMImpl_Fillet1dPoint* theLeft,
 332                                       GEOMImpl_Fillet1dPoint* theRight)
 333 {
 334   Standard_Integer a;
 335   // check the left: if this is solution store it and remove it from the list of researching points of theLeft
 336   a = theLeft->HasSolution(myRadius);
 337   if (a)
 338   {
 339     if (theLeft->IsValid(a))
 340     {
 341       myResultParams.Append(theLeft->GetParam());
 342       myResultOrientation.Append(myStartSide);
 343     }
 344     return;
 345   }
 346
 347   Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
 348   if (aDX < gp::Resolution()/*Precision::Confusion() / 1000000.*/)
 349   {
 350     a = theRight->HasSolution(myRadius);
 351     if (a)
 352       if (theRight->IsValid(a))
 353       {
 354         myResultParams.Append(theRight->GetParam());
 355         myResultOrientation.Append(myStartSide);
 356       }
 357     return;
 358   }
 359
 360   for(a = 1; a <= theLeft->GetNBValues(); a++)
 361   {
 362     Standard_Integer aNear = theLeft->GetNear(a);
 363
 364       Standard_Real aA = (theRight->GetDiff(aNear) - theLeft->GetDiff(a)) / aDX;
 365     Standard_Real aB = theLeft->GetDiff(a) - aA * theLeft->GetParam();
 366     Standard_Real aC = theLeft->GetValue(a) - theLeft->GetDiff(a) * theLeft->GetParam() +
 367                        aA * theLeft->GetParam() * theLeft->GetParam() / 2.0;
 368     Standard_Real aDet = aB * aB - 2.0 * aA * aC;
 369
 370     if (fabs(aA) < Precision::Confusion())
 371     { // linear case
 372       if (fabs(aB) > 10e-20)
 373       {
 374         Standard_Real aX0 = - aC / aB; // use extremum
 375         if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
 376           processPoint(theLeft, theRight, aX0);
 377       }
 378       else
 379       {
 380         processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 381       }
 382     }
 383     else
 384     {
 385       if (fabs(aB) > fabs(aDet * 1000000.))
 386       {  // possible floating point operations accurancy errors
 387         processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 388       }
 389       else
 390       {
 391         if (aDet > 0)
 392         { // two solutions
 393           aDet = sqrt(aDet);
 394           Standard_Boolean aRes = processPoint(theLeft, theRight, (- aB + aDet) / aA);
 395           if (!aRes)
 396             aRes = processPoint(theLeft, theRight, (- aB - aDet) / aA);
 397           if (!aRes)
 398             processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 399         }
 400         else
 401         {
 402           Standard_Real aX0 = - aB / aA; // use extremum
 403           if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
 404             processPoint(theLeft, theRight, aX0);
 405           else
 406             processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 407         }
 408       }
 409     }
 410   }
 411 }
 412
 413 //=======================================================================
 414 //function : Result
 415 //purpose  :
 416 //=======================================================================
 417 TopoDS_Edge GEOMImpl_Fillet1d::Result(const gp_Pnt& thePoint,
 418                                       TopoDS_Edge& theEdge1,
 419                                       TopoDS_Edge& theEdge2)
 420 {
 421   TopoDS_Edge aResult;
 422   gp_Pnt2d aTargetPoint2d;
 423   Standard_Real aX, aY;
 424   ElSLib::PlaneParameters(myPlane->Pln().Position(), thePoint, aX, aY);
 425   aTargetPoint2d.SetCoord(aX, aY);
 426
 427   // choose the nearest circle
 428   Standard_Real aDistance, aP;
 429   GEOMImpl_Fillet1dPoint *aNearest;
 430   Standard_Integer a;
 431   TColStd_ListIteratorOfListOfReal anIter(myResultParams);
 432   for(aNearest = NULL, a = 1; anIter.More(); anIter.Next(), a++)
 433   {
 434     myStartSide = (myResultOrientation.Value(a)) ? Standard_True : Standard_False;
 435     GEOMImpl_Fillet1dPoint *aPoint = new GEOMImpl_Fillet1dPoint(anIter.Value());
 436     fillPoint(aPoint);
 437     if (!aPoint->HasSolution(myRadius))
 438       continue;
 439     aP = fabs(aPoint->GetCenter().Distance(aTargetPoint2d) - myRadius);
 440     if (!aNearest || aP < aDistance)
 441     {
 442       aNearest = aPoint;
 443       aDistance = aP;
 444     }
 445     else
 446     {
 447       delete aPoint;
 448     }
 449    }
 450
 451   if (!aNearest)
 452      return aResult;
 453
 454   // create circle edge
 455   gp_Pnt aCenter = ElSLib::PlaneValue(aNearest->GetCenter().X(),
 456                                       aNearest->GetCenter().Y(),
 457                                       myPlane->Pln().Position());
 458   Handle(Geom_Circle) aCircle =
 459     new Geom_Circle(gp_Ax2(aCenter, myPlane->Pln().Axis().Direction()), myRadius);
 460   gp_Pnt2d aPoint2d1, aPoint2d2;
 461   myCurve1->D0(aNearest->GetParam(), aPoint2d1);
 462   myCurve2->D0(aNearest->GetParam2(), aPoint2d2);
 463   gp_Pnt aPoint1 = ElSLib::PlaneValue(aPoint2d1.X(), aPoint2d1.Y(), myPlane->Pln().Position());
 464   gp_Pnt aPoint2 = ElSLib::PlaneValue(aPoint2d2.X(), aPoint2d2.Y(), myPlane->Pln().Position());
 465
 466   GeomAPI_ProjectPointOnCurve aProj(thePoint, aCircle);
 467   Standard_Real aTarGetParam = aProj.LowerDistanceParameter();
 468   gp_Pnt aPointOnCircle = aProj.NearestPoint();
 469
 470   // Check extrema point manually, because there is a bug in Open CASCADE
 471   //  in calculation of nearest point to a circle near the parameter 0.0
 472   gp_Pnt p0 = ElCLib::Value(0.0, aCircle->Circ());
 473   if (p0.Distance(thePoint) < aPointOnCircle.Distance(thePoint))
 474   {
 475      aTarGetParam = 0.0;
 476      aPointOnCircle = p0;
 477   }
 478
 479   aProj.Perform(aPoint1);
 480   Standard_Real aParam1 = aProj.LowerDistanceParameter();
 481     aProj.Perform(aPoint2);
 482   Standard_Real aParam2 = aProj.LowerDistanceParameter();
 483   Standard_Boolean aIsOut = ((aParam1 < aTarGetParam && aParam2 < aTarGetParam) ||
 484                              (aParam1 > aTarGetParam && aParam2 > aTarGetParam));
 485   if (aParam1 > aParam2)
 486     aIsOut = !aIsOut;
 487   BRepBuilderAPI_MakeEdge aBuilder(aCircle->Circ(),
 488                                    aIsOut ? aParam2 : aParam1,
 489                                    aIsOut? aParam1 : aParam2);
 490   aResult = aBuilder.Edge();
 491
 492   // divide edges
 493   Standard_Real aStart, anEnd;
 494   Handle(Geom_Curve) aCurve = BRep_Tool::Curve(myEdge1, aStart, anEnd);
 495   gp_Vec aDir;
 496   aCurve->D1(aNearest->GetParam(), aPoint1, aDir);
 497
 498   gp_Vec aCircleDir;
 499   aCircle->D1(aParam1, aPoint1, aCircleDir);
 500   if ((aCircleDir.Angle(aDir) > PI / 2.0) ^ aIsOut)
 501     aStart = aNearest->GetParam();
 502   else
 503     anEnd = aNearest->GetParam();
 504
 505   //Check the case when start and end are identical. This happens
 506   //when the edge decreases to size 0. Old ww5 allows such
 507   //cases. So we are again bug compatible
 508   if (fabs(aStart - anEnd) < Precision::PConfusion()/*gp::Resolution()*/)
 509     anEnd = 0.01;
 510   //Divide edge
 511   BRepBuilderAPI_MakeEdge aDivider1(aCurve, aStart, anEnd);
 512   if (myEdgesExchnged)
 513     theEdge2 = aDivider1.Edge();
 514   else
 515     theEdge1 = aDivider1.Edge();
 516
 517
 518   aCurve = BRep_Tool::Curve(myEdge2, aStart, anEnd);
 519   aCurve->D1(aNearest->GetParam2(), aPoint2, aDir);
 520
 521   aCircle->D1(aParam2, aPoint2, aCircleDir);
 522   if ((aCircleDir.Angle(aDir) > PI / 2.0) ^ (!aIsOut))
 523     aStart = aNearest->GetParam2();
 524   else
 525     anEnd = aNearest->GetParam2();
 526
 527   //Check the case when start and end are identical. This happens
 528   //when the edge decreases to size 0. Old ww5 allows such
 529   //cases. So we are again bug compatible
 530   if (fabs(aStart - anEnd) < Precision::PConfusion()/*gp::Resolution()*/)
 531     anEnd = 0.01;
 532   BRepBuilderAPI_MakeEdge aDivider2(aCurve, aStart, anEnd);
 533   if (myEdgesExchnged)
 534     theEdge1 = aDivider2.Edge();
 535   else
 536     theEdge2 = aDivider2.Edge();
 537
 538   delete aNearest;
 539   return aResult;
 540 }
 541
 542 //=======================================================================
 543 //function : AddValue
 544 //purpose  :
 545 //=======================================================================
 546 void GEOMImpl_Fillet1dPoint::AddValue(Standard_Real theValue, Standard_Boolean theValid)
 547 {
 548   Standard_Integer a;
 549   for(a = 1; a <= myV.Length(); a++)
 550   {
 551     if (theValue < myV.Value(a))
 552     {
 553       myV.InsertBefore(a, theValue);
 554       myValid.InsertBefore(a, (Standard_Integer)theValid);
 555       return;
 556     }
 557   }
 558   myV.Append(theValue);
 559   myValid.Append((Standard_Integer)theValid);
 560 }
 561
 562 //=======================================================================
 563 //function : ComputeDifference
 564 //purpose  :
 565 //=======================================================================
 566 Standard_Boolean GEOMImpl_Fillet1dPoint::ComputeDifference(GEOMImpl_Fillet1dPoint* thePoint)
 567 {
 568   Standard_Integer a;
 569   Standard_Boolean aDiffsSet = (myD.Length() != 0);
 570   Standard_Real aDX = thePoint->GetParam() - myParam, aDY;
 571   if (thePoint->myV.Length() == myV.Length())
 572   { // absolutely the same points
 573     for(a = 1; a <= myV.Length(); a++)
 574     {
 575       aDY = thePoint->myV.Value(a) - myV.Value(a);
 576       if (aDiffsSet)
 577         myD.SetValue(a, aDY / aDX);
 578       else
 579         myD.Append(aDY / aDX);
 580     }
 581     return Standard_True;
 582   }
 583   // between the diffeerent points searching for nearest analogs
 584   Standard_Integer b;
 585   for(a = 1; a <= myV.Length(); a++)
 586   {
 587     for(b = 1; b <= thePoint->myV.Length(); b++)
 588     {
 589       if (b == 1 || fabs(thePoint->myV.Value(b) - myV.Value(a)) < fabs(aDY))
 590         aDY = thePoint->myV.Value(b) - myV.Value(a);
 591     }
 592     if (aDiffsSet)
 593     {
 594       if (fabs(aDY / aDX) < fabs(myD.Value(a)))
 595         myD.SetValue(a, aDY / aDX);
 596     }
 597     else
 598     {
 599       myD.Append(aDY / aDX);
 600     }
 601   }
 602
 603   return Standard_False;
 604 }
 605
 606 //=======================================================================
 607 //function : FilterPoints
 608 //purpose  :
 609 //=======================================================================
 610 void GEOMImpl_Fillet1dPoint::FilterPoints(GEOMImpl_Fillet1dPoint* thePoint)
 611 {
 612   Standard_Integer a, b;
 613   TColStd_SequenceOfReal aDiffs;
 614   Standard_Real aY, aY2, aDX = thePoint->GetParam() - myParam;
 615   for(a = 1; a <= myV.Length(); a++)
 616   {
 617     // searching for near point from thePoint
 618     Standard_Integer aNear = 0;
 619     Standard_Real aDiff = aDX * 10000.;
 620     aY = myV.Value(a) + myD.Value(a) * aDX;
 621     for(b = 1; b <= thePoint->myV.Length(); b++)
 622     {
 623       // calculate hypothesis value of the Y2 with the constant first and second derivative
 624       aY2 = aY + aDX * (thePoint->myD.Value(b) - myD.Value(a)) / 2.0;
 625       if (aNear == 0 || fabs(aY2 - thePoint->myV.Value(b)) < fabs(aDiff))
 626       {
 627         aNear = b;
 628         aDiff = aY2 - thePoint->myV.Value(b);
 629       }
 630     }//for b...
 631
 632     if (aNear)
 633     {
 634       if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
 635       {// the same sign at the same sides of the interval
 636         if (myV.Value(a) * myD.Value(a) > 0)
 637         {
 638           if (fabs(myD.Value(a)) > Precision::Confusion())
 639             aNear = 0;
 640         }
 641         else
 642         {
 643           if (fabs(myV.Value(a)) > fabs(thePoint->myV.Value(aNear)))
 644             if (thePoint->myV.Value(aNear) * thePoint->myD.Value(aNear) < 0 &&
 645                 fabs(thePoint->myD.Value(aNear)) > Precision::Confusion())
 646             {
 647               aNear = 0;
 648             }
 649         }
 650       }
 651     }
 652
 653     if (aNear)
 654     {
 655       if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
 656       {
 657         if ((myV.Value(a) + myD.Value(a) * aDX) * myV.Value(a) > Precision::Confusion() &&
 658         (thePoint->myV.Value(aNear) + thePoint->myD.Value(aNear) * aDX) * thePoint->myV.Value(aNear) > Precision::Confusion())
 659         {
 660           aNear = 0;
 661         }
 662       }
 663     }
 664
 665     if (aNear)
 666     {
 667       if (fabs(aDiff / aDX) > 1.e+7)
 668       {
 669         aNear = 0;
 670       }
 671     }
 672
 673     if (aNear == 0)
 674     {  // there is no near: remove it from the list
 675       myV.Remove(a);
 676       myD.Remove(a);
 677       myValid.Remove(a);
 678       a--;
 679     }
 680     else
 681     {
 682       Standard_Boolean aFound = Standard_False;
 683       for(b = 1; b <= myNear.Length(); b++)
 684       {
 685         if (myNear.Value(b) == aNear)
 686         {
 687           if (fabs(aDiffs.Value(b)) < fabs(aDiff))
 688           { // return this 'near'
 689             aFound = Standard_True;
 690             myV.Remove(a);
 691             myD.Remove(a);
 692             myValid.Remove(a);
 693             a--;
 694             break;
 695           }
 696           else
 697           { // remove the old 'near'
 698             myV.Remove(b);
 699             myD.Remove(b);
 700             myValid.Remove(b);
 701             myNear.Remove(b);
 702             aDiffs.Remove(b);
 703             a--;
 704             break;
 705           }
 706         }
 707       }//for b...
 708       if (!aFound)
 709       {
 710         myNear.Append(aNear);
 711         aDiffs.Append(aDiff);
 712       }
 713     }
 714   }//for a...
 715 }
 716
 717 //=======================================================================
 718 //function : Copy
 719 //purpose  :
 720 //=======================================================================
 721 GEOMImpl_Fillet1dPoint* GEOMImpl_Fillet1dPoint::Copy()
 722 {
 723   GEOMImpl_Fillet1dPoint* aCopy = new GEOMImpl_Fillet1dPoint(myParam);
 724   Standard_Integer a;
 725   for(a = 1; a <= myV.Length(); a++)
 726   {
 727     aCopy->myV.Append(myV.Value(a));
 728     aCopy->myD.Append(myD.Value(a));
 729     aCopy->myValid.Append(myValid.Value(a));
 730   }
 731   return aCopy;
 732 }
 733
 734 //=======================================================================
 735 //function : HasSolution
 736 //purpose  :
 737 //=======================================================================
 738 Standard_Integer GEOMImpl_Fillet1dPoint::HasSolution(const Standard_Real theRadius)
 739 {
 740   Standard_Integer a;
 741   for(a = 1; a <= myV.Length(); a++)
 742   {
 743     if (fabs(sqrt(fabs(fabs(myV.Value(a)) + theRadius * theRadius)) - theRadius) < Precision::Confusion() / 10.)
 744       return a;
 745   }
 746   return 0;
 747 }
 748
 749 //=======================================================================
 750 //function : RemoveSolution
 751 //purpose  :
 752 //=======================================================================
 753 void GEOMImpl_Fillet1dPoint::RemoveSolution(Standard_Integer theIndex)
 754 {
 755   myV.Remove(theIndex);
 756   myD.Remove(theIndex);
 757   myValid.Remove(theIndex);
 758   myNear.Remove(theIndex);
 759 }