%%
% Optical Imaging and Spectroscopy
%
% David J. Brady
% Duke University
% www.opticalimaging.org
%
%
% Convolution backprojection example, figure 2.36
%
%
xMax=2;
nP=256;
xStep=2*xMax/(nP-1);
uMax=1/(2*xStep);uStep=2*uMax/(nP-1);
urange=linspace(-uMax,uMax,nP);
xrange=linspace(-xMax,xMax,nP);
thetaRange=0:.5:179.5;
[x,y]=meshgrid(xrange,xrange);
figure(1);
set(gcf,'color','white'); colormap gray;
%I=[zeros(64,256);zeros(128,64),phantom(128),zeros(128,64);zeros(64,256)];
I=f3(x,y);
subplot(3,2,1);imagesc(xrange,xrange,I);title('Object'); axis square; xlabel('x');ylabel('y');
rTf1=radon(I,thetaRange);
uP=size(rTf1,1);
xRadonRange=linspace(-uP*xMax/nP,uP*xMax/nP,uP);
uRadonRange=linspace(-uP*uMax/nP,uP*uMax/nP,uP);
subplot(3,2,2);imagesc(thetaRange,xRadonRange,rTf1);title('h(l,\theta )');axis square; 
xlabel('\theta');ylabel('l');
ftR=fftshift(fft(fftshift(rTf1,1)),1);
subplot(3,2,3);imagesc(thetaRange,uRadonRange,abs(ftR));title('|g(u_l,\theta )|'); axis square;
xlabel('\theta');ylabel('u_l');
subplot(3,2,4);imagesc(urange,urange,abs(fftshift(fft2(I))));title('|f(u,v)|');axis square;
xlabel('u');ylabel('v');
uP=size(ftR,1);
Q=kron(abs((-(uP-1)*uStep/2):uStep:((uP-1)*uStep/2))',ones(1,360)).*ftR;
Q=fftshift(ifft(fftshift(Q,1)),1);
subplot(3,2,5);imagesc(thetaRange,xRadonRange,abs(Q));title('|Q(l,\theta)|');axis square;
ylabel('l');xlabel('\theta');
lRange=linspace(-uP*xMax/nP,uP*xMax/nP,uP);
fEst=zeros(nP);
for pop=0:.5:179.5
    fEst=fEst+interp1(lRange,Q(:,2*pop+1),x.*cos(pop*pi/180)+y.*sin(pop*pi/180));
end
subplot(3,2,6);imagesc(xrange,xrange,abs(fEst));title('Reconstructed object');axis square;
xlabel('x');ylabel('y');