[MatLab] Filter a tomographic image in Fourier’s space

MatLab has built-in functions to simulate acquisition and elaboration of tomographic (as in CT) images. When it comes to filtering the image prior to back-projecting it, it is possible to do it yourself without relying on the (good as in Shepp-Logan) filters MatLab has.

 

We will filter our image in Fourier’s space. For each column of the sinogram, we:

• apply the Fourier transform
• filter it
• apply the Fourier antitransform

When we reconstruct our image WITHOUT having MatLab apply any filter, we’ll see the image, filtered with our filter, as a result.

%filter an image in Fourier's space

%clear matlab environment
clear all;
close all;
theta = [0:179]; %projection angle 180 degrees
F = phantom(256); %create new test phantom 256x256 pixels aka alien
figure, imshow(F,[]),title('Alien vs Predator'); %show alien
R = radon(F,theta); %apply Radon aka scan the alien to create a sinogram
figure, imshow(R,[]),title('Sinoalien'); %show sinogram
%get Shepp-Logan filter, you can use any filter you want
Fi = phantom(128);
Ri = radon(Fi,theta);
[I, Filter]=iradon(Ri,theta, 'Shepp-Logan');
%add image noise to the sinogram
maximum=max(max(R));
minimum=min(min(R));
R=(R-minimum)/(maximum-minimum);
%normalize between 0 and 1
fact=1001;%set the number of X-rays, the higher the better (and deadlier)
R=(fact/10^12)*R;
Noised=imnoise(R,'poisson');%add Poissonian noise
figure, imshow(Noised,[]),title('Dirty alien');%show noisy sinogram
%add zero-padding
Padded(1:512, 1:180)=0;
Padded(1:367, :) = Noised;
figure, imshow(Padded,[]),title('Padded alien');%show padded sinogram
%filter the noise with our filter in Fourier's space for each column of the sinogram
for k=1:180
	%FPadded = filtered and padded
    FPadded(:, k)=fft(Padded(:, k));%Fourier transform
    FPadded(:, k)=FPadded(:, k).*Filter;%Filter in Fourier
    Padded(:, k)=real(ifft(FPadded(:, k)));%Fourier antitransform of the filtered and still padded sinogram column
end
%remove padding
R=Padded(1:367,:);
figure, imshow(R,[]),title('Unpadded filtered alien');%show final sinogram
I=iradon(R, theta, 'None'); %back project without filtering to show final result
figure, imshow(I,[]),title('Final alien');