Shepp-Logan filtering (spatial domain)

Hey there,

I got a school assignment where I have to re-construct a CT image, using filtered backprojection.

It is specified that the filtering has to be done in the spatial domain (probably because it's more challenging). I have succeeded (I think?) creating the filter in the frequency domain, but attempting to convert it to the spatial domain, fails miserably. Here is the code I got:

% Create Shepp-Logan filter in frequency domain
frqFit  = linspace(-2*pi*B, 2*pi*B, 256);
frqFilt = abs(frqFilt);
frqFilt = frqFilt.* (sin(frqFilt/(4*B)) ./ (frqFilt/(4*B)) ) ;
filtWin = hanning(length(frqFilt));
frqFilt = win'.*frqFilt;

% Convert filter to spatial domain
sptFilt = ifftshift(frqFilt,1);
sptFilt = ifft(sptFilt);
sptFilt = sptFilt(1:round(end/2));
sptFilt = real(sptFilt);

%Use filter
result = conv(test_line, sptFilt);   % test_line is one of the projections

This is not the full code (since that is kinda messy), but it should be the important part. What the duck am I doing wrong?

Answer 1

[Matt J]

I doubt you're meant to be using FFTs to build the filter. There are direct formulas for the ramp filter in the space domain and more than likely, you're meant to be using those.

Apart from that, the following line will not create a frequency axis that includes DC (frqFit=0) as one of its samples.

frqFit  = linspace(-2*pi*B, 2*pi*B, 256);

[Jakob Sørensen]

That is actually a good point. Is that really what is causing me that much problems? Also, as I understood the professor, it should be possible to create the filter in the frequency domain, and then convert it to the spatial domain...

[Matt J]

One can only guess what's causing you problems. You haven't actually described the problems you're experiencing.

If you create the filter in the frequency domain and then ifft to the spatial domain, you will introduce aliasing and other discrete sampling errors into the filter. (Remember, the ifft is only a discrete approximation to the continuous IFT).

If your prof. told you to perform the filtering in the non-Fourier domain, it's hard to imagine why he wouldn't want you to construct the filter there as well. There is nothing to be learned purely from performing the convolution in the space domain. It produces the same result, but less efficiently, than fft-based convolution.