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