% Try out ft.m on lightning image. This time with realistic shot noise and readout noise.

test_setup

readout_noise = 5; % electrons (= counts, photons) RMS


figure(1)
hold off
imagesc(image.^contrast)
colormap(gray)
title('Incident intensity (expectation value)')
print('1-incident-intensity.pdf','-dpdfwrite')

% Now produce and display the smeared version, with noise.
smeared = ft(image, T, f);
smeared = poissrnd(smeared) + readout_noise*randn(Ny,Nx); % Add noise.

figure(2)
hold off
imagesc(real(smeared.^contrast))
colormap(gray)
title('Smeared CCD readout (with noise)')
print('2-smeared-ft-ccd.pdf','-dpdfwrite')

% Attempt to invert, and reconstruct the original image
reconstructed = invert_ft(smeared, T, f);
max_diff_reconstructed_smeared = max(max(abs(reconstructed-smeared)))
image_std = std(image(:))
reconstructed_error = reconstructed - image*T;
reconstructed_max_error = max(max(abs(reconstructed_error)))
reconstructed_std_error = std(reconstructed_error(:))
reconstructed_SNR = image_std/reconstructed_std_error

figure(3)
hold off
imagesc(real(reconstructed.^contrast))
colormap(gray)
title('Reconstructed image')
print('3-reconstructed.pdf','-dpdfwrite')

% How would this compare to a shuttered exposure?
shuttered = poissrnd(image*T) + readout_noise*randn(Ny,Nx); % Add noise.
shuttered_error = shuttered - image*T;
shuttered_std_error = std(shuttered_error(:))
shuttered_SNR = image_std/shuttered_std_error

figure(4)
hold off
imagesc(real(shuttered.^contrast))
colormap(gray)
title('Shuttered full frame CCD readout')
print('4-fullframe-shuttered-ccd.pdf','-dpdfwrite')

figure(5)
hold off
plot(std(reconstructed_error,0,2),'g', 'linewidth', 2)
hold on
plot(std(shuttered_error,0,2),'b')
xlabel('Image row')
ylabel('Standard deviation')
legend('Reconstructed from frame transfer CCD','Shuttered CCD')
print('5-compare-noise-fullframe.pdf','-dpdfwrite')



