title: 1D image deblurring using line spread function and Wiener Filter
body: <p>I have been reading up on using the wiener filter to deblur or "restore" a 3D image. The difference between my approach and that used in the MATLAB function "deconvwnr" is that I am not interested in deblurring the image in all 3D dimensions (x,y,z) only in the z-dimension. My theory is that if I know the line spread function (LSF) in the z-dimension (LSF_z) then I can use the wiener filter and the LSF_z to deblur the image in the z-dimension.</p>
  
  <h2>My stab at the code is as follows:</h2>
  
  <pre><code>close all
  % Define LSF (i.e. gaussian with mean = 0.5 &amp; FWHM = 0.2)
      space=linspace(0,1,50)';
      space=fftshift(space);
      FWHM=.2;
      sigma=FWHM/2.3548;
      c=0.5;
      h = gaussmf(space,[sigma c]); % LSF 
  % Define degradation function frequency cutoff
      Hcut=1e-7; % 
      Flpcut=1e-7;
  % Define Transfer Function
      H=fftshift(fft(h)');
  
      [nx ny nz]=size(I);
      Ideconv=zeros(size(I)); % allocate memory for "true" image
  
      eps=1e6; % regularization parameter (i.e. constant noise to signal ratio)
  
  % Loop over all voxels in 2D space
  for j=1:ny
    for i=1:nx
      % Flp: FFT of line profile in z-dimension 
          lp=squeeze(I(i,j,:));
          Flp=fftshift(fft(lp));
  
      % cut-off all frequences below Flpcut
      Flp=(Flp&gt;Flpcut).*Flp; 
  
  % Define power spectrum of degradation function and cut-off all
  % frequences below Hcut
      magsquared=(H&gt;Hcut).*abs(H).^2; 
  
  % Inverse FFT of (Wiener Filter x Flp) 
      Ideconv_1D_z=ifft(ifftshift((H2&gt;Hcut.*magsquared./(H.*(magsquared+eps))).*Flp)); 
  
  % Place "deblurred" line profile into true image matricx
      Ideconv(i,j,:)=Ideconv_1D_z; 
  
    end
  end
  </code></pre>
  
  <p>Data sample:</p>
  
  <pre><code> a=squeeze(I(250,500,:));
   a=[70
      64
      64
      57
      57
      54
      58
      59
      55
      55
      51
      52
      51
      52
      55
      57
      58
      66
      64
      66
      69
      68
      65
      66
      66
      62
      56
      55
      54
      54
      55
      58
      63
      69
      65
      52
      52
      53
      55
      58
      58
      58
      59
      60
      60
      56
      56
      55
      55
      57
      60
      58
      61
      61];
  </code></pre>
  
  <p>Pseudo-code for plotting results:</p>
  
  <p>NOTE: For the sake of simplicity,this is only plotting the results for a single line profile (i=400 and j=600) </p>
  
  <pre><code>subplot(2,2,1); 
       plot(real(Flp),'b')
       plot(imag(Flp),'r')
       title('FT of line profile')
       legend('real(Flp)','imag(Flp')
  subplot(2,2,2); 
       plot(real(H),'b')
       plot(imag(H),'r')
       title('FT of LSF')
       legend('real(H)','imag(H)')
  subplot(2,2,3)
       plot(magsquared,'b');
       title('Power Spectrum: magsquared')
  subplot(2,2,4)
       plot(lp,'r')
       plot(I,'b');
       legend('original','deconvolved')
       title(sprintf('Line Profile results with FWHM=%1.2f',FWHM))
  </code></pre>
  
  <p>I believe my issue is with indexing and/or the frequency cut-off. When I plot the results they are not what I am looking for. Specifically, they are not nearly as dependent on the FWHM of the LSF as they should be.</p>
  
  <p>Furthermore, I am surprised that the FT of a gaussian has an imaginary component since a gaussian is a real and even fxn which should produce a real and even FT.</p>
  
  <p><img src="https://i.stack.imgur.com/eDpnZ.jpg" alt="Output plots"></p>