point-to-point distance algorithms survey

distance.cc

// point-to-point distance algorithms survey
// - rlyeh, public domain.

// possible output {
// manhattan_distance: 0.608632s. (result: 6)
// octagonal_distance: 0.651783s. (result: 3.52249)
// baptista_distance: 1.10458s. (result: 4.79492)
// baptista_distance_b: 1.09981s. (result: 4.64063)
// euclidean_distance: 1.25s. (result: 4.47214)
// euclidean_distance_fast: 1.03207s. (result: 4.5)  
// }

#include <iostream>

struct point {
    float x, y;
};

float baptista_distance( point a, point b );
float baptista_distance_b( point a, point b );
float euclidean_distance( point a, point b );
float euclidean_distance_fast( point a, point b );
float manhattan_distance( point a, point b );
float octagonal_distance( point a, point b );

#include <omp.h>

int main() {
    point a { 0, 0 }, b { 4, 2 };

    #define TEST(x) do { \
    auto result = x(a,b); \
    auto time = -omp_get_wtime(); \
    for( int i = 0; i < 100000000; ++i) x(a,b); \
    time += omp_get_wtime(); \
    std::cout << #x ": " << time << "s. (result: " << result << ")" << std::endl; \
    } while(0)

    TEST( manhattan_distance );
    TEST( octagonal_distance );
    TEST( baptista_distance );
    TEST( baptista_distance_b );
    TEST( euclidean_distance );
    TEST( euclidean_distance_fast );
}


#include <math.h>

float manhattan_distance( point a, point b ) {
    float dx = a.x > b.x ? a.x - b.x : b.x - a.x;
    float dy = a.y > b.y ? a.y - b.y : b.y - a.y;
    return dx + dy;
}

float octagonal_distance( point a, point b ) {
    // [ref] https://en.wikibooks.org/wiki/Algorithms/Distance_approximations
    float dx = a.x > b.x ? a.x - b.x : b.x - a.x;
    float dy = a.y > b.y ? a.y - b.y : b.y - a.y;
    return dx >= dy ? 0.41f * dx + 0.941246f * dy : 0.41f * dy + 0.941246f * dx;
}

float baptista_distance( point a, point b ) {
   // [ref] http://www.flipcode.com/archives/Fast_Approximate_Distance_Functions.shtml
   unsigned int dx = a.x > b.x ? a.x - b.x : b.x - a.x;
   unsigned int dy = a.y > b.y ? a.y - b.y : b.y - a.y;

   unsigned int min, max;

   if ( dx < dy ) {
      min = dx;
      max = dy;
   } else {
      min = dy;
      max = dx;
   }

   return ( ( max * 1007 ) + ( min * 441 ) ) / 1024.f;
} 

float baptista_distance_b( point a, point b ) {
   // [ref] http://www.flipcode.com/archives/Fast_Approximate_Distance_Functions.shtml
   unsigned int dx = a.x > b.x ? a.x - b.x : b.x - a.x;
   unsigned int dy = a.y > b.y ? a.y - b.y : b.y - a.y;

   unsigned int min, max;

   if ( dx < dy ) {
      min = dx;
      max = dy;
   } else {
      min = dy;
      max = dx;
   }

   return ( ( 123 * max ) + ( 51 * min ) ) / 128.f;
} 

float euclidean_distance( point a, point b ) {
    float dx = b.x - a.x;
    float dy = b.y - a.y;
    return sqrt( dx * dx + dy * dy );
}

float euclidean_distance_fast( point a, point b ) {
    // [ref] https://en.wikibooks.org/wiki/Algorithms/Distance_approximations
    /* Assumes that float is in the IEEE 754 single precision floating point format and that int is 32 bits. */
    auto isqrt = [](float x) -> float {
        float xhalf = 0.5f*x;
        union {
            float x;
            int i;
        } u;
        u.x = x;
        u.i = 0x5f3759df - (u.i >> 1);
        /* The next line can be repeated any number of times to increase accuracy */
        u.x = u.x * (1.5f - xhalf * u.x * u.x);
        return u.x;
    };
    auto ssqrt = [](float x) -> float {
      unsigned int i = *(unsigned int*) &x;
      // adjust bias
      i  += 127 << 23;
      // approximation of square root
      i >>= 1;
      return *(float*) &i;
    };
    float dx = b.x - a.x;
    float dy = b.y - a.y;
    return ssqrt( dx * dx + dy * dy );     // less accurate
    return 1 / isqrt( dx * dx + dy * dy ); // more accurate; extra division though.
}