% This code generates all FIGURES presented in the paper:
%--------------------------------------------------------------------------
% [1]  Joy A, Paul J S. "Multi-channel MR Image Reconstruction using
%      Statistically Optimized Non-Linear Diffusion", Magnetic Resonance in 
%      Medicine, June 2017, doi:10.1002/mrm.26774.
% Email: ajin.joy@iiitmk.ac.in, j.paul@iiitmk.ac.in
%
% Coded by: Ajin Joy on March 2017
% 
% These codes and data are shared to help other researchers understand
% and reproduce our work.
%
% This implimentation uses data avaliable for download at
% http://www.quxiaobo.org/ (included)
%--------------------------------------------------------------------------

clc; clear all; close all;
disp('Time to run this code:');
disp('~20 minutes on a quad core(2.4GHz) machine running MATLAB 2015a'); 
params.gamma  = .1;  % moderate diffusion strength. (0<gamma<.25]
params.c      = 1;   % maximum consistency. [0<c<1]
params.numitr = 200; % maximum number of iterations
load datasets;


% COMPUTE RESULTS
%--------------------------------------------------------------------------
disp('::::running PM diffusion');
DF =@(dI,alpha_sq) 1./(1 + (dI.*conj(dI))/(alpha_sq));
disp('  Figure 3..');
[PMmad]                   = generate_fig3(dataset,DF,params);
disp('  Figure 3 completed');
disp('  Figure 4..');
I                         = phantom(128); % Figure 4 in paper used phantom(256)
addn                      =@(I,std) I+std*(wgn(128,128,0)+1i*wgn(128,128,0));
gamma                     = linspace(.0001,.25,100);
[PMgam]                   = generate_fig4(dataset,DF,params,I,addn);
disp('  Figure 5..');
[OutPM11,OutPM12,OutPM13] = generate_fig5and6(dataset{4},DF,params);
PMerr11                   = OutPM11.Idif; 
PMerr12                   = OutPM12.Idif;
PMerr13                   = OutPM13.Idif;
disp('  Figure 6..');
[OutPM21,OutPM22,OutPM23] = generate_fig5and6(dataset{5},DF,params);
PMerr21                   = OutPM21.Idif; 
PMerr22                   = OutPM22.Idif;
PMerr23                   = OutPM23.Idif;
%--------------------------------------------------------------------------
disp('::::running TV');
DF =@(dI,dummy) 1./((dI.*conj(dI) + 1e-10).^(1/2));
disp('  Figure 4..');
[TVgam]                   = generate_fig4(dataset,DF,params,I,addn);
disp('  Figure 4 completed');
disp('  Figure 5..');
params.gamma              = 4.6057e-07;
[OutTV11,OutTV12,OutTV13] = generate_fig5and6(dataset{4},DF,params);
TVerr11                   = OutTV11.Idif; 
TVerr12                   = OutTV12.Idif;
TVerr13                   = OutTV13.Idif;
disp('  Figure 5 completed');
disp('  Figure 6..');
[OutTV21,OutTV22,OutTV23] = generate_fig5and6(dataset{5},DF,params);
TVerr21                   = OutTV21.Idif; 
TVerr22                   = OutTV22.Idif;
TVerr23                   = OutTV23.Idif;
disp('  Figure 6 completed');
%--------------------------------------------------------------------------


% DISPLAY FIGURES
%--------------------------------------------------------------------------
figure(1);
for i=1:5
subplot(2,3,i); plot(PMmad{i}{1},'LineWidth',2,'Color',[0 0.4470 0.7410]);
title(sprintf('dataset %d',i),'fontweight','bold');xlabel('iterations');ylabel('RLNE');drawnow
hold on; plot(PMmad{i}{2},'LineWidth',2,'Color',[0.8500 0.3250 0.0980]);
hold on; plot(PMmad{i}{3},'LineWidth',2,'Color',[0.9290 0.6940 0.1250]);
hold on; plot(PMmad{i}{4},'LineWidth',2,'Color',[0.4940 0.1840 0.5560]);
hold on; plot(PMmad{i}{5},'LineWidth',2,'Color',[0.4660 0.6740 0.1880]);
hold on; plot(PMmad{i}{6},'LineWidth',2,'Color',[0.3010 0.7450 0.9330]);
legend('20% MAD','30% MAD','50% MAD','70% MAD','100% MAD','150% MAD');
end
suptitle('Figure 3');
figure(2);
for i=1:4
subplot(3,2,i); plot(gamma,PMgam{i},'LineWidth',2,'Color',[0 0.4470 0.7410]);
xlabel('iterations');ylabel('RLNE');drawnow
hold on; plot(gamma,TVgam{i},'LineWidth',2,'Color',[0.8500 0.3250 0.0980]);
legend('PM','TV');
end
subplot(3,2,5); 
plot(gamma,PMgam{1},'LineWidth',2,'Color',[0 0.4470 0.7410]);axis([0,.25,0,.5]);
xlabel('iterations');ylabel('RLNE');drawnow
hold on; plot(gamma,PMgam{2},'LineWidth',2,'Color',[0.8500 0.3250 0.0980]);axis([0,.25,0,.5]);
hold on; plot(gamma,PMgam{3},'LineWidth',2,'Color',[0.9290 0.6940 0.1250]);axis([0,.25,0,.5]);
hold on; plot(gamma,PMgam{4},'LineWidth',2,'Color',[0.4940 0.1840 0.5560]);axis([0,.25,0,.5]);
legend('\sigma = .005','\sigma = .010','\sigma = .015','\sigma = .020');
subplot(3,2,6); 
plot(gamma,TVgam{1},'LineWidth',2,'Color',[0 0.4470 0.7410]);axis([0,.03,0,.5]);
xlabel('iterations');ylabel('RLNE');drawnow
hold on; plot(gamma,TVgam{2},'LineWidth',2,'Color',[0.8500 0.3250 0.0980]);axis([0,.03,0,.5]);
hold on; plot(gamma,TVgam{3},'LineWidth',2,'Color',[0.9290 0.6940 0.1250]);axis([0,.03,0,.5]);
hold on; plot(gamma,TVgam{4},'LineWidth',2,'Color',[0.4940 0.1840 0.5560]);axis([0,.03,0,.5]);
legend('\sigma = .005','\sigma = .010','\sigma = .015','\sigma = .020');
suptitle('Figure 4');
figure(3);
cl1 = cat(1,abs(OutPM11.Izf),abs(OutPM12.Izf),abs(OutPM13.Izf));
cl2 = cat(1,abs(OutPM11.Irecon),abs(OutPM12.Irecon),abs(OutPM13.Irecon));
cl3 = cat(1,abs(OutTV11.Irecon),abs(OutTV12.Irecon),abs(OutTV13.Irecon));
subplot(1,3,1); imagesc(cl1);colormap gray;axis off
title(sprintf('(Zero-Filled)'));drawnow
subplot(1,3,2); imagesc(cl2);colormap gray;axis off
title(sprintf('(PM diffusion)'));drawnow
subplot(1,3,3); imagesc(cl3);colormap gray;axis off
title(sprintf('(TV)'));drawnow
suptitle('Figure 5: rows - x3,x4,x5 undersampling');
figure(4);
cl1 = cat(1,abs(OutPM21.Izf),abs(OutPM22.Izf),abs(OutPM23.Izf));
cl2 = cat(1,abs(OutPM21.Irecon),abs(OutPM22.Irecon),abs(OutPM23.Irecon));
cl3 = cat(1,abs(OutTV21.Irecon),abs(OutTV22.Irecon),abs(OutTV23.Irecon));
subplot(1,3,1); imagesc(cl1);colormap gray;axis off
title(sprintf('(Zero-Filled)'));drawnow
subplot(1,3,2); imagesc(cl2);colormap gray;axis off
title(sprintf('(PM diffusion)'));drawnow
subplot(1,3,3); imagesc(cl3);colormap gray;axis off
title(sprintf('(TV)'));drawnow
suptitle('Figure 6: rows - x3,x4,x5 undersampling');
figure(5);
maxerr=max(abs([PMerr11(:);PMerr12(:);PMerr13(:);TVerr11(:);TVerr12(:);TVerr13(:)]));
cl1 = cat(1,PMerr11/maxerr,PMerr12/maxerr,PMerr13/maxerr);
cl2 = cat(1,TVerr11/maxerr,TVerr12/maxerr,TVerr13/maxerr);
maxerr=max(abs([PMerr21(:);PMerr22(:);PMerr23(:);TVerr21(:);TVerr22(:);TVerr23(:)]));
cl3 = cat(1,PMerr21/maxerr,PMerr22/maxerr,PMerr23/maxerr);
cl4 = cat(1,TVerr21/maxerr,TVerr22/maxerr,TVerr23/maxerr);
subplot(1,4,1); imagesc(cl1);colormap hsv; axis off;caxis([-1 1]);
title(sprintf('(Error-PM)'));drawnow
subplot(1,4,2); imagesc(cl2);colormap hsv; axis off;caxis([-1 1]);
title(sprintf('(Error-TV)'));drawnow
subplot(1,4,3); imagesc(cl3);colormap hsv; axis off;caxis([-1 1]);
title(sprintf('(Error-PM)'));drawnow
subplot(1,4,4); imagesc(cl4);colormap hsv; axis off;caxis([-1 1]);
title(sprintf('(Error-TV)'));drawnow
suptitle('Figure 6: Error images of dataset IV and V');
%--------------------------------------------------------------------------