function imRc = myLambdaIRadon(RT, thetas, n, bool)
%MYLAMBDAIRADON Lambda Tomography
%   The idea is to recover \Lambda f, with \Lambda the root of negative
%   laplacian.
%   INPUT ARGS: RT     -  Radon data to invert for
%               thetas - angles corresponding to the projections
%               n      - size of the reconstructed square image
%               bool   - if switched to 1, computes the second s
%               derivative before backprojecting the data. 
%   OUTPUT ARGS: imRc  - reconstructed image
%   See also MYRADON, ROTATEIMAGE

%   Author: Francois Monard
%   Date: 06/24/2011

%--- compute the second derivative of each projection
if bool == 1,
    h = 2/n;
    RT = ([RT(end,:); RT(1:end-1,:)] + [RT(2:end,:); RT(1,:)] - 2*RT)/h^2;
end

imRc = zeros(n); %--- allocate memory for the image

% define the reconstruction grid
x = -1+1/n:2/n:1-1/n;
x = repmat(x, n, 1);

y = (1-1/n:-2/n:-1+1/n).';
y = repmat(y, 1, n);

costhetas = cos(thetas); sinthetas = sin(thetas);

dtheta = thetas(2)-thetas(1);

nr = size(RT,1); 

%--- Backprojection  - vectorized in (x,y), looping over theta
% count = 1;
for i=1:length(thetas)
    proj = RT(:,i);
    t = x.*costhetas(i) + y.*sinthetas(i);
    t = (t-(-nr/n+1/n))*n/2+1;
    a = floor(t);
    
    temp = (t-a).*proj(a+1) + (a+1-t).*proj(a);
    
%     figure(100);
%     if count<=16 && floor(count)==count,
%         subplot(4,4,count); dispmat(temp);
%     end
%     count = count + 0.25;
    imRc = imRc + temp; 
end

imRc = (dtheta/2)*imRc;

%--- take out the part outside the unit disk
imRc = (n/2)*imRc.*double(x.^2+y.^2<=1);

end
