OpenCascade/occt
1
0
Fork 0
mirror of https://git.dev.opencascade.org/repos/occt.git synced 2025-04-10 18:51:21 +03:00
Code

0031661: Modeling Data - Exception when projecting parabola or hyperbola to plane

Browse Source 
ProjLib/ProjLib_ProjectOnPlane.cxx - formatting

0031661: Modeling Data - Algorithm crashes when projecting parabola or hyperbola to plane

ProjLib/ProjLib_ProjectOnPlane.cxx - building of analytical parabola and hyperbola is added
bugs/moddata_3/bug31661_* - new test cases are added
This commit is contained in:
ifv 2022-03-03 15:05:23 +03:00 committed by smoskvin
parent 1f000e5974
commit 1f37f1d50a
4 changed files with 1020 additions and 682 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

395
src/ProjLib/ProjLib_ProjectOnPlane.cxx
View File

@ -47,8 +47,14 @@
#include <GeomLib_Tool.hxx>
#include <math_Jacobi.hxx>
#include <math_Matrix.hxx>
#include <gce_MakeParab.hxx>
#include <gce_MakeDir.hxx>
#include <LProp3d_CLProps.hxx>
#include <math_Function.hxx>
#include <math_BrentMinimum.hxx>
const Standard_Real aParabolaLimit = 20000.;
const Standard_Real aHyperbolaLimit = 10.;
//=======================================================================
//function : OnPlane_Value
@ -269,10 +275,14 @@ public :
}
Standard_Real FirstParameter() const
{return myCurve->FirstParameter();}
{
return myCurve->FirstParameter();
}
Standard_Real LastParameter() const
{return myCurve->LastParameter();}
{
return myCurve->LastParameter();
}
Standard_Boolean Value(const Standard_Real theT,
NCollection_Array1<gp_Pnt2d>& /*thePnt2d*/,
@ -292,6 +302,35 @@ public :
};
//=======================================================================
// class : ProjLib_MaxCurvature
//purpose : Use to search apex of parabola or hyperbola, which is its projection
// on a plane. Apex is point with maximal curvature
//=======================================================================
class ProjLib_MaxCurvature : public math_Function
{
public:
ProjLib_MaxCurvature(LProp3d_CLProps& theProps):
myProps(&theProps)
{
}
virtual Standard_Boolean Value(const Standard_Real X, Standard_Real& F)
{
myProps->SetParameter(X);
F = -myProps->Curvature();
return Standard_True;
}
private:
LProp3d_CLProps* myProps;
};
//=====================================================================//
@ -319,10 +358,20 @@ static void PerformApprox (const Handle(Adaptor3d_Curve)& C,
Standard_Integer Deg1, Deg2;
Deg1 = 8; Deg2 = 8;
if (C->GetType() == GeomAbs_Parabola)
{
Deg1 = 2; Deg2 = 2;
}
Standard_Integer aNbSegm = 100;
if (C->GetType() == GeomAbs_Hyperbola)
{
Deg1 = 14;
Deg2 = 14;
aNbSegm = 1000;
}
Approx_FitAndDivide Fit(Deg1, Deg2, Precision::Approximation(),
Precision::PApproximation(), Standard_True);
Fit.SetMaxSegments(100);
Fit.SetMaxSegments(aNbSegm);
Fit.Perform(F);
if (!Fit.IsAllApproximated())
{
@ -348,9 +397,11 @@ static void PerformApprox (const Handle(Adaptor3d_Curve)& C,
TColStd_Array1OfReal Knots(1, NbCurves + 1); //Noeuds de la BSpline
Standard_Integer Compt = 1;
Standard_Real anErrMax = 0., anErr3d, anErr2d;
for (i = 1; i <= Fit.NbMultiCurves(); i++) {
Fit.Parameters(i, Knots(i), Knots(i + 1));
Fit.Error(i, anErr3d, anErr2d);
anErrMax = Max(anErrMax, anErr3d);
AppParCurves_MultiCurve MC = Fit.Value(i); //Charge la Ieme Curve
TColgp_Array1OfPnt LocalPoles(1, MC.Degree() + 1);//Recupere les poles
MC.Curve(1, LocalPoles);
@ -389,6 +440,17 @@ static void PerformApprox (const Handle(Adaptor3d_Curve)& C,
Mults.SetValue(NbKnots, MaxDeg + 1);
BSplineCurvePtr =
new Geom_BSplineCurve(Poles, Knots, Mults, MaxDeg, Standard_False);
//Try to smooth
Standard_Integer m1 = MaxDeg - 1;
for (i = 2; i < NbKnots; ++i)
{
if (BSplineCurvePtr->Multiplicity(i) == MaxDeg)
{
BSplineCurvePtr->RemoveKnot(i, m1, anErrMax);
}
}
}
@ -537,6 +599,7 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_Curve)& C,
Handle(Geom_Circle) GeomCirclePtr;
Handle(Geom_Ellipse) GeomEllipsePtr;
Handle(Geom_Hyperbola) GeomHyperbolaPtr;
Handle(Geom_Parabola) GeomParabolaPtr;
gp_Lin aLine;
gp_Elips Elips;
@ -806,7 +869,6 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_Curve)& C,
break;
case GeomAbs_Parabola:
{
myKeepParam = Standard_True;
// P(u) = O + (u*u)/(4*f) * Xc + u * Yc
// ==> Q(u) = f(P(u))
// = f(O) + (u*u)/(4*f) * f(Xc) + u * f(Yc)
@ -815,50 +877,51 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_Curve)& C,
gp_Ax2 AxeRef = Parab.Position();
gp_Vec Xc = ProjectVec(myPlane, myDirection, gp_Vec(AxeRef.XDirection()));
gp_Vec Yc = ProjectVec(myPlane, myDirection, gp_Vec(AxeRef.YDirection()));
// fix for case when no one action is done. 28.03.2002
Standard_Boolean alocalIsDone = Standard_False;
if ( Abs( Yc.Magnitude() - 1.) < Precision::Confusion()) {
gp_Pnt P = ProjectPnt(myPlane, myDirection, AxeRef.Location());
if ( Xc.Magnitude() < Precision::Confusion()) {
myIsApprox = Standard_False;
if ((Abs(Yc.Magnitude() - 1.) < Precision::Confusion()) &&
(Xc.Magnitude() < Precision::Confusion()))
{
myType = GeomAbs_Line;
aLine = gp_Lin(P, gp_Dir(Yc));
GeomLinePtr = new Geom_Line(aLine);
// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
GeomAdaptor_Curve aGACurve(GeomLinePtr);
myResult = new GeomAdaptor_Curve(aGACurve);
// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
alocalIsDone = Standard_True;
}
else if (Xc.IsNormal(Yc, Precision::Angular())) {
myType = GeomAbs_Parabola;
Standard_Real F = Parab.Focal() / Xc.Magnitude();
gp_Parab aParab = gp_Parab( gp_Ax2(P,Xc^Yc,Xc), F);
Handle(Geom_Parabola) GeomParabolaPtr =
new Geom_Parabola(aParab) ;
// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
GeomAdaptor_Curve aGACurve(GeomParabolaPtr);
myResult = new GeomAdaptor_Curve(aGACurve);
// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
alocalIsDone = Standard_True;
gp_Parab aProjParab = gp_Parab(gp_Ax2(P, Xc^Yc, Xc), F);
GeomParabolaPtr =
new Geom_Parabola(aProjParab);
}
}
if (!alocalIsDone)/*else*/ {
else if (Yc.Magnitude() < Precision::Confusion() ||
Yc.IsParallel(Xc, Precision::Angular()))
{
myIsApprox = Standard_True;
myType = GeomAbs_BSplineCurve;
PerformApprox(myCurve,myPlane,myDirection,ApproxCurve);
// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
GeomAdaptor_Curve aGACurve(ApproxCurve);
myResult = new GeomAdaptor_Curve(aGACurve);
// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
}
else if(!myKeepParam)
{
// Try building parabola with help of apex position
myIsApprox = !BuildParabolaByApex(GeomParabolaPtr);
}
else
{
myIsApprox = Standard_True;
}
if (!myIsApprox)
{
GetTrimmedResult(GeomParabolaPtr);
}
else
{
BuildByApprox(aParabolaLimit);
}
}
break;
case GeomAbs_Hyperbola:
{
myKeepParam = Standard_True;
// P(u) = O + R1 * Cosh(u) * Xc + R2 * Sinh(u) * Yc
// ==> Q(u) = f(P(u))
// = f(O) + R1 * Cosh(u) * f(Xc) + R2 * Sinh(u) * f(Yc)
@ -871,6 +934,7 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_Curve)& C,
Standard_Real aR1 = Hypr.MajorRadius();
Standard_Real aR2 = Hypr.MinorRadius();
gp_Dir Z = myPlane.Direction();
myIsApprox = Standard_False;
if (Xc.Magnitude() < Precision::Confusion()) {
myType = GeomAbs_Hyperbola;
@ -878,10 +942,6 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_Curve)& C,
Hypr = gp_Hypr(gp_Ax2(P, Z, X), 0., aR2 * Yc.Magnitude());
GeomHyperbolaPtr =
new Geom_Hyperbola(Hypr);
// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
GeomAdaptor_Curve aGACurve(GeomHyperbolaPtr);
myResult = new GeomAdaptor_Curve(aGACurve);
// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
}
else if (Yc.Magnitude() < Precision::Confusion()) {
myType = GeomAbs_Hyperbola;
@ -889,10 +949,6 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_Curve)& C,
gp_Hypr(gp_Ax2(P, Z, gp_Dir(Xc)), aR1 * Xc.Magnitude(), 0.);
GeomHyperbolaPtr =
new Geom_Hyperbola(Hypr);
// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
GeomAdaptor_Curve aGACurve(GeomHyperbolaPtr);
myResult = new GeomAdaptor_Curve(aGACurve);
// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
}
else if (Xc.IsNormal(Yc, Precision::Angular())) {
myType = GeomAbs_Hyperbola;
@ -900,19 +956,27 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_Curve)& C,
aR1 * Xc.Magnitude(), aR2 * Yc.Magnitude());
GeomHyperbolaPtr =
new Geom_Hyperbola(Hypr);
// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
GeomAdaptor_Curve aGACurve(GeomHyperbolaPtr);
myResult = new GeomAdaptor_Curve(aGACurve);
// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
}
else {
else if (Yc.Magnitude() < Precision::Confusion() ||
Yc.IsParallel(Xc, Precision::Angular()))
{
myIsApprox = Standard_True;
myType = GeomAbs_BSplineCurve;
PerformApprox(myCurve,myPlane,myDirection,ApproxCurve);
// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
GeomAdaptor_Curve aGACurve(ApproxCurve);
myResult = new GeomAdaptor_Curve(aGACurve);
// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
}
else if(!myKeepParam)
{
myIsApprox = !BuildHyperbolaByApex(GeomHyperbolaPtr);
}
else
{
myIsApprox = Standard_True;
}
if ( !myIsApprox )
{
GetTrimmedResult(GeomHyperbolaPtr);
}
else
{
BuildByApprox(aHyperbolaLimit);
}
}
break;
@ -1458,7 +1522,7 @@ Handle(Geom_BezierCurve) ProjLib_ProjectOnPlane::Bezier() const
}
//=======================================================================
//function : Bezier
//function : BSpline
//purpose :
//=======================================================================
@ -1470,3 +1534,222 @@ Handle(Geom_BSplineCurve) ProjLib_ProjectOnPlane::BSpline() const
return myResult->BSpline();
}
//=======================================================================
//function : GetTrimmedResult
//purpose :
//=======================================================================
void ProjLib_ProjectOnPlane::GetTrimmedResult(const Handle(Geom_Curve)& theProjCurve)
{
gp_Lin aLin;
gp_Parab aParab;
gp_Hypr aHypr;
if (myType == GeomAbs_Line)
{
aLin = Handle(Geom_Line)::DownCast(theProjCurve)->Lin();
}
else if (myType == GeomAbs_Parabola)
{
aParab = Handle(Geom_Parabola)::DownCast(theProjCurve)->Parab();
}
else if (myType == GeomAbs_Hyperbola)
{
aHypr = Handle(Geom_Hyperbola)::DownCast(theProjCurve)->Hypr();
}
myFirstPar = theProjCurve->FirstParameter();
myLastPar = theProjCurve->LastParameter();
if (!Precision::IsInfinite(myCurve->FirstParameter()))
{
gp_Pnt aP = myCurve->Value(myCurve->FirstParameter());
aP = ProjectPnt(myPlane, myDirection, aP);
if (myType == GeomAbs_Line)
{
myFirstPar = ElCLib::Parameter(aLin, aP);
}
else if (myType == GeomAbs_Parabola)
{
myFirstPar = ElCLib::Parameter(aParab, aP);
}
else if (myType == GeomAbs_Hyperbola)
{
myFirstPar = ElCLib::Parameter(aHypr, aP);
}
else
{
GeomLib_Tool::Parameter(theProjCurve, aP, Precision::Confusion(), myFirstPar);
}
}
if (!Precision::IsInfinite(myCurve->LastParameter()))
{
gp_Pnt aP = myCurve->Value(myCurve->LastParameter());
aP = ProjectPnt(myPlane, myDirection, aP);
if (myType == GeomAbs_Line)
{
myLastPar = ElCLib::Parameter(aLin, aP);
}
else if (myType == GeomAbs_Parabola)
{
myLastPar = ElCLib::Parameter(aParab, aP);
}
else if (myType == GeomAbs_Hyperbola)
{
myLastPar = ElCLib::Parameter(aHypr, aP);
}
else
{
GeomLib_Tool::Parameter(theProjCurve, aP, Precision::Confusion(), myLastPar);
}
}
myResult = new GeomAdaptor_Curve(theProjCurve, myFirstPar, myLastPar);
}
//=======================================================================
//function : BuildParabolaByApex
//purpose :
//=======================================================================
Standard_Boolean ProjLib_ProjectOnPlane::BuildParabolaByApex(Handle(Geom_Curve)& theGeomParabolaPtr)
{
//
//Searching parabola apex as point with maximal curvature
Standard_Real aF = myCurve->Parabola().Focal();
GeomAbs_CurveType aCurType = myType;
myType = GeomAbs_OtherCurve; //To provide correct calculation of derivativesb by projection for
//copy of instance;
Handle(Adaptor3d_Curve) aProjCrv = ShallowCopy();
myType = aCurType;
LProp3d_CLProps aProps(aProjCrv, 2, Precision::Confusion());
ProjLib_MaxCurvature aMaxCur(aProps);
math_BrentMinimum aSolver(Precision::PConfusion());
aSolver.Perform(aMaxCur, -10.*aF, 0., 10.*aF);
if (!aSolver.IsDone())
{
return Standard_False;
}
Standard_Real aT;
aT = aSolver.Location();
aProps.SetParameter(aT);
gp_Pnt aP0 = aProps.Value();
gp_Vec aDY = aProps.D1();
gp_Dir anYDir(aDY);
gp_Dir anXDir;
Standard_Real aCurv = aProps.Curvature();
if (Precision::IsInfinite(aCurv) || aCurv < Precision::Confusion())
{
return Standard_False;
}
aProps.Normal(anXDir);
//
gp_Lin anXLine(aP0, anXDir);
gp_Pnt aP1 = Value(aT + 10.*aF);
//
Standard_Real anX = ElCLib::LineParameter(anXLine.Position(), aP1);
Standard_Real anY = anXLine.Distance(aP1);
Standard_Real aNewF = anY * anY / 4. / anX;
gp_Dir anN = anXDir^anYDir;
gp_Ax2 anA2(aP0, anN, anXDir);
gce_MakeParab aMkParab(anA2, aNewF);
if (!aMkParab.IsDone())
{
return Standard_False;
}
gp_Parab aProjParab = aMkParab.Value();
myType = GeomAbs_Parabola;
theGeomParabolaPtr = new Geom_Parabola(aProjParab);
//GetTrimmedResult(theGeomParabolaPtr);
return Standard_True;
}
//=======================================================================
//function : BuildHyperbolaByApex
//purpose :
//=======================================================================
Standard_Boolean ProjLib_ProjectOnPlane::BuildHyperbolaByApex(Handle(Geom_Curve)& theGeomHyperbolaPtr)
{
//Try to build hyperbola with help of apex position
GeomAbs_CurveType aCurType = myType;
myType = GeomAbs_OtherCurve; //To provide correct calculation of derivativesb by projection for
//copy of instance;
Handle(Adaptor3d_Curve) aProjCrv = ShallowCopy();
myType = aCurType;
//Searching hyperbola apex as point with maximal curvature
LProp3d_CLProps aProps(aProjCrv, 2, Precision::Confusion());
ProjLib_MaxCurvature aMaxCur(aProps);
math_BrentMinimum aSolver(Precision::PConfusion());
aSolver.Perform(aMaxCur, -5., 0., 5.);
if (aSolver.IsDone())
{
Standard_Real aT;
aT = aSolver.Location();
aProps.SetParameter(aT);
Standard_Real aCurv = aProps.Curvature();
if (Precision::IsInfinite(aCurv) || aCurv < Precision::Confusion())
{
return Standard_False;
}
else
{
gp_Hypr Hypr = myCurve->Hyperbola();
gp_Ax2 AxeRef = Hypr.Position();
gp_Pnt P = ProjectPnt(myPlane, myDirection, AxeRef.Location());
gp_Dir Z = myPlane.Direction();
gp_Pnt aP0 = aProps.Value();
gp_Dir anXDir = gce_MakeDir(P, aP0);
gp_Dir anYDir = gce_MakeDir(aProps.D1());
//
Standard_Real aMajRad = P.Distance(aP0);
gp_Pnt aP1 = Value(aT + 1.);
gp_Vec aV(P, aP1);
Standard_Real anX = aV * anXDir;
Standard_Real anY = aV * anYDir;
Standard_Real aMinRad = anY / Sqrt(anX * anX / aMajRad / aMajRad - 1.);
gp_Ax2 anA2(P, Z, anXDir);
gp_Hypr anHypr(anA2, aMajRad, aMinRad);
theGeomHyperbolaPtr =
new Geom_Hyperbola(anHypr);
myType = GeomAbs_Hyperbola;
}
}
else
{
return Standard_False;
}
return Standard_True;
}
//=======================================================================
//function : BuilByApprox
//purpose :
//=======================================================================
void ProjLib_ProjectOnPlane::BuildByApprox(const Standard_Real theLimitParameter)
{
myType = GeomAbs_BSplineCurve;
Handle(Geom_BSplineCurve) anApproxCurve;
if (Precision::IsInfinite(myCurve->FirstParameter()) ||
Precision::IsInfinite(myCurve->LastParameter()))
{
//To avoid exception in approximation
Standard_Real f = Max(-theLimitParameter, myCurve->FirstParameter());
Standard_Real l = Min(theLimitParameter, myCurve->LastParameter());
Handle(Adaptor3d_Curve) aTrimCurve = myCurve->Trim(f, l, Precision::Confusion());
PerformApprox(aTrimCurve, myPlane, myDirection, anApproxCurve);
}
else
{
PerformApprox(myCurve, myPlane, myDirection, anApproxCurve);
}
myFirstPar = anApproxCurve->FirstParameter();
myLastPar = anApproxCurve->LastParameter();
GeomAdaptor_Curve aGACurve(anApproxCurve);
myResult = new GeomAdaptor_Curve(aGACurve);
}

5
src/ProjLib/ProjLib_ProjectOnPlane.hxx
View File

@ -192,9 +192,12 @@ public:
protected:
void GetTrimmedResult(const Handle(Geom_Curve)& theProjCurve);
Standard_Boolean BuildParabolaByApex(Handle(Geom_Curve)& theGeomParabolaPtr);
Standard_Boolean BuildHyperbolaByApex(Handle(Geom_Curve)& theGeomParabolaPtr);
void BuildByApprox(const Standard_Real theLimitParameter);
private:

26
tests/bugs/moddata_3/bug31661_1 Normal file
View File

@ -0,0 +1,26 @@
puts ""
puts "=========================================================================="
puts "OCC31661: Modeling Data - Algorithm crashes when projecting parabola or hyperbola to plane"
puts "=========================================================================="
puts ""
parabola p 0 0 0 1 1 1 2 0 -2 10
plane pln 0 0 0 0 0 1
projonplane r p pln 0
if {![regexp {Parabola} [dump r]]} {
puts "ERROR: Projected curve is not a parabola"
}
trim p p -100 100
projonplane rp p pln 0
if {![regexp {Parabola} [dump rp]]} {
puts "ERROR: Projected curve is not a parabola"
}
checklength rp -l 408.40363195229503
bounds rp t1 t2
checkreal t1 [dval t1] -91.077748943768597 1.e-7 1.e-7
checkreal t2 [dval t2] 72.221567418462357 1.e-7 1.e-7

26
tests/bugs/moddata_3/bug31661_2 Normal file
View File

@ -0,0 +1,26 @@
puts ""
puts "=========================================================================="
puts "OCC31661: Modeling Data - Algorithm crashes when projecting parabola or hyperbola to plane"
puts "=========================================================================="
puts ""
hyperbola h 0 0 0 1 1 1 2 0 -2 10 10
plane pln 0 0 0 0 0 1
projonplane rh h pln 0
if {![regexp {Hyperbola} [dump rh]]} {
puts "ERROR: Projected curve is not a hyperbola"
}
trim h h -5 5
projonplane rh h pln 0
if {![regexp {Hyperbola} [dump rh]]} {
puts "ERROR: Projected curve is not a hyperbola"
}
checklength rh -l 1664.3732976598988
bounds rh t1 t2
checkreal t1 [dval t1] -5.23179933356147 1.e-7 1.e-7
checkreal t2 [dval t2] 4.76820064934972 1.e-7 1.e-7