1 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 // This file is part of Open CASCADE Technology software library.
5 // This library is free software; you can redistribute it and/or modify it under
6 // the terms of the GNU Lesser General Public License version 2.1 as published
7 // by the Free Software Foundation, with special exception defined in the file
8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
9 // distribution for complete text of the license and disclaimer of any warranty.
11 // Alternatively, this file may be used under the terms of Open CASCADE
12 // commercial license or contractual agreement.
14 #include <ChFi2d_FilletAlgo.hxx>
16 #include <GeomProjLib.hxx>
17 #include <BRep_Tool.hxx>
18 #include <Precision.hxx>
22 #include <Geom2dAPI_ProjectPointOnCurve.hxx>
23 #include <GeomAPI_ProjectPointOnCurve.hxx>
24 #include <Geom2dAPI_InterCurveCurve.hxx>
27 #include <TopoDS_Iterator.hxx>
28 #include <TColStd_ListIteratorOfListOfReal.hxx>
30 #include <Geom_Circle.hxx>
31 #include <Geom2d_Line.hxx>
33 #include <BRepBuilderAPI_MakeEdge.hxx>
34 #include <BRepAdaptor_Curve.hxx>
36 ChFi2d_FilletAlgo::ChFi2d_FilletAlgo()
42 myStartSide (Standard_False),
43 myEdgesExchnged(Standard_False),
44 myDegreeOfRecursion(0)
48 ChFi2d_FilletAlgo::ChFi2d_FilletAlgo(const TopoDS_Wire& theWire, const gp_Pln& thePlane)
54 myStartSide (Standard_False),
55 myEdgesExchnged(Standard_False),
56 myDegreeOfRecursion(0)
58 Init(theWire, thePlane);
61 ChFi2d_FilletAlgo::ChFi2d_FilletAlgo(const TopoDS_Edge& theEdge1,
62 const TopoDS_Edge& theEdge2,
63 const gp_Pln& thePlane)
71 myStartSide (Standard_False),
72 myEdgesExchnged(Standard_False),
73 myDegreeOfRecursion(0)
75 Init(theEdge1, theEdge2, thePlane);
78 void ChFi2d_FilletAlgo::Init(const TopoDS_Wire& theWire, const gp_Pln& thePlane)
80 TopoDS_Edge theEdge1, theEdge2;
81 TopoDS_Iterator itr(theWire);
82 for (; itr.More(); itr.Next())
84 if (theEdge1.IsNull())
85 theEdge1 = TopoDS::Edge(itr.Value());
86 else if (theEdge2.IsNull())
87 theEdge2 = TopoDS::Edge(itr.Value());
91 if (theEdge1.IsNull() || theEdge2.IsNull())
92 throw Standard_ConstructionError("The fillet algorithms expects a wire consisting of two edges.");
93 Init(theEdge1, theEdge2, thePlane);
96 void ChFi2d_FilletAlgo::Init(const TopoDS_Edge& theEdge1,
97 const TopoDS_Edge& theEdge2,
98 const gp_Pln& thePlane)
100 myPlane = new Geom_Plane(thePlane);
102 myEdgesExchnged = Standard_False;
104 BRepAdaptor_Curve aBAC1(theEdge1);
105 BRepAdaptor_Curve aBAC2(theEdge2);
106 if (aBAC1.GetType() < aBAC2.GetType())
107 { // first curve must be more complicated
110 myEdgesExchnged = Standard_True;
118 Handle(Geom_Curve) aCurve1 = BRep_Tool::Curve(myEdge1, myStart1, myEnd1);
119 Handle(Geom_Curve) aCurve2 = BRep_Tool::Curve(myEdge2, myStart2, myEnd2);
121 myCurve1 = GeomProjLib::Curve2d(aCurve1, myStart1, myEnd1, myPlane);
122 myCurve2 = GeomProjLib::Curve2d(aCurve2, myStart2, myEnd2, myPlane);
124 while (myCurve1->IsPeriodic() && myStart1 >= myEnd1)
125 myEnd1 += myCurve1->Period();
126 while (myCurve2->IsPeriodic() && myStart2 >= myEnd2)
127 myEnd2 += myCurve2->Period();
129 if (aBAC1.GetType() == aBAC2.GetType())
131 if (myEnd2 - myStart2 < myEnd1 - myStart1)
132 { // first curve must be parametrically shorter
133 TopoDS_Edge anEdge = myEdge1;
136 Handle(Geom2d_Curve) aCurve = myCurve1;
139 Standard_Real a = myStart1;
145 myEdgesExchnged = Standard_True;
150 //! This function returns true if linear segment from start point of the
151 //! fillet arc to the end point is intersected by the first or second
152 //! curve: in this case fillet is invalid.
153 static Standard_Boolean IsRadiusIntersected(const Handle(Geom2d_Curve)& theCurve, const Standard_Real theCurveMin, const double theCurveMax,
154 const gp_Pnt2d theStart, const gp_Pnt2d theEnd, const Standard_Boolean theStartConnected)
156 //Check the given start and end if they are identical. If yes
158 if (theStart.SquareDistance(theEnd) < Precision::SquareConfusion())
160 return Standard_False;
162 Handle(Geom2d_Line) line = new Geom2d_Line(theStart, gp_Dir2d(gp_Vec2d(theStart, theEnd)));
163 Geom2dAPI_InterCurveCurve anInter(theCurve, line, Precision::Confusion());
166 for(a = anInter.NbPoints(); a > 0; a--)
168 aPoint = anInter.Point(a);
169 Geom2dAPI_ProjectPointOnCurve aProjInt(aPoint, theCurve, theCurveMin, theCurveMax);
170 if (aProjInt.NbPoints() < 1 || aProjInt.LowerDistanceParameter() > Precision::Confusion())
171 continue; // point is not on edge
173 if (aPoint.Distance(theStart) < Precision::Confusion())
175 if (!theStartConnected)
176 return Standard_True;
178 if (aPoint.Distance(theEnd) < Precision::Confusion())
179 return Standard_True;
180 if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), Precision::Angular()))
181 return Standard_True;
183 const Handle(Geom2d_Curve)& aCurve = theCurve;
184 for(a = anInter.NbSegments(); a > 0; a--)
186 //anInter.Segment(a, aCurve); //not implemented (bug in OCC)
187 aPoint = aCurve->Value(aCurve->FirstParameter());
189 Geom2dAPI_ProjectPointOnCurve aProjInt(aPoint, theCurve, theCurveMin, theCurveMax);
190 if (aProjInt.NbPoints() && aProjInt.LowerDistanceParameter() < Precision::Confusion())
191 { // point is on edge
192 if (aPoint.Distance(theStart) < Precision::Confusion())
193 if (!theStartConnected)
194 return Standard_True;
195 if (aPoint.Distance(theEnd) < Precision::Confusion())
196 return Standard_True;
197 if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), Precision::Angular()))
198 return Standard_True;
200 aPoint = aCurve->Value(aCurve->LastParameter());
202 aProjInt.Init(aPoint, theCurve, theCurveMin, theCurveMax);
203 if (aProjInt.NbPoints() && aProjInt.LowerDistanceParameter() < Precision::Confusion())
204 { // point is on edge
205 if (aPoint.Distance(theStart) < Precision::Confusion())
206 if (!theStartConnected)
207 return Standard_True;
208 if (aPoint.Distance(theEnd) < Precision::Confusion())
209 return Standard_True;
210 if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), Precision::Angular()))
211 return Standard_True;
214 return Standard_False;
217 void ChFi2d_FilletAlgo::FillPoint(FilletPoint* thePoint, const Standard_Real theLimit)
220 // on the intersection point
221 Standard_Boolean aValid = Standard_False;
222 Standard_Real aStep = Precision::Confusion();
223 gp_Pnt2d aCenter, aPoint; // center of fillet and point on curve1
224 Standard_Real aParam = thePoint->getParam();
225 if (theLimit < aParam) aStep = -aStep;
226 for(aValid = Standard_False; !aValid; aParam += aStep)
228 if ((aParam - aStep - theLimit) * (aParam - theLimit) <= 0)
229 break; // limit was exceeded
232 myCurve1->D1(aParam, aPoint, aVec);
233 if (aVec.SquareMagnitude() < Precision::Confusion())
236 gp_Vec2d aPerp(((myStartSide)?-1:1) * aVec.Y(), ((myStartSide)?1:-1) * aVec.X());
238 aPerp.Multiply(myRadius);
239 aCenter = aPoint.Translated(aPerp);
242 Geom2dAPI_ProjectPointOnCurve aProjInt(aPoint, myCurve2, myStart2, myEnd2);
243 if (aProjInt.NbPoints() == 0 || aPoint.Distance(aProjInt.NearestPoint()) > Precision::Confusion())
245 aValid = Standard_True;
251 thePoint->setParam(aParam);
252 thePoint->setCenter(aCenter);
253 aValid = !IsRadiusIntersected(myCurve2, myStart2, myEnd2, aPoint, aCenter, Standard_True);
256 Geom2dAPI_ProjectPointOnCurve aProj(aCenter, myCurve2);
257 int a, aNB = aProj.NbPoints();
258 for(a = aNB; a > 0; a--)
260 if (aPoint.SquareDistance(aProj.Point(a)) < Precision::Confusion())
263 Standard_Boolean aValid2 = aValid;
265 aValid2 = !IsRadiusIntersected(myCurve1, myStart1, myEnd1, aCenter, aProj.Point(a), Standard_False);
267 // checking the right parameter
268 Standard_Real aParamProj = aProj.Parameter(a);
269 while(myCurve2->IsPeriodic() && aParamProj < myStart2)
270 aParamProj += myCurve2->Period();
272 const Standard_Real d = aProj.Distance(a);
273 thePoint->appendValue(d * d - myRadius * myRadius, (aParamProj >= myStart2 && aParamProj <= myEnd2 && aValid2));
274 if (Abs(d - myRadius) < Precision::Confusion())
275 thePoint->setParam2(aParamProj);
279 void ChFi2d_FilletAlgo::FillDiff(FilletPoint* thePoint, Standard_Real theDiffStep, Standard_Boolean theFront)
281 Standard_Real aDelta = theFront?(theDiffStep):(-theDiffStep);
282 FilletPoint* aDiff = new FilletPoint(thePoint->getParam() + aDelta);
283 FillPoint(aDiff, aDelta * 999.);
284 if (!thePoint->calculateDiff(aDiff))
286 aDiff->setParam(thePoint->getParam() - aDelta);
287 FillPoint(aDiff, - aDelta * 999);
288 thePoint->calculateDiff(aDiff);
293 // returns true, if at least one result was found
294 Standard_Boolean ChFi2d_FilletAlgo::Perform(const Standard_Real theRadius)
296 myDegreeOfRecursion = 0;
297 myResultParams.Clear();
298 myResultOrientation.Clear();
300 Standard_Real aNBSteps;
301 Geom2dAdaptor_Curve aGAC(myCurve1);
302 switch (aGAC.GetType())
310 case GeomAbs_Ellipse:
313 case GeomAbs_BSplineCurve:
314 aNBSteps = 2 + aGAC.Degree() * aGAC.NbPoles();
316 default: // unknown: maximum
319 //std::cout<<"aNBSteps = "<<aNBSteps<<std::endl;
321 myRadius = theRadius;
322 Standard_Real aParam, aStep, aDStep;
323 aStep = (myEnd1 - myStart1) / aNBSteps;
324 aDStep = 1.e-4 * aStep;
327 for(aCycle = 2, myStartSide = Standard_False; aCycle; myStartSide = !myStartSide, aCycle--)
329 FilletPoint *aLeft = NULL, *aRight;
331 for(aParam = myStart1 + aStep; aParam < myEnd1 || Abs(myEnd1 - aParam) < Precision::Confusion(); aParam += aStep)
335 aLeft = new FilletPoint(aParam - aStep);
336 FillPoint(aLeft, aParam);
337 FillDiff(aLeft, aDStep, Standard_True);
340 aRight = new FilletPoint(aParam);
341 FillPoint(aRight, aParam - aStep);
342 FillDiff(aRight, aDStep, Standard_False);
344 aLeft->FilterPoints(aRight);
345 PerformNewton(aLeft, aRight);
353 return !myResultParams.IsEmpty();
356 Standard_Boolean ChFi2d_FilletAlgo::ProcessPoint(FilletPoint* theLeft, FilletPoint* theRight, Standard_Real theParameter)
358 if (theParameter >= theLeft->getParam() && theParameter < theRight->getParam())
360 Standard_Real aDX = (theRight->getParam() - theLeft->getParam());
361 if (theParameter - theLeft->getParam() < aDX/100.0)
363 theParameter = theLeft->getParam() + aDX/100.0;
365 if (theRight->getParam() - theParameter < aDX/100.0)
367 theParameter = theRight->getParam() - aDX/100.0;
370 // Protection on infinite loops.
371 myDegreeOfRecursion++;
372 Standard_Real diffx = 0.001 * aDX;
373 if (myDegreeOfRecursion > 100)
376 if (myDegreeOfRecursion > 1000)
379 if (myDegreeOfRecursion > 3000)
381 return Standard_True;
386 FilletPoint* aPoint1 = theLeft->Copy();
387 FilletPoint* aPoint2 = new FilletPoint(theParameter);
388 FillPoint(aPoint2, aPoint1->getParam());
389 FillDiff(aPoint2, diffx, Standard_True);
391 aPoint1->FilterPoints(aPoint2);
392 PerformNewton(aPoint1, aPoint2);
393 aPoint2->FilterPoints(theRight);
394 PerformNewton(aPoint2, theRight);
398 return Standard_True;
401 return Standard_False;
404 void ChFi2d_FilletAlgo::PerformNewton(FilletPoint* theLeft, FilletPoint* theRight)
407 // check the left: if this is solution store it and remove it from the list of researching points of theLeft
408 a = theLeft->hasSolution(myRadius);
411 if (theLeft->isValid(a))
413 myResultParams.Append(theLeft->getParam());
414 myResultOrientation.Append(myStartSide);
419 Standard_Real aDX = theRight->getParam() - theLeft->getParam();
420 if (aDX < 1.e-6 * Precision::Confusion())
422 a = theRight->hasSolution(myRadius);
423 if (a && theRight->isValid(a))
425 myResultParams.Append(theRight->getParam());
426 myResultOrientation.Append(myStartSide);
430 for(a = 1; a <= theLeft->getNBValues(); a++)
432 int aNear = theLeft->getNear(a);
434 Standard_Real aA = (theRight->getDiff(aNear) - theLeft->getDiff(a)) / aDX;
435 Standard_Real aB = theLeft->getDiff(a) - aA * theLeft->getParam();
436 Standard_Real aC = theLeft->getValue(a) - theLeft->getDiff(a) * theLeft->getParam() + aA * theLeft->getParam() * theLeft->getParam() / 2.0;
437 Standard_Real aDet = aB * aB - 2.0 * aA * aC;
439 if (Abs(aA) < Precision::Confusion())
441 //std::cout<<"###"<<std::endl;
442 if (Abs(aB) > 10e-20)
444 Standard_Real aX0 = - aC / aB; // use extremum
445 if (aX0 > theLeft->getParam() && aX0 < theRight->getParam())
446 ProcessPoint(theLeft, theRight, aX0);
450 ProcessPoint(theLeft, theRight, theLeft->getParam() + aDX / 2.0); // linear division otherwise
455 if (Abs(aB) > Abs(aDet * 1000000.))
456 { // possible floating point operations accurancy errors
458 ProcessPoint(theLeft, theRight, theLeft->getParam() + aDX / 2.0); // linear division otherwise
465 Standard_Boolean aRes = ProcessPoint(theLeft, theRight, (- aB + aDet) / aA);
467 aRes = ProcessPoint(theLeft, theRight, (- aB - aDet) / aA);
469 ProcessPoint(theLeft, theRight, theLeft->getParam() + aDX / 2.0); // linear division otherwise
473 //std::cout<<"%%%"<<std::endl;
474 Standard_Real aX0 = - aB / aA; // use extremum
475 if (aX0 > theLeft->getParam() && aX0 < theRight->getParam())
476 ProcessPoint(theLeft, theRight, aX0);
478 ProcessPoint(theLeft, theRight, theLeft->getParam() + aDX / 2.0); // linear division otherwise
485 // returns number of possible solutions.
486 int ChFi2d_FilletAlgo::NbResults(const gp_Pnt& thePoint)
488 Standard_Real aX, aY;
489 gp_Pnt2d aTargetPoint2d;
490 ElSLib::PlaneParameters(myPlane->Pln().Position(), thePoint, aX, aY);
491 aTargetPoint2d.SetCoord(aX, aY);
493 //iterate through all possible solutions.
495 TColStd_ListIteratorOfListOfReal anIter(myResultParams);
496 for(; anIter.More(); anIter.Next(), i++)
498 myStartSide = (myResultOrientation.Value(i)) ? Standard_True : Standard_False;
499 FilletPoint *aPoint = new FilletPoint(anIter.Value());
500 FillPoint(aPoint, anIter.Value() + 1.);
501 if (aPoint->hasSolution(myRadius))
509 // returns result (fillet edge, modified edge1, modified edge2), nearest to the given point <thePoint>
510 TopoDS_Edge ChFi2d_FilletAlgo::Result(const gp_Pnt& thePoint, TopoDS_Edge& theEdge1, TopoDS_Edge& theEdge2, const int iSolution)
513 gp_Pnt2d aTargetPoint2d;
514 Standard_Real aX, aY;
515 ElSLib::PlaneParameters(myPlane->Pln().Position(), thePoint, aX, aY);
516 aTargetPoint2d.SetCoord(aX, aY);
518 // choose the nearest circle
519 Standard_Real aDistance = 0.0, aP;
520 FilletPoint *aNearest;
522 TColStd_ListIteratorOfListOfReal anIter(myResultParams);
523 for(aNearest = NULL, a = 1; anIter.More(); anIter.Next(), a++)
525 myStartSide = (myResultOrientation.Value(a))?Standard_True:Standard_False;
526 FilletPoint *aPoint = new FilletPoint(anIter.Value());
527 FillPoint(aPoint, anIter.Value() + 1.);
528 if (!aPoint->hasSolution(myRadius))
536 aP = Abs(aPoint->getCenter().Distance(aTargetPoint2d) - myRadius);
538 else if (iSolution == iSol)
542 if (!aNearest || aP < aDistance)
551 if (iSolution == iSol)
559 // create circle edge
560 gp_Pnt aCenter = ElSLib::PlaneValue(aNearest->getCenter().X(), aNearest->getCenter().Y(), myPlane->Pln().Position());
561 Handle(Geom_Circle) aCircle = new Geom_Circle(gp_Ax2(aCenter, myPlane->Pln().Axis().Direction()), myRadius);
562 gp_Pnt2d aPoint2d1, aPoint2d2;
563 myCurve1->D0(aNearest->getParam(), aPoint2d1);
564 myCurve2->D0(aNearest->getParam2(), aPoint2d2);
565 gp_Pnt aPoint1 = ElSLib::PlaneValue(aPoint2d1.X(), aPoint2d1.Y(), myPlane->Pln().Position());
566 gp_Pnt aPoint2 = ElSLib::PlaneValue(aPoint2d2.X(), aPoint2d2.Y(), myPlane->Pln().Position());
568 GeomAPI_ProjectPointOnCurve aProj(thePoint, aCircle);
569 Standard_Real aTargetParam = aProj.LowerDistanceParameter();
570 gp_Pnt aPointOnCircle = aProj.NearestPoint();
572 // There is a bug in Open CASCADE in calculation of nearest point to a circle near the parameter 0.0
573 // Therefore I check this extrema point manually:
574 gp_Pnt p0 = ElCLib::Value(0.0, aCircle->Circ());
575 if (p0.Distance(thePoint) < aPointOnCircle.Distance(thePoint))
581 aProj.Perform(aPoint1);
582 Standard_Real aParam1 = aProj.LowerDistanceParameter();
583 aProj.Perform(aPoint2);
584 Standard_Real aParam2 = aProj.LowerDistanceParameter();
585 Standard_Boolean aIsOut = ((aParam1 < aTargetParam && aParam2 < aTargetParam) || (aParam1 > aTargetParam && aParam2 > aTargetParam));
586 if (aParam1 > aParam2)
588 BRepBuilderAPI_MakeEdge aBuilder(aCircle->Circ(), aIsOut ? aParam2 : aParam1, aIsOut? aParam1 : aParam2);
589 aResult = aBuilder.Edge();
592 Standard_Real aStart, anEnd;
593 Handle(Geom_Curve) aCurve = BRep_Tool::Curve(myEdge1, aStart, anEnd);
595 aCurve->D1(aNearest->getParam(), aPoint1, aDir);
598 aCircle->D1(aParam1, aPoint1, aCircleDir);
600 if ((aCircleDir.Angle(aDir) > M_PI / 2.0) ? !aIsOut : aIsOut)
601 aStart = aNearest->getParam();
603 anEnd = aNearest->getParam();
605 //Check the case when start and end are identical. This happens
606 //when the edge decreases to size 0. Old ww5 allows such
607 //cases. So we are again bug compatible
608 if (fabs(aStart - anEnd) < Precision::Confusion())
609 anEnd = aStart + Precision::Confusion();
611 BRepBuilderAPI_MakeEdge aDivider1(aCurve, aStart, anEnd);
613 theEdge2 = aDivider1.Edge();
615 theEdge1 = aDivider1.Edge();
617 aCurve = BRep_Tool::Curve(myEdge2, aStart, anEnd);
618 aCurve->D1(aNearest->getParam2(), aPoint2, aDir);
620 aCircle->D1(aParam2, aPoint2, aCircleDir);
622 if ((aCircleDir.Angle(aDir) > M_PI / 2.0) ? aIsOut : !aIsOut)
623 aStart = aNearest->getParam2();
625 anEnd = aNearest->getParam2();
627 //Check the case when start and end are identical. This happens
628 //when the edge decreases to size 0. Old ww5 allows such
629 //cases. So we are again bug compatible
630 if (fabs(aStart - anEnd) < Precision::Confusion())
631 anEnd = aStart + Precision::Confusion();
632 BRepBuilderAPI_MakeEdge aDivider2(aCurve, aStart, anEnd);
634 theEdge1 = aDivider2.Edge();
636 theEdge2 = aDivider2.Edge();
642 FilletPoint::FilletPoint(const Standard_Real theParam)
643 : myParam (theParam),
648 void FilletPoint::appendValue(Standard_Real theValue, Standard_Boolean theValid)
651 for(a = 1; a <= myV.Length(); a++)
653 if (theValue < myV.Value(a))
655 myV.InsertBefore(a, theValue);
656 myValid.InsertBefore(a, theValid);
660 myV.Append(theValue);
661 myValid.Append(theValid);
664 Standard_Boolean FilletPoint::calculateDiff(FilletPoint* thePoint)
667 Standard_Boolean aDiffsSet = (myD.Length() != 0);
668 Standard_Real aDX = thePoint->getParam() - myParam, aDY = 0.0;
669 if (thePoint->myV.Length() == myV.Length())
670 { // absolutely the same points
671 for(a = 1; a <= myV.Length(); a++)
673 aDY = thePoint->myV.Value(a) - myV.Value(a);
675 myD.SetValue(a, aDY / aDX);
677 myD.Append(aDY / aDX);
679 return Standard_True;
681 // between the diffeerent points searching for nearest analogs
683 for(a = 1; a <= myV.Length(); a++)
685 for(b = 1; b <= thePoint->myV.Length(); b++)
687 if (b == 1 || Abs(thePoint->myV.Value(b) - myV.Value(a)) < Abs(aDY))
688 aDY = thePoint->myV.Value(b) - myV.Value(a);
692 if (Abs(aDY / aDX) < Abs(myD.Value(a)))
693 myD.SetValue(a, aDY / aDX);
697 myD.Append(aDY / aDX);
701 return Standard_False;
704 void FilletPoint::FilterPoints(FilletPoint* thePoint)
706 Standard_Integer a, b;
707 TColStd_SequenceOfReal aDiffs;
708 Standard_Real aY, aY2, aDX = thePoint->getParam() - myParam;
709 for(a = 1; a <= myV.Length(); a++)
711 // searching for near point from thePoint
712 Standard_Integer aNear = 0;
713 Standard_Real aDiff = aDX * 10000.;
714 aY = myV.Value(a) + myD.Value(a) * aDX;
715 for(b = 1; b <= thePoint->myV.Length(); b++)
717 // calculate hypothesis value of the Y2 with the constant first and second derivative
718 aY2 = aY + aDX * (thePoint->myD.Value(b) - myD.Value(a)) / 2.0;
719 if (aNear == 0 || Abs(aY2 - thePoint->myV.Value(b)) < Abs(aDiff))
722 aDiff = aY2 - thePoint->myV.Value(b);
728 if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
729 {// the same sign at the same sides of the interval
730 if (myV.Value(a) * myD.Value(a) > 0)
732 if (Abs(myD.Value(a)) > Precision::Confusion())
737 if (Abs(myV.Value(a)) > Abs(thePoint->myV.Value(aNear)))
738 if (thePoint->myV.Value(aNear) * thePoint->myD.Value(aNear) < 0 && Abs(thePoint->myD.Value(aNear)) > Precision::Confusion())
748 if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
750 if ((myV.Value(a) + myD.Value(a) * aDX) * myV.Value(a) > Precision::Confusion() &&
751 (thePoint->myV.Value(aNear) + thePoint->myD.Value(aNear) * aDX) * thePoint->myV.Value(aNear) > Precision::Confusion())
760 if (Abs(aDiff / aDX) > 1.e+7)
767 { // there is no near: remove it from the list
775 Standard_Boolean aFound = Standard_False;
776 for(b = 1; b <= myNear.Length(); b++)
778 if (myNear.Value(b) == aNear)
780 if (Abs(aDiffs.Value(b)) < Abs(aDiff))
781 { // return this 'near'
782 aFound = Standard_True;
790 { // remove the old 'near'
803 myNear.Append(aNear);
804 aDiffs.Append(aDiff);
810 FilletPoint* FilletPoint::Copy()
812 FilletPoint* aCopy = new FilletPoint(myParam);
814 for(a = 1; a <= myV.Length(); a++)
816 aCopy->myV.Append(myV.Value(a));
817 aCopy->myD.Append(myD.Value(a));
818 aCopy->myValid.Append(myValid.Value(a));
823 int FilletPoint::hasSolution(const Standard_Real theRadius)
826 for(a = 1; a <= myV.Length(); a++)
828 if (Abs(sqrt(Abs(Abs(myV.Value(a)) + theRadius * theRadius)) - theRadius) < Precision::Confusion())
834 void FilletPoint::remove(int theIndex)
836 myV.Remove(theIndex);
837 myD.Remove(theIndex);
838 myValid.Remove(theIndex);
839 myNear.Remove(theIndex);