function plotlamest(lamest, M)
%Usage: plotlamest(lamest, M)
%Given output "lamest" from the routine backproject.m, plot the output
%image as a grayscale image, on (2*M+1)^2 resolution grid.

[m,n] = size(lamest);

%Storage for actual image values on grid
lamp = zeros(2*M+1);

%Loop over grid points, get lambda estimate at each
for j = 0:2*M
    x = (j-M)/M;
    for k = 0:2*M
        y = (k-M)/M;
        r = sqrt(x.^2 + y.^2);
        if r<1.0
            alpha = atan2(y,x);
            if y < 0
                alpha = alpha + 2*pi;
            end
            tk = floor(alpha*m/2/pi)+1;
            rk = ceil(n*r);
            if r==0 %Pesky case
                rk = 1;
            end
            
            %Estimate lambda by using tk, rk. Crazy indexing is a pain
            lamp(j+1,2*M-k+1) = lamest(tk,rk);
        else
            lamp(j+1,2*M-k+1) = 0.0;
        end
    end
end

miny = min(min(lamp));
maxy = max(max(lamp));

figure('Position',[200 200 650 650]);
graymap(255);
%colormap('default')

X = (-M:M)/M;
image(X,X,255-255*(lamp'-miny)/(maxy-miny));
axis off;
title('CT Reconstruction');