%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function v = psf(sigma,x,y)
% returns the value of the point spread function
% We're using a Gaussian PSF (truncated far away from
% its center -- the origin).

if x >= -8 & x <= 8 & y >= -8 & y <= 8
   v = exp(-sigma*(x.^2+y.^2));
else
   v = 0;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function g = applyGaussianPSF(sigma, f, x)
% Returns the blurred graph of f
% if x = 0 the mesh plot is generated,
% otherwise it is not graphed.
% sigma determines how much the light is spread
% (How bad the blur is.)

% Get size
[m n] = size(f);

% create Toeplitz blocks
for i = -8:8
   for k = 1:n
      c(k) = psf(sigma, i, k-1);
      r(k) = psf(sigma, i, 1-k);
   end;
   
   block(i+9,1:n,1:n) = toeplitz(c,r);
end;

% apply the convolution product row by row
for i = 1:m
   g(i,1:n) = zeros(1,n);

   % apply convolution product
   % to the graph on row i
   for j = -8:8
      if i-j >= 1 & i-j <= m
         b(:,:) = block(j+9,:,:);
         g(i,1:n) = g(i,1:n) + f(i-j,1:n) * b;
      end;
   end;
end;

if x == 0, mesh(g); end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function n = inner_product(x,y)
% turn x and y into columns by transposing each row
% into a column and concatenating together.
% then find regular inner product
% (x and y are sz x ?? matrices)

sz = size(x,1);

n = 0;
for i = 1:sz
   n = n + x(i,:) * y(i,:)' ;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [g, error, niter, flag] = CGmethod(sigma, g, f, maxit, tol)
% Conjugate Gradient method for solving linear systems.
%
% sigma is a parameter for the Gaussian psf
%
% We are solving "Ag = f":
% g as input is the initial guess at a solution and then
% g as output is the solution to the system.
% f is the blurred image.
%
% maxit is the maximum number of iterations to attempt
% niter is the number of iterations actually preformed
% error is the difference between the exact answer and our answer
% tol is the allowable error
% flag is set to 1 if the function returns before error <= tol

flag = 0;
niter = 0;

% r is the residual (error).
r = f - applyGaussianPSF(sigma,g,1);

% Initial search direction.
d = r;

% norm of f (as a vector: turn rows into cols and concatenate)
norm_f = sqrt(inner_product(f,f));

if(norm_f == 0.0), norm_f = 1.0; end;

error = sqrt(inner_product(r,r)) / norm_f;

for niter = 1:maxit
   if(error < tol), return, end; % bail out if the error is small enough.
   
   % multiply our search direction by the blur matrix.
   Ad = applyGaussianPSF(sigma,d,1);
   
   % the square of the blur matrix norm of d. 
   inp_Ad_d = inner_product(Ad,d);
   
   % Step size (how far to move along the gradient direction).
   alpha = inner_product(r,d) / inp_Ad_d;
   
   % Improve the 
   g = g + alpha * d;
   
   % update the residual (error).
   r = f - applyGaussianPSF(sigma,g,1);
   
   % the conjugate factor to force our search directions
   % to be orthogonal (in some sense).
   beta = -inner_product(r,Ad) / inp_Ad_d;
   
   % New search direction.
   d = r + beta * d;
   
   % measure the relative error.
   error = sqrt(inner_product(r,r)) / norm_f;
   
end;

if(error > tol) flag = 1; end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function B = blur(pic, gry, flag)
% Blurs "pic" (a jpeg) and returns the blurred image.
% The image is converted to grayscale if gry = 1 and 
% the image is displayed if flag = 1.

% Load up the picture. A is a m by n matrix if grayscale
% and an m by n by 3 tensor if in color (think 3 matrices
% stacked on top of each other -- one each for RGB values.
A = loadImage(pic, gry, 0);

if gry == 0  % Color
    % Blur Red, Greens, and then Blues.
    B(:,:,1) = applyGaussianPSF(0.1, A(:,:,1), 1);
    B(:,:,2) = applyGaussianPSF(0.1, A(:,:,2), 1);
    B(:,:,3) = applyGaussianPSF(0.1, A(:,:,3), 1);
else % Grayscale
    B = applyGaussianPSF(0.1, A, 1);
end;

if flag == 1  % Display the image
    figure; 
    if gry == 1
        colormap(gray(256)); 
        
        % rescale to make image show up.
        x = max([max(max(B)),0.000001]);
        scaleB = (1/x)*B;
    else 
        % rescale to make image show up.
        x = max([max(max(max(B))),0.000001]);
        scaleB = (1/x)*B;
    end;
    %image(myImg);
    imshow(scaleB);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = deblur(B, n, gry, flag)
% Deblurs the image stored in the matrix/tensor B 
% and returns the "unblurred" image.
% The image is considered to be grayscale if gry = 1 
% and the image is displayed if flag = 1.


% Our initial "guess" to feed into the CG method is the zero image.
A = zeros(size(B));

if gry == 0 % Color

    % Run the CG method on RGB values. 
    % Each call to the CG method is iterated n times.
    [A(:,:,1), error, niter, flag] = CGmethod(0.1, A(:,:,1), B(:,:,1), n, 0);
    [A(:,:,2), error, niter, flag] = CGmethod(0.1, A(:,:,2), B(:,:,2), n, 0);
    [A(:,:,3), error, niter, flag] = CGmethod(0.1, A(:,:,3), B(:,:,3), n, 0);
    
else % Grayscale

    % Each call to the CG method is iterated n times.
    [A, error, niter, flag] = CGmethod(0.1, A, B, n, 0);

end;


% Display image when flag is 1.
if flag == 1
    figure; 
    if gry == 1
        colormap(gray(256)); 
        
        % rescale to make image show up.
        x = max([max(max(A)),0.000001]);
        scaleA = (1/x)*A;
    else 
        % rescale to make image show up.
        x = max([max(max(max(A))),0.000001]);
        scaleA = (1/x)*A;
    end;
    %image(myImg);
    imshow(scaleA);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function myImg = loadImage(pic, gry, flag)
% This function loads the file "pic" (should be a jpeg).
% If gry = 0, it's loaded as a color picture.
% If gry = 1, it's turned into a grayscale image.
% flag determines whether the image is displayed.
% flag = 0 not shown / flag = 1 shown

myImg = imread(pic);

if gry == 1 
    tmp = .2989*myImg(:,:,1) + .5870*myImg(:,:,2) + .1140*myImg(:,:,3);
    myImg = tmp;
end;

if flag == 1
    figure; 
    if gry == 1
            colormap(gray(256)); 
    end;
    %image(myImg);
    imshow(myImg);
end;

% convert to double so the image can be manipulated.
myImg = double(myImg);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f = cityScape(x)
% Returns the graph of the skyline
% if x = 0 then mesh plot is generated,
% otherwise it is not graphed.

f = zeros(160,100);
%f = zeros(64);
f(21:44,20:30) = 1;
f(16:49,40:45) = 2;
f(11:55,50:55) = 3;
f(30:40,50:55) = 2.2;

if x == 0
    mesh(f); 
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Z = GaussianPSF(sigma,x)
% Returns the n by n graph of the point spread function
% if x = 0 the mesh plot is generated, 
% otherwise, it is not graphed.
% sigma determines how much the light is spread.

[X,Y] = meshgrid(-8:8,-8:8);
Z = exp(-sigma*(X.^2+Y.^2));

if x == 0, mesh(-8:8,-8:8,Z); end;