MATLAB Answers

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

8 views (last 30 days)
Gianni Schena
Gianni Schena on 7 Nov 2016
Edited: Matt J on 7 Nov 2016
%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)])

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!