1
0
mirror of https://git.dev.opencascade.org/repos/occt.git synced 2025-04-04 18:06:22 +03:00

0026525: Wrong result obtained by curve / surface intersection algorithm.

Added protection from double overflow caused by untrimmed parameters space.
Obsolete code deleted.

Draw command OCC26525 added.

Test cases for issue CR26525
This commit is contained in:
aml 2015-09-24 13:53:09 +03:00 committed by kgv
parent 657ccd95ab
commit 81b47143f1
6 changed files with 351 additions and 183 deletions

View File

@ -306,204 +306,65 @@ void IntCurveSurface_Inter::DoNewBounds(
//purpose : Decompose la surface si besoin est
//=======================================================================
void IntCurveSurface_Inter::Perform(const TheCurve& curve,
const TheSurface& surface) {
const TheSurface& surface)
{
ResetFields();
done = Standard_True;
Standard_Integer NbUOnS = TheSurfaceTool::NbUIntervals(surface,GeomAbs_C2);
Standard_Integer NbVOnS = TheSurfaceTool::NbVIntervals(surface,GeomAbs_C2);
Standard_Real U0,U1,V0,V1;
if(NbUOnS > 1) {
if(NbUOnS > 1)
{
TColStd_Array1OfReal TabU(1,NbUOnS+1);
TheSurfaceTool::UIntervals(surface,TabU,GeomAbs_C2);
for(Standard_Integer iu = 1;iu <= NbUOnS; iu++) {
for(Standard_Integer iu = 1;iu <= NbUOnS; iu++)
{
U0 = TabU.Value(iu);
U1 = TabU.Value(iu+1);
if(NbVOnS > 1) {
TColStd_Array1OfReal TabV(1,NbVOnS+1);
TheSurfaceTool::VIntervals(surface,TabV,GeomAbs_C2);
for(Standard_Integer iv = 1;iv <= NbVOnS; iv++) {
V0 = TabV.Value(iv);
V1 = TabV.Value(iv+1);
Perform(curve,surface,U0,V0,U1,V1);
}
if(NbVOnS > 1)
{
TColStd_Array1OfReal TabV(1,NbVOnS+1);
TheSurfaceTool::VIntervals(surface,TabV,GeomAbs_C2);
for(Standard_Integer iv = 1;iv <= NbVOnS; iv++)
{
// More than one interval on U and V param space.
V0 = TabV.Value(iv);
V1 = TabV.Value(iv+1);
Perform(curve,surface,U0,V0,U1,V1);
}
}
else { //------ composite en U non composite en V
V0 = TheSurfaceTool::FirstVParameter(surface);
V1 = TheSurfaceTool::LastVParameter(surface);
Perform(curve,surface,U0,V0,U1,V1);
else
{
// More than one interval only on U param space.
V0 = TheSurfaceTool::FirstVParameter(surface);
V1 = TheSurfaceTool::LastVParameter(surface);
Perform(curve,surface,U0,V0,U1,V1);
}
}
}
else if(NbVOnS > 1) { //---------- non composite en U composite en V
else if(NbVOnS > 1)
{
// More than one interval only on V param space.
U0 = TheSurfaceTool::FirstUParameter(surface);
U1 = TheSurfaceTool::LastUParameter(surface);
TColStd_Array1OfReal TabV(1,NbVOnS+1);
TheSurfaceTool::VIntervals(surface,TabV,GeomAbs_C2);
for(Standard_Integer iv = 1;iv <= NbVOnS; iv++) {
for(Standard_Integer iv = 1;iv <= NbVOnS; iv++)
{
V0 = TabV.Value(iv);
V1 = TabV.Value(iv+1);
Perform(curve,surface,U0,V0,U1,V1);
}
}
else {
}
else
{
// One interval on U and V param space.
V0 = TheSurfaceTool::FirstVParameter(surface);
V1 = TheSurfaceTool::LastVParameter(surface);
U0 = TheSurfaceTool::FirstUParameter(surface);
U1 = TheSurfaceTool::LastUParameter(surface);
U1 = TheSurfaceTool::LastUParameter(surface);
#if 0
//-- jgv patch (from)
//Computing of local bounds
Standard_Real LocalU0 = U0, LocalU1 = U1, LocalV0 = V0, LocalV1 = V1;
Standard_Real Umin = RealLast(), Vmin = RealLast(), Umax = RealFirst(), Vmax = RealFirst();
Bnd_Box CurveBox;
Standard_Integer i, j, k;
//Making GeomAdaptor_Curve
Standard_Real f = TheCurveTool::FirstParameter(curve);
Standard_Real l = TheCurveTool::LastParameter(curve);
GeomAbs_CurveType CurveType = TheCurveTool::GetType(curve);
Handle(Geom_Curve) theCurve;
switch (CurveType)
{
case GeomAbs_Line:
theCurve = new Geom_Line( TheCurveTool::Line(curve) );
break;
case GeomAbs_Circle:
theCurve = new Geom_Circle( TheCurveTool::Circle(curve) );
break;
case GeomAbs_Ellipse:
theCurve = new Geom_Ellipse( TheCurveTool::Ellipse(curve) );
break;
case GeomAbs_Hyperbola:
theCurve = new Geom_Hyperbola( TheCurveTool::Hyperbola(curve) );
break;
case GeomAbs_Parabola:
theCurve = new Geom_Parabola( TheCurveTool::Parabola(curve) );
break;
case GeomAbs_BezierCurve:
theCurve = TheCurveTool::Bezier(curve);
break;
case GeomAbs_BSplineCurve:
theCurve = TheCurveTool::BSpline(curve);
break;
}
if (!theCurve.IsNull())
{
GeomAdaptor_Curve GACurve( theCurve, f, l );
BndLib_Add3dCurve::Add( GACurve, Precision::Confusion(), CurveBox );
}
else
{
Standard_Integer nbp = TheCurveTool::NbSamples(curve,f,l);
Standard_Real delta = (f-l)/nbp;
for (i = 0; i <= nbp; i++)
{
gp_Pnt aPnt = TheCurveTool::Value(curve, f + i*delta);
CurveBox.Add(aPnt);
}
}
Standard_Real X[2], Y[2], Z[2];
CurveBox.Get( X[0], Y[0], Z[0], X[1], Y[1], Z[1] );
Standard_Real Ures = TheSurfaceTool::UResolution(surface, Precision::Confusion());
Standard_Real Vres = TheSurfaceTool::VResolution(surface, Precision::Confusion());
//Making GeomAdaptor_Surface
GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
Handle(Geom_Surface) theSurface;
switch (SurfaceType)
{
case GeomAbs_Plane:
theSurface = new Geom_Plane( TheSurfaceTool::Plane(surface) );
break;
case GeomAbs_Cylinder:
theSurface = new Geom_CylindricalSurface( TheSurfaceTool::Cylinder(surface) );
break;
case GeomAbs_Cone:
theSurface = new Geom_ConicalSurface( TheSurfaceTool::Cone(surface) );
break;
case GeomAbs_Torus:
theSurface = new Geom_ToroidalSurface( TheSurfaceTool::Torus(surface) );
break;
case GeomAbs_Sphere:
theSurface = new Geom_SphericalSurface( TheSurfaceTool::Sphere(surface) );
break;
case GeomAbs_BezierSurface:
theSurface = TheSurfaceTool::Bezier(surface);
break;
case GeomAbs_BSplineSurface:
theSurface = TheSurfaceTool::BSpline(surface);
break;
case GeomAbs_SurfaceOfRevolution:
{
gp_Ax1 Axis = TheSurfaceTool::AxeOfRevolution(surface);
Handle(Adaptor3d_HCurve) AdBC = TheSurfaceTool::BasisCurve(surface);
Handle(Geom_Curve) BC = GetCurve(AdBC);
if (!BC.IsNull())
theSurface = new Geom_SurfaceOfRevolution( BC, Axis );
break;
}
case GeomAbs_SurfaceOfExtrusion:
{
gp_Dir Direction = TheSurfaceTool::Direction(surface);
Handle(Adaptor3d_HCurve) AdBC = TheSurfaceTool::BasisCurve(surface);
Handle(Geom_Curve) BC = GetCurve(AdBC);
if (!BC.IsNull())
theSurface = new Geom_SurfaceOfLinearExtrusion( BC, Direction );
break;
}
case GeomAbs_OffsetSurface:
{
Standard_Real OffsetValue = TheSurfaceTool::OffsetValue(surface);
Handle(Adaptor3d_HSurface) AdBS = TheSurfaceTool::BasisSurface(surface);
Handle(Geom_Surface) BS = GetSurface(AdBS);
if (!BS.IsNull())
theSurface = new Geom_OffsetSurface( BS, OffsetValue );
break;
}
}
if (!theSurface.IsNull())
{
GeomAdaptor_Surface GASurface( theSurface );
Extrema_ExtPS Projector;
Projector.Initialize( GASurface, U0, U1, V0, V1, Ures, Vres );
for (i = 0; i <= 1; i++)
for (j = 0; j <= 1; j++)
for (k = 0; k <= 1; k++)
{
gp_Pnt aPoint( X[i], Y[j], Z[k] );
Projector.Perform( aPoint );
if (Projector.IsDone() && Projector.NbExt() > 0)
{
Standard_Real mindist = RealLast();
Standard_Integer minind, ind;
for (ind = 1; ind <= Projector.NbExt(); ind++)
{
Standard_Real dist = Projector.Value(ind);
if (dist < mindist)
{
mindist = dist;
minind = ind;
}
}
Extrema_POnSurf pons = Projector.Point(minind);
Standard_Real U, V;
pons.Parameter( U, V );
if (U < Umin) Umin = U;
if (U > Umax) Umax = U;
if (V < Vmin) Vmin = V;
if (V > Vmax) Vmax = V;
}
}
Umin -= Ures; Umax += Ures; Vmin -= Vres; Vmax += Vres;
if (Umin > U0 && Umin <= U1) LocalU0 = Umin;
if (Umax < U1 && Umax >= U0) LocalU1 = Umax;
if (Vmin > V0 && Vmin <= V1) LocalV0 = Vmin;
if (Vmax < V1 && Vmax >= V0) LocalV1 = Vmax;
U0 = LocalU0; U1 = LocalU1; V0 = LocalV0; V1 = LocalV1;
}
//-- jgv patch (to)
#endif
Perform(curve,surface,U0,V0,U1,V1);
}
}
@ -512,36 +373,49 @@ void IntCurveSurface_Inter::Perform(const TheCurve& curve,
//purpose :
//=======================================================================
void IntCurveSurface_Inter::Perform(const TheCurve& curve,
const TheSurface& surface,
const Standard_Real U1,const Standard_Real V1,
const Standard_Real U2,const Standard_Real V2) {
const TheSurface& surface,
const Standard_Real U1,const Standard_Real V1,
const Standard_Real U2,const Standard_Real V2)
{
// Protection from double type overflow.
// This may happen inside square magnitude computation based on normal,
// which was computed on bound parameteres (bug26525).
Standard_Real UU1 = U1, UU2 = U2, VV1 = V1, VV2 = V2;
if (U1 < -1.0e50)
UU1 = -1.0e50;
if (U2 > 1.0e50)
UU2 = 1.0e50;
if (V1 < -1.0e50)
VV1 = -1.0e50;
if (V2 > 1.0e50)
VV2 = 1.0e50;
GeomAbs_CurveType CurveType = TheCurveTool::GetType(curve);
switch(CurveType) {
case GeomAbs_Line:
{
PerformConicSurf(TheCurveTool::Line(curve),curve,surface,U1,V1,U2,V2);
PerformConicSurf(TheCurveTool::Line(curve),curve,surface,UU1,VV1,UU2,VV2);
break;
}
case GeomAbs_Circle:
{
PerformConicSurf(TheCurveTool::Circle(curve),curve,surface,U1,V1,U2,V2);
PerformConicSurf(TheCurveTool::Circle(curve),curve,surface,UU1,VV1,UU2,VV2);
break;
}
case GeomAbs_Ellipse:
{
PerformConicSurf(TheCurveTool::Ellipse(curve),curve,surface,U1,V1,U2,V2);
PerformConicSurf(TheCurveTool::Ellipse(curve),curve,surface,UU1,VV1,UU2,VV2);
break;
}
case GeomAbs_Parabola:
{
PerformConicSurf(TheCurveTool::Parabola(curve),curve,surface,U1,V1,U2,V2);
PerformConicSurf(TheCurveTool::Parabola(curve),curve,surface,UU1,VV1,UU2,VV2);
break;
}
case GeomAbs_Hyperbola:
{
PerformConicSurf(TheCurveTool::Hyperbola(curve),curve,surface,U1,V1,U2,V2);
PerformConicSurf(TheCurveTool::Hyperbola(curve),curve,surface,UU1,VV1,UU2,VV2);
break;
}
default:
@ -568,7 +442,7 @@ void IntCurveSurface_Inter::Perform(const TheCurve& curve,
// IntCurveSurface_ThePolygon polygon(curve,u1,u2,TheCurveTool::NbSamples(curve,u1,u2));
IntCurveSurface_ThePolygon polygon(curve, aPars->Array1());
InternalPerform(curve,polygon,surface,U1,V1,U2,V2);
InternalPerform(curve,polygon,surface,UU1,VV1,UU2,VV2);
}
}
else {
@ -583,7 +457,7 @@ void IntCurveSurface_Inter::Perform(const TheCurve& curve,
// IntCurveSurface_ThePolygon polygon(curve,TheCurveTool::NbSamples(curve,u1,u2));
IntCurveSurface_ThePolygon polygon(curve, aPars->Array1());
InternalPerform(curve,polygon,surface,U1,V1,U2,V2);
InternalPerform(curve,polygon,surface,UU1,VV1,UU2,VV2);
}
}
else { //-- la surface est une quadrique

View File

@ -4133,6 +4133,90 @@ static Standard_Integer OCC26313(Draw_Interpretor& di,Standard_Integer n,const c
return 0;
}
//=======================================================================
//function : OCC26525
//purpose : check number of intersection points
//=======================================================================
#include <BRepAdaptor_Curve.hxx>
#include <BRepAdaptor_Surface.hxx>
#include <BRepAdaptor_HCurve.hxx>
#include <BRepAdaptor_HSurface.hxx>
#include <IntCurveSurface_HInter.hxx>
Standard_Integer OCC26525 (Draw_Interpretor& di,
Standard_Integer n,
const char** a)
{
TopoDS_Shape aS1, aS2;
TopoDS_Edge aE;
TopoDS_Face aF;
if (n<4)
{
di << " use OCC26525 r edge face \n";
return 1;
}
aS1 = DBRep::Get(a[2]);
aS2 = DBRep::Get(a[3]);
if (aS1.IsNull() || aS2.IsNull())
{
di << " Null shapes are not allowed \n";
return 0;
}
if (aS1.ShapeType()!=TopAbs_EDGE)
{
di << " Shape" << a[2] << " should be of type EDGE\n";
return 0;
}
if (aS2.ShapeType()!=TopAbs_FACE)
{
di << " Shape" << a[3] << " should be of type FACE\n";
return 0;
}
aE=TopoDS::Edge(aS1);
aF=TopoDS::Face(aS2);
char buf[128];
Standard_Boolean bIsDone;
Standard_Integer i, aNbPoints;
Standard_Real aU, aV, aT;
gp_Pnt aP;
BRepAdaptor_Curve aBAC;
BRepAdaptor_Surface aBAS;
IntCurveSurface_TransitionOnCurve aTC;
IntCurveSurface_HInter aHInter;
aBAC.Initialize(aE);
aBAS.Initialize(aF);
Handle(BRepAdaptor_HCurve) aHBAC=new BRepAdaptor_HCurve(aBAC);
Handle(BRepAdaptor_HSurface) aHBAS = new BRepAdaptor_HSurface(aBAS);
aHInter.Perform(aHBAC, aHBAS);
bIsDone=aHInter.IsDone();
if (!bIsDone)
{
di << " intersection is not done\n";
return 0;
}
aNbPoints=aHInter.NbPoints();
sprintf (buf, " Number of intersection points found: %d", aNbPoints);
di << buf << "\n";
for (i=1; i<=aNbPoints; ++i)
{
const IntCurveSurface_IntersectionPoint& aIP=aHInter.Point(i);
aIP.Values(aP, aU, aV, aT, aTC);
//
sprintf (buf, "point %s_%d %lg %lg %lg ", a[1], i, aP.X(), aP.Y(), aP.Z());
di << buf << "\n";
}
return 0;
}
void QABugs::Commands_19(Draw_Interpretor& theCommands) {
const char *group = "QABugs";
@ -4222,6 +4306,7 @@ void QABugs::Commands_19(Draw_Interpretor& theCommands) {
__FILE__, OCC26462, group);
theCommands.Add ("OCC26313", "OCC26313 result shape", __FILE__, OCC26313, group);
theCommands.Add ("OCC26525", "OCC26525 result edge face ", __FILE__, OCC26525, group);
return;
}

View File

@ -0,0 +1,52 @@
puts "================"
puts "OCC26525"
puts "================"
puts ""
#######################################################################
# Wrong result obtained by curve / surface intersection algorithm.
#######################################################################
pload QAcommands
restore [locate_data_file bug26525_a.brep] b1
restore [locate_data_file bug26525_b.brep] b2
mksurface sb1 b1
trimv sb1t sb1 -30000 30000
mkface b1t sb1t
explode b2 e
# Case 1. The curve is from the edge b2_3
# 1.1 The face b1 is based on untrimmed surface
set log [OCC26525 p b2_3 b1]
regexp {([-0-9]+)} $log full Number
if { $Number == 2} {
set tol_abs 0.0001
set tol_rel 0.01
regexp {point p_1 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x1 y1 z1
set expected_x1 48.4205
set expected_y1 -22.5336
set expected_z1 82.7431
checkreal "x1" ${x1} ${expected_x1} ${tol_abs} ${tol_rel}
checkreal "y1" ${y1} ${expected_y1} ${tol_abs} ${tol_rel}
checkreal "z1" ${z1} ${expected_z1} ${tol_abs} ${tol_rel}
regexp {point p_2 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x2 y2 z2
set expected_x2 32.5621
set expected_y2 -5.89907
set expected_z2 82.7431
checkreal "x2" ${x2} ${expected_x2} ${tol_abs} ${tol_rel}
checkreal "y2" ${y2} ${expected_y2} ${tol_abs} ${tol_rel}
checkreal "z2" ${z2} ${expected_z2} ${tol_abs} ${tol_rel}
} else {
puts "Error: Bad Number of intersection points"
}

View File

@ -0,0 +1,52 @@
puts "================"
puts "OCC26525"
puts "================"
puts ""
#######################################################################
# Wrong result obtained by curve / surface intersection algorithm.
#######################################################################
pload QAcommands
restore [locate_data_file bug26525_a.brep] b1
restore [locate_data_file bug26525_b.brep] b2
mksurface sb1 b1
trimv sb1t sb1 -30000 30000
mkface b1t sb1t
explode b2 e
# Case 1. The curve is from the edge b2_3
# 1.2 The face b1 is based on trimmed surface
set log [OCC26525 pt b2_3 b1t]
regexp {([-0-9]+)} $log full Number
if { $Number == 2} {
set tol_abs 0.0001
set tol_rel 0.01
regexp {point pt_1 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x1 y1 z1
set expected_x1 48.4205
set expected_y1 -22.5336
set expected_z1 82.7431
checkreal "x1" ${x1} ${expected_x1} ${tol_abs} ${tol_rel}
checkreal "y1" ${y1} ${expected_y1} ${tol_abs} ${tol_rel}
checkreal "z1" ${z1} ${expected_z1} ${tol_abs} ${tol_rel}
regexp {point pt_2 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x2 y2 z2
set expected_x2 32.5621
set expected_y2 -5.89907
set expected_z2 82.7431
checkreal "x2" ${x2} ${expected_x2} ${tol_abs} ${tol_rel}
checkreal "y2" ${y2} ${expected_y2} ${tol_abs} ${tol_rel}
checkreal "z2" ${z2} ${expected_z2} ${tol_abs} ${tol_rel}
} else {
puts "Error: Bad Number of intersection points"
}

View File

@ -0,0 +1,52 @@
puts "================"
puts "OCC26525"
puts "================"
puts ""
#######################################################################
# Wrong result obtained by curve / surface intersection algorithm.
#######################################################################
pload QAcommands
restore [locate_data_file bug26525_a.brep] b1
restore [locate_data_file bug26525_b.brep] b2
mksurface sb1 b1
trimv sb1t sb1 -30000 30000
mkface b1t sb1t
explode b2 e
# Case 2. The curve is from the edge: b2_1
# 2.1 The face b1 is based on untrimmed surface
set log [OCC26525 p b2_1 b1]
regexp {([-0-9]+)} $log full Number
if { $Number == 2} {
set tol_abs 0.0001
set tol_rel 0.01
regexp {point p_1 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x1 y1 z1
set expected_x1 39.0504
set expected_y1 -12.8696
set expected_z1 82.6099
checkreal "x1" ${x1} ${expected_x1} ${tol_abs} ${tol_rel}
checkreal "y1" ${y1} ${expected_y1} ${tol_abs} ${tol_rel}
checkreal "z1" ${z1} ${expected_z1} ${tol_abs} ${tol_rel}
regexp {point p_2 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x2 y2 z2
set expected_x2 43.1172
set expected_y2 -17.16
set expected_z2 82.6048
checkreal "x2" ${x2} ${expected_x2} ${tol_abs} ${tol_rel}
checkreal "y2" ${y2} ${expected_y2} ${tol_abs} ${tol_rel}
checkreal "z2" ${z2} ${expected_z2} ${tol_abs} ${tol_rel}
} else {
puts "Error: Bad Number of intersection points"
}

View File

@ -0,0 +1,53 @@
puts "================"
puts "OCC26525"
puts "================"
puts ""
#######################################################################
# Wrong result obtained by curve / surface intersection algorithm.
#######################################################################
pload QAcommands
restore [locate_data_file bug26525_a.brep] b1
restore [locate_data_file bug26525_b.brep] b2
mksurface sb1 b1
trimv sb1t sb1 -30000 30000
mkface b1t sb1t
explode b2 e
# Case 2. The curve is from the edge: b2_1
#
# 2.2 The face b1 is based on trimmed surface
set log [OCC26525 pt b2_1 b1t]
regexp {([-0-9]+)} $log full Number
if { $Number == 2} {
set tol_abs 0.0001
set tol_rel 0.01
regexp {point pt_1 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x1 y1 z1
set expected_x1 39.0504
set expected_y1 -12.8696
set expected_z1 82.6099
checkreal "x1" ${x1} ${expected_x1} ${tol_abs} ${tol_rel}
checkreal "y1" ${y1} ${expected_y1} ${tol_abs} ${tol_rel}
checkreal "z1" ${z1} ${expected_z1} ${tol_abs} ${tol_rel}
regexp {point pt_2 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x2 y2 z2
set expected_x2 43.1172
set expected_y2 -17.16
set expected_z2 82.6048
checkreal "x2" ${x2} ${expected_x2} ${tol_abs} ${tol_rel}
checkreal "y2" ${y2} ${expected_y2} ${tol_abs} ${tol_rel}
checkreal "z2" ${z2} ${expected_z2} ${tol_abs} ${tol_rel}
} else {
puts "Error: Bad Number of intersection points"
}