%% This script shows Filtering of Sinogram and backprojection explicitly.
N = 256;
SheppLogan = false; % true;
if SheppLogan
    im = phantom(N);
else
    im = zeros(N);
    im = im + MakeCircle(192,192,32,N);
    im(32:131,32:131) = 1;

end
theta3 =  [0:1:179];
%theta3 =  [0:1:59];
%theta3 =  [0:15:179];

g3 = radon(im,theta3);
reconUBP3 = iradon(g3,theta3,'linear','none',1,256);
[np,ntht] = size(g3);


%% add noise

%g3 = g3*(1+0.45*randn(size(g3)));

%% 
figure(1); clf;
subplot(2,3,1); imagesc(im);title('Original Image');
subplot(2,3,4); imagesc(g3); title('Radon Transform');

%% filter the sinogram
fg3 = fft(g3,512,1);
ak = [ 0:255, [256:-1:1]]; % absolute k in shifted domain.
aks = real(ifft(ak));
ffg3 = fg3.*repmat(ak',1,ntht);
g3_filt = ifft(ffg3,[],1);
g3_filt = g3_filt(1:np,:);


figure(2);clf;
subplot(2,3,1); imagesc(log(abs(fftshift(fg3))));title('log Fourier of SinoGram');
subplot(2,3,2); plot(fftshift(ak)); title('Abs(k) in Fourier')
subplot(2,3,3); imagesc(log(abs(fftshift(ffg3))));title('log Filtered Fourier of SinoGram');
subplot(2,3,5); plot(ifftshift(aks)); title('Abs(k) in Spatial domain')

%%
figure(1);  
subplot(2,3,5); imagesc(g3_filt); title('filtered Radon Transform');

%%
reconFBP3 = iradon(g3_filt,theta3,'linear','none',1,256);
figure(1); 
subplot(2,3,2); imagesc(reconUBP3); title('unfiltered Back projection');
subplot(2,3,3); imagesc(reconFBP3); title('filtered Back projection');

%%
fimff = fft2(reconFBP3);
fimuf = fft2(reconUBP3);
figure(2);
subplot(2,3,4); imagesc(log(abs(fftshift(fimuf))));title('log Fourier Unfiltered BackProjection');
subplot(2,3,6); imagesc(log(abs(fftshift(fimff))));title('log Fourier Filtered BackProjection');
%%
function im = MakeCircle(x0,y0,R,N)
im = zeros(N);
for i = 1:N
    for j = 1:N
        if( (i-x0)^2 + (j - y0)^2 < R^2)
            im(i,j) = 1;
        end
    end
end
end