% reblur.m    Created by:  Dion Monstavicius
%
% This function takes an input image, and a blur map
% of the same size as the input image.  A blur map
% is a matrix where each entry indicates the std. of
% the gaussian filter to be applied to that pixel in
% the input image.
% 
% This function returns a blurred image the same size
% as the inputs.
%
% This function corrects for boundary conditions by
% reflecting the image at the boundaries before
% performing the isotropic gaussian filter.
% This eliminates the "image darkening" at the edges.

function [output1] = reblur(image, blur)

stepsize = .3;
quantizedblur = stepsize*round(blur/stepsize);

lowerlimit = min(min(quantizedblur));
upperlimit = max(max(quantizedblur));

w = 50;
x = -w:w;

[A B] = size(image);
%tempimage = zeros(A+2*w+1, B+2*w+1);
tempimage(w+1:A+w, w+1:B+w) = image;

% this part reflects the image at the boundary so we don't loose energy!
for i = 1:w,
 tempimage(w-i+1,w+1:B+w) = image(i,:);
 tempimage(w+1:A+w,w-i+1) = image(:, i);
 tempimage(A+w+i,w+1:w+B) = image(A-i,:);
 tempimage(w+1:w+A, B+w+i) = image(:, B-i);
end
inputimage = image;
blur2 = zeros(A+2*w, B+2*w);
blur2(w+1:A+w, w+1:B+w) = quantizedblur;
quantizedblur = blur2;
image = tempimage;

[A B] = size(image);

output1 = zeros(A,B);

for sigma = lowerlimit:stepsize:upperlimit,
 if sigma >= w/2,
  break;
 end
 if sigma ~= 0;
  sigma
  G = exp(-x.^2/(2*sigma^2))/(2*pi*sigma^2);
  G = G/sum(G);
  filteredimage = conv2(G,G,image);
  filteredimage = filteredimage(w+1:A+w, w+1:B+w);
  [C D] = find(quantizedblur == sigma);
  map = find(quantizedblur == sigma);
  for i = 1:length(C)
   output1(C(i),D(i)) = filteredimage(C(i),D(i));
  end
 else
  sigma
  [C D] =  find(quantizedblur == sigma);
  for i = 1:length(C)
   output1(C(i),D(i)) = image(C(i),D(i));
  end
 end
end

temp = output1;
clear output1;
[A B] = size(inputimage);

output1 = temp(w+1:w+A, w+1: w+B);