Back to index

salome-geom  6.5.0
GEOMImpl_Fillet1d.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 //  File   : GEOMImpl_Fillet1d.cxx
00021 //  Module : GEOMImpl
00022 
00023 #include "GEOMImpl_Fillet1d.hxx"
00024 
00025 #include <BRep_Tool.hxx>
00026 #include <BRepAdaptor_Curve.hxx>
00027 #include <BRepBuilderAPI_MakeEdge.hxx>
00028 #include <ElCLib.hxx>
00029 #include <ElSLib.hxx>
00030 #include <gp_Circ.hxx>
00031 #include <Geom2d_Line.hxx>
00032 #include <Geom2dAPI_ProjectPointOnCurve.hxx>
00033 #include <Geom2dAPI_InterCurveCurve.hxx>
00034 #include <GeomAPI_ProjectPointOnCurve.hxx>
00035 #include <GeomProjLib.hxx>
00036 #include <Geom_Circle.hxx>
00037 #include <Precision.hxx>
00038 #include <TColStd_ListIteratorOfListOfReal.hxx>
00039 
00044 //=======================================================================
00045 //function : Constructor
00046 //purpose  :
00047 //=======================================================================
00048 GEOMImpl_Fillet1d::GEOMImpl_Fillet1d(const TopoDS_Edge& theEdge1,
00049                                      const TopoDS_Edge& theEdge2,
00050                                      const gp_Pln& thePlane)
00051 : myEdgesExchnged( Standard_False )
00052 {
00053   myPlane = new Geom_Plane(thePlane);
00054 
00055   BRepAdaptor_Curve aBAC1(theEdge1);
00056   BRepAdaptor_Curve aBAC2(theEdge2);
00057   if (aBAC1.GetType() < aBAC2.GetType())
00058   { // first curve must be more complicated
00059     myEdge1 = theEdge2;
00060     myEdge2 = theEdge1;
00061     myEdgesExchnged = Standard_True;
00062   }
00063   else
00064   {
00065     myEdge1 = theEdge1;
00066     myEdge2 = theEdge2;
00067   }
00068 
00069   Handle(Geom_Curve) aCurve1 = BRep_Tool::Curve(myEdge1, myStart1, myEnd1);
00070   Handle(Geom_Curve) aCurve2 = BRep_Tool::Curve(myEdge2, myStart2, myEnd2);
00071 
00072   myCurve1 = GeomProjLib::Curve2d(aCurve1, myStart1, myEnd1, myPlane);
00073   myCurve2 = GeomProjLib::Curve2d(aCurve2, myStart2, myEnd2, myPlane);
00074 
00075   while (myCurve1->IsPeriodic() && myStart1 >= myEnd1)
00076     myEnd1 += myCurve1->Period();
00077   while (myCurve2->IsPeriodic() && myStart2 >= myEnd2)
00078     myEnd2 += myCurve2->Period();
00079 
00080   if (aBAC1.GetType() == aBAC2.GetType())
00081   {
00082     if (myEnd2 - myStart2 < myEnd1 - myStart1)
00083     { // first curve must be parametrically shorter
00084       TopoDS_Edge anEdge = myEdge1;
00085       myEdge1 = myEdge2;
00086       myEdge2 = anEdge;
00087       Handle(Geom2d_Curve) aCurve = myCurve1;
00088       myCurve1 = myCurve2;
00089       myCurve2 = aCurve;
00090       Standard_Real a = myStart1;
00091       myStart1 = myStart2;
00092       myStart2 = a;
00093       a = myEnd1;
00094       myEnd1 = myEnd2;
00095       myEnd2 = a;
00096       myEdgesExchnged = Standard_True;
00097     }
00098   }
00099 }
00100 
00101 //=======================================================================
00102 //function : isRadiusIntersected
00103 //purpose  : local function
00104 //=======================================================================
00105 static Standard_Boolean isRadiusIntersected(const Handle(Geom2d_Curve)& theCurve,
00106                                             const gp_Pnt2d theStart,
00107                                             const gp_Pnt2d theEnd,
00108                                             const Standard_Boolean theStartConnected)
00109 {
00110   const Standard_Real aTol = Precision::Confusion();
00111   const Standard_Real anAngTol = Precision::Angular();
00112   Geom2dAPI_InterCurveCurve anInter(theCurve, new Geom2d_Line(theStart,
00113     gp_Dir2d(gp_Vec2d(theStart, theEnd))), aTol);
00114   Standard_Integer a;
00115   gp_Pnt2d aPoint;
00116   for(a = anInter.NbPoints(); a > 0; a--)
00117   {
00118     aPoint = anInter.Point(a);
00119     if ( aPoint.Distance(theStart) < aTol && !theStartConnected )
00120       return Standard_True;
00121     if (aPoint.Distance(theEnd) < aTol * 200)
00122       return Standard_True;
00123     if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
00124       return Standard_True;
00125   }
00126   Handle(Geom2d_Curve) aCurve;
00127   for(a = anInter.NbSegments(); a > 0; a--)
00128   {
00129     anInter.Segment(a, aCurve);
00130     aPoint = aCurve->Value(aCurve->FirstParameter());
00131     if (aPoint.Distance(theStart) < aTol)
00132       if (!theStartConnected)
00133         return Standard_True;
00134     if (aPoint.Distance(theEnd) < aTol)
00135       return Standard_True;
00136     if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
00137       return Standard_True;
00138     aPoint = aCurve->Value(aCurve->LastParameter());
00139     if (aPoint.Distance(theStart) < aTol)
00140       if (!theStartConnected)
00141         return Standard_True;
00142     if (aPoint.Distance(theEnd) < aTol)
00143       return Standard_True;
00144     if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
00145       return Standard_True;
00146   }
00147   return Standard_False;
00148 }
00149 
00150 
00151 //=======================================================================
00152 //function : fillPoint
00153 //purpose  :
00154 //=======================================================================
00155 void GEOMImpl_Fillet1d::fillPoint(GEOMImpl_Fillet1dPoint* thePoint)
00156 {
00157   gp_Pnt2d aPoint;
00158   gp_Vec2d aVec;
00159   const Standard_Real aTol = Precision::Confusion();
00160   myCurve1->D1(thePoint->GetParam(), aPoint, aVec);
00161   if (aVec.SquareMagnitude() < aTol)
00162     return;
00163 
00164   gp_Vec2d aPerp(((myStartSide)?-1:1) * aVec.Y(), ((myStartSide)?1:-1) * aVec.X());
00165   aPerp.Normalize();
00166   aPerp.Multiply(myRadius);
00167   gp_Pnt2d aCenter = aPoint.Translated(aPerp);
00168   thePoint->SetCenter(aCenter);
00169 
00170   // on the intersection point
00171   Standard_Boolean aValid = Standard_True;
00172   Geom2dAPI_ProjectPointOnCurve aProjInt(aPoint, myCurve2);
00173   if (aProjInt.NbPoints() && aPoint.Distance(aProjInt.NearestPoint()) < aTol)
00174     aValid = Standard_False;
00175   else
00176     aValid = !isRadiusIntersected(myCurve2, aPoint, aCenter, Standard_True);
00177 
00178   Geom2dAPI_ProjectPointOnCurve aProj(aCenter, myCurve2);
00179   Standard_Integer a, aNB = aProj.NbPoints();
00180   for(a = aNB; a > 0; a--)
00181   {
00182     if (aPoint.Distance(aProj.Point(a)) < aTol)
00183       continue;
00184 
00185     Standard_Boolean aValid2 = aValid;
00186     if (aValid2)
00187       aValid2 = !isRadiusIntersected(myCurve1, aCenter, aProj.Point(a), Standard_False);
00188 
00189     // checking the right parameter
00190     Standard_Real aParam = aProj.Parameter(a);
00191     while(myCurve2->IsPeriodic() && aParam < myStart2)
00192       aParam += myCurve2->Period();
00193 
00194     thePoint->AddValue(aProj.Distance(a) * aProj.Distance(a) - myRadius * myRadius,
00195                        (aParam >= myStart2 && aParam <= myEnd2 && aValid2));
00196     if (fabs(fabs(aProj.Distance(a)) - myRadius) < aTol)
00197       thePoint->SetParam2(aParam);
00198   }
00199 }
00200 
00201 //=======================================================================
00202 //function : fillDiff
00203 //purpose  :
00204 //=======================================================================
00205 void GEOMImpl_Fillet1d::fillDiff(GEOMImpl_Fillet1dPoint* thePoint, Standard_Real theDiffStep, Standard_Boolean theFront)
00206 {
00207   GEOMImpl_Fillet1dPoint* aDiff =
00208     new GEOMImpl_Fillet1dPoint(thePoint->GetParam() + (theFront?(theDiffStep):(-theDiffStep)));
00209   fillPoint(aDiff);
00210   if (!thePoint->ComputeDifference(aDiff))
00211   {
00212     aDiff->SetParam(thePoint->GetParam() + (theFront?(-theDiffStep):(theDiffStep)));
00213     fillPoint(aDiff);
00214     thePoint->ComputeDifference(aDiff);
00215   }
00216   delete aDiff;
00217 }
00218 
00219 //=======================================================================
00220 //function : Perform
00221 //purpose  :
00222 //=======================================================================
00223 Standard_Boolean GEOMImpl_Fillet1d::Perform(const Standard_Real theRadius)
00224 {
00225   myDegreeOfRecursion = 0;
00226   myResultParams.Clear();
00227   myResultOrientation.Clear();
00228 
00229   Standard_Real aNBSteps = 100;
00230   Geom2dAdaptor_Curve aGAC(myCurve1);
00231   switch (aGAC.GetType())
00232   {
00233     case GeomAbs_Line:
00234       aNBSteps = 2;
00235       break;
00236     case GeomAbs_Circle:
00237       aNBSteps = 4;
00238       break;
00239     case GeomAbs_Ellipse:
00240       aNBSteps = 5;
00241       break;
00242     case GeomAbs_BezierCurve:
00243     case GeomAbs_BSplineCurve:
00244       aNBSteps = 2 + aGAC.Degree() * aGAC.NbPoles();
00245       break;
00246     default: // unknown: maximum
00247       aNBSteps = 100;
00248   }
00249 
00250   myRadius = theRadius;
00251   Standard_Real aParam, aStep, aDStep;
00252   aStep = (myEnd1 - myStart1) / aNBSteps;
00253   aDStep = aStep/1000.;
00254 
00255   Standard_Integer aCycle;
00256   for(aCycle = 2, myStartSide = Standard_False; aCycle; myStartSide = !myStartSide, aCycle--)
00257   {
00258     GEOMImpl_Fillet1dPoint *aLeft = NULL, *aRight = NULL;
00259 
00260     for(aParam = myStart1 + aStep; aParam < myEnd1 || fabs(myEnd1 - aParam) < Precision::Confusion(); aParam += aStep)
00261     {
00262       if (!aLeft)
00263       {
00264         aLeft = new GEOMImpl_Fillet1dPoint(aParam - aStep);
00265         fillPoint(aLeft);
00266         fillDiff(aLeft, aDStep, Standard_True);
00267       }
00268 
00269       aRight = new GEOMImpl_Fillet1dPoint(aParam);
00270       fillPoint(aRight);
00271       fillDiff(aRight, aDStep, Standard_False);
00272 
00273       aLeft->FilterPoints(aRight);
00274       performNewton(aLeft, aRight);
00275 
00276       delete aLeft;
00277       aLeft = aRight;
00278     }
00279     delete aLeft;
00280   }
00281 
00282   if (myResultParams.Extent())
00283     return Standard_True;
00284 
00285   return Standard_False;
00286 }
00287 
00288 //=======================================================================
00289 //function : processPoint
00290 //purpose  :
00291 //=======================================================================
00292 Standard_Boolean GEOMImpl_Fillet1d::processPoint(GEOMImpl_Fillet1dPoint* theLeft,
00293                                                  GEOMImpl_Fillet1dPoint* theRight,
00294                                                  Standard_Real           theParameter)
00295 {
00296   if (theParameter >= theLeft->GetParam() && theParameter < theRight->GetParam())
00297   {
00298     Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
00299     if (theParameter - theLeft->GetParam() < aDX / 100.)
00300     {
00301       theParameter = theLeft->GetParam() + aDX / 100.;
00302     }
00303     if (theRight->GetParam() - theParameter < aDX / 100.)
00304     {
00305       theParameter = theRight->GetParam() - aDX / 100.;
00306     }
00307 
00308     // Protection on infinite loop.
00309     myDegreeOfRecursion++;
00310     Standard_Real diffx = 0.001 * aDX;
00311     if (myDegreeOfRecursion > 1000)
00312     {
00313         diffx *= 10.0;
00314         if (myDegreeOfRecursion > 10000)
00315         {
00316             diffx *= 10.0;
00317             if (myDegreeOfRecursion > 100000)
00318             {
00319                 return Standard_True;
00320             }
00321         }
00322     }
00323 
00324     GEOMImpl_Fillet1dPoint* aPoint1 = theLeft->Copy();
00325     GEOMImpl_Fillet1dPoint* aPoint2 = new GEOMImpl_Fillet1dPoint(theParameter);
00326     fillPoint(aPoint2);
00327     fillDiff(aPoint2, diffx, Standard_True);
00328 
00329     aPoint1->FilterPoints(aPoint2);
00330     performNewton(aPoint1, aPoint2);
00331     aPoint2->FilterPoints(theRight);
00332     performNewton(aPoint2, theRight);
00333 
00334     delete aPoint1;
00335     delete aPoint2;
00336     return Standard_True;
00337   }
00338 
00339   return Standard_False;
00340 }
00341 
00342 //=======================================================================
00343 //function : performNewton
00344 //purpose  :
00345 //=======================================================================
00346 void GEOMImpl_Fillet1d::performNewton(GEOMImpl_Fillet1dPoint* theLeft,
00347                                       GEOMImpl_Fillet1dPoint* theRight)
00348 {
00349   Standard_Integer a;
00350   // check the left: if this is solution store it and remove it from the list of researching points of theLeft
00351   a = theLeft->HasSolution(myRadius);
00352   if (a)
00353   {
00354     if (theLeft->IsValid(a))
00355     {
00356       myResultParams.Append(theLeft->GetParam());
00357       myResultOrientation.Append(myStartSide);
00358     }
00359     return;
00360   }
00361 
00362   Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
00363   if ( aDX < Precision::Confusion() / 1000000.)
00364   {
00365     a = theRight->HasSolution(myRadius);
00366     if (a)
00367       if (theRight->IsValid(a))
00368       {
00369         myResultParams.Append(theRight->GetParam());
00370         myResultOrientation.Append(myStartSide);
00371       }
00372     return;
00373   }
00374 
00375   for(a = 1; a <= theLeft->GetNBValues(); a++)
00376   {
00377     Standard_Integer aNear = theLeft->GetNear(a);
00378 
00379     Standard_Real aA = (theRight->GetDiff(aNear) - theLeft->GetDiff(a)) / aDX;
00380     Standard_Real aB = theLeft->GetDiff(a) - aA * theLeft->GetParam();
00381     Standard_Real aC = theLeft->GetValue(a) - theLeft->GetDiff(a) * theLeft->GetParam() +
00382                        aA * theLeft->GetParam() * theLeft->GetParam() / 2.0;
00383     Standard_Real aDet = aB * aB - 2.0 * aA * aC;
00384 
00385     if ( fabs(aDet) < gp::Resolution() )
00386       continue;
00387 
00388     if (fabs(aA) < Precision::Confusion())
00389     { // linear case
00390       if (fabs(aB) > 10e-20)
00391       {
00392         Standard_Real aX0 = - aC / aB; // use extremum
00393         if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
00394           processPoint(theLeft, theRight, aX0);
00395       }
00396       else
00397       {
00398         processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
00399       }
00400     }
00401     else
00402     {
00403       if (fabs(aB) > fabs(aDet * 1000000.))
00404       {  // possible floating point operations accurancy errors
00405         processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
00406       }
00407       else
00408       {
00409         if (aDet > 0)
00410         { // two solutions
00411           aDet = sqrt(aDet);
00412           Standard_Boolean aRes = processPoint(theLeft, theRight, (- aB + aDet) / aA);
00413           if (!aRes)
00414             aRes = processPoint(theLeft, theRight, (- aB - aDet) / aA);
00415           if (!aRes)
00416             processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
00417         }
00418         else
00419         {
00420           Standard_Real aX0 = - aB / aA; // use extremum
00421           if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
00422             processPoint(theLeft, theRight, aX0);
00423           else
00424             processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
00425         }
00426       }
00427     }
00428   }
00429 }
00430 
00431 //=======================================================================
00432 //function : Result
00433 //purpose  :
00434 //=======================================================================
00435 TopoDS_Edge GEOMImpl_Fillet1d::Result(const gp_Pnt& thePoint,
00436                                       TopoDS_Edge& theEdge1,
00437                                       TopoDS_Edge& theEdge2)
00438 {
00439   TopoDS_Edge aResult;
00440   gp_Pnt2d aTargetPoint2d;
00441   Standard_Real aX, aY;
00442   ElSLib::PlaneParameters(myPlane->Pln().Position(), thePoint, aX, aY);
00443   aTargetPoint2d.SetCoord(aX, aY);
00444 
00445   // choose the nearest circle
00446   Standard_Real aDistance, aP;
00447   GEOMImpl_Fillet1dPoint *aNearest;
00448   Standard_Integer a;
00449   TColStd_ListIteratorOfListOfReal anIter(myResultParams);
00450   for(aNearest = NULL, a = 1; anIter.More(); anIter.Next(), a++)
00451   {
00452     myStartSide = (myResultOrientation.Value(a)) ? Standard_True : Standard_False;
00453     GEOMImpl_Fillet1dPoint *aPoint = new GEOMImpl_Fillet1dPoint(anIter.Value());
00454     fillPoint(aPoint);
00455     if (!aPoint->HasSolution(myRadius))
00456       continue;
00457     aP = fabs(aPoint->GetCenter().Distance(aTargetPoint2d) - myRadius);
00458     if (!aNearest || aP < aDistance)
00459     {
00460       aNearest = aPoint;
00461       aDistance = aP;
00462     }
00463     else
00464     {
00465       delete aPoint;
00466     }
00467    }
00468 
00469   if (!aNearest)
00470      return aResult;
00471 
00472   // create circle edge
00473   gp_Pnt aCenter = ElSLib::PlaneValue(aNearest->GetCenter().X(),
00474                                       aNearest->GetCenter().Y(),
00475                                       myPlane->Pln().Position());
00476   Handle(Geom_Circle) aCircle =
00477     new Geom_Circle(gp_Ax2(aCenter, myPlane->Pln().Axis().Direction()), myRadius);
00478   gp_Pnt2d aPoint2d1, aPoint2d2;
00479   myCurve1->D0(aNearest->GetParam(), aPoint2d1);
00480   myCurve2->D0(aNearest->GetParam2(), aPoint2d2);
00481   gp_Pnt aPoint1 = ElSLib::PlaneValue(aPoint2d1.X(), aPoint2d1.Y(), myPlane->Pln().Position());
00482   gp_Pnt aPoint2 = ElSLib::PlaneValue(aPoint2d2.X(), aPoint2d2.Y(), myPlane->Pln().Position());
00483 
00484   GeomAPI_ProjectPointOnCurve aProj(thePoint, aCircle);
00485   Standard_Real aTarGetParam = aProj.LowerDistanceParameter();
00486   gp_Pnt aPointOnCircle = aProj.NearestPoint();
00487 
00488   // Check extrema point manually, because there is a bug in Open CASCADE
00489   //  in calculation of nearest point to a circle near the parameter 0.0
00490   gp_Pnt p0 = ElCLib::Value(0.0, aCircle->Circ());
00491   if (p0.Distance(thePoint) < aPointOnCircle.Distance(thePoint))
00492   {
00493      aTarGetParam = 0.0;
00494      aPointOnCircle = p0;
00495   }
00496 
00497   aProj.Perform(aPoint1);
00498   Standard_Real aParam1 = aProj.LowerDistanceParameter();
00499     aProj.Perform(aPoint2);
00500   Standard_Real aParam2 = aProj.LowerDistanceParameter();
00501   Standard_Boolean aIsOut = ((aParam1 < aTarGetParam && aParam2 < aTarGetParam) ||
00502                              (aParam1 > aTarGetParam && aParam2 > aTarGetParam));
00503   if (aParam1 > aParam2)
00504     aIsOut = !aIsOut;
00505   BRepBuilderAPI_MakeEdge aBuilder(aCircle->Circ(),
00506                                    aIsOut ? aParam2 : aParam1,
00507                                    aIsOut? aParam1 : aParam2);
00508   aResult = aBuilder.Edge();
00509 
00510   // divide edges
00511   Standard_Real aStart, anEnd;
00512   Handle(Geom_Curve) aCurve = BRep_Tool::Curve(myEdge1, aStart, anEnd);
00513   gp_Vec aDir;
00514   aCurve->D1(aNearest->GetParam(), aPoint1, aDir);
00515 
00516   gp_Vec aCircleDir;
00517   aCircle->D1(aParam1, aPoint1, aCircleDir);
00518   if ((aCircleDir.Angle(aDir) > M_PI / 2.0) ^ aIsOut)
00519     aStart = aNearest->GetParam();
00520   else
00521     anEnd = aNearest->GetParam();
00522 
00523   if (fabs(aStart - anEnd) > Precision::Confusion())
00524   {
00525       //Divide edge
00526       BRepBuilderAPI_MakeEdge aDivider1(aCurve, aStart, anEnd);
00527       if (myEdgesExchnged)
00528         theEdge2 = aDivider1.Edge();
00529       else
00530         theEdge1 = aDivider1.Edge();
00531   }
00532 
00533   aCurve = BRep_Tool::Curve(myEdge2, aStart, anEnd);
00534   aCurve->D1(aNearest->GetParam2(), aPoint2, aDir);
00535 
00536   aCircle->D1(aParam2, aPoint2, aCircleDir);
00537   if ((aCircleDir.Angle(aDir) > M_PI / 2.0) ^ (!aIsOut))
00538     aStart = aNearest->GetParam2();
00539   else
00540     anEnd = aNearest->GetParam2();
00541 
00542   if (fabs(aStart - anEnd) > Precision::Confusion())
00543   {
00544       BRepBuilderAPI_MakeEdge aDivider2(aCurve, aStart, anEnd);
00545       if (myEdgesExchnged)
00546         theEdge1 = aDivider2.Edge();
00547       else
00548         theEdge2 = aDivider2.Edge();
00549   }
00550 
00551   delete aNearest;
00552   return aResult;
00553 }
00554 
00555 //=======================================================================
00556 //function : AddValue
00557 //purpose  :
00558 //=======================================================================
00559 void GEOMImpl_Fillet1dPoint::AddValue(Standard_Real theValue, Standard_Boolean theValid)
00560 {
00561   Standard_Integer a;
00562   for(a = 1; a <= myV.Length(); a++)
00563   {
00564     if (theValue < myV.Value(a))
00565     {
00566       myV.InsertBefore(a, theValue);
00567       myValid.InsertBefore(a, (Standard_Integer)theValid);
00568       return;
00569     }
00570   }
00571   myV.Append(theValue);
00572   myValid.Append((Standard_Integer)theValid);
00573 }
00574 
00575 //=======================================================================
00576 //function : ComputeDifference
00577 //purpose  :
00578 //=======================================================================
00579 Standard_Boolean GEOMImpl_Fillet1dPoint::ComputeDifference(GEOMImpl_Fillet1dPoint* thePoint)
00580 {
00581   Standard_Integer a;
00582   Standard_Boolean aDiffsSet = (myD.Length() != 0);
00583   Standard_Real aDX = thePoint->GetParam() - myParam, aDY;
00584   if (thePoint->myV.Length() == myV.Length())
00585   { // absolutely the same points
00586     for(a = 1; a <= myV.Length(); a++)
00587     {
00588       aDY = thePoint->myV.Value(a) - myV.Value(a);
00589       if ( aDiffsSet )
00590         myD.SetValue(a, fabs(aDX) > gp::Resolution() ? (aDY/aDX) : 0);
00591       else
00592         myD.Append( fabs(aDX) > gp::Resolution() ? (aDY/aDX) : 0);
00593     }
00594     return Standard_True;
00595   }
00596   // between the diffeerent points searching for nearest analogs
00597   Standard_Integer b;
00598   for(a = 1; a <= myV.Length(); a++)
00599   {
00600     for(b = 1; b <= thePoint->myV.Length(); b++)
00601     {
00602       if (b == 1 || fabs(thePoint->myV.Value(b) - myV.Value(a)) < fabs(aDY))
00603         aDY = thePoint->myV.Value(b) - myV.Value(a);
00604     }
00605     if (aDiffsSet)
00606     {
00607       if ( fabs(aDX) > gp::Resolution() && fabs(aDY / aDX) < fabs(myD.Value(a)))
00608         myD.SetValue(a, aDY / aDX);
00609       else
00610         myD.SetValue(a, 0);
00611     }
00612     else
00613     {
00614       myD.Append( fabs(aDX) > gp::Resolution() ? aDY/aDX : 0);
00615     }
00616   }
00617 
00618   return Standard_False;
00619 }
00620 
00621 //=======================================================================
00622 //function : FilterPoints
00623 //purpose  :
00624 //=======================================================================
00625 void GEOMImpl_Fillet1dPoint::FilterPoints(GEOMImpl_Fillet1dPoint* thePoint)
00626 {
00627   Standard_Integer a, b;
00628   TColStd_SequenceOfReal aDiffs;
00629   Standard_Real aY, aY2, aDX = thePoint->GetParam() - myParam;
00630   for(a = 1; a <= myV.Length(); a++)
00631   {
00632     // searching for near point from thePoint
00633     Standard_Integer aNear = 0;
00634     Standard_Real aDiff = aDX * 10000.;
00635     aY = myV.Value(a) + myD.Value(a) * aDX;
00636     for(b = 1; b <= thePoint->myV.Length(); b++)
00637     {
00638       // calculate hypothesis value of the Y2 with the constant first and second derivative
00639       aY2 = aY + aDX * (thePoint->myD.Value(b) - myD.Value(a)) / 2.0;
00640       if (aNear == 0 || fabs(aY2 - thePoint->myV.Value(b)) < fabs(aDiff))
00641       {
00642         aNear = b;
00643         aDiff = aY2 - thePoint->myV.Value(b);
00644       }
00645     }//for b...
00646 
00647     if (aNear)
00648     {
00649       if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
00650       {// the same sign at the same sides of the interval
00651         if (myV.Value(a) * myD.Value(a) > 0)
00652         {
00653           if (fabs(myD.Value(a)) > Precision::Confusion())
00654             aNear = 0;
00655         }
00656         else
00657         {
00658           if (fabs(myV.Value(a)) > fabs(thePoint->myV.Value(aNear)))
00659             if (thePoint->myV.Value(aNear) * thePoint->myD.Value(aNear) < 0 &&
00660                 fabs(thePoint->myD.Value(aNear)) > Precision::Confusion())
00661             {
00662               aNear = 0;
00663             }
00664         }
00665       }
00666     }
00667 
00668     if (aNear)
00669     {
00670       if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
00671       {
00672         if ((myV.Value(a) + myD.Value(a) * aDX) * myV.Value(a) > Precision::Confusion() &&
00673         (thePoint->myV.Value(aNear) + thePoint->myD.Value(aNear) * aDX) * thePoint->myV.Value(aNear) > Precision::Confusion())
00674         {
00675           aNear = 0;
00676         }
00677       }
00678     }
00679 
00680     if (aNear)
00681     {
00682       if (  fabs(aDX) < gp::Resolution() || fabs(aDiff / aDX) > 1.e+7)
00683       {
00684         aNear = 0;
00685       }
00686     }
00687 
00688     if (aNear == 0)
00689     {  // there is no near: remove it from the list
00690       myV.Remove(a);
00691       myD.Remove(a);
00692       myValid.Remove(a);
00693       a--;
00694     }
00695     else
00696     {
00697       Standard_Boolean aFound = Standard_False;
00698       for(b = 1; b <= myNear.Length(); b++)
00699       {
00700         if (myNear.Value(b) == aNear)
00701         {
00702           if (fabs(aDiffs.Value(b)) < fabs(aDiff))
00703           { // return this 'near'
00704             aFound = Standard_True;
00705             myV.Remove(a);
00706             myD.Remove(a);
00707             myValid.Remove(a);
00708             a--;
00709             break;
00710           }
00711           else
00712           { // remove the old 'near'
00713             myV.Remove(b);
00714             myD.Remove(b);
00715             myValid.Remove(b);
00716             myNear.Remove(b);
00717             aDiffs.Remove(b);
00718             a--;
00719             break;
00720           }
00721         }
00722       }//for b...
00723       if (!aFound)
00724       {
00725         myNear.Append(aNear);
00726         aDiffs.Append(aDiff);
00727       }
00728     }
00729   }//for a...
00730 }
00731 
00732 //=======================================================================
00733 //function : Copy
00734 //purpose  :
00735 //=======================================================================
00736 GEOMImpl_Fillet1dPoint* GEOMImpl_Fillet1dPoint::Copy()
00737 {
00738   GEOMImpl_Fillet1dPoint* aCopy = new GEOMImpl_Fillet1dPoint(myParam);
00739   Standard_Integer a;
00740   for(a = 1; a <= myV.Length(); a++)
00741   {
00742     aCopy->myV.Append(myV.Value(a));
00743     aCopy->myD.Append(myD.Value(a));
00744     aCopy->myValid.Append(myValid.Value(a));
00745   }
00746   return aCopy;
00747 }
00748 
00749 //=======================================================================
00750 //function : HasSolution
00751 //purpose  :
00752 //=======================================================================
00753 Standard_Integer GEOMImpl_Fillet1dPoint::HasSolution(const Standard_Real theRadius)
00754 {
00755   Standard_Integer a;
00756   for(a = 1; a <= myV.Length(); a++)
00757   {
00758     if (fabs(sqrt(fabs(fabs(myV.Value(a)) + theRadius * theRadius)) - theRadius) < Precision::Confusion() / 10.)
00759       return a;
00760   }
00761   return 0;
00762 }
00763 
00764 //=======================================================================
00765 //function : RemoveSolution
00766 //purpose  :
00767 //=======================================================================
00768 void GEOMImpl_Fillet1dPoint::RemoveSolution(Standard_Integer theIndex)
00769 {
00770   myV.Remove(theIndex);
00771   myD.Remove(theIndex);
00772   myValid.Remove(theIndex);
00773   myNear.Remove(theIndex);
00774 }