pifft

 PIFFT - Homodyne reconstruciton of partial MRI k-space (Spatial Frequency) data

  IMG = PIFFT(KSP)
  The input KSP is assumed to be a rectangular matrix with the number of rows
  equaling the number of phase encodes (lines) of the image's k-space.  The 
  first rows of the matrix should be the highest frequency phase encodes 
  acquired so that the center of k-space is near the bottom of the matrix.
  The image is reconstructed for the next highest power of 2 greater than the
  number of phase encodes in the partial k-space data set.

  Example:

  load mri
  img = double( squeeze(D(:,:,1,16)) );
  ksp = fftshift( fft2( img ) );
  pksp = ksp(1:68,:);
  m=pifft(pksp);

  figure;
  subplot(1,3,1); imagesc(img); axis image; title('128 of 128 of k-space lines');
  subplot(1,3,2); imagesc(m); axis image; title('68 of 128 of k-space lines');
  subplot(1,3,3); imagesc(abs(img-m)); axis image; title('abs. difference');
  colormap(gray);

  Reference:

  Noll DC, Nishimura DG, Macovski A. Homodyne Detection in Magnetic 
  Resonance Imaging.  IEEE Trans. on Medical Imaging 1991; 10(2):154-163

0001 function [m] = pifft(pksp)
0002 
0003 % PIFFT - Homodyne reconstruciton of partial MRI k-space (Spatial Frequency) data
0004 %
0005 %  IMG = PIFFT(KSP)
0006 %  The input KSP is assumed to be a rectangular matrix with the number of rows
0007 %  equaling the number of phase encodes (lines) of the image's k-space.  The
0008 %  first rows of the matrix should be the highest frequency phase encodes
0009 %  acquired so that the center of k-space is near the bottom of the matrix.
0010 %  The image is reconstructed for the next highest power of 2 greater than the
0011 %  number of phase encodes in the partial k-space data set.
0012 %
0013 %  Example:
0014 %
0015 %  load mri
0016 %  img = double( squeeze(D(:,:,1,16)) );
0017 %  ksp = fftshift( fft2( img ) );
0018 %  pksp = ksp(1:68,:);
0019 %  m=pifft(pksp);
0020 %
0021 %  figure;
0022 %  subplot(1,3,1); imagesc(img); axis image; title('128 of 128 of k-space lines');
0023 %  subplot(1,3,2); imagesc(m); axis image; title('68 of 128 of k-space lines');
0024 %  subplot(1,3,3); imagesc(abs(img-m)); axis image; title('abs. difference');
0025 %  colormap(gray);
0026 %
0027 %  Reference:
0028 %
0029 %  Noll DC, Nishimura DG, Macovski A. Homodyne Detection in Magnetic
0030 %  Resonance Imaging.  IEEE Trans. on Medical Imaging 1991; 10(2):154-163
0031 %
0032 %
0033 
0034 %
0035 %    Written by Edward Brian Welch (edwardbrianwelch@yahoo.com)
0036 %    MRI Research Lab, Mayo Graduate School, November 2001
0037 %
0038 
0039 % ERROR CHECKING
0040 if ndims(pksp)~=2 | ~isnumeric(pksp),
0041     error('PFFT operates on two-dimensional numerical data');
0042 end
0043 
0044 % Convert data into hybrid space by inverse FFT along FE dir.
0045 % for better numerical accuracy.
0046 pksp = ifft(pksp,[],2);
0047 
0048 % Frequency Encode direction is assumed to be the longer direction
0049 [Npp Nf] = size(pksp);
0050 
0051 % Number of phase encodes in the full k-space
0052 Np = 2^nextpow2(Npp+1);
0053 
0054 % Number of high frequency phase encodes
0055 NH = Np-Npp;
0056 
0057 % NUmber of low frequency phase encodes
0058 NL = Npp - NH;
0059 
0060 % Row indices of high and low frequency lines
0061 HFind = 1:NH;
0062 LFind = (NH+1):(NH+NL);
0063 
0064 % Create weighting window
0065 %
0066 % It multiplies all un-partnered high frequency lines by 2.0
0067 % Low frequency lines are multiplied by a ramp that approaches
0068 % zero at the zero-padded edge of k-space.  The weighting
0069 % factors of two partnered low frequency lines should have a
0070 % sum of 2.0
0071 w = zeros(Np,1);
0072 w(HFind) = 2;
0073 rstep = 2/(NL+1);
0074 w(LFind) = [(2-rstep):-rstep:rstep];
0075 
0076 % Create weighted partial k-space
0077 HFksp = zeros(Np,Nf);
0078 HFksp(1:Npp,:)=pksp;
0079 HFksp = HFksp.*repmat(w,1,Nf);
0080 
0081 % Create Low Frequency k-space
0082 LFksp = zeros(Np,Nf);
0083 LFksp(LFind,:) = pksp(LFind,:);
0084 
0085 % Low frequency image
0086 Rc = ifft(LFksp,[],1);
0087 
0088 % Unsynchronous image
0089 Ic = ifft(HFksp,[],1);
0090 
0091 % Synchronous image
0092 Is = Ic.*exp(-i*angle(Rc));
0093 
0094 % Demodulated image
0095 m = real(Is);
0096 
0097 return