5
\$\begingroup\$

I want to compute the shortest distance between a position (x,y) and a rectangular box defined by (x_min, y_min) and (x_max, y_max). See the picture below with some examples. I know that in two dimensions there are 9 cases I have to take into account.

enter image description here

Since this part of my code is used very often, I would like to know what is the fastest possible way to do that? Here is what I have done so far:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define sqr(x) ((x)*(x))

double distanceBoxParticle2D(double x, double y, double x_min, double y_min, double x_max, double y_max);

int main(){
    double x = -6;
    double y = -6;
    double x_min = -5;
    double x_max = 5;
    double y_min = -5;
    double y_max = 5;

    double distance = distanceBoxParticle2D(x, y, x_min, y_min, x_max, y_max);
    printf("distance = %f\n", distance);

}

double distanceBoxParticle2D(double x, double y, double x_min, double y_min, double x_max, double y_max){
    double r = 0;

    if(x < x_min){
        if(y < y_min){
            return r = sqrt( sqr(x_min-x) + sqr(y_min-y) );
        }else if(y > y_max){
            return r = sqrt( sqr(x_min-x) + sqr(y_max-y) );
        }else if((y >= y_min) && (y <= y_max)){
            return r = x_min - x;
        }
    }else{
        if((x > x_min) && (x < x_max)){
            if(y < y_min){
                return r = y_min - y;
            }else if(y > y_max){
                return r = y - y_max;
            }
        }else{
            if(x > x_max){
                if(y < y_min){
                    return r = sqrt( sqr(x_max-x) + sqr(y_min-y) );
                }else if(y > y_max){
                    return r = sqrt( sqr(x_max-x) + sqr(y_max-y) );
                }else if((y >= y_min) && (y <= y_max)){
                    return r = x - x_max;
                }
            }
        }
    }

    return r = 0; // if (x,y) is inside of the box
}

The example is only for two dimensions. Later I want to do the same in 3D. There I have to consider 27 cases, which is a lot. Is there any smart and fast solution for this case or do I have to hard-code every case?

\$\endgroup\$
  • \$\begingroup\$ @chux corrected \$\endgroup\$ – Samuel Sep 13 '17 at 17:01
  • \$\begingroup\$ What are you doing with this code? Is this possibly an XY problem? Would a boolean function that checked for collisions, or containment, be enough? \$\endgroup\$ – Austin Hastings Sep 13 '17 at 22:22
  • \$\begingroup\$ If you are trying to find the particles inside the box, or close to the box an AABB Tree (azurefromthetrenches.com/…) or a grid could help \$\endgroup\$ – Harald Scheirich Sep 14 '17 at 12:06
2
\$\begingroup\$

Bug

If either x == x_min or x == x_max, your function will return 0 even if the point is outside of the rectangle. This is because none of your if statements included the case where x == x_min or x == x_max, so the flow falls through to the end of the function where you return 0.

Extraneous variable

Your function has a variable r that is set but never used. You should remove that variable so that lines like this:

        return r = x_min - x;

can turn into this:

        return x_min - x;

Extraneous third comparison

Consider this chain of three if statements:

    if(y < y_min){
        return r = sqrt( sqr(x_min-x) + sqr(y_min-y) );
    }else if(y > y_max){
        return r = sqrt( sqr(x_min-x) + sqr(y_max-y) );
    }else if((y >= y_min) && (y <= y_max)){
        return r = x_min - x;
    }

The point can either be above the rectangle, below the rectangle, or within the rectangle. If the point is neither above nor below the rectangle, then it must be within the rectangle. Therefore, the third comparison can be removed, like this:

    if(y < y_min){
        return r = sqrt( sqr(x_min-x) + sqr(y_min-y) );
    }else if(y > y_max){
        return r = sqrt( sqr(x_min-x) + sqr(y_max-y) );
    }else{
        return r = x_min - x;
    }

Rewrite

In this rewrite, I defined a HYPOT macro to simplify the expression that you use 4 times. I also modified the indentation so that each case would be at the same indentation level. Lastly, I reordered the if statements so the order would be "below, inside, above" instead of "below, above, inside", although the order really doesn't matter. This rewrite was a tiny bit faster (2%) than the original code, possibly due to removing the extraneous third comparisons mentioned above.

#define HYPOT(x, y) sqrt(sqr(x) + sqr(y))

double distanceBoxParticle2D(double x, double y, double x_min, double y_min,
        double x_max, double y_max)
{
    if (x < x_min) {
        if (y < y_min)
            return HYPOT(x_min-x, y_min-y);
        else if (y <= y_max)
            return x_min - x;
        else
            return HYPOT(x_min-x, y_max-y);
    } else if (x <= x_max) {
        if (y < y_min)
            return y_min - y;
        else if (y <= y_max)
            return 0;
        else
            return y - y_max;
    } else {
        if (y < y_min)
            return HYPOT(x_max-x, y_min-y);
        else if (y <= y_max)
            return x - x_max;
        else
            return HYPOT(x_max-x, y_max-y);
    }
}

Here is a more concise variant you might prefer:

double distanceBoxParticle2D(double x, double y, double x_min, double y_min,
        double x_max, double y_max)
{
    if (x < x_min) {
        if (y <  y_min) return HYPOT(x_min-x, y_min-y);
        if (y <= y_max) return x_min - x;
                        return HYPOT(x_min-x, y_max-y);
    } else if (x <= x_max) {
        if (y <  y_min) return y_min - y;
        if (y <= y_max) return 0;
                        return y - y_max;
    } else {
        if (y <  y_min) return HYPOT(x_max-x, y_min-y);
        if (y <= y_max) return x - x_max;
                        return HYPOT(x_max-x, y_max-y);
    }
}
\$\endgroup\$
  • \$\begingroup\$ Your answer really helped me! Is there an easy way to do the same calculations for the 3D case? It is kind of error-prone to do the same for 27 cases in 3D. \$\endgroup\$ – Samuel Oct 20 '17 at 12:41
  • \$\begingroup\$ @Samuel You can write the 3d code in a clean way, but it won't be the fastest. Or you can write it with the ugly 27 cases, but it will be fast. To write it cleanly, you can compute the dx, dy, and dz separately and just take the hypot of all 3. \$\endgroup\$ – JS1 Oct 21 '17 at 8:48
2
\$\begingroup\$

tl;dr: Your code is super fast. You could improve the formatting a bit, but that's about it.

Simpler code

Looking at the picture, I see that when \$x_{min} \le x \le x_{max}\$, the distance is the same as if \$x\$ were \$x_{min}\$.

Furthermore, if \$x > x_{max}\$, you can mirror it so that it becomes smaller than \$x_{min}\$.

These two observations reduce the number of cases you have to consider.

By the way, since you are surely not the first to solve this task, just search for the solution. Googling for "distance point rectangle" gave me this: https://gamedev.stackexchange.com/questions/44483/how-do-i-calculate-distance-between-a-point-and-an-axis-aligned-rectangle.

To calculate the distance, first convert the rectangle to the center-width-height form. Then, computing the distance is easy. In Java, the code looks like this:

static double sqrDistance(double xmin, double ymin, double xmax, double ymax, double px, double py) {
    double rx = (xmin + xmax) / 2;
    double ry = (ymin + ymax) / 2;
    double rwidth = xmax - xmin;
    double rheight = ymax - ymin;

    double dx = Math.max(Math.abs(px - rx) - rwidth / 2, 0);
    double dy = Math.max(Math.abs(py - ry) - rheight / 2, 0);
    return dx * dx + dy * dy;
}

public static void main(String[] args) {
    System.out.println(sqrDistance(0, 0, 2, 2,   1, 1)); // point inside rectangle
    System.out.println(sqrDistance(0, 0, 2, 2,   2, 2)); // point at corner of rectangle
    System.out.println(sqrDistance(0, 0, 2, 2,   3, 3)); // point outside of rectangle, away
    System.out.println(sqrDistance(0, 0, 2, 2,   3, 1)); // point outside of rectangle, y is facing the rectangle
}

It is trivial to translate this code to C and to extend it to 3D as well.

Note that the code calculates the square of the distance, to avoid the costly square root calculation.

Performance

Performance-wise, my above code is bad. For the worst-case invocation, it performs 3 add, 8 sub, 6 mul, 1 sqrt, 4 cmp. No wonder it is slow.

The original code performs much better: 1 add, 2 sub, 2 mul, 1 sqrt, 3 cmp, all in the worst case. Two of the comparisons can even be merged with the subtractions.

So even if it is more to write, your code is really good.

Benchmark

Here is the code I used to benchmark several implementations. I played around with the use_system constants, and especially the use_system_hypot made the code slow down from 1 second to 8 seconds. This is on Windows, G++ 6.4.0.

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <ctime>

namespace Math {

    constexpr bool use_system_fabs = true;
    constexpr bool use_system_fmax = false;
    constexpr bool use_system_hypot = false;
    constexpr bool use_system_fdim = true;

    static double abs(double x) {
        return use_system_fabs ? std::fabs(x) : x < 0.0 ? -x : x;
    }
    static double max(double a, double b) {
        return use_system_fmax ? std::fmax(a, b) : a > b ? a : b;
    }
    static double hypot(double x, double y) {
        return use_system_hypot ? std::hypot(x, y) : std::sqrt(x * x + y * y);
    }
    static double hypotsqr(double x, double y) {
        return x * x + y * y;
    }
    static double dim(double x, double y) {
        return use_system_fdim ? std::fdim(x, y) : max(x - y, 0.0);
    }
}

static double distance_samuel(
    double x, double y,
    double x_min, double y_min, double x_max, double y_max) {
  if (x < x_min) {
    if (y < y_min) return Math::hypot(x_min - x, y_min - y);
    if (y > y_max) return Math::hypot(x_min - x, y_max - y);
    return x_min - x;
  }
  if (x > x_max) {
    if (y < y_min) return Math::hypot(x - x_max, y_min - y);
    if (y > y_max) return Math::hypot(x - x_max, y_max - y);
    return x - x_max;
  }
  if (y < y_min) return y_min - y;
  if (y > y_max) return y - y_max;
  return 0.0;
}

static double distance_samuel_sqr(
    double x, double y,
    double x_min, double y_min, double x_max, double y_max) {
  if (x < x_min) {
    if (y < y_min) return Math::hypotsqr(x_min - x, y_min - y);
    if (y > y_max) return Math::hypotsqr(x_min - x, y_max - y);
    return x_min - x;
  }
  if (x > x_max) {
    if (y < y_min) return Math::hypotsqr(x - x_max, y_min - y);
    if (y > y_max) return Math::hypotsqr(x - x_max, y_max - y);
    return x - x_max;
  }
  if (y < y_min) return y_min - y;
  if (y > y_max) return y - y_max;
  return 0.0;
}

static double distance_rillig(
    double px, double py,
    double xmin, double ymin, double xmax, double ymax) {
  double rx = (xmin + xmax) / 2;
  double ry = (ymin + ymax) / 2;
  double rwidth = xmax - xmin;
  double rheight = ymax - ymin;
  double dx = Math::max(Math::abs(px - rx) - rwidth / 2, 0);
  double dy = Math::max(Math::abs(py - ry) - rheight / 2, 0);
  return Math::hypot(dx, dy);
}

static double distance_rillig_sqr(
    double px, double py,
    double xmin, double ymin, double xmax, double ymax) {
  double rx = (xmin + xmax) / 2;
  double ry = (ymin + ymax) / 2;
  double rwidth = xmax - xmin;
  double rheight = ymax - ymin;
  double dx = Math::max(Math::abs(px - rx) - rwidth / 2, 0);
  double dy = Math::max(Math::abs(py - ry) - rheight / 2, 0);
  return Math::hypotsqr(dx, dy);
}

static double distance_rillig_fdim(
    double px, double py,
    double xmin, double ymin, double xmax, double ymax) {
  double rx = (xmin + xmax) / 2;
  double ry = (ymin + ymax) / 2;
  double rwidth = xmax - xmin;
  double rheight = ymax - ymin;
  double dx = Math::dim(Math::abs(px - rx), rwidth / 2);
  double dy = Math::dim(Math::abs(py - ry), rheight / 2);
  return Math::hypot(dx, dy);
}

#define N 40000000
double x[N];
double y[N];

static void benchmark(
    const char *name, 
    double (*func)(double, double, double, double, double, double)) {
  double x_min = RAND_MAX / 3;
  double x_max = RAND_MAX / 3 * 2;
  double y_min = RAND_MAX / 3;
  double y_max = RAND_MAX / 3 * 2;

  clock_t t0 = std::clock();
  for (unsigned i = N; --i;) {
    (void) func(x[i], y[i], x_min, y_min, x_max, y_max);
  }
  clock_t t1 = std::clock();

  double diff = double(t1 - t0) / CLOCKS_PER_SEC;
  std::printf("%30s: %.3f\n", name, diff);
}

int main() {
  for (unsigned i = N; --i;) {
    x[i] = std::rand();
    y[i] = std::rand();
  }

#define B(fn) benchmark(#fn, fn)

  B(distance_samuel);
  B(distance_rillig);
  B(distance_samuel_sqr);
  B(distance_rillig_sqr);
  B(distance_rillig_fdim);
  std::printf("\n");
  B(distance_samuel);
  B(distance_rillig);
  B(distance_samuel_sqr);
  B(distance_rillig_sqr);
  B(distance_rillig_fdim);
}
\$\endgroup\$
  • \$\begingroup\$ I simulated this against OP's code and surprisingly found it around 1.7x slower. Even with removing the redundant rx,ry,rwidth,rheight \$\endgroup\$ – chux Sep 13 '17 at 14:49
  • \$\begingroup\$ Oh, that's bad. But the code looks so nice and harmless. :-/ \$\endgroup\$ – Roland Illig Sep 13 '17 at 14:54
  • \$\begingroup\$ It is nice tight code. Perhaps try the test harness with your machine and see if you receive similar results. \$\endgroup\$ – chux Sep 13 '17 at 15:33
  • \$\begingroup\$ Your Benchmark is not valid, you are doing the conversion from max/min to center and extent for every test, this should only be done once for every square. By using the same signature for both checks you are siding the results. Even if there are more than one square you and the representation is kept to the max/min form one would still only convert once. \$\endgroup\$ – Harald Scheirich Sep 14 '17 at 10:10
2
\$\begingroup\$

"what is the fastest possible " and "hard code every case" often go together.


Tried stream-lining OP's' code and only managed a marginally consistent 2% improvement.

#define HYPOT(a,b) (sqrt((a)*(a) + (b)*(b)))

double distanceBoxParticle2D_4(double x, double y, double x_min, double y_min,
    double x_max, double y_max) {
  if (x <= x_min) {
    if (y <= y_min) return HYPOT(x_min - x, y_min - y);
    if (y_max <= y) return HYPOT(x_min - x, y - y_max);
    return x_min - x;
  }
  if (x >= x_max) {
    if (y <= y_min) return HYPOT(x - x_max, y_min - y);
    if (y_max <= y) return HYPOT(x - x_max, y - y_max);
    return x - x_max;
  }
  if (y <= y_min) return y_min - y;
  if (y_max <= y) return y - y_max;
  return 0.0; // if (x,y) is inside of the box
}

From Roland Illig nice tight solution: about 70% slower.

double distanceBoxParticle2D_3(double px, double py, double xmin, double ymin,
    double xmax, double ymax) {
  static bool first = 1;
  static double rx, ry, rwidth, rheight;
  if (first) {
    first = 0;
    rx = (xmin + xmax) / 2;
    ry = (ymin + ymax) / 2;
    rwidth = xmax - xmin;
    rheight = ymax - ymin;
  }
  double dx = fmax(fabs(px - rx) - rwidth / 2, 0);
  double dy = fmax(fabs(py - ry) - rheight / 2, 0);
  return sqrt(dx * dx + dy * dy);  // square root needed?
}

Test harness: Random points about 1/3 left, middle and right each and 1/3 above, middle, below each. Variations tried - one is below.

double randd() {
  return (double) rand();
}

#define N 100000000
double x[N];
double y[N];

int main() {
  double x_min = RAND_MAX / 3.0;
  double x_max = RAND_MAX / 3.0 * 2;
  double y_min = RAND_MAX / 3.0;
  double y_max = RAND_MAX / 3.0 * 2;
  for (unsigned i = N; --i;) {
    x[i] = randd();
    y[i] = randd();
  }
  for (unsigned i = N; --i;) {  // exercise cache
    (void) distanceBoxParticle2D(x[i], y[i], x_min, y_min, x_max, y_max);
  }

  clock_t t0, t1;
  t0 = clock();
  for (unsigned i = N; --i;) {
    (void) distanceBoxParticle2D(x[i], y[i], x_min, y_min, x_max, y_max);
  }
  t1 = clock();
  double diff = 1.0 * (t1 - t0) / CLOCKS_PER_SEC;
  printf("orig: %f\n", diff);

  t0 = clock();
  for (unsigned i = N; --i;) {
    (void) distanceBoxParticle2D_4(x[i], y[i], x_min, y_min, x_max, y_max);
  }
  t1 = clock();
  double diff2 = 1.0 * (t1 - t0) / CLOCKS_PER_SEC;
  printf("test: %f\n", diff2);
  printf("%%: %.3f\n", diff2 / diff * 100);

  return 0;
}

Sample output. YMMV

orig: 1.498000
test: 1.466000
%: 97.864

Using Roland Illig suggestion of returning the square of the distance was significantly faster (3-4x) in all cases. If code can use the square of the distance rather than the distance, certainly that will make for the best performance increase.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.