Most accurate line intersection crosses ordinate with float?

Question

I calculate the ordinate y of the point on the line for a given abscissa x. The line is defined by the two coordinates of the end points (x0, y0) (x1, y1). Endpoint coordinates are floats, and the calculation must be accurate to float for use in the GPU.

Mathematics, and thus the naive realization, is trivial.

Let t = (x - x0) / (x1 - x0), then y = (1 - t) * y0 + t * y1 = y0 + t * (y1 - y0).

The problem is that x1 - x0 is small. As a result, a cancellation error will be introduced. In combination with one of x - x0 in division, I expect a significant error in t.

The question is, is there another way to determine y with better accuracy?

i.e. should you first calculate (x - x0) * (y1 - y0) and divide by (x1 - x0) after?

The difference y1 - y0 will always be large.

Answer 1

To a large extent, your underlying problem is fundamental. When (x1-x0) is small, this means that in the mantissa x1 and x0 there are only a few bits that differ. And besides, there is only a limited number of floats between x0 and x1. For example. if only the lower 4 bits of the mantissa differ, no more than 14 values ​​between them.

t . , x0 x1 4 , t 16 . . 3E0/14E0 3E-12/14E-12, 3/14.

: y0 <= y <= y1, 0 <= t <= 1

( , , , "(x1-x0) " " x1 x0". 1E-1 x0 = 1E3, , x0 = 1E-6)

Answer 2

Qt "QLine" ( ) ; , " " ( , EDonkey ), , , ( , ).

Answer 3

, , abs (x1-x0) < (1-0). (x1-x0) abs (y1-y0), x y y x.

. , . , .

// Input is X
xmin = min(x0,x1);
xmax = max(x0,x1);
ymin = min(y0,y1);
ymax = max(y0,y1);
for (int i=0;i<20;i++) // get 20 bits in result
{
  xmid = (xmin+xmax)*0.5;
  ymid = (ymin+ymax)*0.5;
  if ( x < xmid ) { xmax = xmid; ymax = ymid; } // first half
  else { xmin = xmid; ymin = ymid; } // second half
}
// Output is some value in [ymin,ymax]
Y = ymin;

Answer 4

.

y, , y, .

:

inline double getYDbl( double x, double x0, double y0, double x1, double y1 )
{
    double const t = (x - x0)/(x1 - x0);
    return y0 + t*(y1 - y0);
} 

inline float getYFlt1( float x, float x0, float y0, float x1, float y1 )
{
    double const t = (x - x0)/(x1 - x0);
    return y0 + t*(y1 - y0);
} 

inline float getYFlt2( float x, float x0, float y0, float x1, float y1 )
{
    double const t = (x - x0)*(y1 - y0);
    return y0 + t/(x1 - x0);
} 

inline float getYFlt3( float x, float x0, float y0, float x1, float y1 )
{
    double const t = (y1 - y0)/(x1 - x0);
    return y0 + t*(x - x0);
} 

inline float getYFlt4( float x, float x0, float y0, float x1, float y1 )
{
    double const t = (x1 - x0)/(y1 - y0);
    return y0 + (x - x0)/t;
} 

stdDev .

, 1000 10K . icc , g++.

, isnan() . , .

, .

, . ( ).

Answer 5

, .

, , . 2- .

1: . , , y, x ( ).

2: , . , , , . ? ? ? , ? , . , .

, , - , , 2, 1 ( , librray). . ,

Answer 6

, x0 x1 , fabs (x1 - x0) < EPS. y , .. y x. y, .

Answer 7

- :

t = sign * power2 ( sqrt (abs(x - x0))/ sqrt (abs(x1 - x0)))

, (x1-x0) . ( , , , )

Answer 8

MSalters, .

/ , ( ).

canot . , , , - .

---

:
, ((x0, y0), (x1, y1)) (x0, y0, , ). , , ... .

Of course, this will not work if you need an endpoint often and you have so many lines that you cannot store additional data, I have no idea. But maybe there is another view that works well for your needs.

In most situations, doubles the resolution, but it doubles the working set.