Introduction

Hi, it's Kevin again. I'm back to continue our exploration of building a software Hollerith card reader using a document scanner and image processing. In this episode, I'll be covering a technique called the fast discrete Radon transformation, which overcomes some of the limitations of the Hough transformation for identifying the straight lines in images. In the last couple of episodes - there's a link somewhere nearby - we scanned a card, used k-means clustering to separate the paper from the green screen, did some elementary image processing to clean up the segmentation, and found straight lines using the Hough transform. I confessed at the very end that I told a little white lie about using the Hough transform. In this episode, let's take a look at what I actually did.

Review: the Hough transform

Recall that the input to the Hough transform is a binary image of edges - that is, it's an array of bits, which has ones in the positions where a pixel is at a location that we detected as being on an edge, and zeroes elsewhere. The output of the Hough transform is a 2-d array that counts up the number of pixels that lie on each possible straight line through the image. The array indices correspond to two parameters that define a line. One parameter, that we name theta, is the angle that the line makes with one of the axes. This angle will vary from zero to pi. A given value of theta selects a particular pencil of parallel lines. The other parameter, which we label rho, is an offset from the origin that selects a particular line from the pencil. For the Hough transform, it's customary to place the origin at the image center, and rho will range between plus and minus the image diagonal. To compute the Hough transform, we go through every pixel in the image that we identified as being on an edge. For each of those pixels, we go through all the possible angles and compute the offset from the origin to a line through the given pixel at the given angle. We increment the element in the array corresponding to that angle and offset. At the end, the most prominent local maxima correspond to the longest straight lines in the image.

Problems with Hough transform

There are a number of problems with the basic Hough transform. First, it's quite sensitive to speckle noise. If there are a significant number of false detections, it will start placing false lines all over the place. Next, it's extremely sensitive to position errors. If a line appears curved by even a few pixels - as in the punched card, where the paper is slightly warped - it will come through as multiple short lines that approximate the curve, and may not have a clear enough local maximum for detection. Finally, it's excruciatingly slow. For every edge pixel, it needs to do the offset calculation for every possible angle. For a high-resolution image, this is a lot of trips through a loop. It gets even worse if the count array is larger than the processor cache, because the counts being incremented don't lie in contiguous memory. We can address the first two issues by using fuzzy logic. We assign a value to pixels that expresses a relative likelihood of their being on a line, and replace pixel counts with the sums of the relative likelihoods. But now we've made it even slower. Rather than considering only edge pixels, we have to consider any pixel with a nonzero likelihood - and that's usually a big fraction of the image. We wind up doing big O of N cubed operations on an N by N image - and practically every trip through the inner loop is a cache miss. The computation slows to a glacial pace.

Discrete Radon transform

We can solve all of this by looking at a closely-related technique called the Radon transformation. It was first published in 1917 bu the Austrian mathematician Johann Radon, who used it to approach various problems in differential geometry. Its name has nothing to do with the chemical element radon - and I've inwardly seethed when I've heard colleagues mispronounce it. Some have even speculated about whether the chemical element relates to why the name was chosen - because we'll eventually see that the Radon transformation has many uses in nuclear medicine. The Radon transform is defined with this rather intimidating looking integral, which will make even an applied mathematician like me reach for a notebook and pencil and start sketching the path of integration and trying to figure out just what is being computed here. What it's actually doing, though, is pretty familiar by now. If we define a line by the angle theta and the distance to the origin rho, then it's integrating the value of the function on that line. If we approximate the integral by summing over a discrete grid of pixels, and the pixel values are all ones and zeroes, this transformation and the Hough transform are the same thing except for perhaps some constant factors. But because the Radon transform views both the coordinates and the pixel values in the continuous domain, it has many important uses that the Hough transform cannot model. In addition, the discrete Radon transformation has a fast computation - reducing the number of operations from big-O of N cubed to big-O of N squared log N. Even more important, the discrete Radon transformation has an exact inverse that can be computed fast - big-O of N squared polylog N - but that's a topic for another day. All of this material is covered in a 2006 paper by numerical-analysis rockstar William Press, who's one of the authors of the "Numerical Recipes" textbook. There's a link in the video description. The basic ideas were discovered by several researchers independently in the late 1990's, but as far as I know, Press's paper was the first to tie them all together. It deserves to be more widely known than it appears to be.

Why is discrete Radon slow?

To understand how to approach the performance problems, let's have a look at what makes the discrete Radon transform so slow when we do it in a naïve way. Look at the pixels that contribute to the integral along one line in the image. Now let's choose another nearby line - that is, another line with nearly the same angle and displacement. The process of calculating the sums has had the effect of calculating the total of over half the pixels repeatedly, so it's done half again the work necessary. If we consider more nearby lines, we'll see that there are many groups of pixels that have been summed over and over. The key to making the algorithm fast will be to carry out the repeated calculations only once.

Fast DRT - Setup

We need to work out a way to do that. In order to set up the problem, let's make a change to how we model the lines. Rather than using a form defined by angle and displacement, we'll use slope-and-intercept form. Here, I use u and v for the horizontal and vertical pixel coordinates. A little simple trigonometry converts one form to the other. We want to avoid singularities in secant theta, so let's for now consider only the slopes between 0 and 1, or the angles between zero and a quarter pi. We'll see later how to generalize back to the entire range of angles. We'll divide that range of slopes into N equal steps, ranging from 0 to 1 inclusive of the endpoints. Here's a grid overlaid upon an actual image The slopes of lines range from 0 to 1, in N equal steps. The array index in the output array is the numerator of a fraction defining the slope, with N-1 in the denominator. Lines with slope 0 are horizontal. Their intercept can range from 0 to N-1. As the slope gets steeper, the maximum intercept remains N-1, but the minimum goes more and more negative. At a slope of 1, the intercept ranges between plus and minus N-1. Given the slope, it's simple to compute the range of intercepts. This model gives us a packing for quantized angles and intercepts in a 2N by N array. One corner of the array will be zero, because it corresponds to lines that do not intersect the image.

Fast DRT - Divide and conquer

Now we come to the key aspect that makes the algorithm fast. We take a divide-and-conquer approach. Suppose that we've already computed the solution for two images, each having k columns and N rows. In a single pass over the pixels, we can compute the solution for the combined image that has 2k columns and N rows. Choose a slope m and intercept b in the left image, and let p be the column number corresponding to the slope. That choice defines a set of pixels in the left half of the combined image. We are going to draw a line with the same slope in the right-hand image, shifted downward. We can draw it in two places. One of the two places has been shifted down by p pixels relative to the left-hand line The other line is shifted down by p+1 pixels. It also has the same slope. With two addition operations between elements of the left- and right-hand arrays, we can assemble discretized lines for the starting intercept, and two values of the slope. Repeating this operation for every slope-intercept pair in the input will fill the output array completely. And this gives us our N2 log N solution. We start with the input array, and columns of width 1. The sums are trivial - the pixel value is the same thing as the integral along a 1-pixel horizontal line through the pixel. We run the doubling repeatedly, combining into strips of 2 columns, 4 columns, and so on, until we have a single strip covering the entire width of the image. This gives us log (to the base 2) of N passes, each of which executes its inner loop N2 times. And that's the fundamental algorithm. But wait!

Fast DRT - Error analysis

I put one over on you! Did you see what I did there? If we draw the actual best discrete approximation to a line through the image, here shown in green crosses, and then look at the line that was produced by our computation, you'll see that they're different at a few spots. Our computation has drawn the line a pixel too high or too low. These errors, in fact, continue to accumulate with each pass of the algorithm. Fortunately, a numerical analyst much smarter or more patient than me computed the maximum error, way back in 1965. The calculation fills a 32-page technical report, so I'll present only the conclusion: for an N by N image, the error is bounded by 1/6 log-to-the-base 2 N plus big O of 1. For the 4k by 4k images that I used on the punched cards, this is a deviation of less than three pixels This was only one of many sources of uncertainty in the imaging. I had to contend with distortion of the image, noise in the photodetector, pixel misclassification, and warpage of the physical punched card, among other factors. Fortunately, they all added up to only a handful of pixels, and the solution was simply to be tolerant of that much uncertainty in the results. The blurring of the edges that we saw in the image preparation step is now explained. I was simply spreading out the likelihood function to cover all these uncertainties.

Fast DRT - Symmetry

We're almost done! We're still considering only slopes between 0 and 1. Rather than trying to expand the algorithm to cover other slopes, there's a simpler approach: use symmetry, and manipulate the input array. To get slopes from 1 to infinity, or angles from a quarter pi to half pi, all that we have to do is transpose the input image. We can get negative slopes up to a quarter pi by flipping the image vertically. And we can get negative slopes from a quarter pi to a half pi by doing both, which is equivalent to rotating the image by a quarter turn. Each of these four orientations gives us another 2N by N array. Some simple trigonometry, which I won't trouble to derive, gives us the mapping between the cells of these arrays and the angles and displacements. I actually made the picture that I called the Hough transform in the previous video by sampling the arrays according to that mapping. It's interesting that although the textbooks all seem to show both the Hough and Radon transforms in the angle-and-displacement form, Hough's original patent also used slope and intercept. There's language in the patent that shows that Hough was well aware of the symmetries, which he applied by simply rotating the subject in front of the camera.

Fast DRT - Topology

If all we need is a likelihood estimate of individual lines, read out in either slope-intercept or angle-displacement form, we're done. Recall, though, that I mentioned needing to find local maxima. In order to find those, we have to make sure that a pixel right on the edge of one of the four arrays doesn't have a higher neighbour in one of the adjacent arrays. To a mathematician, that's a question of topology. Although the popular literature often describes topology as "geometry on a rubber sheet", it's more accurate to say that topology is the study of what points are connected to what other points. Label the array corresponding to the upright image as A, and the one corresponding to the transposed images as B. A slope of 1 in the array A corresponds to a slope of 1 in the array B, but the intercepts change sign. The rightmost column of A therefore corresponds to the rightmost column of B, rotated about the long axis. Now let C be the vertically flipped image. A slope of zero in C corresponds to a slope of zero in A, but an intercept of 0 in one corresponds to an intercept of N-1 in the other. Therefore, the leftmost column of C corresponds to the leftmost column of A, again rotated about the long axis. The rotated image is D. The diagonal of the rotated image is the same as the diagonal of the vertically flipped image, again with a change of sign on the intercept. We can pair D off against C along the long axes. Now we have the full range of angles. It's natural to guess that the pattern of flipping at each junction will continue, with B and D adjacent. But that's not how it works. Look at the angle represented by the column as we go through the spliced images. As we go through the first image, it advances from minus half pi to minus a quarter pi With each array, it spins anther quarter pi, to zero to a quarter pi to a half pi At a half pi, the direction that the column represents is the same direction in terms of the equations of lines, but its sense is reversed. Flipping the array has also reversed the intercept, so the intercepts wind up having the same sign. We have to reverse the vertical axis to line things up, mirroring the array. And of course, the bookkeeping for adjacencies at the other end reverses correspondingly. This means that the arrays make a half twist when they join. This sort of connection is what a topologist means by a Möbius strip. While it's fun to model a Möbius strip using a twisted strip of paper, the real use of the Möbius strip is in mathematical structures like the one we're analyzing here. Here, it gives us the necessary insight to work out which elements need to be compared to find local maxima. There's more to how I found local maxima, but I'll save that for another video! As an aside, I was once working on an image analysis problem that involved this structure, and wound up using a representation where the intercepts would wrap around - an intercept that went beyond the border of an array would come in from the other side. (The more mathematical will probably recognize that I was doing a Fourier transform along that axis.) In working out the adjacencies, I realized that what I was dealing with was the topological figure called a Klein bottle. I annoyed a whole bunch of colleagues because I just had to share that it was the first time I'd spotted a wild Klein bottle in its natural habitat! While I saw topological spheres and projective planes all the time, actually seeing a Klein bottle other than as a curiosity in a topology class or popular literature was a new one for me! While I wasn't quite as maniacal about it as Cliff Stoll (there's a link to his video somewhere nearby) I surely came close to it. My colleagues uniformly replied, "yes, Kevin, that's nice." and obviously had no idea that a wild Klein bottle might be an exciting find. I'm sure they're glad I've retired. Go ahead and hit the "dislike" button now! See if I care!

Conclusions

Klein bottles are cool! With that, I'm pretty much done with the machinery for scanning the punched cards. I've still, though, not covered even half the great material in the paper by William Press. The rest, though, needs a little more background. Before coming back to that - which I plan to do - I have some other ideas that I want to cover.

Coming attractions

I'd like to look a little more into the history of these techniques. I got fascinated by how early they emerged. That led me down a rabbit hole with passages that led to 1950s- and 1960s-era particle physics, to very early exploration in video games - borrowing techniques from state-of-the-art radar of the day, and to the birth of numerically-controlled machine tools, upon which the modern factory depends. I want to explore some of the uses of the Radon transformation, from how it can help drivers stay in the lane, to creating architectural drawings from sets of photographs. I want to discuss the watershed algorithm, which models the prominence of mountain peaks, but also turns up in many image processing algorithms, such as separating the areas of an image corresponding to things like parts on a conveyor. But, as I said, most of all, I want to come back and show how the Radon transform is at the foundation of a lot of today's medical imaging technology. So there's a lot of material to come! Until next time, stay safe, stay healthy, and keep calculating!