Search code examples
optimizationmathematical-optimization

Is there an iterative way to calculate radii along a scanline?


I am processing a series of points which all have the same Y value, but different X values. I go through the points by incrementing X by one. For example, I might have Y = 50 and X is the integers from -30 to 30. Part of my algorithm involves finding the distance to the origin from each point and then doing further processing.

After profiling, I've found that the sqrt call in the distance calculation is taking a significant amount of my time. Is there an iterative way to calculate the distance?

In other words:

I want to efficiently calculate: r[n] = sqrt(x[n]*x[n] + y*y)). I can save information from the previous iteration. Each iteration changes by incrementing x, so x[n] = x[n-1] + 1. I can not use sqrt or trig functions because they are too slow except at the beginning of each scanline.

I can use approximations as long as they are good enough (less than 0.l% error) and the errors introduced are smooth (I can't bin to a pre-calculated table of approximations).

Additional information: x and y are always integers between -150 and 150

I'm going to try a couple ideas out tomorrow and mark the best answer based on which is fastest.

Results

I did some timings

  • Distance formula: 16 ms / iteration
  • Pete's interperlating solution: 8 ms / iteration
  • wrang-wrang pre-calculation solution: 8ms / iteration

I was hoping the test would decide between the two, because I like both answers. I'm going to go with Pete's because it uses less memory.


Solution

  • Just to get a feel for it, for your range y = 50, x = 0 gives r = 50 and y = 50, x = +/- 30 gives r ~= 58.3. You want an approximation good for +/- 0.1%, or +/- 0.05 absolute. That's a lot lower accuracy than most library sqrts do.

    Two approximate approaches - you calculate r based on interpolating from the previous value, or use a few terms of a suitable series.

    Interpolating from previous r

    r = ( x2 + y2 ) 1/2

    dr/dx = 1/2 . 2x . ( x2 + y2 ) -1/2 = x/r

        double r = 50;
        
        for ( int x = 0; x <= 30; ++x ) {
            
            double r_true = Math.sqrt ( 50*50 + x*x );
            
            System.out.printf ( "x: %d r_true: %f r_approx: %f error: %f%%\n", x, r, r_true, 100 * Math.abs ( r_true - r ) / r );
            
            r = r + ( x + 0.5 ) / r; 
        }
    

    Gives:

    x: 0 r_true: 50.000000 r_approx: 50.000000 error: 0.000000%
    x: 1 r_true: 50.010000 r_approx: 50.009999 error: 0.000002%
    ....
    x: 29 r_true: 57.825065 r_approx: 57.801384 error: 0.040953%
    x: 30 r_true: 58.335225 r_approx: 58.309519 error: 0.044065%
    

    which seems to meet the requirement of 0.1% error, so I didn't bother coding the next one, as it would require quite a bit more calculation steps.

    Truncated Series

    The taylor series for sqrt ( 1 + x ) for x near zero is

    sqrt ( 1 + x ) = 1 + 1/2 x - 1/8 x2 ... + ( - 1 / 2 )n+1 xn

    Using r = y sqrt ( 1 + (x/y)2 ) then you're looking for a term t = ( - 1 / 2 )n+1 0.36n with magnitude less that a 0.001, log ( 0.002 ) > n log ( 0.18 ) or n > 3.6, so taking terms to x^4 should be Ok.