Linear Time Maximum Clearance Computation on a Grid Graph?

Question

I have a uniform NxN grid with a non-empty subset of vertices marked as obstacles. My goal is to compute, for each non-obstacle vertex, the "maximum clearance" from the obstacle set. In other words, the radius of the largest circle center at the vertex not overlapping with any obstacle. Or the nearest "gas station" for each vertex. [sorry if there is some standard terminology I should have used]

I would like a linear-time algorithm to compute this. Can anyone help?

My gut feeling is that it is possible. I have searched but didn't get a straight answer, but I know the following:

(1) if the distance metric is Manhattan, this can be solved by a simple flooding.

(2) I know a little bit about Level Set, and it may help, but is there something simpler?

(3) I understand Fortune's elegant sweeping algorithm to calculate Voronoi diagram in O(n log n), where n is the number of points. For my particular application, since the wavefront never contains more than N vertices, the complexity is O(N log N). The O(log N) is due to updates to a binary tree. I hope this can be reduced from O(N log N) to O(N) because of the nature of my problem?

I am willing to do this in multiple passes in order to get linear time. For example, one pass each for N, E, S, W, which hopefully somewhat simplifies Fortune's algorithm. Also I am not interested in the Voronoi vertices but rather the distances of the grid vertices from the obstacles.

Is this a known, solved problem? If so, source code or (more detailed) pseudo code?

If not, so far what I have in mind is to do this in two passes (Eastbound sweep and Westbound sweep):

Consider the East-bound sweep, I update one column of N vertices at a time. There is a set of exactly N "back" vertices (each is assigned an X coordinate, or -infinity if they are out of scope (OOS)) such that a partition is defined by the bisecting lines between immediate non-OOS neighbors (similar to Fortune). In other words, each non-OOS vertex is given a range [ys, ye] which together partitions [0, N - 1]. I don't maintain a binary tree of non-OOS vertices. Instead, I always keep them in an array of size N.

When we move to the next column, the partition changes. The partition is recalculated, and if a vertex's range becomes negative (ys > ye), the vertex is marked as OOS accordingly. Hopefully, this can be done in O(N) time and hence the total runtime is O(NxN).

Am I going in the right direction? Actually by spelling it out it already helps a bit :D

Answer 1

All the sites are themselves grid points too, right? So I think it should work to compute the Delaunay triangulation of the sites in linear time using the word-model algorithm of Buchin and Mulzer (FOCS 2009), dualize to get the Voronoi diagram, and then fill out the grid neighbors within each Voronoi cell. The Delaunay part should be linear in the number of sites (its bottleneck is integer sorting, but for the magnitude of the integers in your problem sorting is easily done in time linear in the grid size).

To fill out the grid points within a single Voronoi cell, partition the cell along the lines connecting the cell's site to its neighbors and then find the grid points within each of the resulting quadrilaterals. That way it requires only a constant number of comparisons to determine whether any particular grid point belongs to the quadrilateral. Some quadrilaterals may be empty of grid cells but the total number of quadrilaterals formed in this way is linear so there isn't too much wasted work in the empty ones.

Answer 2

This sounds like a perfect problem for a GPU. For each site (a,b), render the paraboloid z^2 = (x-a)^2 + (y-b)^2 into an nxn Z-buffer. At the end of the loop, the z-buffer contains the results you want.

Answer 3

First of all thanks David for the pointer. I believe his citation addresses the original formulation, but requires (1) a grid based on integer (or linear-time sortable) coordinates, and (2) a way to take the (non-gridded) Voronoi result and convert to the distance assignment for the grid vertices.

I feel uneasy about that, and in parallel I've been working on the approach I suggested at the end of my original question (it was more of a thinking-out-loud). It does solve the problem in linear time (in the size of the grid), and does not require integer or linear-time sortable coordinates, so I mark this as the answer.

I'll run this in two passes: a East-bound pass and a West-bound pass. The coordinates are $x_1, x_2, \dots, x_N$, and $y_1, y_2, \dots, y_N$.

WLOG in the East-bound pass, I'll process column by column but will only look at obstacle vertices on or to the left of the column.

I maintain an array $A$ of obstacle vertices of size $N$. Suppose now I am processing column $c$ ($x = x_c$). The $k$th element in $A$ stores the obstacle vertex with the largest $x$-coordinate such that $y = y_k$ and $x \le x_c$, or null if there is no such obstacle vertex. Note that if there are multiple obstacle vertices with the same $y$ coordinates, all but the East-most one are guaranteed not to matter.

Now, I want to partition the vertical line $x = x_c$ into ranges closest to the same obstacle vertex. Each partition $p_i$ is associated with an index $i$ such that $A[i]$ is the obstacle vertex closest to this partition. Similar to the argument behind Fortune's Voronoi algorithm, it is easy to see that the $y_i$'s are strictly ordered. Furthermore, the first and last non-null vertices in A (could be the same) both have their own associated ranges, since the first owns $[-\infty, \cdots]$ and the last owns $[\cdots, +\infty]$.

So we go through the array A and process all non-null vertices in order. We use a stack to represent the current ranges and process the vertices one by one. The stack represents the current partitioning of the line using all vertices considered so far. Each vertex $v_i$ in the stack owns a range $[LB_i, UB_i]$. Initially, the first vertex owns the range $[-\infty, +\infty]$.

Let's say we consider the next vertex $v_j$, and the current top-of-stack vertex is $v_k$ (note that $v_k$ has the range [LBk, +infinity]). We compute the bisector between $v_j$ and $v_k$. This bisector will intersect with the line $x = x_c$ at y (guaranteed to hit, because $v_j$ and $v_k$ have different y coordinates). If $y < LB_k$, we know that $v_k$ is subsumed by its neighbors -- meaning every point in the line $x = x_c$ is closer to some other vertices than $v_k$. In this case, we say $v_k$ is "subsumed" and we simply remove $v_k$ from the stack.

Once $v_k$ is removed, we have to repeat this exercise, because it is possible that new top-of-stack vertex can be proven to be subsumed. We repeat this until the top-of-stack vertex $v_k$ is not "subsumed" (and we're guaranteed to find such vertex, because the bottom-most vertex has a range $[-\infty, \cdots]$), and we simply update the upper bound of the range of $v_k$ to y, and then push $v_j$ onto the stack, with a range of $[y, +\infty]$.

The stack may grow and shrink, but the runtime is $O(N)$ per column by amortization. At each point within the column processing, the stack is guaranteed to not contain any "subsumed" vertices considering all the vertices seen so far. Therefore, at the end of the processing, no vertices are subsumed (for the purist this may require a more rigorous proof by contradiction). Once the column is process, we simply revisit the ranges on the stack and assign the distance accordingly.

Last but not least, I talked in terms of an $N\times N$ grid. The argument holds just the same for an $M\times N$ grid, namely there is an $O(M\times N)$-time algorithm to compute the maximum clearance for each grid point. Moreover, the grid can also be non-uniform (in edge lengths).