central slice theorem experiments : comparing fft1 of the sinogram to the correspondent section of the fft2 of the phantom image

%This simple script was written to compare the fft1 of one column of the sinogram of one phantom to the 
%correspondent (according the Fourier central slice theorem) section of the fft2 (of the phantom image).
%The two transforms are overlapped but it comes out that the overlapping is not really satisfying. 
%I make %use of imrotate but at list for phi = 0 and phi =90 I cannot blame the interpolation of imrotate; the %overlapping should be perfect. 
%Any idea ? 3 alternative phantom are available in the first rows of the 
%testing script.
%%%%%%%%%%%%%%%%%%%%%%%%%CODE%%%%%%%%%%%%%%%%%%%%%%%%%
% 3 alternative phantoms
n=128; phan=phantom(n);
%phan=imread('cameraman.tif'); n=size(phan,1);
%n=20; phan=zeros(n,n); c=round(n/2); phan(c-5:c+5,c-2:c+3)=1;

phan=double(phan); phan=phan./max(phan(:));
phan=imresize(phan,[128 128]);

if (1) %phantom  image in a disk , 0 outside 
    disk=logical(0*phan); mask=disk;% + rand(n,n).*0.1;
    szd=size(disk);
    disk(round(szd(1)./2),round(szd(2)./2))=1;
    disk=bwdist(logical(disk));
    mask(disk<=szd(1)./2)=1;
    phan=phan.*double(mask);
   figure, imshow(phan,[]), title('starting phantom in a disk');
end

%phan=padarray(phan,[20,20],0,'both');
%%%%%%%%%%%%%%%%%%%%%
prj_angles=0:1:179;
proj_sino=radon(phan,prj_angles);
% phir changes angle phi only in case of different origin 
% of imrotate & radon 

phi= 0; % take one proj
% to extract the fft1 data from the fft2 you should counter-rotate 
phir=-phi; % rotation fft2(phantom) 

%%%%fft - iradon (phi=0 in the first column) 

proj_phi=proj_sino(:,phi+1); % proj( :,  phi degrees)
% trick of no help :  mean 2 close projections

if 0
    % only with delta amgle very small 
% proj_phi=0.5.*(proj_phi+proj_sino(:,min(phi+2,size(proj_sino,2)))); %
%
end

% dx is to crop to make the length od the sino equal to the side
% of the image , dx=0 id you do not like this option
dx=(( length(proj_sino(:,1))-size(phan,1) )); 

%dx=0; % !! do nothing 
dx_dx=round(dx./2);
dx_sx=dx-dx_dx; %nc=dx-dx_sx-dx_dx;
nx=length(proj_phi(dx_sx+1:end-dx_dx));
proj_phi=proj_phi(dx_dx+1:end-dx_sx); % !! crop
% filter the 'phi' projection taken from sinogram 
f_p=ifftshift(proj_phi(:));
f_p=fft(f_p,n); % fourier transf for radon
f_p=fftshift(f_p);

%%%%%%%%%%%%%%%%fft2 image :central slice 
ft2_phan=ifftshift(phan);
ft2_phan=fft2(ft2_phan,n,n); %fft2 phantom
ft2_phan=fftshift(ft2_phan);

interp='bicubic';
interp='bilinear';
%interp='nearest';

if(phir==0)
f_s_r=(real(ft2_phan));
f_s_i=(imag(ft2_phan));
f_s_m=(abs(ft2_phan));
end

if(phir==-90)
    k=-1;
f_s_r=rot90(real(ft2_phan),k);
f_s_i=rot90(imag(ft2_phan),k);
f_s_m=rot90(abs(ft2_phan),k);
end

if (phi~=0 & phi~=90)
   % phi
 % rotate 1 component of the complex numbers at time
%box='loose';
box='crop'; % use crop for phantom is in a disk
f_s_r=imrotate(real(ft2_phan),phir,box,interp);
f_s_i=imrotate(imag(ft2_phan),phir,box,interp);
f_s_m=imrotate(abs(ft2_phan), phir,box,interp);

end

f_s_c=f_s_r+1i.*f_s_i; %reconst complex fft2 rotated 

% central slice / central row 
nr=size(f_s_m,1);
% 2nd trick  : it helps
irw=ceil(nr./2); irw1=irw+1;
c_r_f_s_r=mean(f_s_r(irw:irw1,:),1); % real 
c_r_f_s_m=mean(f_s_m(irw:irw1,:),1); % modulo 
c_r_f_s_c=mean(f_s_c(irw:irw1,:),1); % complex

%%%%%%%%%%
axi_m=[0 300 0 +20];  % module
axi_r=[0 50 0 +10]; % real
% show fft iradon 
figure,
subplot(2,2,1)
plot(log(1+real(f_p)),'.'); grid on,
title([ 'real(fft(radon(Img))) , phi= ' num2str(phi)] )
axis(axi_r);

subplot(2,2,2)
%plot(abs(f_p));
plot(log(1+abs(f_p))); grid on
title( 'abs(fft(radon(Img)))')
axis(axi_m);

%%%%%one row of fft2 
subplot(2,2,3)
plot(log(1+real(c_r_f_s_r)),'.r'); grid on,
title('real( fft2(Img)) , central row')
axis(axi_r);

error=(log(1+abs(f_p(:)))-log(1+(c_r_f_s_m(:))));

subplot(2,2,4)
%plot(abs(c_r_f_s));
plot(log(1+(c_r_f_s_m)),'r'); grid on,
%plot(log(1+abs(c_r_f_s_c)),'r'); grid on,
hold on % overlaping
plot(log(1+abs(f_p))); grid on
title('abs(fft2(Img)) , central row')
axis(axi_m);

saer=sum(abs(error))/(n.^2);
%figure,plot(error); title(['error in magnitude ' num2str(saer)])