From b8349e90244664558e61d12272112cc339362711 Mon Sep 17 00:00:00 2001 From: Agah Date: Wed, 9 Oct 2019 08:52:05 -0400 Subject: [PATCH 1/6] Add data folder to gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index cd2946a..c641214 100644 --- a/.gitignore +++ b/.gitignore @@ -45,3 +45,5 @@ $RECYCLE.BIN/ Network Trash Folder Temporary Items .apdisk + +data From fb4a2c3a4857e7c8362858d80851e730afd46c42 Mon Sep 17 00:00:00 2001 From: Agah Date: Wed, 9 Oct 2019 08:56:29 -0400 Subject: [PATCH 2/6] Fix check for nargin count * First varargin is scoped when nargin==6 * Set default invEFF to 0.96 (based on Demos) * Minor syntax clean up --- func/T1B1correctpackageTFL.m | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/func/T1B1correctpackageTFL.m b/func/T1B1correctpackageTFL.m index 860ba44..ec1143f 100644 --- a/func/T1B1correctpackageTFL.m +++ b/func/T1B1correctpackageTFL.m @@ -1,4 +1,4 @@ -function [ T1temp MP2RAGEcorrected] = T1B1correctpackageTFL( B1img,MP2RAGEimg,T1,MP2RAGE,brain,varargin) +function [ T1temp,MP2RAGEcorrected] = T1B1correctpackageTFL( B1img,MP2RAGEimg,T1,MP2RAGE,brain,varargin) % usage % % [ T1corr MP2RAGEcorr] = T1B1correctpackage(B1,MP2RAGEimg,T1,MP2RAGE,brain,varargin) @@ -43,19 +43,19 @@ % % please cite: % Marques, J.P., Gruetter, R., 2013. New Developments and Applications of the MP2RAGE Sequence - Focusing the Contrast and High Spatial Resolution R1 Mapping. PLoS ONE 8. doi:10.1371/journal.pone.0069294 -% Marques, J.P., Kober, T., Krueger, G., van der Zwaag, W., Van de Moortele, P.-F., Gruetter, R., 2010a. MP2RAGE, a self bias-field corrected sequence for improved segmentation and T1-mapping at high field. NeuroImage 49, 1271–1281. doi:10.1016/j.neuroimage.2009.10.002 +% Marques, J.P., Kober, T., Krueger, G., van der Zwaag, W., Van de Moortele, P.-F., Gruetter, R., 2010a. MP2RAGE, a self bias-field corrected sequence for improved segmentation and T1-mapping at high field. NeuroImage 49, 1271–1281. doi:10.1016/j.neuroimage.2009.10.002 % -if nargin==8 +if nargin==6 invEFF=varargin{1} else - invEFF=0.99; + invEFF=0.96; -end; +end if isempty(brain) @@ -71,9 +71,9 @@ brain.img=ones(size(brain.img)); - end; + end -end; +end %% sanity check to see how B1 sensitive your sequence was @@ -85,7 +85,7 @@ for B1val=0.6:0.2:1.4 - [MP2RAGEamp T1vector IntensityBeforeComb]=MP2RAGE_lookuptable(2,MP2RAGE.TR,MP2RAGE.TIs,B1val*MP2RAGE.FlipDegrees,MP2RAGE.NZslices,MP2RAGE.TRFLASH,'normal'); + [MP2RAGEamp,T1vector,IntensityBeforeComb]=MP2RAGE_lookuptable(2,MP2RAGE.TR,MP2RAGE.TIs,B1val*MP2RAGE.FlipDegrees,MP2RAGE.NZslices,MP2RAGE.TRFLASH,'normal'); plot(MP2RAGEamp,T1vector,'color',[0.5 0.5 0.5]*B1val,'Linewidth',2) @@ -125,9 +125,9 @@ T1CSF=3.5; - end; + end -end; +end plot([-0.5 0.5],[T1CSF T1CSF;T1GM T1GM;T1WM T1WM]','Linewidth',2) @@ -161,7 +161,7 @@ MP2RAGEimg.img=double(MP2RAGEimg.img)/4095-0.5; -end; +end %% now the fun starts @@ -178,7 +178,7 @@ MP2RAGEmatrix(k,:)=interp1(T1vector,Intensity,T1_vector); -end; +end k=0; @@ -212,9 +212,9 @@ - end; + end -end; +end %% correcting the estimates of T1 and B1 iteratively @@ -295,7 +295,7 @@ -end; +end end From 2c809cd7b9f67eb8d94950a3394a033cce1a7260 Mon Sep 17 00:00:00 2001 From: Agah Date: Wed, 9 Oct 2019 09:00:14 -0400 Subject: [PATCH 3/6] Reflect changes from fb4a2c3 --- func/T1B1correctpackageTFL.asv | 204 ++++++++++++++++----------------- 1 file changed, 100 insertions(+), 104 deletions(-) diff --git a/func/T1B1correctpackageTFL.asv b/func/T1B1correctpackageTFL.asv index b683bc9..034c792 100644 --- a/func/T1B1correctpackageTFL.asv +++ b/func/T1B1correctpackageTFL.asv @@ -1,5 +1,5 @@ -function [ T1temp MP2RAGEcorrected] = T1B1correctpackageTFL( B1img,MP2RAGEimg,T1,MP2RAGE,brain,varargin) -% usage +function [ T1temp,MP2RAGEcorrected] = T1B1correctpackageTFL( B1img,MP2RAGEimg,T1,MP2RAGE,brain,varargin) +% usage % % [ T1corr MP2RAGEcorr] = T1B1correctpackage(B1,MP2RAGEimg,T1,MP2RAGE,brain,varargin) % @@ -7,7 +7,7 @@ function [ T1temp MP2RAGEcorrected] = T1B1correctpackageTFL( B1img,MP2RAGEimg,T1 % the MP2RAGE (MP2RAGEimg, T1) and the result of some B1 mapping technique % with load_nii or load_untouch_nii % -% the variable B1 is compulsory +% the variable B1 is compulsory % % Only MP2RAGEimg or the T1 have to be loaded (I usually use the MP2RAGEimg) % @@ -20,7 +20,7 @@ function [ T1temp MP2RAGEcorrected] = T1B1correctpackageTFL( B1img,MP2RAGEimg,T1 % from 0-2 % % MP2RAGE.B0=7; % in Tesla -% MP2RAGE.TR=6; % MP2RAGE TR in seconds +% MP2RAGE.TR=6; % MP2RAGE TR in seconds % MP2RAGE.TRFLASH=6.7e-3; % TR of the GRE readout % MP2RAGE.TIs=[800e-3 2700e-3];% inversion times - time between middle of refocusing pulse and excitatoin of the k-space center encoding % MP2RAGE.NZslices=[40 80];% Slices Per Slab * [PartialFourierInSlice-0.5 0.5] @@ -32,48 +32,48 @@ function [ T1temp MP2RAGEcorrected] = T1B1correctpackageTFL( B1img,MP2RAGEimg,T1 % everywhere % % additionally the inversion efficiency of the adiabatic inversion can be -% set as a last optional variable. Ideally it should be 1. -% In the first implementation of the MP2RAGE the inversino efficiency was +% set as a last optional variable. Ideally it should be 1. +% In the first implementation of the MP2RAGE the inversino efficiency was % measured to be ~0.96 % % outputs are: % T1corr - T1map corrected for B1 bias -% MP2RAGEcorr - MP2RAGE image corrected for B1 bias -% -% +% MP2RAGEcorr - MP2RAGE image corrected for B1 bias +% +% % please cite: % Marques, J.P., Gruetter, R., 2013. New Developments and Applications of the MP2RAGE Sequence - Focusing the Contrast and High Spatial Resolution R1 Mapping. PLoS ONE 8. doi:10.1371/journal.pone.0069294 -% Marques, J.P., Kober, T., Krueger, G., van der Zwaag, W., Van de Moortele, P.-F., Gruetter, R., 2010a. MP2RAGE, a self bias-field corrected sequence for improved segmentation and T1-mapping at high field. NeuroImage 49, 1271–1281. doi:10.1016/j.neuroimage.2009.10.002 +% Marques, J.P., Kober, T., Krueger, G., van der Zwaag, W., Van de Moortele, P.-F., Gruetter, R., 2010a. MP2RAGE, a self bias-field corrected sequence for improved segmentation and T1-mapping at high field. NeuroImage 49, 1271–1281. doi:10.1016/j.neuroimage.2009.10.002 % -if nargin==8 - +if nargin==6 + invEFF=varargin{1} - + else - - invEFF=0.99; - -end; + + invEFF=0.96; + +end if isempty(brain) - + if isempty(MP2RAGEimg) - + brain=T1; - + brain.img=ones(size(brain.img)); - + else - + brain=MP2RAGEimg; - + brain.img=ones(size(brain.img)); - - end; - -end; + + end + +end %% sanity check to see how B1 sensitive your sequence was @@ -84,13 +84,13 @@ set(gcf,'Color',[1 1 1]); hold off for B1val=0.6:0.2:1.4 - - [MP2RAGEamp T1vector IntensityBeforeComb]=MP2RAGE_lookuptable(2,MP2RAGE.TR,MP2RAGE.TIs,B1val*MP2RAGE.FlipDegrees,MP2RAGE.NZslices,MP2RAGE.TRFLASH,'normal'); - + + [MP2RAGEamp,T1vector,IntensityBeforeComb]=MP2RAGE_lookuptable(2,MP2RAGE.TR,MP2RAGE.TIs,B1val*MP2RAGE.FlipDegrees,MP2RAGE.NZslices,MP2RAGE.TRFLASH,'normal'); + plot(MP2RAGEamp,T1vector,'color',[0.5 0.5 0.5]*B1val,'Linewidth',2) - + hold on - + end legend('B1=0.6','B1=0.8','B1=1','B1=1.2','B1=1.4') @@ -98,36 +98,36 @@ legend('B1=0.6','B1=0.8','B1=1','B1=1.2','B1=1.4') % examples of T1 values at 3T if ~isfield(MP2RAGE,'B0') - + T1WM=1.1; - + T1GM=1.85; - + T1CSF=3.5; - + else - + if MP2RAGE.B0==3 - + T1WM=0.85; - + T1GM=1.35; - + T1CSF=2.8; - + else - + % examples of T1 values at 7T - + T1WM=1.1; - + T1GM=1.85; - + T1CSF=3.5; - - end; - -end; + + end + +end plot([-0.5 0.5],[T1CSF T1CSF;T1GM T1GM;T1WM T1WM]','Linewidth',2) @@ -150,18 +150,18 @@ T1_vector=0.5:0.05:5.2; [MP2RAGE.Intensity MP2RAGE.T1vector ]=MP2RAGE_lookuptable(2,MP2RAGE.TR,MP2RAGE.TIs,MP2RAGE.FlipDegrees,MP2RAGE.NZslices,MP2RAGE.TRFLASH,'normal',invEFF) if isempty(MP2RAGEimg) - + T1.img=double(T1.img)/1000; - + MP2RAGEimg.img=reshape(interp1(MP2RAGE.T1vector,MP2RAGE.Intensity,T1.img(:)),size(B1img.img)); - + MP2RAGEimg.img(isnan(MP2RAGEimg.img))=-0.5 - + else - + MP2RAGEimg.img=double(MP2RAGEimg.img)/4095-0.5; - -end; + +end %% now the fun starts @@ -171,14 +171,14 @@ end; k=0; for b1val=B1_vector - + k=k+1; - + [Intensity T1vector ]=MP2RAGE_lookuptable(2,MP2RAGE.TR,MP2RAGE.TIs,b1val*MP2RAGE.FlipDegrees,MP2RAGE.NZslices,MP2RAGE.TRFLASH,'normal',invEFF); - + MP2RAGEmatrix(k,:)=interp1(T1vector,Intensity,T1_vector); - -end; + +end k=0; @@ -193,36 +193,34 @@ MP2RAGE_vector=linspace(-0.5,0.5,npoints); k=0; for b1val=B1_vector - + k=k+1; - + try - + T1matrix(k,:)=interp1(MP2RAGEmatrix(k,:),T1_vector,MP2RAGE_vector,'pchirp'); - + catch - + temp=MP2RAGEmatrix(k,:); - - temp(isnan(temp))=linspace(-0.5,-1,length(find(isnan(temp)==1))); - + + temp(isnan(temp))=linspace(-0.5-eps,-1,length(find(isnan(temp)==1))); + temp=interp1(temp,T1_vector,MP2RAGE_vector); - + T1matrix(k,:)=temp; - - - end; - -end; + + + end + +end %% correcting the estimates of T1 and B1 iteratively T1temp=MP2RAGEimg; -B1temp=B1; - -brain.img(B1.img==0)=0; +brain.img(B1img.img==0)=0; brain.img(MP2RAGEimg.img==0)=0; @@ -230,37 +228,36 @@ T1temp.img(brain.img==0)=0; T1temp.img(brain.img==1)=0; -B1temp.img(brain.img==0)=0; +B1img.img(brain.img==0)=0; - btemp2=squeeze(B1temp.img(:,end/2,:)); - temp=squeeze(T1temp.img(:,end/2,:)); +temp=squeeze(T1temp.img(:,end/2,:)); - - T1temp.img(brain.img~=0)=interp2(MP2RAGE_vector,B1_vector,T1matrix,MP2RAGEimg.img(brain.img~=0),B1temp.img(brain.img~=0)) - T1temp.img(isnan(T1temp.img))=4; +T1temp.img(brain.img~=0)=interp2(MP2RAGE_vector,B1_vector,T1matrix,MP2RAGEimg.img(brain.img~=0),B1img.img(brain.img~=0)) - +T1temp.img(isnan(T1temp.img))=4; - temp2=squeeze(T1temp.img(:,end/2,:)); - - gcf=figure(1); +temp2=squeeze(T1temp.img(:,end/2,:)); + - set(gcf,'Color',[1 1 1]); - imagesc(temp2-temp);colorbar +gcf=figure(1); + +set(gcf,'Color',[1 1 1]); - title(['T1 correction ']); +imagesc(temp2-temp);colorbar + +title(['T1 correction ']); + +colormap(gray) - colormap(gray) - %% creates an MP2RAGEcorrected image and puts both the B1 and T1 in the ms scale [MP2RAGE.Intensity MP2RAGE.T1vector ]=MP2RAGE_lookuptable(2,MP2RAGE.TR,MP2RAGE.TIs,MP2RAGE.FlipDegrees,MP2RAGE.NZslices,MP2RAGE.TRFLASH,'normal',invEFF) @@ -275,31 +272,30 @@ MP2RAGEcorrected.img=round(4095*(MP2RAGEcorrected.img+0.5)); T1temp.img=(T1temp.img)*1000; -B1temp.img=(B1temp.img)*1000; %% showimages=0 if showimages==1; - + imagesc(MP2RAGE_vector,B1_vector,T1matrix,[0.4 5]);colorbar; - + xlabel ('MP2RAGE','FontSize',12,'FontWeight','bold') - + ylabel ('B_1','FontSize',12,'FontWeight','bold') - + title('T_1 look-up table','FontSize',12,'FontWeight','bold') - + H=gca - + set(H,'FontSize',12,'LineWidth',2) - + axes(H) - - -end; + + +end end From ed5b12d4d24de04b1a572cd5f34e0b4cd6f5b603 Mon Sep 17 00:00:00 2001 From: Agah Date: Wed, 9 Oct 2019 09:01:09 -0400 Subject: [PATCH 4/6] Minor syntax cedit --- func/T1M0estimateMP2RAGE.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/func/T1M0estimateMP2RAGE.m b/func/T1M0estimateMP2RAGE.m index 518ac3d..28cbc7f 100644 --- a/func/T1M0estimateMP2RAGE.m +++ b/func/T1M0estimateMP2RAGE.m @@ -35,9 +35,9 @@ invEFF=varargin{1}; else invEFF=0.96; -end; +end -[Intensity T1vector IntensityUncomb]=MP2RAGE_lookuptable(2,MP2RAGE.TR,MP2RAGE.TIs,MP2RAGE.FlipDegrees,MP2RAGE.NZslices,MP2RAGE.TRFLASH,'normal',invEFF); +[Intensity,T1vector,IntensityUncomb]=MP2RAGE_lookuptable(2,MP2RAGE.TR,MP2RAGE.TIs,MP2RAGE.FlipDegrees,MP2RAGE.NZslices,MP2RAGE.TRFLASH,'normal',invEFF); % in a first instance the T1 map is computed. @@ -47,7 +47,7 @@ T1=interp1(Intensity,T1vector,-0.5+1/4095*double(MP2RAGEnii.img(:))); else T1=interp1(Intensity,T1vector,MP2RAGEnii.img(:)); -end; +end T1(isnan(T1))=0; % copies the header from the MP2RAGEnii From 4030b13cc3f463146b4f7891c6fcde229590f412 Mon Sep 17 00:00:00 2001 From: Agah Date: Wed, 9 Oct 2019 09:01:46 -0400 Subject: [PATCH 5/6] Minor syntax edit --- func/B1mappingSa2RAGElookuptable.m | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/func/B1mappingSa2RAGElookuptable.m b/func/B1mappingSa2RAGElookuptable.m index e2828e1..cc0e649 100644 --- a/func/B1mappingSa2RAGElookuptable.m +++ b/func/B1mappingSa2RAGElookuptable.m @@ -1,4 +1,4 @@ -function [B1vector Intensity Signal]=B1mappingSa2RAGElookuptable(nimage,MPRAGE_tr,invtimesAB,flipangleABdegree,nZslices,FLASH_tr,varargin) +function [B1vector,Intensity,Signal]=B1mappingSa2RAGElookuptable(nimage,MPRAGE_tr,invtimesAB,flipangleABdegree,nZslices,FLASH_tr,varargin) % usage % [B1map]=B1mappingSa2RAGE(Rationii,nimage,MPRAGE_tr,invtimesAB,flipangleABdegree,nZslices,FLASH_tr,varargin) % varargin{1} is a phase image @@ -9,7 +9,7 @@ T1average=1.5; else T1average=varargin{2}; -end; +end T1average; B1vector=0.005:0.005:2.5; @@ -24,7 +24,7 @@ nZ_bef=nZslices/2; nZ_aft=nZslices/2; nZslices2=(nZslices); -end; +end @@ -35,8 +35,8 @@ Signal(m,1:2)=1*MPRAGEfunc(nimage,MPRAGE_tr,invtimesAB,nZslices2,FLASH_tr,B1*[flipangleABdegree],'normal',T1average,-cos((B1*pi/2))); else Signal(m,1:2)=0; - end; -end; + end +end Intensity=squeeze(real(Signal(:,1))./(real(Signal(:,2)))); B1vector=squeeze(B1vector); From 41279d81956daaa9462f406442ead50a19c63e52 Mon Sep 17 00:00:00 2001 From: Agah Date: Wed, 9 Oct 2019 09:02:38 -0400 Subject: [PATCH 6/6] Minor syntax edit --- func/MP2RAGE_lookuptable.m | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/func/MP2RAGE_lookuptable.m b/func/MP2RAGE_lookuptable.m index a562be3..deca6a8 100644 --- a/func/MP2RAGE_lookuptable.m +++ b/func/MP2RAGE_lookuptable.m @@ -1,5 +1,5 @@ -function [Intensity T1vector IntensityBeforeComb]=MP2RAGE_lookuptable(nimages,MPRAGE_tr,invtimesAB,flipangleABdegree,nZslices,FLASH_tr,sequence,varargin) +function [Intensity,T1vector,IntensityBeforeComb]=MP2RAGE_lookuptable(nimages,MPRAGE_tr,invtimesAB,flipangleABdegree,nZslices,FLASH_tr,sequence,varargin) % first extra parameter is the inversion efficiency % second extra parameter is the alldata % if ==1 all data is shown @@ -31,7 +31,7 @@ nZ_bef=nZslices/2; nZ_aft=nZslices/2; nZslices2=(nZslices); -end; +end @@ -55,16 +55,16 @@ else Signal(j,m,1:2)=1*MPRAGEfunc(nimages,MPRAGEtr,inversiontimes2,nZslices2,FLASH_tr,B1*[flipanglea flipangleb],sequence,T1); - end; - end; + end + end else Signal(j,m,1:2)=0; - end; - end; - end; - end; - end; -end; + end + end + end + end + end +end Intensity=squeeze(real(Signal(:,:,1).*conj(Signal(:,:,2)))./(abs(Signal(:,:,1)).^2+abs(Signal(:,:,2)).^2)); T1vector=squeeze(T1vector); if alldata==0 @@ -78,4 +78,4 @@ Intensity=squeeze(real(Signal(:,:,1).*conj(Signal(:,:,2)))./(abs(Signal(:,:,1)).^2+abs(Signal(:,:,2)).^2)); T1vector=squeeze(T1vector); IntensityBeforeComb=squeeze(Signal(:,1,:)); -end; +end