Algorithms/Distance approximations

Calculating distances is common in spatial and other search algorithms, as well as in computer game physics engines. However, the common Euclidean distance requires calculating square roots, which is often a relatively heavy operation on a CPU.

Contents

Given (x1, y1) and (x2, y2), which is closer to the origin by Euclidean distance? You might be tempted to calculate the two Euclidean distances, and compare them:

d1 = sqrt(x1^2 + y1^2)
d2 = sqrt(x2^2 + y2^2)
return d1 > d2

But those square roots are often heavy to compute, and what's more, you don't need to compute them at all. Do this instead:

dd1 = x1^2 + y1^2
dd2 = x2^2 + y2^2
return dd1 > dd2

The result is exactly the same (because the positive square root is a strictly monotonic function). This only works for comparing distances though, not for calculating individual values, which is sometimes what you need. So we look at approximations.

Note that you can also use it as a "first pass" since it's never lower than the Euclidean distance. You could check if data points are within a particular bounding box, as a first pass for checking if they are within the bounding sphere that you're really interested in. In fact, if you take this idea further, you end up with an efficient spatial data structure such as a w:Kd-tree.

However, be warned that taxicab distance is not w:isotropic - if you're in a Euclidean space, taxicab distances change a lot depending on which way your "grid" is aligned. This can lead to big discrepancies if you use it as a drop-in replacement for Euclidean distance. Octagonal distance approximations help to knock some of the problematic corners off, giving better isotropy:

A fast approximation of 2D distance based on an octagonal boundary can be computed as follows.

Given two points and , Let (w:absolute value) and . If , approximated distance is .

Some years ago I developed a similar distance approximation algorithm using three terms, instead of just 2, which is much more accurate, and because it uses power of 2 denominators for the coefficients can be implemented without using division hardware. The formula is: 1007/1024 max(|x|,|y|) + 441/1024 min(|x|,|y|) - if ( max(|x|.|y|)<16min(|x|,|y|), 40/1024 max(|x|,|y|), 0 ). Also it is possible to implement a distance approximation without using either multiplication or division when you have very limited hardware: ((( max << 8 ) + ( max << 3 ) - ( max << 4 ) - ( max << 1 ) + ( min << 7 ) - ( min << 5 ) + ( min << 3 ) - ( min << 1 )) >> 8 ). This is just like the 2 coefficient min max algorithm presented earlier, but with the coefficients 123/128 and 51/128. I have an article about it at http://web.oroboro.com:90/rafael/docserv.php/index/programming/article/distance --Rafael