I went to the library looking for the The Radon Transform and some of its applications book, but it wasn't on the shelf. The librarian (who had eyed me suspiciously and pointed at the sign that says "Summer students must have their advisors check out books for them" when I walked in) took up my query and was ired to discover that the book's last borrower had left the Institute without returning his library books. There was one book, 75 Years of the Radon Transform or something like that that contained an amusing article of the The Radon Transform and Me nature, in which the author explained that, as a physicist in Capetown, he had been asked to serve as a radiologist at the understaffed hospital, and there he came up with a rudementary yet prescient CAT-scan scheme out of curiosity.
Asmahan and I went over Ed's script for particle-picking using cross-correlation. The script works by taking an image containing particles, and a smaller reference image of a single particle, and looks for locations in the big image that match the smaller image. It does this by padding the reference image (increasing its dimensions to that of the larger image) and then looking for peaks in the cross-correlation between these two images. It does this once for 45 rotations of the reference image. Finally, the script does some nearest-neighbor elimination.
I wrote a C program (spi_conv.c) that reads in a SPIDER image file and dumps it as space-deliminated ASCII. This was pretty easy; the only trickery is that SPIDER uses floating point numbers for everything (even those things that are integral, such as the image dimensions in pixels) and big-endian byte order (the files were created on SGI IRIX). My program works and I was able to load its output into MATLAB with dlmread without problem. I think it will be quite easy to reimplement Ed's script in Matlab in a very concise and readable way.
Ed's scheme is well suited to parallelization along "embarassingly parallel" lines, namely that each rotation of the reference image is an independent task. That provides a clean separation into approx. 45 equally difficult yet independent tasks. It would be fun to play with parallelizing the linear algebra / fourier transform tasks, but that's probably not necessary. Maybe there's legitimate work to be done in parallelizing the Radon transform. The annoying part about reimplementing the algorithm in C is that library functions for image rotation with interpolation (etc) might not be available/convenient.
I'd like to learn more about correlation functions; right now I just take it for granted that they work. The correlation function we're using is simply convolution, so it can be done easily by transforming to fourier space, multiplying elementwise, and then taking the inverse transform. (Although, unless there's some other stage that's easy to do in the fourier domain, this might be more work than simply doing a coordinate-space convolution.)
Ate lunch at Santa Maria (?) with Asmahan and her friend, and Ilija (the yugoslavian) showed up too.
The reason it is better to do convolution in fourier space is that, in one dimension, direct convolution is a Θ(n2) operation, while the fast fourier transform is only Θ(n log n). Once in fourier space, convolution is simply an elementwise multiplication, Θ(n). Thus it's asymptotically 'cheaper' to go into fourier space, do the multiplcation, and take the inverse transform, than it is to just do the convolution directly. (c.f. CLR ยง30).
At 16:30 we had our meeting ("party") with Greta at Charlie's. The other students who have arrived so far and were at the meeting were: Jool (Joel?) from Canada, Adam from USA, Timur from Russia, Toni from Germany, Ira from India, Dogan from Turkey, Marija from Yugoslavia/Serbia, Taylan from Turkey, Arielle from Canada, Yaakov from USA, me, and Ilija from Yugoslavia/Serbia. It seems a little like we all came in pairs from our respective countries. Apparently we've all received the same directions from our parents; officially nearly all of us are restricted to Rehovot or even sequestered to the Institute campus.
I went into the lab to work on getting what will now be called the pattern matching by cross-correlation algorithm working. It turns out that it really was working all along, but the results looked wrong because I was viewing the absolute value of the results instead of the real components. This should be fine, because the imaginary components are nearly zero, but our pixel values are negative, so this results in an inverted image. In any case, it works now. Now I have to work on a scheme to extract the particle locations from the cross-correlation image.
There are several possibilities for optimization. The small reference image should be rotated before the padding is added, otherwise much time will be wasted on interpolating values known to be zero. Also, I suspect that rotation and fourier transform operations may commute, which could be useful if the padding operation could be done in fourier space as well. Actually, I wouldn't be surprised if the padding operation commuted with the fourier transform as well. On the implementation side, we could use the fftw fourier transform library which is multithreaded, to take advantage of the multiple CPUs available to us.
I figured out exactly what the Radon Transform does today, and implemented it in Matlab. My MATLAB implementation is based on a description I found on a web page for a DSP course somewhere; the essential suggestion was to rotate the image itself (of which we are taking the radon transform) rather than trying to reposition the projection surface. It's the sort of thing that seems completely obvious once you know it, and it makes the implementation completely trivial: use imrotate to rotate the image into position, and then use sum to take the sum along the first dimension. That's all! Weighted back-projection (which is reconstruction in real-space, as opposed to fourier-space reconstruction via the central-line theorem) can also be implemented in this way. Nifty! I'm still not sure how the Radon Transform is useful to us; it seems to me that it simply describes the CAT-scan-like process that's actually done by the microscope.
I modified my pattern-matching-by-cross-correlation program to return not all of the cross-correlation plots, but simply (1) the maximum of all of the cross-correlation plots for each pixel and (2) the rotation angle of the reference image which resulted in that maximum. Some simple thresholding and multiplication results in a nice picture that shows where the helicies lie and their detected orientation. However, I still haven't implemented anything to actually extract the position of helical segments. Along each helix there is a nice chain of little peaks in the maximum-cross-correlation plot (resulting from the periodicity with which the reference should match the helix, I think), and I think there should be a good way to lock onto these little peaks and return them.
Yes, rotation does commute with the Fourier transform. The essential part of the proof is that the Fourier transform includes the inner product <x,k>. If we apply a rotation matrix R to the coordinate position vector x, then we can move it over to the wave-vector as R-1k. That's really the most important part. Unfortunately it doesn't work in our sampled space, so this elegant property is useless to us. Padding with zeros does some curious things I don't fully understand yet.
Shabat
Copyright © 2002 by Tobin Fricke