JE_main.m

% Description: 
%   - Simulates MRI data according to the model
%
%     s[m]=\sum_n f[n] e^{-(1/T2_star[n]+i*w[n])(m1*d1+m2*d2)} e^{-i(2*pi*m1*delta1*n1/N1 + 2*pi*m2*delta*n2)}
%     where d1,d2 - time sampling rate (sec)
%           Delta1, Delta2 - sampling intervals in k-space (should be 1)
%           f[n] - imaged object
%           w[n] - field map (rad/s)
%           T2_star[n] - relaxation T2* map (sec)
%
%   - Apply Harmonic Retrieval algorithm to jointly estimate the image, field map, and T2* map
%   - Display and compute the PSNR (in dB) of the estimates
%
% For further information, see Nguyen et al., Joint estimation and
% correction of geometric distortions for EPI functional MRI

% Hien M. Nguyen, ECE, University of Illinois at Urbana-Champaign, Aug 2008
%--------------------------------------------------------------------------
clear all;close all;

% Constants for simulation-------------------------------------------------
N1=64;
N2=N1/2;
M1=N1;
M2=2*N2;
d1=2.7e-6;
d2=5e-4;
delta1=1;
delta2=1;
L1=(N1-1)/2;
L2=(N2-1)/2;

% Load object - ground truth-----------------------------------------------
obj=load('mimb.mat');                             
obj=obj.mimb;
obj(find(obj==0))=10;       % for 0-background image

% Load field map ----------------------------------------------------------
we=load('we_b.mat');
we=we.we_b/(1.7);   

% Resize image and field map ----------------------------------------------
we=imresize(we,[N2 N1],'bilinear');
obj=imresize(obj,[N2 N1],'bilinear');
mask_im = abs(obj)>0.1*max(abs(obj(:)));       % to mask out the background

% Defining T2star map------------------------------------------------------
%Usually T2star is in the range from 10ms to 100ms; less than 10ms the
%signal will have a severe decay (2nd half of the data)
T2star=10*10^(-3)*ones(N2,N1);
R2star=1./T2star;

% Simulated 64x64 evenly spaced data --------------------------------------
[x1,obj_hat]=simdata_T2star(obj,we,N1,N2,M1,M2,d1,d2,delta1,delta2, T2star);   % T2* included
[x2,obj_hat]=simdata_noT2star(obj,we,N1,N2,M1,M2,d1,d2,delta1,delta2);         % no T2* included
figure, imagesc(abs(x1)),title('Simulated MRI data with field map and T2*'),colorbar;

% DFT reconstruction-------------------------------------------------------
a=ifft2(x2(1:N2,:));    % for a fair comparison, perform DFT on the data set without T2* effect

% Run Harmonic Restoration-------------------------------------------------
s_hat=ifft(x1,[],2);        
figure, imagesc(abs(s_hat)),title('Simulated data in half k-space domain'),colorbar;
est_we1=[];
est_image1=[];
lambda=exp(-d2*R2star(1,1));    % note: in the paper - alpha 
betta=1;                        % note: in the paper - lambda (due to history...)
                                % these are important parameters to be tuned
                                
for ii=1:M1 
    [z(:,ii),a_found,R2star_found(:,ii),z_noswap(:,ii),ind_z(1,ii)]= ...
    HR(s_hat(:,ii),N2,M2-1,d2,lambda,betta);
 
    u=angle(z(:,ii));
    
    n2=[0:N2-1];  
    a_found(n2+1)=a_found(n2+1)./(exp(j*pi*(L1*delta1*ii/N1+L2*delta2*n2(:)/N2)));    % Need this, but doesn't affect magnitude of the image
    est_image1=[est_image1 a_found];
    
    %Phase unwarping ------------------------------------------------------
    % can improve phase unwarping using other techniques...
    l=2:1:N2;
    delta_anglez=angle(z(l,ii))-angle(z(l-1,ii));
    g=delta_anglez./(2*pi)+repmat(delta2/N2,N2-1,1);
    delta_k=round(g);    
    tg=0;
    for l=2:N2
        k(l,1)=tg+delta_k(l-1,1);
        tg=k(l,1);
    end;
   
    wrapped_linear = angle(exp(-j * (2 * pi * delta2 * 1/N2*(0:1:N2-1).'-2*pi*k.*(0:1:N2-1).')));
    wrapped_linear =- (2 * pi * delta2 * 1/N2*(0:1:N2-1).'-2*pi*k);
    est_field_map=(-u+wrapped_linear)./d2;
    est_we1=[est_we1 est_field_map];

end;        % loop of ii

% Display outputs----------------------------------------------------------

figure, imagesc(interp2(obj)), title('True image'), colormap('gray');
figure, imagesc(interp2(abs(est_image1))), title('HR recovered image'),colormap('gray');
figure,imagesc(interp2(abs(a)),[0 11000]), title('DFT recovered image'),colormap('gray');
figure, imagesc(interp2(we./(2*pi)),[-75 10]),title('True field map'),colormap('gray'),colorbar;
figure, imagesc(interp2(est_we1.*mask_im./(2*pi)),[-75 10]),title('HR recovered field map'),colormap('gray'),colorbar;
figure, imagesc((1./R2star).*mask_im,[0.004 0.022]),colorbar,colormap('gray'),title('True R2star map');
figure, imagesc(interp2((1./R2star_found).*mask_im),[0.004 0.022]),colorbar,colormap('gray'),title('Recovered T2star map');

% Conventional method estimating field map---------------------------------
lr1 = fftshift(fft(fftshift(s_hat(1:M2/2,:),1),[],1),1);
lr2 = fftshift(fft(fftshift(s_hat(M2/2+1:M2,:),1),[],1),1);
FM1 = (angle(lr2) - angle(lr1));
FM2 = angle(exp(-j*angle(lr2))./exp(-j*angle(lr1)));
FM2=FM2/(2*pi*d2);  
FM2=flipud(fftshift(FM2,1));
FM2=FM2.*mask_im;
figure,imagesc(FM2),colormap('gray'),colorbar;

% Computing PSNRs (dB)-----------------------------------------------------
 SNR(we,est_we1.*mask_im)
 SNR(obj, abs(a))
 SNR(obj, abs(est_image1))
 
% To see there is really a shift btw DFT recon img and orig img
figure, imagesc(interp2(obj-abs(a))), colormap('gray'), title(' Difference between original image and DFT reconstruction'), colorbar

% check monotonicity condition
ii=1:N2-1;
diff=we(ii+1,:)-we(ii,:);
bound=2*pi*delta2/(N2*d2);
[m,n]=find(diff>bound);     %[m,n] should be empty