6

I'm working on a platform that has only integer arithmetic. The application uses geographic information, and I'm representing points by (x, y) coordinates where x and y are distances measured in meters. As an approximation, I want to compute the Euclidean distance between two points. But to do this I have to square distances, and with 32-bit integers, the largest distance I can represent is 32 kilometers. Not good. My needs are more on the order of 1000 kilometers. But I'd like to be able to resolve distances on a scale smaller than 30 meters.

Hence my question: how can I compute Euclidean distance, using only integer arithmetic, without overflow, on distances whose squares don't fit in a single word?

ETA: I would like to be able to compute distances, but I might settle for being able to compare them.

Norman Ramsey
  • 198,648
  • 61
  • 360
  • 533
  • 2
    You actually want to *really* compute the distances, or compare distances? Hope you won't kill a kitten when you see my comment. :) – gsamaras Feb 18 '16 at 23:57
  • 1
    @gsamaras Prefer to compute, but will settle for comparison. All kittens alive and happy! – Norman Ramsey Feb 19 '16 at 01:33
  • Is there a memory limit ? I've done something similar by precomputing values in a grid and storing those then doing an approximate bilinear interpolation. – yhenon Feb 19 '16 at 03:59
  • There are some clever algorithms in [this SO answer](http://stackoverflow.com/questions/3506404/fast-hypotenuse-algorithm-for-embedded-processor/3506633). – Raymond Chen Feb 19 '16 at 05:25

3 Answers3

2

Perhaps comparing the octagonal distance approximation would be sufficient?

Slightly more up to date is this article on fast approximate distance functions.

whybird
  • 1,108
  • 8
  • 19
0

I would leave the square root out of play, so that I can approximate the Euclidean distance. However, when comparing distances, this approach gives you 100% accuracy, since the comparison would be the same if you squared the distances.

I am pretty sure about that, since I had use that approach when searching for nearest neighbours in high dimensional spaces. You can check my code and the theory in kd-GeRaF.

gsamaras
  • 71,951
  • 46
  • 188
  • 305
  • 1
    Did you meant 'leave the _square root_ out of the play'? The distance I found [here](https://github.com/gsamaras/kd_GeRaF/blob/6e24152c28e258ed68c304489e9ee8608517268f/source/Division_Euclidean_space.h#L221) calculates the square. – minus one Dec 22 '21 at 10:03
  • @minusone yes, good catch! Answer updated, thanks! – gsamaras Dec 23 '21 at 20:45
0

I would recommend to use fixed point calculation using integers and then the distance approximation is already not too complicated.

First step is to choose some fixed point representation for our needs:

For example in case we need a number range for 1000km with 1m resolution we can use 20bits that would be 2^20 = 1,048,576. So we have around 10bits for fractions.

Then we need to implement the approximation we choose:

For example in case we select the following approximation:

h ≈ b (1 + 0.337 (a/b)) = b + 0.337 a AND assuming 0 ≤ a ≤ b

We will implement as follows:

  int32_t dx = (x1 > x2 ? x1 - x2 : x2 - x1);
  int32_t dy = (y1 > y2 ? y1 - y2 : y2 - y1);
  int32_t a = dx > dy ? dy : dx;
  int32_t b = dx > dy ? dx : dy;
  int32_t h = b + (345 * a >> 10);      /* 345.088 = 0.337 * 2^10 */

About overflow:

Adding two <+20.0> positive numbers will result a maximum of <+21.0> number. That is Ok.

The multiplication is also safe while we use numbers in a range of -1..1. In this case the result will also remain in the same range. In our case <+20.0> * <+0.10> will result <+20.10> numbers. That we convert back to <+20.0>.

There is one step here we need to pay attention. During the multiplication we will get temporary a number with <+20.10> that is already near to our 32bits limit.

Exact calculation

We can also calculate the exact distance using the following consideration:

h = b sqrt(1 + (a/b)^2) AND assuming 0 < b ≤ a

In tis case we also need to calculate the square root:

In case the a/b still significantly larger than one or too large to calculate the square of it, we can simplify the calculation to:

h = a

See the implementation here

minus one
  • 642
  • 1
  • 7
  • 28