N = 256;
M = 50;
d = 10

for sigma = 5:5:25   
    
    sigma
    
    h = compute_psf(N, sigma);
   
    f = zeros(N, N);
    f(N/2 - d - M : N/2 - d, :) = 200;
    f(N/2 + d : N/2 + d + M, :) = 200;

    figure(1);
    colormap(gray(256));
    image(f);
    
    ffile = sprintf('f_%d.jpg', sigma);
    imwrite(uint8(f), ffile);
    
    g = filter2(h, f);
    g = 200*g/max(max(g));

    figure(3);
    colormap(gray(256));
    image(g);  
    
    gfile = sprintf('g_%d.jpg', sigma); 
    imwrite(uint8(g), gfile);

    % now let's look at a single line profile as a 1D plot
    figure(4);
    plot(g(:, N/2));
    hold on
    plot([N/2 - d - M, N/2 - d - M], [0, 200], 'r');
    plot([N/2 + d + M, N/2 + d + M], [0, 200], 'r');
    plot([N/2 - d, N/2 - d], [0, 200], 'r');
    plot([N/2 + d, N/2 + d], [0, 200], 'r');
    axis([0 300 0 250])
    hold off
    
    pfile = sprintf('p_%d.jpg', sigma); 
    print('-djpeg', '-r72', pfile);

end