我正在计算给定横坐标 x 线上的点的纵坐标 y。该线由其两个端点坐标 (x0,y0)(x1,y1) 定义。端点坐标是浮点数，必须以浮点精度进行计算才能在 GPU 中使用。
数学以及因此幼稚的实现是微不足道的。
让 t = (x - x0)/(x1 - x0)，然后 y = (1 - t) * y0 + t * y1 = y0 + t * (y1 - y0)。
问题是当 x1 - x0 很小时。结果将引入取消错误。当与 x - x0 之一结合时，在除法中，我预计 t 会出现重大错误。
问题是是否存在另一种更准确地确定 y 的方法？
即我应该先计算 (x - x0)*(y1 - y0)，然后除以 (x1 - x0) 吗？
y1 - y0 的差异总是很大。
To a large degree, your underlying problem is fundamental. When (x1-x0) is small, it means there are only a few bits in the mantissa of x1 and x0 which differ. And by extension, there are only a limted number of floats between x0 and x1. E.g. if only the lower 4 bits of the mantissa differ, there are at most 14 values between them.
In your best algorithm, the t term represents these lower bits. And to continue or example, if x0 and x1 differ by 4 bits, then t can take on only 16 values either. The calculation of these possible values is fairly robust. Whether you're calculating 3E0/14E0 or 3E-12/14E-12, the result is going to be close to the mathematical value of 3/14.
Your formula has the additional advantage of having y0 <= y <= y1, since 0 <= t <= 1
(I'm assuming that you know enough about float representations, and therefore "(x1-x0) is small" really means "small, relative to the values of x1 and x0 themselves". A difference of 1E-1 is small when x0=1E3 but large if x0=1E-6 )
You may have a look at Qt's "QLine" (if I remember it right) sources; they have implemented an intersection determination algorithm taken from one the "Graphics Gems" books (the reference must be in the code comments, the book was on EDonkey a couple of years ago), which, in turn, has some guarantees on applicability for a given screen resolution when calculations are performed with given bit-width (they use fixed-point arithmetics if I'm not wrong).
如果您有可能这样做，您可以在计算中引入两种情况，具体取决于 abs(x1-x0) < abs(y1-y0)。在垂直情况下 abs(x1-x0) < abs(y1-y0)，从 y 计算 x 而不是从 x 计算 y。
编辑。另一种可能性是使用二分搜索的变体逐位获得结果。这会更慢，但在极端情况下可能会改善结果。
// 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;
I have implemented a benchmark program to compare the effect of the different expression.
I computed y using double precision and then compute y using single precision with different expressions.
Here are the expression tested:
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;
}
I computed the average and stdDev of the difference between the double precision result and single precision result.
The result is that there is none on the average over 1000 and 10K random value sets. I used icc compiler with and without optimization as well as g++.
Note that I had to use the isnan() function to filter out bogus values. I suspect these result from underflow in the difference or division.
I don't know if the compilers rearrange the expression.
Anyway, the conclusion from this test is that the above rearrangements of the expression have no effect on the computation precision. The error remains the same (on average).
如果您的源数据已经是浮点数，那么您已经存在根本的不准确性。
为了进一步解释，想象一下您是否以图形方式执行此操作。您有一张 2D 方格纸，并标记了 2 个点。
案例1：这些点非常准确，并且已经用非常锋利的铅笔做了标记。很容易画出连接它们的线，然后很容易在给定 x 的情况下得到 y（反之亦然）。
Case 2: These point have been marked with a big fat felt tip pen, like a bingo marker. Clearly the line you draw will be less accurate. Do you go through the centre of the spots? The top edge? The bottom edge? Top of one, bottom of the other? Clearly there are many different options. If the two dots are close to each other then the variation will be even greater.
Floats have a certain level of inaccuracy inherent in them, due to the way they represent numbers, ergo they correspond more to case 2 than case 1 (which one could suggest is the equivalent of using an arbitrary precision librray). No algorithm in the world can compensate for that. Imprecise data in, Imprecise data out
Check if the distance between x0 and x1 is small, i.e. fabs(x1 - x0) < eps. Then the line is parallell to the y axis of the coordinate system, i.e. you can't calculuate the y values of that line depending on x. You have infinite many y values and therefore you have to treat this case differently.
How about computing something like:
t = sign * power2 ( sqrt (abs(x - x0))/ sqrt (abs(x1 - x0)))
The idea is to use a mathematical equivalent formula in which low (x1-x0) has less effect. (not sure if the one I wrote matches this criteria)
As MSalters said, the problem is already in the original data.
Interpolation / extrapolation requires the slope, which already has low accuracy in the given conditions (worst for very short line segments far away from the origin).
Choice of algorithm canot regain this accuracy loss. My gut feeling is that the different evaluation order will not change things, as the error is introduced by the subtractions, not the devision.
Idea:
If you have more accurate data when the lines are generated, you can change the representation from ((x0, y0), (x1, y1)) to (x0,y0, angle, length). You could store angle or slope, slope has a pole, but angle requires trig functions... ugly.
Of course that won't work if you need the end point frequently, and you have so many lines that you can't store additional data, I have no idea. But maybe there is another representation that works well for your needs.
doubles have enough resolution in most situations, but that would double the working set too.