1 // Copyright (C) 2007-2015 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, or (at your option) any later version.
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
23 #include "GEOMImpl_Fillet1d.hxx"
25 #include <BRep_Tool.hxx>
26 #include <BRepAdaptor_Curve.hxx>
27 #include <BRepBuilderAPI_MakeEdge.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 #include <IntRes2d_IntersectionSegment.hxx>
44 * This function returns Standard_True if it is possible to divide edge, i.e.
45 * if one parameter either start or end one is inside the edge. This function
46 * is used in the method GEOMImpl_Fillet1d::Result.
48 * \param theEdge the edge
49 * \param theStart the start parameter
50 * \param theEnd the end parameter
51 * \return Standard_True if it is possible to split edge;
52 * Standard_False otherwise.
54 static Standard_Boolean IsDivideEdge(const TopoDS_Edge &theEdge,
55 const Standard_Real theStart,
56 const Standard_Real theEnd)
60 Handle(Geom_Curve) aCurve = BRep_Tool::Curve(theEdge, aFirst, aLast);
61 gp_Pnt aPStart = aCurve->Value(theStart);
62 gp_Pnt aPEnd = aCurve->Value(theEnd);
63 TopoDS_Vertex aVFirst = TopExp::FirstVertex(theEdge);
64 TopoDS_Vertex aVLast = TopExp::FirstVertex(theEdge);
65 Standard_Real aTolFirst = BRep_Tool::Tolerance(aVFirst);
66 Standard_Real aTolLast = BRep_Tool::Tolerance(aVLast);
67 Standard_Real aTolConf = Precision::Confusion();
68 gp_Pnt aPFirst = BRep_Tool::Pnt(aVFirst);
69 gp_Pnt aPLast = BRep_Tool::Pnt(aVLast);
70 Standard_Real aDistSF = aPStart.Distance(aPFirst);
71 Standard_Real aDistSL = aPStart.Distance(aPLast);
72 Standard_Real aDistEF = aPEnd.Distance(aPFirst);
73 Standard_Real aDistEL = aPEnd.Distance(aPLast);
74 Standard_Boolean isSplit = Standard_True;
76 if (aDistSF <= aTolFirst + aTolConf ||
77 aDistSL <= aTolLast + aTolConf) {
78 if (aDistEF <= aTolFirst + aTolConf ||
79 aDistEL <= aTolLast + aTolConf) {
80 isSplit = Standard_False;
88 * class GEOMImpl_Fillet1d
91 //=======================================================================
92 //function : Constructor
94 //=======================================================================
95 GEOMImpl_Fillet1d::GEOMImpl_Fillet1d(const TopoDS_Edge& theEdge1,
96 const TopoDS_Edge& theEdge2,
97 const gp_Pln& thePlane)
98 : myEdgesExchnged( Standard_False )
100 myPlane = new Geom_Plane(thePlane);
102 BRepAdaptor_Curve aBAC1(theEdge1);
103 BRepAdaptor_Curve aBAC2(theEdge2);
104 if (aBAC1.GetType() < aBAC2.GetType())
105 { // first curve must be more complicated
108 myEdgesExchnged = Standard_True;
116 Handle(Geom_Curve) aCurve1 = BRep_Tool::Curve(myEdge1, myStart1, myEnd1);
117 Handle(Geom_Curve) aCurve2 = BRep_Tool::Curve(myEdge2, myStart2, myEnd2);
119 myCurve1 = GeomProjLib::Curve2d(aCurve1, myStart1, myEnd1, myPlane);
120 myCurve2 = GeomProjLib::Curve2d(aCurve2, myStart2, myEnd2, myPlane);
122 while (myCurve1->IsPeriodic() && myStart1 >= myEnd1)
123 myEnd1 += myCurve1->Period();
124 while (myCurve2->IsPeriodic() && myStart2 >= myEnd2)
125 myEnd2 += myCurve2->Period();
127 if (aBAC1.GetType() == aBAC2.GetType())
129 if (myEnd2 - myStart2 < myEnd1 - myStart1)
130 { // first curve must be parametrically shorter
131 TopoDS_Edge anEdge = myEdge1;
134 Handle(Geom2d_Curve) aCurve = myCurve1;
137 Standard_Real a = myStart1;
143 myEdgesExchnged = Standard_True;
148 //=======================================================================
149 //function : isRadiusIntersected
150 //purpose : local function
151 //=======================================================================
152 static Standard_Boolean isRadiusIntersected(const Handle(Geom2d_Curve)& theCurve,
153 const gp_Pnt2d theStart,
154 const gp_Pnt2d theEnd,
155 const Standard_Boolean theStartConnected)
157 const Standard_Real aTol = Precision::Confusion();
158 const Standard_Real anAngTol = Precision::Angular();
159 Geom2dAPI_InterCurveCurve anInter(theCurve, new Geom2d_Line(theStart,
160 gp_Dir2d(gp_Vec2d(theStart, theEnd))), aTol);
163 for(a = anInter.NbPoints(); a > 0; a--)
165 aPoint = anInter.Point(a);
166 if ( aPoint.Distance(theStart) < aTol && !theStartConnected )
167 return Standard_True;
168 if (aPoint.Distance(theEnd) < aTol * 200)
169 return Standard_True;
170 if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
171 return Standard_True;
173 Handle(Geom2d_Curve) aCurve;
174 for(a = anInter.NbSegments(); a > 0; a--)
176 anInter.Segment(a, aCurve);
177 aPoint = aCurve->Value(aCurve->FirstParameter());
178 if (aPoint.Distance(theStart) < aTol)
179 if (!theStartConnected)
180 return Standard_True;
181 if (aPoint.Distance(theEnd) < aTol)
182 return Standard_True;
183 if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
184 return Standard_True;
185 aPoint = aCurve->Value(aCurve->LastParameter());
186 if (aPoint.Distance(theStart) < aTol)
187 if (!theStartConnected)
188 return Standard_True;
189 if (aPoint.Distance(theEnd) < aTol)
190 return Standard_True;
191 if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
192 return Standard_True;
194 return Standard_False;
198 //=======================================================================
199 //function : fillPoint
201 //=======================================================================
202 void GEOMImpl_Fillet1d::fillPoint(GEOMImpl_Fillet1dPoint* thePoint)
206 const Standard_Real aTol = Precision::Confusion();
207 myCurve1->D1(thePoint->GetParam(), aPoint, aVec);
208 if (aVec.SquareMagnitude() < aTol)
211 gp_Vec2d aPerp(((myStartSide)?-1:1) * aVec.Y(), ((myStartSide)?1:-1) * aVec.X());
213 aPerp.Multiply(myRadius);
214 gp_Pnt2d aCenter = aPoint.Translated(aPerp);
215 thePoint->SetCenter(aCenter);
217 // on the intersection point
218 Standard_Boolean aValid = Standard_True;
219 Geom2dAPI_ProjectPointOnCurve aProjInt(aPoint, myCurve2);
220 if (aProjInt.NbPoints() && aPoint.Distance(aProjInt.NearestPoint()) < aTol)
221 aValid = Standard_False;
223 aValid = !isRadiusIntersected(myCurve2, aPoint, aCenter, Standard_True);
225 Geom2dAPI_ProjectPointOnCurve aProj(aCenter, myCurve2);
226 Standard_Integer a, aNB = aProj.NbPoints();
227 for(a = aNB; a > 0; a--)
229 if (aPoint.Distance(aProj.Point(a)) < aTol)
232 Standard_Boolean aValid2 = aValid;
234 aValid2 = !isRadiusIntersected(myCurve1, aCenter, aProj.Point(a), Standard_False);
236 // checking the right parameter
237 Standard_Real aParam = aProj.Parameter(a);
238 while(myCurve2->IsPeriodic() && aParam < myStart2)
239 aParam += myCurve2->Period();
241 thePoint->AddValue(aProj.Distance(a) * aProj.Distance(a) - myRadius * myRadius,
242 (aParam >= myStart2 && aParam <= myEnd2 && aValid2));
243 if (fabs(fabs(aProj.Distance(a)) - myRadius) < aTol)
244 thePoint->SetParam2(aParam);
248 //=======================================================================
249 //function : fillDiff
251 //=======================================================================
252 void GEOMImpl_Fillet1d::fillDiff(GEOMImpl_Fillet1dPoint* thePoint, Standard_Real theDiffStep, Standard_Boolean theFront)
254 GEOMImpl_Fillet1dPoint* aDiff =
255 new GEOMImpl_Fillet1dPoint(thePoint->GetParam() + (theFront?(theDiffStep):(-theDiffStep)));
257 if (!thePoint->ComputeDifference(aDiff))
259 aDiff->SetParam(thePoint->GetParam() + (theFront?(-theDiffStep):(theDiffStep)));
261 thePoint->ComputeDifference(aDiff);
266 //=======================================================================
269 //=======================================================================
270 Standard_Boolean GEOMImpl_Fillet1d::Perform(const Standard_Real theRadius)
272 myDegreeOfRecursion = 0;
273 myResultParams.Clear();
274 myResultOrientation.Clear();
276 Standard_Real aNBSteps = 100;
277 Geom2dAdaptor_Curve aGAC(myCurve1);
278 switch (aGAC.GetType())
286 case GeomAbs_Ellipse:
289 case GeomAbs_BezierCurve:
290 case GeomAbs_BSplineCurve:
291 aNBSteps = 2 + aGAC.Degree() * aGAC.NbPoles();
293 default: // unknown: maximum
297 myRadius = theRadius;
299 // Compute the intervals.
300 const Standard_Real aTol = Precision::Confusion();
301 Geom2dAPI_InterCurveCurve anAPIInter(myCurve1, myCurve2, aTol);
302 const Geom2dInt_GInter &anInter = anAPIInter.Intersector();
303 Standard_Integer aNb = anInter.NbPoints();
305 TColStd_ListOfReal aParams;
306 TColStd_ListIteratorOfListOfReal anIter;
308 // Treat intersection points.
309 for(i = 1; i <= aNb; i++) {
310 const IntRes2d_IntersectionPoint &aPoint = anInter.Point(i);
311 Standard_Real aParam = aPoint.ParamOnFirst();
313 // Adjust parameter on periodic curve.
314 if (myCurve1->IsPeriodic()) {
315 aParam = ElCLib::InPeriod
316 (aParam, myStart1, myStart1 + myCurve1->Period());
319 if (aParam > myStart1 + aTol && aParam < myEnd1 - aTol) {
320 // Add the point in the list in increasing order.
321 for(anIter.Initialize(aParams); anIter.More(); anIter.Next()) {
322 if (anIter.Value() > aParam) {
323 aParams.InsertBefore(aParam, anIter);
328 if (!anIter.More()) {
329 aParams.Append(aParam);
334 // Treat intersection segments.
335 aNb = anInter.NbSegments();
337 for(i = 1; i <= aNb; i++) {
338 const IntRes2d_IntersectionSegment &aSegment = anInter.Segment(i);
340 if (aSegment.HasFirstPoint() && aSegment.HasLastPoint()) {
341 Standard_Real aParam1 = aSegment.FirstPoint().ParamOnFirst();
342 Standard_Real aParam2 = aSegment.LastPoint().ParamOnFirst();
344 // Adjust parameters on periodic curve.
345 if (myCurve1->IsPeriodic()) {
346 ElCLib::AdjustPeriodic(myStart1, myStart1 + myCurve1->Period(),
347 aTol, aParam1, aParam2);
350 if (aParam1 > myStart1 + aTol && aParam1 < myEnd1 - aTol &&
351 aParam2 > myStart1 + aTol && aParam2 < myEnd1 - aTol) {
352 // Add the point in the list in increasing order.
353 const Standard_Real aParam = 0.5*(aParam1 + aParam2);
355 for(anIter.Initialize(aParams); anIter.More(); anIter.Next()) {
356 if (anIter.Value() > aParam) {
357 aParams.InsertBefore(aParam, anIter);
362 if (!anIter.More()) {
363 aParams.Append(aParam);
369 // Add start and end parameters to the list.
370 aParams.Prepend(myStart1);
371 aParams.Append(myEnd1);
372 anIter.Initialize(aParams);
374 // Perform each interval.
375 Standard_Real aStart = anIter.Value();
377 for (anIter.Next(); anIter.More(); anIter.Next()) {
378 const Standard_Real anEnd = anIter.Value();
380 // Perform the interval.
381 performInterval(aStart, anEnd, aNBSteps);
385 if (myResultParams.Extent())
386 return Standard_True;
388 return Standard_False;
391 //=======================================================================
392 //function : performInterval
394 //=======================================================================
395 void GEOMImpl_Fillet1d::performInterval(const Standard_Real theStart,
396 const Standard_Real theEnd,
397 const Standard_Integer theNBSteps)
399 Standard_Real aParam, aStep, aDStep;
400 aStep = (theEnd - theStart) / theNBSteps;
401 aDStep = aStep/1000.;
403 Standard_Integer aCycle;
404 for(aCycle = 2, myStartSide = Standard_False; aCycle; myStartSide = !myStartSide, aCycle--)
406 GEOMImpl_Fillet1dPoint *aLeft = NULL, *aRight = NULL;
408 for(aParam = theStart + aStep; aParam < theEnd || fabs(theEnd - aParam) < Precision::Confusion(); aParam += aStep)
412 aLeft = new GEOMImpl_Fillet1dPoint(aParam - aStep);
414 fillDiff(aLeft, aDStep, Standard_True);
417 aRight = new GEOMImpl_Fillet1dPoint(aParam);
419 fillDiff(aRight, aDStep, Standard_False);
421 aLeft->FilterPoints(aRight);
422 performNewton(aLeft, aRight);
431 //=======================================================================
432 //function : processPoint
434 //=======================================================================
435 Standard_Boolean GEOMImpl_Fillet1d::processPoint(GEOMImpl_Fillet1dPoint* theLeft,
436 GEOMImpl_Fillet1dPoint* theRight,
437 Standard_Real theParameter)
439 if (theParameter >= theLeft->GetParam() && theParameter < theRight->GetParam())
441 Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
442 if (theParameter - theLeft->GetParam() < aDX / 100.)
444 theParameter = theLeft->GetParam() + aDX / 100.;
446 if (theRight->GetParam() - theParameter < aDX / 100.)
448 theParameter = theRight->GetParam() - aDX / 100.;
451 // Protection on infinite loop.
452 myDegreeOfRecursion++;
453 Standard_Real diffx = 0.001 * aDX;
454 if (myDegreeOfRecursion > 1000)
457 if (myDegreeOfRecursion > 10000)
460 if (myDegreeOfRecursion > 100000)
462 return Standard_True;
467 GEOMImpl_Fillet1dPoint* aPoint1 = theLeft->Copy();
468 GEOMImpl_Fillet1dPoint* aPoint2 = new GEOMImpl_Fillet1dPoint(theParameter);
470 fillDiff(aPoint2, diffx, Standard_True);
472 aPoint1->FilterPoints(aPoint2);
473 performNewton(aPoint1, aPoint2);
474 aPoint2->FilterPoints(theRight);
475 performNewton(aPoint2, theRight);
479 return Standard_True;
482 return Standard_False;
485 //=======================================================================
486 //function : performNewton
488 //=======================================================================
489 void GEOMImpl_Fillet1d::performNewton(GEOMImpl_Fillet1dPoint* theLeft,
490 GEOMImpl_Fillet1dPoint* theRight)
493 // check the left: if this is solution store it and remove it from the list of researching points of theLeft
494 a = theLeft->HasSolution(myRadius);
497 if (theLeft->IsValid(a))
499 myResultParams.Append(theLeft->GetParam());
500 myResultOrientation.Append(myStartSide);
505 Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
506 if ( aDX < Precision::Confusion() / 1000000.)
508 a = theRight->HasSolution(myRadius);
510 if (theRight->IsValid(a))
512 myResultParams.Append(theRight->GetParam());
513 myResultOrientation.Append(myStartSide);
518 for(a = 1; a <= theLeft->GetNBValues(); a++)
520 Standard_Integer aNear = theLeft->GetNear(a);
522 Standard_Real aA = (theRight->GetDiff(aNear) - theLeft->GetDiff(a)) / aDX;
523 Standard_Real aB = theLeft->GetDiff(a) - aA * theLeft->GetParam();
524 Standard_Real aC = theLeft->GetValue(a) - theLeft->GetDiff(a) * theLeft->GetParam() +
525 aA * theLeft->GetParam() * theLeft->GetParam() / 2.0;
526 Standard_Real aDet = aB * aB - 2.0 * aA * aC;
528 if ( fabs(aDet) < gp::Resolution() )
531 if (fabs(aA) < Precision::Confusion())
533 if (fabs(aB) > 10e-20)
535 Standard_Real aX0 = - aC / aB; // use extremum
536 if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
537 processPoint(theLeft, theRight, aX0);
541 processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
546 if (fabs(aB) > fabs(aDet * 1000000.))
547 { // possible floating point operations accuracy errors
548 processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
555 Standard_Boolean aRes = processPoint(theLeft, theRight, (- aB + aDet) / aA);
557 aRes = processPoint(theLeft, theRight, (- aB - aDet) / aA);
559 processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
563 Standard_Real aX0 = - aB / aA; // use extremum
564 if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
565 processPoint(theLeft, theRight, aX0);
567 processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
574 //=======================================================================
577 //=======================================================================
578 TopoDS_Edge GEOMImpl_Fillet1d::Result(const gp_Pnt& thePoint,
579 TopoDS_Edge& theEdge1,
580 TopoDS_Edge& theEdge2)
583 gp_Pnt2d aTargetPoint2d;
584 Standard_Real aX, aY;
585 ElSLib::PlaneParameters(myPlane->Pln().Position(), thePoint, aX, aY);
586 aTargetPoint2d.SetCoord(aX, aY);
588 // choose the nearest circle
589 Standard_Real aDistance, aP;
590 GEOMImpl_Fillet1dPoint *aNearest;
592 TColStd_ListIteratorOfListOfReal anIter(myResultParams);
593 for(aNearest = NULL, a = 1; anIter.More(); anIter.Next(), a++)
595 myStartSide = (myResultOrientation.Value(a)) ? Standard_True : Standard_False;
596 GEOMImpl_Fillet1dPoint *aPoint = new GEOMImpl_Fillet1dPoint(anIter.Value());
598 if (!aPoint->HasSolution(myRadius))
600 aP = fabs(aPoint->GetCenter().Distance(aTargetPoint2d) - myRadius);
601 if (!aNearest || aP < aDistance)
615 // create circle edge
616 gp_Pnt aCenter = ElSLib::PlaneValue(aNearest->GetCenter().X(),
617 aNearest->GetCenter().Y(),
618 myPlane->Pln().Position());
619 Handle(Geom_Circle) aCircle =
620 new Geom_Circle(gp_Ax2(aCenter, myPlane->Pln().Axis().Direction()), myRadius);
621 gp_Pnt2d aPoint2d1, aPoint2d2;
622 myCurve1->D0(aNearest->GetParam(), aPoint2d1);
623 myCurve2->D0(aNearest->GetParam2(), aPoint2d2);
624 gp_Pnt aPoint1 = ElSLib::PlaneValue(aPoint2d1.X(), aPoint2d1.Y(), myPlane->Pln().Position());
625 gp_Pnt aPoint2 = ElSLib::PlaneValue(aPoint2d2.X(), aPoint2d2.Y(), myPlane->Pln().Position());
627 GeomAPI_ProjectPointOnCurve aProj(thePoint, aCircle);
628 Standard_Real aTarGetParam = aProj.LowerDistanceParameter();
629 gp_Pnt aPointOnCircle = aProj.NearestPoint();
631 // Check extrema point manually, because there is a bug in Open CASCADE
632 // in calculation of nearest point to a circle near the parameter 0.0
633 gp_Pnt p0 = ElCLib::Value(0.0, aCircle->Circ());
634 if (p0.Distance(thePoint) < aPointOnCircle.Distance(thePoint))
640 aProj.Perform(aPoint1);
641 Standard_Real aParam1 = aProj.LowerDistanceParameter();
642 aProj.Perform(aPoint2);
643 Standard_Real aParam2 = aProj.LowerDistanceParameter();
644 Standard_Boolean aIsOut = ((aParam1 < aTarGetParam && aParam2 < aTarGetParam) ||
645 (aParam1 > aTarGetParam && aParam2 > aTarGetParam));
646 if (aParam1 > aParam2)
648 BRepBuilderAPI_MakeEdge aBuilder(aCircle->Circ(),
649 aIsOut ? aParam2 : aParam1,
650 aIsOut? aParam1 : aParam2);
651 aResult = aBuilder.Edge();
654 Standard_Real aStart, anEnd;
655 Handle(Geom_Curve) aCurve = BRep_Tool::Curve(myEdge1, aStart, anEnd);
657 aCurve->D1(aNearest->GetParam(), aPoint1, aDir);
660 aCircle->D1(aParam1, aPoint1, aCircleDir);
661 if ((aCircleDir.Angle(aDir) > M_PI / 2.0) ^ aIsOut)
662 aStart = aNearest->GetParam();
664 anEnd = aNearest->GetParam();
666 if (IsDivideEdge(myEdge1, aStart, anEnd))
669 BRepBuilderAPI_MakeEdge aDivider1(aCurve, aStart, anEnd);
671 theEdge2 = aDivider1.Edge();
673 theEdge1 = aDivider1.Edge();
676 aCurve = BRep_Tool::Curve(myEdge2, aStart, anEnd);
677 aCurve->D1(aNearest->GetParam2(), aPoint2, aDir);
679 aCircle->D1(aParam2, aPoint2, aCircleDir);
680 if ((aCircleDir.Angle(aDir) > M_PI / 2.0) ^ (!aIsOut))
681 aStart = aNearest->GetParam2();
683 anEnd = aNearest->GetParam2();
685 if (IsDivideEdge(myEdge2, aStart, anEnd))
687 BRepBuilderAPI_MakeEdge aDivider2(aCurve, aStart, anEnd);
689 theEdge1 = aDivider2.Edge();
691 theEdge2 = aDivider2.Edge();
698 //=======================================================================
699 //function : AddValue
701 //=======================================================================
702 void GEOMImpl_Fillet1dPoint::AddValue(Standard_Real theValue, Standard_Boolean theValid)
705 for(a = 1; a <= myV.Length(); a++)
707 if (theValue < myV.Value(a))
709 myV.InsertBefore(a, theValue);
710 myValid.InsertBefore(a, (Standard_Integer)theValid);
714 myV.Append(theValue);
715 myValid.Append((Standard_Integer)theValid);
718 //=======================================================================
719 //function : ComputeDifference
721 //=======================================================================
722 Standard_Boolean GEOMImpl_Fillet1dPoint::ComputeDifference(GEOMImpl_Fillet1dPoint* thePoint)
725 Standard_Boolean aDiffsSet = (myD.Length() != 0);
726 Standard_Real aDX = thePoint->GetParam() - myParam, aDY;
727 if (thePoint->myV.Length() == myV.Length())
728 { // absolutely the same points
729 for(a = 1; a <= myV.Length(); a++)
731 aDY = thePoint->myV.Value(a) - myV.Value(a);
733 myD.SetValue(a, fabs(aDX) > gp::Resolution() ? (aDY/aDX) : 0);
735 myD.Append( fabs(aDX) > gp::Resolution() ? (aDY/aDX) : 0);
737 return Standard_True;
739 // between the diffeerent points searching for nearest analogs
741 for(a = 1; a <= myV.Length(); a++)
743 for(b = 1; b <= thePoint->myV.Length(); b++)
745 if (b == 1 || fabs(thePoint->myV.Value(b) - myV.Value(a)) < fabs(aDY))
746 aDY = thePoint->myV.Value(b) - myV.Value(a);
750 if ( fabs(aDX) > gp::Resolution() && fabs(aDY / aDX) < fabs(myD.Value(a)))
751 myD.SetValue(a, aDY / aDX);
757 myD.Append( fabs(aDX) > gp::Resolution() ? aDY/aDX : 0);
761 return Standard_False;
764 //=======================================================================
765 //function : FilterPoints
767 //=======================================================================
768 void GEOMImpl_Fillet1dPoint::FilterPoints(GEOMImpl_Fillet1dPoint* thePoint)
770 Standard_Integer a, b;
771 TColStd_SequenceOfReal aDiffs;
772 Standard_Real aY, aY2, aDX = thePoint->GetParam() - myParam;
773 for(a = 1; a <= myV.Length(); a++)
775 // searching for near point from thePoint
776 Standard_Integer aNear = 0;
777 Standard_Real aDiff = aDX * 10000.;
778 aY = myV.Value(a) + myD.Value(a) * aDX;
779 for(b = 1; b <= thePoint->myV.Length(); b++)
781 // calculate hypothesis value of the Y2 with the constant first and second derivative
782 aY2 = aY + aDX * (thePoint->myD.Value(b) - myD.Value(a)) / 2.0;
783 if (aNear == 0 || fabs(aY2 - thePoint->myV.Value(b)) < fabs(aDiff))
786 aDiff = aY2 - thePoint->myV.Value(b);
792 if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
793 {// the same sign at the same sides of the interval
794 if (myV.Value(a) * myD.Value(a) > 0)
796 if (fabs(myD.Value(a)) > Precision::Confusion())
801 if (fabs(myV.Value(a)) > fabs(thePoint->myV.Value(aNear)))
802 if (thePoint->myV.Value(aNear) * thePoint->myD.Value(aNear) < 0 &&
803 fabs(thePoint->myD.Value(aNear)) > Precision::Confusion())
813 if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
815 if ((myV.Value(a) + myD.Value(a) * aDX) * myV.Value(a) > Precision::Confusion() &&
816 (thePoint->myV.Value(aNear) + thePoint->myD.Value(aNear) * aDX) * thePoint->myV.Value(aNear) > Precision::Confusion())
825 if ( fabs(aDX) < gp::Resolution() || fabs(aDiff / aDX) > 1.e+7)
832 { // there is no near: remove it from the list
840 Standard_Boolean aFound = Standard_False;
841 for(b = 1; b <= myNear.Length(); b++)
843 if (myNear.Value(b) == aNear)
845 if (fabs(aDiffs.Value(b)) < fabs(aDiff))
846 { // return this 'near'
847 aFound = Standard_True;
855 { // remove the old 'near'
868 myNear.Append(aNear);
869 aDiffs.Append(aDiff);
875 //=======================================================================
878 //=======================================================================
879 GEOMImpl_Fillet1dPoint* GEOMImpl_Fillet1dPoint::Copy()
881 GEOMImpl_Fillet1dPoint* aCopy = new GEOMImpl_Fillet1dPoint(myParam);
883 for(a = 1; a <= myV.Length(); a++)
885 aCopy->myV.Append(myV.Value(a));
886 aCopy->myD.Append(myD.Value(a));
887 aCopy->myValid.Append(myValid.Value(a));
892 //=======================================================================
893 //function : HasSolution
895 //=======================================================================
896 Standard_Integer GEOMImpl_Fillet1dPoint::HasSolution(const Standard_Real theRadius)
899 for(a = 1; a <= myV.Length(); a++)
901 if (fabs(sqrt(fabs(fabs(myV.Value(a)) + theRadius * theRadius)) - theRadius) < Precision::Confusion() / 10.)
907 //=======================================================================
908 //function : RemoveSolution
910 //=======================================================================
911 void GEOMImpl_Fillet1dPoint::RemoveSolution(Standard_Integer theIndex)
913 myV.Remove(theIndex);
914 myD.Remove(theIndex);
915 myValid.Remove(theIndex);
916 myNear.Remove(theIndex);