5

I have an assignment to implement a Ram-Lak filter, but nearly no information given on it (except look at fft, ifft, fftshift, ifftshift).

I have a sinogram that I have to filter via Ram-Lak. Also the number of projections is given.

I try to use the filter 

                      1/4              if I == 0

(b^2)/(2*pi^2)  *     0                if I even

                      -1/(pi^2 * I^2)  if I odd

b seems to be the cut-off frequency, I has something to do with the sampling rate?

Also it is said that the convolution of two functions is a simple multiplication in Fourier space.

I do not understand how to implement the filter at all, especially with no b given, not told what I is and no idea how to apply this to the sinogram, I hope someone can help me here. I spent 2hrs googling and trying to understand what is needed to do here, but I could not understand how to implement it.

Dennis Jaheruddin
  • 21,208
  • 8
  • 66
  • 122
mmlac
  • 1,091
  • 1
  • 9
  • 24
  • Your question is very unclear. What is I? How do we find b? None of the texts I've looked through use the same notation. Ram Lak filter is just an abs(f) function, like a double ramp. If you tell me what those variables are, I'll be able to help you. – Phonon Jul 15 '11 at 16:08
  • My problem is I do not know. In literature (i.e. Algorithms for Reconstruction with Nondiffracting Sources Page 72 - [link](www.umiacs.umd.edu/~mingyliu/enee631/CTI_Ch03.pdf) they use just a I (or teta there, the sampling interval) Can you help me to implement a simple abs() filter based on the sinogram and the number of projections? – mmlac Jul 15 '11 at 16:09
  • Try this link. http://laskin.mis.hiroshima-u.ac.jp/Kougi/07a/PIP/PIP12pr.pdf It's a good text on it. Figure 2B is a Ram Lak filter – Phonon Jul 15 '11 at 16:10

1 Answers1

11

The formula you listed is an intermediate result if you wanted to do an inverse Radon transform without filtering in the Fourier domain. An alternative is to do the entire filtered back projection algorithm using convolution in the spatial domain, which might have been faster 40 years ago; you would eventually rederive the formula you posted. However, I wouldn't recommended it now, especially not for your first reconstruction; you should really understand the Hilbert transform first.

Anyway, here's some Matlab code which does the obligatory Shepp-Logan phantom filtered back projection reconstruction. I show how you can do your own filtering with the Ram-Lak filter. If I was really motivated, I would replace radon/iradon with some interp2 commands and summations.

phantomData=phantom();

N=size(phantomData,1);

theta = 0:179;
N_theta = length(theta);
[R,xp] = radon(phantomData,theta);

% make a Ram-Lak filter. it's just abs(f).
N1 = length(xp);
freqs=linspace(-1, 1, N1).';
myFilter = abs( freqs );
myFilter = repmat(myFilter, [1 N_theta]);

% do my own FT domain filtering
ft_R = fftshift(fft(R,[],1),1);
filteredProj = ft_R .* myFilter;
filteredProj = ifftshift(filteredProj,1);
ift_R = real(ifft(filteredProj,[],1));

% tell matlab to do inverse FBP without a filter
I1 = iradon(ift_R, theta, 'linear', 'none', 1.0, N);

subplot(1,3,1);imagesc( real(I1) ); title('Manual filtering')
colormap(gray(256)); axis image; axis off

% for comparison, ask matlab to use their Ram-Lak filter implementation
I2 = iradon(R, theta, 'linear', 'Ram-Lak', 1.0, N);

subplot(1,3,2);imagesc( real(I2) ); title('Matlab filtering')
colormap(gray(256)); axis image; axis off

% for fun, redo the filtering wrong on purpose
% exclude high frequencies to create a low-resolution reconstruction
myFilter( myFilter > 0.1 ) = 0;
ift_R = real(ifft(ifftshift(ft_R .* myFilter,1),[],1));
I3 = iradon(ift_R, theta, 'linear', 'none', 1.0, N);

subplot(1,3,3);imagesc( real(I3) ); title('Low resolution filtering')
colormap(gray(256)); axis image; axis off

Demonstration of manual filtering

user244795
  • 698
  • 6
  • 14