% demonstrates the effect of a single projection and backprojection
N = 256;
Shapes = true;
if Shapes
    
    im = zeros(N);
    im = im + MakeCircle(192,192,40,N);
    im(32:131,32:131) = 1;
%    im(64:113,64:113) = 4;
%    im(148:197,32:131)  = 2;
%    im(32:231,148:247)  = 0.5;
    
else
    impath = '../../StandardTestImages/';
    ln = imread([impath 'Lena512.png'],'png'); ln = double(imresize(ln,[N N]));
    hs = imread([impath 'house.png'],'png'); hs = double(imresize(hs,[N N]));
    pt = imread([impath 'Cameraman256.png'],'png'); pt = double(imresize(pt,[N N]));
    ot = imread([impath 'ot.pgm']); ot = double(imresize(ot,[N,N]));
    sl = phantom(N);
    im = ot;
end

figure(1);clf;
subplot(2,3,1); imagesc(im);title('Original Image');
%%
t1 = 135; %;[0:179];
ang = [t1,t1];

g = radon(im,ang);
subplot(2,3,2); plot(g); title('Projection')
bg = iradon(g,ang,'linear','none');
figure(1); subplot(2,3,3);imagesc(bg);title('BackProjection');

%% filter 
ng = size(g,1);
fg = fft(g(:,1));
kmax = floor(ng/2);
ak = [ 0:kmax-1, [ng-kmax:-1:1]]; % absolute k in shifted domain.
% apodise the high frequencies
gk = exp( - ak.^2/(2*(ng/8)^2));
ak = ak.*gk;

ffg = fg.*ak';
g_filt = real(ifft(ffg,[],1));
%g_filt = g_filt(1:np,:);

%
ig = iradon(g,ang);
figure(1); subplot(2,3,5); plot(g_filt); title('Filtered Projection')
figure(1); subplot(2,3,6);imagesc(ig);title('Filtered BackProjection');

%% repeat as a movie over angle
igsum = zeros(size(ig));
for t1 = 0:179
    ang = [t1,t1];
    figure(2); clf;
    subplot(2,3,1); imagesc(im);title('Original Image');
    
    g = radon(im,ang);
    subplot(2,3,2); plot(g); title(['Projection \theta=',num2str(t1)])
    bg = iradon(g,ang,'linear','none');
    subplot(2,3,3);imagesc(bg);title('BackProjection');
    
    % filter
    ng = size(g,1);
    fg = fft(g(:,1));
    kmax = floor(ng/2);
    ak = [ 0:kmax-1, [ng-kmax:-1:1]]; % absolute k in shifted domain
    % apodise the high frequencies
    gk = exp( - ak.^2/(2*(ng/8)^2));
    ak = ak.*gk;
    
    ffg = fg.*ak';
    g_filt = real(ifft(ffg,[],1));
    %g_filt = g_filt(1:np,:);
    
    %
    ig = iradon(g,ang);
    igsum = ig + igsum;
    subplot(2,3,5); plot(g_filt); title('Filtered Projection')
    subplot(2,3,6);imagesc(ig);title('Filtered BackProjection');
    subplot(2,3,4); imagesc(igsum);title('Summed Filtered BackProjection');
    pause(0.1);
end
%%
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