title: Image deblurring using gaussian filter in matlab without additive noise
body: <p>I have to use an inverse filter to remove the blurring from this image</p>
  
  <p><a href="https://i.stack.imgur.com/tCTel.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/tCTel.jpg" alt="Blurred Image"></a>. </p>
  
  <p>Unfortunately, I have to figure out the transfer function <code>H</code> of the imaging
  system used to get these sharper images, It should be Gaussian. So, I should determine the approximate width of the Gaussian by trying different Gaussian widths in an inverse filter and judging which resulting images look the “best”. </p>
  
  <p>The best result will be optimally sharp – i.e., edges will look sharp but will not have visible ringing.</p>
  
  <p>I tried by using 3 approaches: </p>
  
  <ol>
  <li>I created a transfer function with <code>N</code> dimensions (odd number, for simplicity), by creating a grid of <code>N</code> dimensions, and then applying the Gaussian function to this grid. After that, we add zeroes to this transfer function in order to get the same size as the original image. However, after applying the filter to the original image, I just see noise (too many artifacts). </li>
  <li>I created the transfer function with size as high as the original image, by creating a grid of the same size as the original image. If <code>sigma</code> is too small, then the <code>PSF FFT</code> magnitude is wide. Otherwise it gets thinner. If <code>sigma</code> is small, then the image is even more blurred, but if we set a very high <code>sigma</code> value then we get the same image (not better at all).</li>
  <li>I used the <code>fspecial</code> function, playing with sizes of <code>sigma</code> and <code>h</code>. But still I do not get anything sharper than the original blurred image.  </li>
  </ol>
  
  <p>Any ideas?</p>
  
  <p>Here is the code used for creating the transfer function in Approach 1:</p>
  
  <pre><code>%Create Gaussian Filter
  function h = transfer_function(N, sigma, I) %N is the dimension of the kernel
  %create a 2D-grid that is the same size as the Gaussian filter matrix
  grid = -floor(N/2) : floor(N/2);
  [x, y] = meshgrid(grid, grid);
  arg = -(x.*x + y.*y)/(2*sigma*sigma);
  h = exp(arg); %gaussian 2D-function
  kernel = h/sum(h(:)); %Normalize so that total weight equals 1
  
  [rows,cols] = size(I);
  add_zeros_w = (rows - N)/2;
  add_zeros_h = (cols - N)/2;
  
  h = padarray(kernel,[add_zeros_w  add_zeros_h],0,'both'); % h = kernel_final_matrix
  
  end 
  </code></pre>
  
  <p>And this is the code for every approach: </p>
  
  <pre><code>I = imread('lena_blur.jpg');
  I1 = rgb2gray(I);
  figure(1),
  I1 = double(I1);
  %---------------Approach 1
  % N = 5; %Dimension Assume is an odd number
  % sigma = 20; %The bigger number, the thinner the PSF in FREQ
  % H = transfer_function(N, sigma, I1);
  %I1=I1(2:end,2:end); %To simplify operations
  imagesc(I1); colormap('gray'); title('Original Blurred Image')
  
  I_fft = fftshift(fft2(I1)); %Shift the image in Fourier domain to let its DC part in the center of the image
  
  
  
  % %FILTER-----------Approach 2---------------
  % N = 5; %Dimension Assume is an odd number
  % sigma = 20; %The bigger number, the thinner the PSF in FREQ
  % 
  % 
  % [x,y] = meshgrid(-size(I,2)/2:size(I,2)/2-1, -size(I,1)/2:size(I,1)/2-1);
  % H = exp(-(x.^2+y.^2)*sigma/2);
  % %// Normalize so that total area (sum of all weights) is 1
  % H = H /sum(H(:));
  % 
  % %Avoid zero freqs
  % for i = 1:size(I,2) %Cols
  %     for j = 1:size(I,1) %Rows
  %         if (H(i,j) == 0)
  %             H(i,j) = 1e-8;
  %         end
  %     end
  % end
  % 
  % [rows columns z] = size(I);
  % G_filter_fft = fft2(H,rows,columns);
  %FILTER---------------------------------
  
  
  %Filter--------- Aproach 3------------
  N = 21; %Dimension Assume is an odd number
  sigma = 1.25; %The bigger number, the thinner the PSF in FREQ
  
  H = fspecial('gaussian',N,sigma)
  [rows columns z] = size(I);
  G_filter_fft = fft2(H,rows,columns);
  
  %Filter--------- Aproach 3------------
  
  
  
  %DISPLAY FFT PSF MAGNITUDE
  figure(2),
  imshow(fftshift(abs(G_filter_fft)),[]); title('FFT PSF magnitude 2D');
  
  
  % Yest = Y_blurred/Gaussian_Filter
  I_restoration_fft = I_fft./G_filter_fft;
  I_restoration = (ifft2(I_restoration_fft));
  I_restoration = abs(I_restoration);
  
  
  
  
  I_fft = abs(I_fft);
  
  % Display of Frequency domain (To compare with the slides) 
  figure(3),
  subplot(1,3,1); 
  imagesc(I_fft);colormap('gray');title('|DFT Blurred Image|')
  subplot(1,3,2)
  imshow(log(fftshift(abs(G_filter_fft))+1),[]) ;title('| Log DFT Point Spread Function + 1|');
  subplot(1,3,3)
  imagesc(abs(I_restoration_fft));colormap('gray'); title('|DFT Deblurred|')
  % imshow(log(I_restoration+1),[])
  
  %Display PSF FFT in 3D
  figure(4)
  hf_abs = abs(G_filter_fft);
  %270x270
  surf([-134:135]/135,[-134:135]/135,fftshift(hf_abs));
  % surf([-134:134]/134,[-134:134]/134,fftshift(hf_abs));
  shading interp, camlight, colormap jet
  xlabel('PSF FFT magnitude')
  
  %Display Result (it should be the de-blurred image)
  figure(5),
  %imshow(fftshift(I_restoration));
  imagesc(I_restoration);colormap('gray'); title('Deblurred Image')
  
  %Pseudo Inverse restoration
  % cam_pinv = real(ifft2((abs(G_filter_fft) &gt; 0.1).*I_fft./G_filter_fft));
  % imshow(fftshift(cam_pinv));
  % xlabel('pseudo-inverse restoration')
  </code></pre>