1 // Copyright (C) 2007-2011 CEA/DEN, EDF R&D, OPEN CASCADE
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 // File : GEOMImpl_Fillet1d.cxx
\r
21 // Module : GEOMImpl
\r
23 #include "GEOMImpl_Fillet1d.hxx"
\r
25 #include <BRep_Tool.hxx>
\r
26 #include <BRepAdaptor_Curve.hxx>
\r
27 #include <BRepBuilderAPI_MakeEdge.hxx>
\r
28 #include <ElCLib.hxx>
\r
29 #include <ElSLib.hxx>
\r
30 #include <gp_Circ.hxx>
\r
31 #include <Geom2d_Line.hxx>
\r
32 #include <Geom2dAPI_ProjectPointOnCurve.hxx>
\r
33 #include <Geom2dAPI_InterCurveCurve.hxx>
\r
34 #include <GeomAPI_ProjectPointOnCurve.hxx>
\r
35 #include <GeomProjLib.hxx>
\r
36 #include <Geom_Circle.hxx>
\r
37 #include <Precision.hxx>
\r
38 #include <TColStd_ListIteratorOfListOfReal.hxx>
\r
41 * class GEOMImpl_Fillet1d
\r
45 //=======================================================================
\r
46 //function : Constructor
\r
48 //=======================================================================
\r
49 GEOMImpl_Fillet1d::GEOMImpl_Fillet1d(const TopoDS_Edge& theEdge1,
\r
50 const TopoDS_Edge& theEdge2,
\r
51 const gp_Pln& thePlane)
\r
52 : myEdgesExchnged( Standard_False )
\r
54 myPlane = new Geom_Plane(thePlane);
\r
56 BRepAdaptor_Curve aBAC1(theEdge1);
\r
57 BRepAdaptor_Curve aBAC2(theEdge2);
\r
58 if (aBAC1.GetType() < aBAC2.GetType())
\r
59 { // first curve must be more complicated
\r
62 myEdgesExchnged = Standard_True;
\r
70 Handle(Geom_Curve) aCurve1 = BRep_Tool::Curve(myEdge1, myStart1, myEnd1);
\r
71 Handle(Geom_Curve) aCurve2 = BRep_Tool::Curve(myEdge2, myStart2, myEnd2);
\r
73 myCurve1 = GeomProjLib::Curve2d(aCurve1, myStart1, myEnd1, myPlane);
\r
74 myCurve2 = GeomProjLib::Curve2d(aCurve2, myStart2, myEnd2, myPlane);
\r
76 while (myCurve1->IsPeriodic() && myStart1 >= myEnd1)
\r
77 myEnd1 += myCurve1->Period();
\r
78 while (myCurve2->IsPeriodic() && myStart2 >= myEnd2)
\r
79 myEnd2 += myCurve2->Period();
\r
81 if (aBAC1.GetType() == aBAC2.GetType())
\r
83 if (myEnd2 - myStart2 < myEnd1 - myStart1)
\r
84 { // first curve must be parametrically shorter
\r
85 TopoDS_Edge anEdge = myEdge1;
\r
88 Handle(Geom2d_Curve) aCurve = myCurve1;
\r
89 myCurve1 = myCurve2;
\r
91 Standard_Real a = myStart1;
\r
92 myStart1 = myStart2;
\r
97 myEdgesExchnged = Standard_True;
\r
102 //=======================================================================
\r
103 //function : isRadiusIntersected
\r
104 //purpose : local function
\r
105 //=======================================================================
\r
106 static Standard_Boolean isRadiusIntersected(const Handle(Geom2d_Curve)& theCurve,
\r
107 const gp_Pnt2d theStart,
\r
108 const gp_Pnt2d theEnd,
\r
109 const Standard_Boolean theStartConnected)
\r
111 const Standard_Real aTol = Precision::Confusion();
\r
112 const Standard_Real anAngTol = Precision::Angular();
\r
113 Geom2dAPI_InterCurveCurve anInter(theCurve, new Geom2d_Line(theStart,
\r
114 gp_Dir2d(gp_Vec2d(theStart, theEnd))), aTol);
\r
115 Standard_Integer a;
\r
117 for(a = anInter.NbPoints(); a > 0; a--)
\r
119 aPoint = anInter.Point(a);
\r
120 if ( aPoint.Distance(theStart) < aTol && !theStartConnected )
\r
121 return Standard_True;
\r
122 if (aPoint.Distance(theEnd) < aTol * 200)
\r
123 return Standard_True;
\r
124 if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
\r
125 return Standard_True;
\r
127 Handle(Geom2d_Curve) aCurve;
\r
128 for(a = anInter.NbSegments(); a > 0; a--)
\r
130 anInter.Segment(a, aCurve);
\r
131 aPoint = aCurve->Value(aCurve->FirstParameter());
\r
132 if (aPoint.Distance(theStart) < aTol)
\r
133 if (!theStartConnected)
\r
134 return Standard_True;
\r
135 if (aPoint.Distance(theEnd) < aTol)
\r
136 return Standard_True;
\r
137 if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
\r
138 return Standard_True;
\r
139 aPoint = aCurve->Value(aCurve->LastParameter());
\r
140 if (aPoint.Distance(theStart) < aTol)
\r
141 if (!theStartConnected)
\r
142 return Standard_True;
\r
143 if (aPoint.Distance(theEnd) < aTol)
\r
144 return Standard_True;
\r
145 if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
\r
146 return Standard_True;
\r
148 return Standard_False;
\r
152 //=======================================================================
\r
153 //function : fillPoint
\r
155 //=======================================================================
\r
156 void GEOMImpl_Fillet1d::fillPoint(GEOMImpl_Fillet1dPoint* thePoint)
\r
160 const Standard_Real aTol = Precision::Confusion();
\r
161 myCurve1->D1(thePoint->GetParam(), aPoint, aVec);
\r
162 if (aVec.SquareMagnitude() < aTol)
\r
165 gp_Vec2d aPerp(((myStartSide)?-1:1) * aVec.Y(), ((myStartSide)?1:-1) * aVec.X());
\r
167 aPerp.Multiply(myRadius);
\r
168 gp_Pnt2d aCenter = aPoint.Translated(aPerp);
\r
169 thePoint->SetCenter(aCenter);
\r
171 // on the intersection point
\r
172 Standard_Boolean aValid = Standard_True;
\r
173 Geom2dAPI_ProjectPointOnCurve aProjInt(aPoint, myCurve2);
\r
174 if (aProjInt.NbPoints() && aPoint.Distance(aProjInt.NearestPoint()) < aTol)
\r
175 aValid = Standard_False;
\r
177 aValid = !isRadiusIntersected(myCurve2, aPoint, aCenter, Standard_True);
\r
179 Geom2dAPI_ProjectPointOnCurve aProj(aCenter, myCurve2);
\r
180 Standard_Integer a, aNB = aProj.NbPoints();
\r
181 for(a = aNB; a > 0; a--)
\r
183 if (aPoint.Distance(aProj.Point(a)) < aTol)
\r
186 Standard_Boolean aValid2 = aValid;
\r
188 aValid2 = !isRadiusIntersected(myCurve1, aCenter, aProj.Point(a), Standard_False);
\r
190 // checking the right parameter
\r
191 Standard_Real aParam = aProj.Parameter(a);
\r
192 while(myCurve2->IsPeriodic() && aParam < myStart2)
\r
193 aParam += myCurve2->Period();
\r
195 thePoint->AddValue(aProj.Distance(a) * aProj.Distance(a) - myRadius * myRadius,
\r
196 (aParam >= myStart2 && aParam <= myEnd2 && aValid2));
\r
197 if (fabs(fabs(aProj.Distance(a)) - myRadius) < aTol)
\r
198 thePoint->SetParam2(aParam);
\r
202 //=======================================================================
\r
203 //function : fillDiff
\r
205 //=======================================================================
\r
206 void GEOMImpl_Fillet1d::fillDiff(GEOMImpl_Fillet1dPoint* thePoint, Standard_Real theDiffStep, Standard_Boolean theFront)
\r
208 GEOMImpl_Fillet1dPoint* aDiff =
\r
209 new GEOMImpl_Fillet1dPoint(thePoint->GetParam() + (theFront?(theDiffStep):(-theDiffStep)));
\r
211 if (!thePoint->ComputeDifference(aDiff))
\r
213 aDiff->SetParam(thePoint->GetParam() + (theFront?(-theDiffStep):(theDiffStep)));
\r
215 thePoint->ComputeDifference(aDiff);
\r
220 //=======================================================================
\r
221 //function : Perform
\r
223 //=======================================================================
\r
224 Standard_Boolean GEOMImpl_Fillet1d::Perform(const Standard_Real theRadius)
\r
226 myDegreeOfRecursion = 0;
\r
227 myResultParams.Clear();
\r
228 myResultOrientation.Clear();
\r
230 Standard_Real aNBSteps = 100;
\r
231 Geom2dAdaptor_Curve aGAC(myCurve1);
\r
232 switch (aGAC.GetType())
\r
237 case GeomAbs_Circle:
\r
240 case GeomAbs_Ellipse:
\r
243 case GeomAbs_BezierCurve:
\r
244 case GeomAbs_BSplineCurve:
\r
245 aNBSteps = 2 + aGAC.Degree() * aGAC.NbPoles();
\r
247 default: // unknown: maximum
\r
251 myRadius = theRadius;
\r
252 Standard_Real aParam, aStep, aDStep;
\r
253 aStep = (myEnd1 - myStart1) / aNBSteps;
\r
254 aDStep = aStep/1000.;
\r
256 Standard_Integer aCycle;
\r
257 for(aCycle = 2, myStartSide = Standard_False; aCycle; myStartSide = !myStartSide, aCycle--)
\r
259 GEOMImpl_Fillet1dPoint *aLeft = NULL, *aRight = NULL;
\r
261 for(aParam = myStart1 + aStep; aParam < myEnd1 || fabs(myEnd1 - aParam) < Precision::Confusion(); aParam += aStep)
\r
265 aLeft = new GEOMImpl_Fillet1dPoint(aParam - aStep);
\r
267 fillDiff(aLeft, aDStep, Standard_True);
\r
270 aRight = new GEOMImpl_Fillet1dPoint(aParam);
\r
272 fillDiff(aRight, aDStep, Standard_False);
\r
274 aLeft->FilterPoints(aRight);
\r
275 performNewton(aLeft, aRight);
\r
283 if (myResultParams.Extent())
\r
284 return Standard_True;
\r
286 return Standard_False;
\r
289 //=======================================================================
\r
290 //function : processPoint
\r
292 //=======================================================================
\r
293 Standard_Boolean GEOMImpl_Fillet1d::processPoint(GEOMImpl_Fillet1dPoint* theLeft,
\r
294 GEOMImpl_Fillet1dPoint* theRight,
\r
295 Standard_Real theParameter)
\r
297 if (theParameter >= theLeft->GetParam() && theParameter < theRight->GetParam())
\r
299 Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
\r
300 if (theParameter - theLeft->GetParam() < aDX / 100.)
\r
302 theParameter = theLeft->GetParam() + aDX / 100.;
\r
304 if (theRight->GetParam() - theParameter < aDX / 100.)
\r
306 theParameter = theRight->GetParam() - aDX / 100.;
\r
309 // Protection on infinite loop.
\r
310 myDegreeOfRecursion++;
\r
311 Standard_Real diffx = 0.001 * aDX;
\r
312 if (myDegreeOfRecursion > 1000)
\r
315 if (myDegreeOfRecursion > 10000)
\r
318 if (myDegreeOfRecursion > 100000)
\r
320 return Standard_True;
\r
325 GEOMImpl_Fillet1dPoint* aPoint1 = theLeft->Copy();
\r
326 GEOMImpl_Fillet1dPoint* aPoint2 = new GEOMImpl_Fillet1dPoint(theParameter);
\r
327 fillPoint(aPoint2);
\r
328 fillDiff(aPoint2, diffx, Standard_True);
\r
330 aPoint1->FilterPoints(aPoint2);
\r
331 performNewton(aPoint1, aPoint2);
\r
332 aPoint2->FilterPoints(theRight);
\r
333 performNewton(aPoint2, theRight);
\r
337 return Standard_True;
\r
340 return Standard_False;
\r
343 //=======================================================================
\r
344 //function : performNewton
\r
346 //=======================================================================
\r
347 void GEOMImpl_Fillet1d::performNewton(GEOMImpl_Fillet1dPoint* theLeft,
\r
348 GEOMImpl_Fillet1dPoint* theRight)
\r
350 Standard_Integer a;
\r
351 // check the left: if this is solution store it and remove it from the list of researching points of theLeft
\r
352 a = theLeft->HasSolution(myRadius);
\r
355 if (theLeft->IsValid(a))
\r
357 myResultParams.Append(theLeft->GetParam());
\r
358 myResultOrientation.Append(myStartSide);
\r
363 Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
\r
364 if ( aDX < Precision::Confusion() / 1000000.)
\r
366 a = theRight->HasSolution(myRadius);
\r
368 if (theRight->IsValid(a))
\r
370 myResultParams.Append(theRight->GetParam());
\r
371 myResultOrientation.Append(myStartSide);
\r
376 for(a = 1; a <= theLeft->GetNBValues(); a++)
\r
378 Standard_Integer aNear = theLeft->GetNear(a);
\r
380 Standard_Real aA = (theRight->GetDiff(aNear) - theLeft->GetDiff(a)) / aDX;
\r
381 Standard_Real aB = theLeft->GetDiff(a) - aA * theLeft->GetParam();
\r
382 Standard_Real aC = theLeft->GetValue(a) - theLeft->GetDiff(a) * theLeft->GetParam() +
\r
383 aA * theLeft->GetParam() * theLeft->GetParam() / 2.0;
\r
384 Standard_Real aDet = aB * aB - 2.0 * aA * aC;
\r
386 if ( fabs(aDet) < gp::Resolution() )
\r
389 if (fabs(aA) < Precision::Confusion())
\r
391 if (fabs(aB) > 10e-20)
\r
393 Standard_Real aX0 = - aC / aB; // use extremum
\r
394 if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
\r
395 processPoint(theLeft, theRight, aX0);
\r
399 processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
\r
404 if (fabs(aB) > fabs(aDet * 1000000.))
\r
405 { // possible floating point operations accurancy errors
\r
406 processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
\r
413 Standard_Boolean aRes = processPoint(theLeft, theRight, (- aB + aDet) / aA);
\r
415 aRes = processPoint(theLeft, theRight, (- aB - aDet) / aA);
\r
417 processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
\r
421 Standard_Real aX0 = - aB / aA; // use extremum
\r
422 if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
\r
423 processPoint(theLeft, theRight, aX0);
\r
425 processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
\r
432 //=======================================================================
\r
433 //function : Result
\r
435 //=======================================================================
\r
436 TopoDS_Edge GEOMImpl_Fillet1d::Result(const gp_Pnt& thePoint,
\r
437 TopoDS_Edge& theEdge1,
\r
438 TopoDS_Edge& theEdge2)
\r
440 TopoDS_Edge aResult;
\r
441 gp_Pnt2d aTargetPoint2d;
\r
442 Standard_Real aX, aY;
\r
443 ElSLib::PlaneParameters(myPlane->Pln().Position(), thePoint, aX, aY);
\r
444 aTargetPoint2d.SetCoord(aX, aY);
\r
446 // choose the nearest circle
\r
447 Standard_Real aDistance, aP;
\r
448 GEOMImpl_Fillet1dPoint *aNearest;
\r
449 Standard_Integer a;
\r
450 TColStd_ListIteratorOfListOfReal anIter(myResultParams);
\r
451 for(aNearest = NULL, a = 1; anIter.More(); anIter.Next(), a++)
\r
453 myStartSide = (myResultOrientation.Value(a)) ? Standard_True : Standard_False;
\r
454 GEOMImpl_Fillet1dPoint *aPoint = new GEOMImpl_Fillet1dPoint(anIter.Value());
\r
456 if (!aPoint->HasSolution(myRadius))
\r
458 aP = fabs(aPoint->GetCenter().Distance(aTargetPoint2d) - myRadius);
\r
459 if (!aNearest || aP < aDistance)
\r
473 // create circle edge
\r
474 gp_Pnt aCenter = ElSLib::PlaneValue(aNearest->GetCenter().X(),
\r
475 aNearest->GetCenter().Y(),
\r
476 myPlane->Pln().Position());
\r
477 Handle(Geom_Circle) aCircle =
\r
478 new Geom_Circle(gp_Ax2(aCenter, myPlane->Pln().Axis().Direction()), myRadius);
\r
479 gp_Pnt2d aPoint2d1, aPoint2d2;
\r
480 myCurve1->D0(aNearest->GetParam(), aPoint2d1);
\r
481 myCurve2->D0(aNearest->GetParam2(), aPoint2d2);
\r
482 gp_Pnt aPoint1 = ElSLib::PlaneValue(aPoint2d1.X(), aPoint2d1.Y(), myPlane->Pln().Position());
\r
483 gp_Pnt aPoint2 = ElSLib::PlaneValue(aPoint2d2.X(), aPoint2d2.Y(), myPlane->Pln().Position());
\r
485 GeomAPI_ProjectPointOnCurve aProj(thePoint, aCircle);
\r
486 Standard_Real aTarGetParam = aProj.LowerDistanceParameter();
\r
487 gp_Pnt aPointOnCircle = aProj.NearestPoint();
\r
489 // Check extrema point manually, because there is a bug in Open CASCADE
\r
490 // in calculation of nearest point to a circle near the parameter 0.0
\r
491 gp_Pnt p0 = ElCLib::Value(0.0, aCircle->Circ());
\r
492 if (p0.Distance(thePoint) < aPointOnCircle.Distance(thePoint))
\r
494 aTarGetParam = 0.0;
\r
495 aPointOnCircle = p0;
\r
498 aProj.Perform(aPoint1);
\r
499 Standard_Real aParam1 = aProj.LowerDistanceParameter();
\r
500 aProj.Perform(aPoint2);
\r
501 Standard_Real aParam2 = aProj.LowerDistanceParameter();
\r
502 Standard_Boolean aIsOut = ((aParam1 < aTarGetParam && aParam2 < aTarGetParam) ||
\r
503 (aParam1 > aTarGetParam && aParam2 > aTarGetParam));
\r
504 if (aParam1 > aParam2)
\r
506 BRepBuilderAPI_MakeEdge aBuilder(aCircle->Circ(),
\r
507 aIsOut ? aParam2 : aParam1,
\r
508 aIsOut? aParam1 : aParam2);
\r
509 aResult = aBuilder.Edge();
\r
512 Standard_Real aStart, anEnd;
\r
513 Handle(Geom_Curve) aCurve = BRep_Tool::Curve(myEdge1, aStart, anEnd);
\r
515 aCurve->D1(aNearest->GetParam(), aPoint1, aDir);
\r
518 aCircle->D1(aParam1, aPoint1, aCircleDir);
\r
519 if ((aCircleDir.Angle(aDir) > PI / 2.0) ^ aIsOut)
\r
520 aStart = aNearest->GetParam();
\r
522 anEnd = aNearest->GetParam();
\r
524 if (fabs(aStart - anEnd) > Precision::Confusion())
\r
527 BRepBuilderAPI_MakeEdge aDivider1(aCurve, aStart, anEnd);
\r
528 if (myEdgesExchnged)
\r
529 theEdge2 = aDivider1.Edge();
\r
531 theEdge1 = aDivider1.Edge();
\r
534 aCurve = BRep_Tool::Curve(myEdge2, aStart, anEnd);
\r
535 aCurve->D1(aNearest->GetParam2(), aPoint2, aDir);
\r
537 aCircle->D1(aParam2, aPoint2, aCircleDir);
\r
538 if ((aCircleDir.Angle(aDir) > PI / 2.0) ^ (!aIsOut))
\r
539 aStart = aNearest->GetParam2();
\r
541 anEnd = aNearest->GetParam2();
\r
543 if (fabs(aStart - anEnd) > Precision::Confusion())
\r
545 BRepBuilderAPI_MakeEdge aDivider2(aCurve, aStart, anEnd);
\r
546 if (myEdgesExchnged)
\r
547 theEdge1 = aDivider2.Edge();
\r
549 theEdge2 = aDivider2.Edge();
\r
556 //=======================================================================
\r
557 //function : AddValue
\r
559 //=======================================================================
\r
560 void GEOMImpl_Fillet1dPoint::AddValue(Standard_Real theValue, Standard_Boolean theValid)
\r
562 Standard_Integer a;
\r
563 for(a = 1; a <= myV.Length(); a++)
\r
565 if (theValue < myV.Value(a))
\r
567 myV.InsertBefore(a, theValue);
\r
568 myValid.InsertBefore(a, (Standard_Integer)theValid);
\r
572 myV.Append(theValue);
\r
573 myValid.Append((Standard_Integer)theValid);
\r
576 //=======================================================================
\r
577 //function : ComputeDifference
\r
579 //=======================================================================
\r
580 Standard_Boolean GEOMImpl_Fillet1dPoint::ComputeDifference(GEOMImpl_Fillet1dPoint* thePoint)
\r
582 Standard_Integer a;
\r
583 Standard_Boolean aDiffsSet = (myD.Length() != 0);
\r
584 Standard_Real aDX = thePoint->GetParam() - myParam, aDY;
\r
585 if (thePoint->myV.Length() == myV.Length())
\r
586 { // absolutely the same points
\r
587 for(a = 1; a <= myV.Length(); a++)
\r
589 aDY = thePoint->myV.Value(a) - myV.Value(a);
\r
591 myD.SetValue(a, fabs(aDX) > gp::Resolution() ? (aDY/aDX) : 0);
\r
593 myD.Append( fabs(aDX) > gp::Resolution() ? (aDY/aDX) : 0);
\r
595 return Standard_True;
\r
597 // between the diffeerent points searching for nearest analogs
\r
598 Standard_Integer b;
\r
599 for(a = 1; a <= myV.Length(); a++)
\r
601 for(b = 1; b <= thePoint->myV.Length(); b++)
\r
603 if (b == 1 || fabs(thePoint->myV.Value(b) - myV.Value(a)) < fabs(aDY))
\r
604 aDY = thePoint->myV.Value(b) - myV.Value(a);
\r
608 if ( fabs(aDX) > gp::Resolution() && fabs(aDY / aDX) < fabs(myD.Value(a)))
\r
609 myD.SetValue(a, aDY / aDX);
\r
611 myD.SetValue(a, 0);
\r
615 myD.Append( fabs(aDX) > gp::Resolution() ? aDY/aDX : 0);
\r
619 return Standard_False;
\r
622 //=======================================================================
\r
623 //function : FilterPoints
\r
625 //=======================================================================
\r
626 void GEOMImpl_Fillet1dPoint::FilterPoints(GEOMImpl_Fillet1dPoint* thePoint)
\r
628 Standard_Integer a, b;
\r
629 TColStd_SequenceOfReal aDiffs;
\r
630 Standard_Real aY, aY2, aDX = thePoint->GetParam() - myParam;
\r
631 for(a = 1; a <= myV.Length(); a++)
\r
633 // searching for near point from thePoint
\r
634 Standard_Integer aNear = 0;
\r
635 Standard_Real aDiff = aDX * 10000.;
\r
636 aY = myV.Value(a) + myD.Value(a) * aDX;
\r
637 for(b = 1; b <= thePoint->myV.Length(); b++)
\r
639 // calculate hypothesis value of the Y2 with the constant first and second derivative
\r
640 aY2 = aY + aDX * (thePoint->myD.Value(b) - myD.Value(a)) / 2.0;
\r
641 if (aNear == 0 || fabs(aY2 - thePoint->myV.Value(b)) < fabs(aDiff))
\r
644 aDiff = aY2 - thePoint->myV.Value(b);
\r
650 if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
\r
651 {// the same sign at the same sides of the interval
\r
652 if (myV.Value(a) * myD.Value(a) > 0)
\r
654 if (fabs(myD.Value(a)) > Precision::Confusion())
\r
659 if (fabs(myV.Value(a)) > fabs(thePoint->myV.Value(aNear)))
\r
660 if (thePoint->myV.Value(aNear) * thePoint->myD.Value(aNear) < 0 &&
\r
661 fabs(thePoint->myD.Value(aNear)) > Precision::Confusion())
\r
671 if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
\r
673 if ((myV.Value(a) + myD.Value(a) * aDX) * myV.Value(a) > Precision::Confusion() &&
\r
674 (thePoint->myV.Value(aNear) + thePoint->myD.Value(aNear) * aDX) * thePoint->myV.Value(aNear) > Precision::Confusion())
\r
683 if ( fabs(aDX) < gp::Resolution() || fabs(aDiff / aDX) > 1.e+7)
\r
690 { // there is no near: remove it from the list
\r
698 Standard_Boolean aFound = Standard_False;
\r
699 for(b = 1; b <= myNear.Length(); b++)
\r
701 if (myNear.Value(b) == aNear)
\r
703 if (fabs(aDiffs.Value(b)) < fabs(aDiff))
\r
704 { // return this 'near'
\r
705 aFound = Standard_True;
\r
713 { // remove the old 'near'
\r
726 myNear.Append(aNear);
\r
727 aDiffs.Append(aDiff);
\r
733 //=======================================================================
\r
736 //=======================================================================
\r
737 GEOMImpl_Fillet1dPoint* GEOMImpl_Fillet1dPoint::Copy()
\r
739 GEOMImpl_Fillet1dPoint* aCopy = new GEOMImpl_Fillet1dPoint(myParam);
\r
740 Standard_Integer a;
\r
741 for(a = 1; a <= myV.Length(); a++)
\r
743 aCopy->myV.Append(myV.Value(a));
\r
744 aCopy->myD.Append(myD.Value(a));
\r
745 aCopy->myValid.Append(myValid.Value(a));
\r
750 //=======================================================================
\r
751 //function : HasSolution
\r
753 //=======================================================================
\r
754 Standard_Integer GEOMImpl_Fillet1dPoint::HasSolution(const Standard_Real theRadius)
\r
756 Standard_Integer a;
\r
757 for(a = 1; a <= myV.Length(); a++)
\r
759 if (fabs(sqrt(fabs(fabs(myV.Value(a)) + theRadius * theRadius)) - theRadius) < Precision::Confusion() / 10.)
\r
765 //=======================================================================
\r
766 //function : RemoveSolution
\r
768 //=======================================================================
\r
769 void GEOMImpl_Fillet1dPoint::RemoveSolution(Standard_Integer theIndex)
\r
771 myV.Remove(theIndex);
\r
772 myD.Remove(theIndex);
\r
773 myValid.Remove(theIndex);
\r
774 myNear.Remove(theIndex);
\r