added epg-x functions

epgx/Test_MET.m

@@ -0,0 +1,29 @@
+% Script to compare to original implementation from
+% https://github.com/mriphysics/EPG-X. Need to add such library to the
+% Matlab path.
+clc;
+clear;
+close all;
+
+f = 0.20; 
+ka = 5e-2; 
+T1 = [1000 500];
+T2 = [100 20];
+esp=5;
+angles = pi*ones(1,50);
+Npulse = length(angles);
+[s,P] = epg_X_CMPG(angles,f,T1,T2,esp,ka,0);
+figure
+subplot(1,2,1)
+plot(abs(s(:,:).'),'linewidth',2)
+title('Own implementation')
+legend('1st compartment','2nd compartment','voxel average','fontsize',13)
+
+a0 = [90 180*ones(1,50)]*pi/180;
+
+[s_gt,Fn_gt] = EPGX_TSE_BM(a0,esp,T1,T2,f,ka,'delta',0);  %function taken from Malik et al.
+subplot(1,2,2)
+plot(abs(s_gt))
+plot(abs([squeeze(Fn_gt(52,:,:)) s_gt.']),'linewidth',2)
+title('Original')
+legend('1st compartment','2nd compartment','voxel average','fontsize',13)

epgx/Test_SPGR.m

@@ -0,0 +1,65 @@
+% Script to compare to original implementation from
+% https://github.com/mriphysics/EPG-X. Need to add such library to the
+% Matlab path.
+
+clc;
+clear;
+close all;
+
+TR = 5;
+alpha = 10;
+phi0 = 117*pi/180;
+T1_MT = [779 779];
+f_MT = 0.117;
+k_MT = 4.3e-3;
+T2_MT = 45;
+
+%%% Relaxation parameters: exchange
+T1x = [1000 500];
+T2x = [100 20];    
+kx = 2e-3;
+fx = 0.2;
+
+
+G = 15.1;% us 
+b1 = 13; % uT
+gam = 267.5221 *1e-3; %< rad /ms /uT
+trf = (pi*alpha/180)/(gam*b1);% ms
+b1sqrdtau = b1^2*trf;
+
+npulse=250;
+phi = RF_phase_cycle(npulse,phi0); %function taken from Malik et al.
+WT = pi*gam^2*b1sqrdtau*G*1e-3; 
+
+[s,P] = epg_X_rfspoil(alpha*pi/180,phi,T1_MT,T2_MT,TR,f_MT,k_MT,WT);
+[s2,P] = epg_X_rfspoil(alpha*pi/180,phi,T1_MT,T2_MT,TR,f_MT*2,k_MT,WT);
+[s3,P] = epg_X_rfspoil(alpha*pi/180,phi,T1_MT,T2_MT,TR,f_MT*3,k_MT,WT);
+
+
+figure
+subplot(2,1,1)
+plot(abs(s(:,:).'),'linewidth',2)
+hold on
+plot(abs(s2(:,:).'),'linewidth',2)
+plot(abs(s3(:,:).'),'linewidth',2)
+legend('f=0.12','f=0.23','f=0.35')
+xlabel('TR number')
+
+title('Own implementation')
+ylim([0 0.17])
+
+%% original
+[smt,Fnmt,Znmt] = EPGX_GRE_MT(d2r(alpha)*ones(npulse,1),phi,b1sqrdtau*ones(npulse,1),TR,T1_MT,T2_MT,f_MT,k_MT,G); %function taken from Malik et al.
+[smt2,Fnmt,Znmt] = EPGX_GRE_MT(d2r(alpha)*ones(npulse,1),phi,b1sqrdtau*ones(npulse,1),TR,T1_MT,T2_MT,f_MT*2,k_MT,G);
+[smt3,Fnmt,Znmt] = EPGX_GRE_MT(d2r(alpha)*ones(npulse,1),phi,b1sqrdtau*ones(npulse,1),TR,T1_MT,T2_MT,f_MT*3,k_MT,G);
+
+
+subplot(2,1,2)
+plot(abs(smt),'linewidth',2)
+hold on
+plot(abs(smt2),'linewidth',2)
+plot(abs(smt3),'linewidth',2)
+title('Original')
+ylim([0 0.17])
+legend('f=0.12','f=0.23','f=0.35')
+xlabel('TR number')

epgx/Test_fitting_MET.m

@@ -0,0 +1,85 @@
+%% File to simulate T2 fitting with different exchange rates
+
+clc;
+clear;
+close all;
+
+
+run('/home/gustavo/Downloads/qMRLab-master/startup.m') %download from https://qmrlab.org/
+warning ('off','all');
+
+T1 = [1000 500];
+T2 = [120 20];    
+list_kxs = [0:1:20]*1e-3; % list of exchange rates to test;
+list_f = [1:40]*1e-2; % list of fractions to test;
+
+for jj=1:length(list_kxs)
+    ka = list_kxs(jj);
+    for ii=1:length(list_f)
+        f = list_f(ii);
+        esp=5;
+        angles = pi*ones(1,50);
+        Npulse = length(angles);
+        [s,P] = epg_X_CMPG(angles,f,T1,T2,esp,ka,0);
+        all_s(ii,jj,:) = abs(s(3,:));
+    end
+end
+
+% load presaved model to use for fitting
+Model = qMRloadObj('mwf.qmrlab.mat');
+Model.Prot.MET2data.Mat = esp*(1:Npulse)';
+data = struct();
+data.MET2data= reshape(all_s,[length(list_f) length(list_kxs) 1 Npulse]);
+data.Mask= ones([length(list_f) length(list_kxs) 1]);
+FitResults = FitData(data,Model,0);
+save('FitResultsMET')
+
+figure
+subplot(1,3,1)
+selected = 1:2:length(list_kxs);
+plot(list_f,list_f);
+hold on
+plot(list_f,FitResults.MWF(:,selected)/100,'*','linewidth',1)
+xlabel('real fraction')
+ylabel('estimated fraction')
+legends={};
+for ii=1:length(selected)
+    legends{ii+1}=['ka=' num2str(list_kxs(selected(ii)))];
+
+end
+legends{1} = 'Ground truth';
+legend(legends);
+xlim([0 0.4])
+ylim([0 0.4])
+title('Fraction estimation')
+
+subplot(1,3,3)
+selected = 1:2:length(list_kxs);
+plot(list_f,T2(2)*ones(size(list_f)));
+hold on
+plot(list_f,FitResults.T2MW(:,selected),'*','linewidth',1)
+xlabel('fraction')
+ylabel('Estimated T2_b (ms)')
+legends={};
+for ii=1:length(selected)
+    legends{ii+1}=['ka=' num2str(list_kxs(selected(ii)))];
+
+end
+legends{1} = 'Ground truth';
+title('T2_b estimation')
+
+subplot(1,3,2)
+selected = 1:2:length(list_kxs);
+plot(list_f,T2(1)*ones(size(list_f)));
+hold on
+plot(list_f,FitResults.T2IEW(:,selected),'*','linewidth',1)
+xlabel('fraction')
+ylabel('Estimated T2_a(ms)')
+legends={};
+for ii=1:length(selected)
+    legends{ii+1}=['ka=' num2str(list_kxs(selected(ii)))];
+
+end
+legends{1} = 'Ground truth';
+% legend(legends);
+title('T2_a estimation')

epgx/epg_X_CMPG.m

@@ -0,0 +1,27 @@
+function [s,P] = epg_X_CMPG(flipangle,f,T1,T2,esp,ka,deltab)
+%Simulate EPG-X CMPG sequence. Only for the T2 case, not the MT case
+% flipangle: flipangle in radians
+% T1: 1x2 vector with the T1 constant of both comparment
+% T2: 1x2 vector or scalar with the T2 constants
+% ka: forward exchange rate
+% f: fraction of second compartment
+
+etl = length(flipangle);
+P = epg_X_m0(f,'full');
+P = [P zeros(6,2*etl)];
+
+% -- 90 excitation
+P = epg_X_rf(P,pi/2,pi/2);	% Do 90 tip.
+s = zeros(3,etl);		% Allocate signal vector to store.
+
+for ech=1:etl
+  P = epg_X_relax(P,esp/2,T1,T2,ka,deltab,f);
+  P = epg_X_grad(P,1);
+  P = epg_X_rf(P,abs(flipangle(ech)),angle(flipangle(ech)));   
+  P = epg_X_grad(P,1);
+  P = epg_X_relax(P,esp/2,T1,T2,ka,deltab,f);
+  s(1:2,ech) = P([1 4],1);  	% Signal is F0 state.
+  s(3,ech) =  sum(P([1 4],1));
+end
+
+end

epgx/epg_X_grad.m

@@ -0,0 +1,24 @@
+function [FpFmZ] = epg_X_grad(FpFmZ,noadd)
+%	Propagate EPG-X states through a "unit" gradient.
+%	
+%	INPUT:
+%		FpFmZ = 4xN or 6xN vector of F+, F- and Z states.
+%		noadd = 1 to NOT add any higher-order states - assume
+%			that they just go to zero.  Be careful - this
+%			speeds up simulations, but may compromise accuracy!
+%
+%	OUTPUT:
+%		Updated FpFmZ state.
+%
+if (nargin < 2) noadd=0; end;	% Add by default.  
+if size(FpFmZ,1) == 4 %MT case
+    temp1 = epg_grad(FpFmZ(1:3,:),noadd);
+    FpFmZ = [temp1; [FpFmZ(4,:) zeros(1,size(temp1,2)-size(FpFmZ,2))]];
+else % full case
+    temp1 = epg_grad(FpFmZ(1:3,:),noadd);
+    temp2 = epg_grad(FpFmZ(4:6,:),noadd);
+    FpFmZ = [temp1;temp2];
+end
+
+end
+

epgx/epg_X_m0.m

@@ -0,0 +1,21 @@
+function [FpFmZ] = epg_X_m0(f,opt)
+%	Creates initial EPG-X magnetizations state.
+%	
+%	INPUT:
+%		f = fraction of the second compartment
+%       opt= put 'MT' to create the reduced MT case (4 components) or
+%       empty for the full state
+%
+%	OUTPUT:
+%		FpFmZ state.
+%
+    if nargin<2
+        opt = 'FULL';
+    end
+    FpFmZ = zeros(6,1);
+    FpFmZ(3)=(1-f);
+    FpFmZ(6)=f;
+    if strcmp(opt,'MT')
+        FpFmZ(4:5)=[];
+    end
+end

epgx/epg_X_relax.m

@@ -0,0 +1,41 @@
+function [FpFmZ] = epg_X_relax(FpFmZ,T,T1,T2,ka,deltab,f)
+%	Relaxation and exchange of  EPG-X.
+%	
+%	INPUT:
+%		FpFmZ = 4xN or 6xN vector of F+, F- and Z states.
+%       T: time for the relaxation
+%       T1: 1x2 vector with the T1 constant of both comparment
+%       T2: 1x2 vector or scalar with the T2 constants
+%       ka: forward exchange rate
+%       deltab: frequency offset of 2nd comparment
+%       f: fraction of second compartment
+%
+%	OUTPUT:
+%		FpFmZ state.
+%
+   M0b = f; 
+   M0a = 1-f; 
+   kb = ka * M0a/M0b; % conservation of magnetization
+
+   R1a=1/T1(1);
+   R1b=1/T1(2); 
+   R2a=1/T2(1);
+   LambdaL = expm([-R1a-ka kb; ka -R1b-kb]*T); %longitudinal relaxation and exchange
+   E1 = exp(-T./T1');
+   Zoff = [M0a;M0b].*(1-E1); % longitudinal recovery 
+   if size(FpFmZ,1) == 4 %MT case
+        E2a = exp(-T*R2a);
+        FpFmZ(1:2,:) = diag([E2a E2a])*FpFmZ(1:2,:);
+        FpFmZ(3:4,:) = LambdaL*FpFmZ(3:4,:);
+        FpFmZ(3:4,1) = FpFmZ(3:4,1)+Zoff;
+   else % full case
+        R2b=1/T2(2); 
+        LambdaT = expm(T*[-R2a-ka 0 kb 0; 0 -R2a-ka 0 kb; ka 0 -R2b-kb-2*pi*1i*deltab 0; 0 ka 0 -R2b-kb+2*pi*1i*deltab]); %xversal relaxation and exchange
+        temp1 = LambdaT*FpFmZ([1;2;4;5],:); 
+        temp2 = LambdaL*FpFmZ([3;6],:);
+        temp2(:,1) = temp2(:,1) + Zoff;
+        FpFmZ = [temp1(1:2,:);temp2(1,:);temp1(3:4,:);temp2(2,:)];
+   end
+
+end
+

epgx/epg_X_rf.m

@@ -0,0 +1,23 @@
+function [FpFmZ,RR] = epg_X_rf(FpFmZ,alpha,phi,WT)
+%	Simulates action of RF pulse on EPG-X states
+%	
+%	INPUT:
+%		FpFmZ = 4xN or 6xN vector of F+, F- and Z states.
+%       alpha: flip angle in radians
+%       phi: phase of the RF pulse (radians)
+%       WT: saturation of the bound pool componen will be exp(-WT), ignore for
+%       full case
+%	OUTPUT:
+%		FpFmZ state.
+%       RR: rotation matrix
+
+    if size(FpFmZ,1) == 4 %MT case      
+        [temp1,RR] = epg_rf(FpFmZ(1:3,:),alpha,phi);
+        FpFmZ = [temp1; exp(-WT)*FpFmZ(4,:)];
+    else % full case
+        [temp1,RR] = epg_rf(FpFmZ(1:3,:),alpha,phi);
+        [temp2,RR] = epg_rf(FpFmZ(4:6,:),alpha,phi);
+        FpFmZ = [temp1;temp2];
+    end
+end
+

epgx/epg_X_rfspoil.m

@@ -0,0 +1,32 @@
+function [s,P] = epg_X_rfspoil(flipangle,rf_phase,T1,T2,TR,f,ka,WT)
+% Simulate RF spoiled sequence, 
+% flipangle: flipangle in radians
+% rf_phase: the vector of phases for the sequence of RF pulses
+% T1: 1x2 vector with the T1 constant of both comparment
+% T2: 1x2 vector or scalar with the T2 constants
+% ka: forward exchange rate
+% f: fraction of second compartment
+% WT: saturation of the bound pool componen will be exp(-WT), ignore for
+% full case
+
+if exist('WT','var')
+    P = epg_X_m0(f,'MT');
+    caso = 'MT';
+else
+    P = epg_X_m0(f,'full');
+    caso = 'full';
+end
+
+P = [P zeros(size(P,1),250)];
+N = length(rf_phase);
+for n=1:N
+    P = epg_X_relax(P,TR,T1,T2,ka,0,f);
+    P = epg_X_grad(P,1);  
+  if strcmp(caso,'MT')
+    P = epg_X_rf(P,flipangle,rf_phase(n),WT);	% RF excitation
+        s(1,n) =  P(1,1);
+  else
+      P = epg_X_rf(P,flipangle,rf_phase(n));	% RF excitation
+        s(1:2,n) = P([1 4],1);  	% Signal is F0 state.
+  end
+end