<p>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.</p>

<p>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.</p>

<p>I'll run this in two passes: a East-bound pass and a West-bound pass.
The coordinates are <span class="math-container">$x_1, x_2, \dots, x_N$</span>, and <span class="math-container">$y_1, y_2, \dots, y_N$</span>.</p>

<p>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.</p>

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

<p>Now, I want to partition the vertical line <span class="math-container">$x = x_c$</span> into ranges closest
to the same obstacle vertex.  Each partition <span class="math-container">$p_i$</span> is associated with an
index <span class="math-container">$i$</span> such that <span class="math-container">$A[i]$</span> is the obstacle vertex closest to this
partition.  Similar to the argument behind Fortune's Voronoi
algorithm, it is easy to see that the <span class="math-container">$y_i$</span>'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 <span class="math-container">$[-\infty, \cdots]$</span> and the
last owns <span class="math-container">$[\cdots, +\infty]$</span>.</p>

<p>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 <span class="math-container">$v_i$</span> in the stack owns a range <span class="math-container">$[LB_i, UB_i]$</span>.  Initially, the
first vertex owns the range <span class="math-container">$[-\infty, +\infty]$</span>.</p>

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

<p>Once <span class="math-container">$v_k$</span> 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 <span class="math-container">$v_k$</span> is not "subsumed" (and we're guaranteed
to find such vertex, because the bottom-most vertex has a range
<span class="math-container">$[-\infty, \cdots]$</span>), and we simply update the upper bound of the range of
<span class="math-container">$v_k$</span> to y, and then push <span class="math-container">$v_j$</span> onto the stack, with a range of <span class="math-container">$[y, +\infty]$</span>.</p>

<p>The stack may grow and shrink, but the runtime is <span class="math-container">$O(N)$</span> 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.</p>

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

CC BY-SA 2.5
2010-09-03T18:24:42.600Z
Yes, I think this works.
David Eppstein
CC BY-SA 2.5
2010-09-05T03:18:10.980Z
Also, this algorithm is massively parallelizable.  You can run the algorithm on each column separately.  Only column-specific input is the array A, which can be precomputed in a parallel fashion (but parallelization will be done on a per-row basis).
Hardy Leung
CC BY-SA 4.0
2010-09-03T18:16:27.267Z
2019-01-25T00:20:21.520Z
2019-01-25T00:20:21.520Z
Hardy Leung
<p>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 (<a href="http://dx.doi.org/10.1109/FOCS.2009.53">FOCS 2009</a>), 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).</p>

<p>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.</p>

CC BY-SA 2.5
2010-09-03T05:22:35.087Z
David, thanks for the reference.  In the original formulation, it is a simple uniform grid graph with unit edge length.  So yes I think the work you cited proves that this takes linear time.  Thanks!

Theory aside, I am a bit hesitant in taking this approach, due to implementation complexity.  I am currently working on what I proposed (seems to work):

That is, given a set S of N points to the left of X=K, already sorted by Y-coordinates, it takes O(N) time to partition the line X=K into ranges s.t. each range is associated with a point in S.  Seems easy to code (unless my logic is flawed).
Hardy Leung
CC BY-SA 2.5
2010-09-03T03:57:04.893Z
2010-09-03T03:57:04.893Z
David Eppstein
<p>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.</p>

CC BY-SA 2.5
2010-09-03T17:48:06.457Z
Not sure if I understand your technique, but there could be up to NxN "sites" (which I assume you mean the obstacle vertices).  Given a 1000x1000 grid, this would take 1M renders which I presume would be expensive even for GPU?
Hardy Leung
CC BY-SA 2.5
2010-09-03T18:14:43.543Z
See, e.g. http://gamma.cs.unc.edu/VORONOI/ for some existing work on GPU computation of rasterized Voronoi diagrams
David Eppstein
CC BY-SA 2.5
2010-09-04T20:57:44.987Z
No, the whole point of GPUs is their massive per-pixel parallelism.  Yes, you'd still need to ship all the paraboloids down the graphics pipeline, but that takes only constant time per site.  Then (effectively) each pixel processes the data in parallel.
Jeffε
CC BY-SA 2.5
2010-09-03T15:53:41.683Z
2010-09-03T15:53:41.683Z
Jeffε
<p>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.  [<em>sorry if there is some standard terminology I should have used</em>]</p>

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

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

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

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

<p>(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?</p>

<p>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.</p>

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

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

<p>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.</p>

<p>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).</p>

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

CC BY-SA 2.5
2010-09-03T02:28:20.940Z
When you say "linear" what do you really mean ? Since there are possible Theta(N^2) locations, even reporting the answer can't be done in O(N) time.
Suresh Venkat
CC BY-SA 2.5
2010-09-03T05:07:35.743Z
Yes by linear I mean in the number of grid points (NxN), not in N.  Sorry about the confusion.
Hardy Leung
<ds.algorithms><cg.comp-geom>
Linear Time Maximum Clearance Computation on a Grid Graph?
Hardy Leung
# Linear Time Maximum Clearance Computation on a Grid Graph?

(Asked by: username_0 on 2010-09-02)

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

**Tags:** ds.algorithms, cg.comp-geom
> **Comments:**
> **username_1**: When you say "linear" what do you really mean ? Since there are possible Theta(N^2) locations, even reporting the answer can't be done in O(N) time.  
> **username_0**: Yes by linear I mean in the number of grid points (NxN), not in N.  Sorry about the confusion.  

---

## 3 answers

### Answer by username_0 on 2010-09-03. Score: 2 | Accepted answer

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).

> **Comments:**
> **username_2**: Yes, I think this works.  
> **username_0**: Also, this algorithm is massively parallelizable.  You can run the algorithm on each column separately.  Only column-specific input is the array A, which can be precomputed in a parallel fashion (but parallelization will be done on a per-row basis).  

---

### Answer by username_2 on 2010-09-03. Score: 8

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](http://dx.doi.org/10.1109/FOCS.2009.53)), 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.

> **Comments:**
> **username_0**: David, thanks for the reference.  In the original formulation, it is a simple uniform grid graph with unit edge length.  So yes I think the work you cited proves that this takes linear time.  Thanks!  Theory aside, I am a bit hesitant in taking this approach, due to implementation complexity.  I am currently working on what I proposed (seems to work):  That is, given a set S of N points to the left of X=K, already sorted by Y-coordinates, it takes O(N) time to partition the line X=K into ranges s.t. each range is associated with a point in S.  Seems easy to code (unless my logic is flawed).  

---

### Answer by username_3 on 2010-09-03. Score: 3

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.

> **Comments:**
> **username_0**: Not sure if I understand your technique, but there could be up to NxN "sites" (which I assume you mean the obstacle vertices).  Given a 1000x1000 grid, this would take 1M renders which I presume would be expensive even for GPU?  
> **username_2**: See, e.g. http://gamma.cs.unc.edu/VORONOI/ for some existing work on GPU computation of rasterized Voronoi diagrams  
> **username_3**: No, the whole point of GPUs is their massive per-pixel parallelism.  Yes, you'd still need to ship all the paraboloids down the graphics pipeline, but that takes only constant time per site.  Then (effectively) each pixel processes the data in parallel.  

---

