function [] = Unblurring(frequency_magnitude_threshold, sigma)
% Read the images to correct in
psf = im2double(imread('stars-psf.png'));
blurred = im2double(imread('stars-blurred.png'));
% Taper the image at the edges to try and reduce ringing
blurred = my_edgetaper(blurred, psf);
% Recover FFTs for the blurred image and its psf
psf_fft = fft2(psf, size(blurred, 1), size(blurred, 2));
blurred_fft = fft2(blurred);
% Thresholding FFT coefficients
frequency_coefficients_to_keep = abs(psf_fft) > frequency_magnitude_threshold;
blurred_fft = frequency_coefficients_to_keep .* blurred_fft;
% Weiner filter on FFT coefficients
sample = im2double(imread('sample-stars.png'));
sample = sample(:, :, 1); % Discard color information
sample_power_spectrum = powerspectrum(sample, size(blurred, 1), size(blurred, 2));
%noise_power_spectrum = powerspectrum(randn(size(sample)), size(blurred, 1), size(blurred, 2)) / 10;
noise_power_spectrum = 2 * sigma ^ 2; % Theoretically spectral density = 2 * sigma ^ 2
weiner = sample_power_spectrum ./ (sample_power_spectrum + noise_power_spectrum);
blurred_fft = weiner .* blurred_fft;
% Finally, unblur the image and display it
unblurred = ifft2(blurred_fft ./ psf_fft);
imshow(real(unblurred) * 255)
end
function [signal_spectrum] = powerspectrum(signal, rows, columns)
signal_fft = fft2(signal, rows, columns);
signal_spectrum = abs(signal_fft) .^ 2;
end
function osf = my_psf2osf(psf, rows, cols)
% This is my reimplementation of the function with the same name from the
% image processing toolbox. It turns a PSF into a frequency domain OSF by
% shifting the PSF so its center is at 1,1 then doing the FFT
padded_psf = psf;
padded_psf(rows, cols) = 0;
shifted_padded_psf = circshift(padded_psf, -floor(size(psf)/2));
osf = fft2(shifted_padded_psf);
end
function output = hamming_window(size)
output = repmat(0, [size 1]);
for j = 1:size,
output(j) = 0.54 - 0.46 * cos(2 * pi * (j / (size - 1)));
end
end
function output = triangle_window(size)
half_size = floor(size / 2);
output = repmat(0, [size 1]);
for j = 1:half_size,
output(j) = (j - 1) / half_size;
end
for j = (half_size + 1):size,
output(j) = (size - j) / half_size;
end
end
function output = flat_window(size)
output = repmat(0, [size 1]);
end
function tapered = my_edgetaper(image, psf)
% This is my reimplementation of the function with the same name from the
% image processing toolbox. The idea is to make the input image look more
% periodic in the time domain (to reduce ringing) by interpolating
% between the image and a blurred version of it. The blurred version
% dominates at the edges and the original everywhere else.
% First we compute the OSF so that we can work in the frequency domain
% and without worrying about the canvas edges
osf = my_psf2osf(psf, size(image, 1), size(image, 2));
% Blur the image by the PSF by using the OSF in the fourier domain
blurred_image = ifft2(fft2(image) .* osf) / 255; % Why a constant factor here??
imshow(blurred_image)
% Build up a weightings matrix
weighting = repmat(1, size(image));
[psf_height, psf_width] = size(psf);
width_window = hamming_window(psf_width);
for j = 1:(psf_width / 2),
weighting(:, j) = weighting(:, j) * width_window(j);
weighting(:, size(weighting, 2) - j + 1) = weighting(:, size(weighting, 2) - j + 1) * width_window(j);
end
height_window = hamming_window(psf_height);
for j = 1:(psf_height / 2),
weighting(j, :) = weighting(j, :) * height_window(j);
weighting(size(weighting, 1) - j + 1, :) = weighting(size(weighting, 1) - j + 1, :) * height_window(j);
end
% Apply the weightings matrix to get the final image
tapered = (image .* weighting) + (blurred_image .* (1 - weighting));
end