a = imread('Shepp_logan.png',  'png');
b = a(:,:,2);

[R,xp] = radon(b,[0:20:180]);
%figure(1)
%clf
%plot(xp,R(:,1)); title('R_{0^o} (x\prime)')
%figure(2)
%clf
%plot(xp,R(:,2)); title('R_{20^o} (x\prime)')

theta = 0:20:180; %0:180;
%I2 = iradon(R,theta);

%[p,theta,filter,d,interp,N] = parse_inputs(varargin{:});

p = R;
theta = pi*theta/180;

% Default values
N = 2*floor( size(p,1)/(2*sqrt(2)) );  % This doesn't always jive with RADON
                        %N = 0;                 % Size of the reconstructed image
d = 1;                 % Defaults to no cropping of filters frequency response
filter = 'ram-lak';    % The ramp filter is the default
interp = 'linear';     % default interpolation is linear
%string_args = {'nearest neighbor', 'linear', 'spline', ...
%      'ram-lak','shepp-logan','cosine','hamming', 'hann'};


% Design the filter
len=size(p,1);   
H = designFilt(filter, len, d);
p(length(H),1)=0;  % Zero pad projections 

% In the code below, I continuously reuse the array p so as to
% save memory.  This makes it harder to read, but the comments
% explain what is going on.

p = fft(p);    % p holds fft of projections

for i = 1:size(p,2)
   p(:,i) = p(:,i).*H; % frequency domain filtering
end

p = real(ifft(p));     % p is the filtered projections
p(len+1:end,:) = [];   % Truncate the filtered projections
img = zeros(N);        % Allocate memory for the image.

% Define the x & y axes for the reconstructed image so that the origin
% (center) is in the spot which RADON would choose.
xax = (1:N)-ceil(N/2);
x = repmat(xax, N, 1);    % x coordinates, the y coordinates are rot90(x)
y = rot90(x);

costheta = cos(theta);
sintheta = sin(theta);
ctrIdx = ceil(len/2);     % index of the center of the projections

% Zero pad the projections to size 1+2*ceil(N/sqrt(2)) if this
% quantity is greater than the length of the projections
imgDiag = 2*ceil(N/sqrt(2))+1;  % largest distance through image.
if size(p,1) < imgDiag 
   rz = imgDiag - size(p,1);  % how many rows of zeros
   p = [zeros(ceil(rz/2),size(p,2)); p; zeros(floor(rz/2),size(p,2))];
   ctrIdx = ctrIdx+ceil(rz/2);
end

% Backprojection - vectorized in (x,y), looping over theta
if strcmp(interp, 'nearest neighbor')
   for i=1:length(theta)   
      proj = p(:,i);
      t = round(x*costheta(i) + y*sintheta(i));
      img = img + proj(t+ctrIdx);
   end
elseif strcmp(interp, 'linear')
   for i=1:length(theta)  
      proj = p(:,i);
      t = x.*costheta(i) + y.*sintheta(i); 
      a = floor(t);  
      img = img + (t-a).*proj(a+1+ctrIdx) + (a+1-t).*proj(a+ctrIdx);
   end
elseif strcmp(interp, 'spline')
   for i=1:length(theta)
      proj = p(:,i);
      taxis = (1:size(p,1)) - ctrIdx;
      t = x.*costheta(i) + y.*sintheta(i);
      projContrib = interp1(taxis,proj,t(:),'*spline');
      img = img + reshape(projContrib,N,N);
   end
end

img = img*pi/(2*length(theta));

figure(2)
clf
imagesc(img)
colormap(gray)
title('ram-lak, 20^o')



