bug in edge-edge intersection

During testing OCC 5.1 my program crashed while performing some simple boolean operations. I traced back the error and found it in the function IntTools_EdgeEdge::FindSimpleRoot(...).
There is an infinite loop ( while(1) {...} ) to perform a bisection alogorithm. The crux is the break condition
if( fabs(b-a) return x0; // = 0.5*(a+b)
where b and a are the endpoints of the actual interval an myEpsT is a small fix number ( with default value 1e-12 ). This is very dangerous, due to the binary representation of floating point numbers. So, if fabs(b-a) >= myEpsT and if there is no representable number in the interval [a,b], x0 = 0.5*(a+b) yields to a or b and the function never returns from the loop.
Ofcourse it is possible to change the value of myEpsT, but this is not practicable. First reason: normally the edge-edge-intersection is performed internally by boolean operations and from this stage the user has no chance to set myEpsT. Second reson: the correct value for myEpsT depends on the input parameters of the function IntTools_EdgeEdge::FindSimpleRoot(...) and cannot determined before.
A solution would be, to choose myEpsT by eps*0.5*(|a|+|b|), where eps is the smallest number, so that 1.0+eps != 1.0 (2^-52 on my computer) and a, b are the initial values for the interval, where bisection is performed ( the parameters tA and tB of the FindSimpleRoot function ).

Here is a simple example to produce the error, intersection of 2 linear edges:

gp_Pnt P11( -3000, 8000, 0 ); // point 1 on line 1
gp_Pnt P12( 2000, -5000, 0 ); // point 2 on line 1
gp_Pnt P21( 3000, 8000, 0 ); // point 1 on line 2
gp_Pnt P22( -2000, -5000, 0 ); // point 2 on line 2
gp_Dir R1( P12.X()-P11.X(), P12.Y()-P11.Y(), 0 ); // vector P11->P12
gp_Dir R2( P22.X()-P21.X(), P22.Y()-P21.Y(), 0 ); // vector P21->P22
gp_Lin L1( P11, R1 );
gp_Lin L2( P21, R2 );

// the ranges for the edges, here the eclidian distances between P11,P12 rsp. P21,P22
double r1 = sqrt( (P12.X()-P11.X())*(P12.X()-P11.X()) + (P12.Y()-P11.Y())*(P12.Y()-P11.Y()) );
double r2 = sqrt( (P22.X()-P21.X())*(P22.X()-P21.X()) + (P22.Y()-P21.Y())*(P22.Y()-P21.Y()) );

TopoDS_Edge E1, E2;
BRepBuilderAPI_MakeEdge EB1( L1, P11, P12 );
if( EB1.IsDone() ) {
E1 = EB1.Edge();
BRepBuilderAPI_MakeEdge EB2( L2, P21, P22 );
if( EB2.IsDone() ) {
E2 = EB2.Edge();

IntTools_EdgeEdge IEE;
IEE.SetEdge1( E1 );
IEE.SetEdge2( E2 );
IEE.SetRange1( 0, r1 );
IEE.SetRange2( 0, r2 );

IEE.Perform(); // never returns

MCV's picture


thank you for that debugging session.

By the way, for numerical reasons, one should never write stuff like " if ( a==b )", while a and b are both floating numbers (float, double or Standard_Real in OCC).
One better writes:
"if ( fabs(a-b)< err )", while err=10e-6 or so, which is the maximum tolerance allowed. Once the difference of a and b is smaller than that, they are regarded as numerically equal.

The same goes for "if (a==0.0)".
Better use "if (a