Introduction

[♪ music ♪] >> KEVIN: Hi, I’m Kevin, and I’m a geek. I really like working on mathematical problems, and explaining them to people. I’ve retired from a job at a big industrial research lab where a lot of my work consisted of just that Now that I’m on my own for a while, I’m looking for whether I can get the same sort of fun out of explaining things in videos. A lot of what I’ve done involves various aspects of geometry, but I’m likely to stray far and wide because the field of mathematics is all connected. >> KEVIN: Today, I’ll get into a question that I was just thinking about a day or two ago - the problem of, given a set of points of interest, finding the distance from every point on the grid to the nearest point of interest. You see, I was hiking in a state forest near me, and looking for a place to set up camp. That particular forest had no official campsites, but there is a general rule that any site is lawful as long as it is more than 150 feet from a road, trail or water source, and below an altitude of 3500 feet. Because I’m an image processing geek, I started musing: given images of a map with lines representing the roads and trails, and areas representing watercourses, and given a digital elevation model, could I produce shading indicating the lawful areas? That’s one instance of the problem that I’m looking at in this video: given a set of points of interest on a grid, compute the distance from every grid point to its nearest point of interest. For this particular problem instance, the points of interest are the points along the trails or covered with water, and I want to compute the distance to the nearest such point.

Voronoi tessellations and Delaunay triangulations

If we had a discrete point set representing the locations of gas stations on the map, this technique could answer the question, “how far is it to the nearest gas station?” over the entire grid. More important, in image processing, there is often a need to answer the question, “how far is it to the nearest boundary of an identified object?” or related questions like, “what is the diameter of the object that we’re examining?” For a concrete example, look at this set of points. For each point, we want to find the grid points that are closer to it than any other, and to compute the distance from each grid point to the nearest point of interest. Those of you who know some computational geometry will recognize this as the problem of finding the “Voronoi tessellation”, named after the Ukrainian mathematician Gyorgy Feodosevich Voronoy. It divides the plane into a polygon for each point of interest, that contains that point’s region of influence - the part of the plane that is closer to that point than to any other point of interest. Each of the edges of the polygons is equidistant between some pair of points of interest, which is to say that it is the perpendicular bisector of the line joining those points. The line joining the pair of points at lower right has a perpendicular bisector defining the edge between the corresponding regions. The line joining the points at upper right again has a perpendicular bisector separating the regions and so on through the entire diagram. The fact that the dividing line is a perpendicular bisector will be important in the computation. It’s important to note that the only lines we bisected are the ones joining points whose regions are adjacent. This collection of line segments forms a mesh of triangles called the "Delaunay triangulation," after the Russian mathematician Boris Nikolaevich Delaunay, who was a student of Voronoy. There are several good computational geometry libraries, such as CGAL and qhull, that have routines to calculate Voronoi tessellations and Delaunay triangulations. Some of them even have n-dimensional versions, useful for cases such as using the Voronoi diagram for statistical clustering of data. I’ll leave links to a couple down in the video description.. One thing to note is that these libraries are often designed to handle relatively sparse discrete point sets, rather than to solve the problem over a dense grid, and they can be relatively slow for the specific problem we’re looking at here. Instead, I’m going to cover an algorithm described by Wang and Tan - again, there should be a link to the paper down in the video description - that’s particularly adapted to image processing. That’s what the rest of this video will be about.

What happens on a grid row?

The algorithm works by moving over the pixels row by row, top-to-bottom and left-to right. To understand it, we need only to look at what happens for a single row, since the processing for every row is the same and doesn’t need to refer to results from other rows. First, note that if we have two points of interest in the same column of the image, we only ever need to look at the closest one. The boundary between the two regions of influence will be horizontal and can never intersect the row. Now examine two points that are not in the same column, and look at their perpendicular bisector. If the slope of the line joining the points is steep enough. the bisector will intersect the row outside the bounds of the image, either to the left or to the right In that case, we can discard the point on that end, the leftmost if the bisector falls to the left and the rightmost if it falls to the right. Again. its region of influence cannot intersect the row. Finally, we look at the bisectors among three points. If we have a pair of points A and B, whose bisector falls somewhere on the row, and add a third point C, B’s region of influence will no longer intersect the row if the bisector of BC intersects to the left of AB’s bisector. This rule holds even if the point C is on the opposite side of the row. We don’t need any special-case computations for points above and below.

Data structure for Voronoi on a grid

These rules suggest a procedure and a data structure. We will be processing a set of feature points from left to right, and we will be maintaining a set of x coordinates where their regions of influence begin on the row. The x coordinates will also be in increasing order, left to right. We keep an array of pairs with the active feature points and their x coordinates. 91 In the algorithm, we add, delete and examine points only at the right end of the array, so the array is actually a stack. Its maximum depth is the width of the image. In the point set shown here, we see in the array that point A’s region of influence starts at the left edge, point C’s starts about a third of the way in, and point D’s is just over half way across.

Algorithm for grid rows

We begin with the stack empty, and the first point we encounter will always have its region of influence tentatively begin at the left end of the row. Let’s move on to the subsequent points. When we come in with a subsequent point, we find the x coordinate where its region of influence separates from the one to the left. If it lies to the right of the coordinate at the top of the stack, we're done and we can push the point and the new coordinate onto the stack. If it lies to the left of the previous coordinate, then the previous point is no longer relevant, like B in the diagram. We pop it off the stack, and loop around to examine the new stack top. When either the coordinates are in sequence or the stack becomes empty, we push the new coordinate, and move on to the next point. We continue this process for all the points in the row. Occasionally, we’ll encounter a point whose intersection is off the right edge of the row. We can simply discard it. That’s easier than stacking and unstacking it. Once we’ve completed the row, a second pass through the pixels can label them with the point of interest or with the distance, by simply reading the information out of the stack.

Performance

Now let’s look at the performance of the algorithm. Every grid cell in the row is processed at most once. That processing step can push at most one element to the stack. An element that is pushed to the stack can be popped only once, and the loop in the cell processing pops an element on every trip through the loop. This means that the algorithm does a constant amount of computation for every grid point that it visits. This is exactly what we need - it will make our performance scale linearly with the size of the image.

Exact arithmetic

There’s still some more to the Yang-Tan paper. They present a clever geometric construction to find the point where the perpendicular bisector of a pair of points intersects a row. Start by constructing the line and bisecting it. The segments AD and BD are equal, and the angles ADC and BDC are both right angles, so are equal as well. Drawing the line segments AC and BC, we can see that the pink and blue triangles are congruent (Side-Angle-Side), making the two line segments equal and verifying that the point C is indeed equidistant from A and B. The squares of the hypotenuses, of course, are also equal. Drop perpendiculars from A and B to the row, with their feet at U and V. Algebraically, this is just subtraction of the y coordinates. This gives us two more right triangles, and we can apply Pythagoras to both of them. Place the point O at the left end of the row. The x coordinates of the feature points are the lengths of the segments OU and OV, and the x coordinate of the intersection is the length of the segment OC. We can subtract these coordinates to give us formulas for UC and VC, and square both sides of the equations. Substitute these formulas into the equation above. The two squares of OC cancel out. Now comes the clever part. Construct the segments OA and OB, and we have right triangles OUA and OVB. We can apply Pythagoras to those triangles, and then subsitute into the last equation, giving a much tidier formula. Let me take a moment to neaten the chalkboard. Rearrange terms to put the OC terms together. Factor the left hand side Divide both sides by a common denominator, and we have solved for OC. This is a really elegant formula, because it allows for exact computation. There's no floating-point required, and no square root operations, just two subtractions and a division. Of course, we have to precompute the squares of OA and OB. Let me tidy up again and keep just the intersection formula.

Putting it all together

We can compute the squares at the same time as we identify A and B as being the closest feature points to the row in their respective columns. There's a simple, linear-time first pass to find the closest points to each row in each column, and fill in the squared distances. The algorithm for processing rows then runs on a second pass to identify the two-dimensional closest points and fill in the intersection points and distances. When the first pass is processing a column, it neither uses nor produces data pertaining to any other column, so it can be broken into as many threads as there are columns, without one thread interfering with another. Similarly, when the second pass is processing a row, it doesn’t use or produce data pertaining to any other row. There can be as many threads as there are rows, again without interfering. Essentially, the algorithm can be parallelized perfectly on a multicore machine. It's therefore an extremely fast computation on that sort of architecture.

Next steps

And that’s really all there is to the algorithm. It’s one of the fundamental building blocks for image processing. Now I want to make further progress to the processing of an image of coins, intending to count them up. We'll use the k-means clustering that we saw in a previous video to separate the coins from the background and then use the algorithm we saw in this video to find distances from the edges of the coins. From that, the next step will be to locate the centers of the coins and find their sizes. For that, we’ll need an algorithm that I mentioned several times in the past: watershed segmentation. I think now that I’ve presented all the tools we need to develop that algorithm, so I’ll begin that next time. [♪ music builds ♪] Until then, stay safe, stay healthy, and keep calculating! [♪ music ♪]