% In: g, 513-by-360: matrix of g(s,theta)
load siemens.mat
[q,p] = size(g);
figure(1); clf;
imagesc(g); colormap(1-hot)

%% fan2para
thetas = 2*pi*(1:360)/360;  % in radians
thetas_deg = 1:360;         % in degrees
a = asin(1/2.87);
new_g = zeros(q,p);
k = -180:179;
for s = 1:q
    dtheta = -((1+q)/2 - s)/((q-1)/2) * a;
    new_g(s,:) = real(ifft(fftshift(exp(-1i*k*dtheta)).*fft(g(s,:))));
end

figure(2); clf
imagesc(new_g); colormap(1-hot)

%% Filtering    
ghat = fftshift(fft(new_g));
k = -(q-1)/2:(q-1)/2;
ghat = ghat .* (abs(k)'*ones(1,p));
new_g = real(ifft(ifftshift(ghat)));

%% Backprojection
N = 200;
L = 2;
h = L/N;
xx = -L/2:h:(L/2-h);
[XX,YY]=ndgrid(xx,xx);
I = zeros(N);
s = linspace(-1,1,513);
for j=1:N
    for t=1:p 
        x = [XX(j,:); YY(j,:)];
        theta = [cos(thetas(t)); sin(thetas(t))];
        new_s = x'*theta;
        G = new_g(:,t);
        I(j,:) = I(j,:) + interp1(s,G,new_s,'linear')';
    end
end

% if (0),
% for j = 1:N
%     for k = 1:N
%         for t = 1:p
%             x = [XX(j,k), YY(j,k)];
%             theta = [cos(2*pi*t/360); sin(2*pi*t/360)];
%             new_s = x*theta;
%             G = new_g(:,t);
%             I(j,k) = I(j,k) + interp1(s,G,new_s);
%         end
%     end
% end
% end

figure(3); clf
imagesc(I); colormap(hot)