function [img,j,d] = sharp(i, sigma_r, sigma_d)

% transform to double if necessary (im2double will determine the necessity)
i = im2double(i);

c_r = 1/(2*sigma_r^2);
c_d = 1/(2*sigma_d^2);
hwin = max(1, ceil(2.1*sigma_d));   % half window size for domain filter
squares = ([-hwin:hwin].^2)'; gaussian = exp(-c_d*squares);
gaussian_d = gaussian*gaussian';

j = zeros(size(i)); d = zeros(size(i));

[rows, cols] = size(i);
for r = (-hwin):hwin, 
  for c = (-hwin):hwin,
    
    g = gaussian_d(r+hwin+1,c+hwin+1);

    s = zeros(size(i));
    is = zeros(size(i));
    
    rs = 1 + max(0, r);
    re = rows + min(0, r);
    cs = 1 + max(0, c);
    ce = cols + min(0, c);
    
    is(rs:re, cs:ce) = i((rs:re)-r, (cs:ce)-c);
    s(rs:re, cs:ce) = exp(-c_r*(i(rs:re, cs:ce) - is(rs:re, cs:ce)).^2) * g;
  
    d = d + s;
    j = j + s.*is;
    
  end 
end

j = j ./ d;

img=j+d;
figure;imagesc(img);colormap(gray);

