How to find number of points with same minimal distances on matrix

Question

I'm trying to find the number of points in a matrix with the same minimal distances.

Start with a MxN matrix, where M and N < 50000. There is given a set of fixed points, with their respective coordinates.
The problem is to find the number of points in the matrix such as the minimal distance from any point in the set can be reached for at least two point in the set. The distance is measured by moving one units a time horizonatally or vertically. You can't move diagonally.

An example would make things clearer.
Let's say we have 3 fixed point: (1,3), (3,1) and (3,6). Here the point (3,3) would be one point since the minimal distance is 2 units and it can be reached for the first two fixed points in the set. However, the point (4,4) doesn't satisfy the requirements. Although the distance from (1,3) and (3,1) is same, 4 units, the minimal distance is 3 units.

My idea was to use any point with integer coordinates that lies on the bisector of a segment joining each par as a point that would fill our requirements. But I find a lot of false positives were a point satisfies the requirements, but doesn't lie on the perpendicular bisector.

Brute-force methods won't work, because we need too check a billion of distance for each point and that one heck of a job to do.

How should I go about solving this problem?

Answer 1

As depicted here, what you are looking for is the 1st-order 2-dimensional Voronoi diagram under the Manhattan, or L1-metric. This is a quite non-trivial problem (to solve efficiently), fortunately with many existing algorithms and software.

You actually want the subset of the Voronoi diagram that coincides with the discrete grid defined by your matrix coordinates, but it is straightforward to find which points of the diagram happen to have integral coordinates.

The "fixed points" are called sites, seeds, generators, or source points in different contexts. Since the number of sites is small, it is more appropriate to use a computational geometry algorithm, that is, find a vector representation of the line segments the diagram consists of, independently of the size of your matrix (for 10 sites, it should be extremely fast).

A highly recognized library is CGAL, which includes a large number of computational geometry algorithms in C++, and in particular 2D segment Delaunay graphs (Voronoi diagrams). The user manual talks only about Euclidean distance, but there is also an L1 Voronoi diagram demo, so I hope you'll find your way.

Another perspective is a raster representation, using algorithms that sweep the entire plane (your matrix) once, in a particular order similar to wave propagation. This is recommended if the number of sites is large, forming arbitrary shapes on the plane. It is related to the distance transform, the medial axis and the topological skeleton. The latter two are generalizations of the Voronoi diagram, in that they consist of curves rather than line segments in general, regardless of the metric.

A very good solution under this perspective is distance transforms of sampled functions, also giving source code. This is limited to computing the distance transform and to separable distance metrics for instance squared Euclidean or L1 (your case).

To find the medial axis given the distance transform, a good choice is my own medial feature detector, with binary code available. This work also uses an alternative method for distance transform computation that has no limitation on the distance metric. It is slower than the previous one but still linear-time on the size of domain (your matrix). Computation of the medial axis itself is extremely fast (linear-time on the number of points of the medial axis).

As an example, you can have a look at this image, where white points are the sites (actually forming a closed curve on the plane), and red-yellowish curves are the (Euclidean) medial axis. Note that the latter have sub-pixel accuracy, i.e. appear as 1-pixel-wide anti-aliased curves, although the input sites do not. This illustrates that you have the required information (exact position) to find out whether a candidate point on the medial axis has integral coordinates or not. To get an idea of speed, this result was generated in around 300ms.

If you care about efficiency, my advice is to look for existing algorithms or code. If you insist on your own implementation, you will have to read a lot first, so that you have a good understanding of the problem.

Answer 2

You can easily see that if a point you are searching for it's at the same distance from two or more points, then they must all be at a distance such that d(a,b) is even.

Given this, you can iterate over all the pairs {a,b} of fixed points, and add to your set of result points all the matrix points that are half the distance d(a,b) from both the points you have chosen. To do it, consider the two "circles" (in your geometry it's a tilted square) of radius d(a,b)/2 and center a and b. All the matrix points that belong to both circles are possibly good points for you.

Now you have to check if every point you found at the previous step are actually good or false positives. To do so, just iterate over your result set and check the distances from all the fixed points.

The algorithm should be quite fast given the fact you have at most 10 fixed points