15

For 3 points in 2D : 

P1(x1,y1), 
P2(x2,y2), 
P3(x3,y3)

I need to find a point P(x,y), such that the maximum of the manhattan distances 

max(dist(P,P1), 
    dist(P,P2), 
    dist(P,P3))

will be minimal.

Any ideas about the algorithm?

I would really prefer an exact algorithm.

moooeeeep
  • 31,622
  • 22
  • 98
  • 187
YAKOVM
  • 9,805
  • 31
  • 116
  • 217

2 Answers2

17

There is an exact, noniterative algorithm for the problem; as Knoothe pointed out, the Manhattan distance is rotationally equivalent to the Chebyshev distance, and P is trivially computable for the Chebyshev distance as the mean of the extreme coordinates.

The points reachable from P within the Manhattan distance x form a diamond around P. Therefore, we need to find the minimum diamond that encloses all points, and its center will be P.

If we rotate the coordinate system by 45 degrees, the diamond is a square. Therefore, the problem can be reduced to finding the smallest enclosing square of the points.

The center of a smallest enclosing square can be found as the center of the smallest enclosing rectangle (which is trivially computed as the max and min of the coordinates). There is an infinite number of smallest enclosing squares, since you can shift the center along the shorter edge of the minimum rectangle and still have a minimal enclosing square. For our purposes, we can simply use the one whose center coincides with the enclosing rectangle.

So, in algorithmic form:

  1. Rotate and scale the coordinate system by assigning x' = x/sqrt(2) - y/sqrt(2), y' = x/sqrt(2) + y/sqrt(2)
  2. Compute x'_c = (max(x'_i) + min(x'_i))/2, y'_c = (max(y'_i) + min(y'_i))/2
  3. Rotate back with x_c = x'_c/sqrt(2) + y'_c/sqrt(2), y_c = - x'_c/sqrt(2) + y'_c/sqrt(2)

Then x_c and y_c give the coordinates of P.

thiton
  • 35,651
  • 4
  • 70
  • 100
  • Won't some of the points be lost, since rotation by 45 degree produces quite a number of non-integral number? – nhahtdh Mar 11 '13 at 20:29
  • So it's just (1) rotate by 45° and (2) _then_ take the mean of the most distant points, separately for x and y, and (3) rotate back? Still can't believe it's so simple, but can't find a flaw either... +1 – tobias_k Mar 11 '13 at 20:29
  • @nhahtdh: Don't do it in integer arithmetics. If an integer result is desired, round the result. – thiton Mar 11 '13 at 20:30
  • 1
    @tobias_k: Don't forget to (3) rotate the mean point by -45. ;-) Thanks for the nice summary, will edit that in. – thiton Mar 11 '13 at 20:31
  • (It seems that I have this incorrect presumption that Manhattan distance only applies to integral coordinates. Well, I learn something here). But how do you enlarge the rectangle to a square? – nhahtdh Mar 11 '13 at 20:37
  • 2
    @nhahtdh: You don't need to explicitly do it. The center of *the* minimum rectangle is also the center of *a* minimum square. In most cases, there are many minimum squares because you can shift the center along the shorter edge of the minimum rectangle. – thiton Mar 11 '13 at 20:39
  • @thiton: That's true. And I think you should add that to your answer. – nhahtdh Mar 11 '13 at 20:45
  • 3
    If you are interested, you have just (re)discovered that in 2D, L_{\infty} and L_{1} norm are related by a 45 degree rotation. I was about to add an answer with this rotation business, but you beat me to it. +1. – Knoothe Mar 11 '13 at 20:54
  • @Knoothe: Thanks for the hint, and yes, I had to rediscover that fact. Shows how math can serve you :-). – thiton Mar 11 '13 at 20:59
  • 1
    I'd avoid using sqrt(2). You can just use x' = x+y; y' = x-y. – tmyklebu Mar 11 '13 at 21:50
  • 1
    why do you use sqrt(2) and not 1/sqrt(2) ?(sqrt(2) is not sin(45) – YAKOVM Mar 11 '13 at 21:58
  • @tmyklebu: You'll need the back-transform. Can be done with 1 and 0.5, but 1/sqrt(2) is more symmetrical. – thiton Mar 11 '13 at 23:06
4

If an approximate solution is okay, you could try a simple optimization algorithm. Here's an example, in Python

import random
def opt(*points):
    best, dist = (0, 0), 99999999
    for i in range(10000):
        new = best[0] + random.gauss(0, .5), best[1] + random.gauss(0, .5)
        dist_new = max(abs(new[0] - qx) + abs(new[1] - qy) for qx, qy in points)
        if dist_new < dist:
            best, dist = new, dist_new
            print new, dist_new
    return best, dist

Explanation: We start with the point (0, 0), or any other random point, and modify it a few thousand times, each time keeping the better of the new and the previously best point. Gradually, this will approximate the optimum.

Note that simply picking the mean or median of the three points, or solving for x and y independently does not work when minimizing the maximum manhattan distance. Counter-example: Consider the points (0,0), (0,20) and (10,10), or (0,0), (0,1) and (0,100). If we pick the mean of the most separated points, this would yield (10,5) for the first example, and if we take the median this would be (0,1) for the second example, which both have a higher maximum manhattan distance than the optimum.

Update: Looks like solving for x and y independently and taking the mean of the most distant points does in fact work, provided that one does some pre- and postprocessing, as pointed out by thiton.

Community
  • 1
  • 1
tobias_k
  • 81,265
  • 12
  • 120
  • 179