diff --git a/License b/License index 86611d8..46d4823 100644 --- a/License +++ b/License @@ -13,5 +13,5 @@ You should have received a copy of the GNU General Public License along with this program. If not, see . -This software is not for clinical use. +This software is for academic research purpose only. Commercial or clinical use is not permitted. Dependencies such as SEPIA toolbox and SPM12 have their own licenses. Please see README.md file for details about the dependencies. diff --git a/MPM_QSM.m b/MPM_QSM.m index 9957cff..275fd93 100644 --- a/MPM_QSM.m +++ b/MPM_QSM.m @@ -35,9 +35,6 @@ % (also known as iterative Tikhonov), % may give more contrast than Star-QSM but is less robust to noise % 'Star-QSM' - is very robust to noise and quick -% calc_mean_qsm : 'yes' or 'no', if 'yes' mean QSM from the three -% MPM acquisitions will be calculated, ATTENTION: currently nor -% coregistration between the scans is implemented %%%% Inputs - directories, parameters and files specific to given contrast % ATTENTION: ensure only niftis you want to use are in that folder, with increasing echo numbering: @@ -72,7 +69,6 @@ % script created by Barbara Dymerska % @ UCL FIL Physics -% last modifications 09/09/2021 function [QSM_V, QSM , QSMinvrot_V, QSMinvrot] = MPM_QSM(para) @@ -82,111 +78,72 @@ ph_fulldir = fullfile(para.in_root_dir, para.ph_dir) ; output_fulldir = fullfile(para.out_root_dir, para.output_dir) ; + if ~exist(output_fulldir, 'dir') mkdir(output_fulldir) end cd(output_fulldir) - TEs = para.TEs/10^3 ; + TEs = para.TEs ; ph_files = spm_select('FPList', ph_fulldir, '^.*\.(nii|img)$'); mag_files = spm_select('FPList', mag_fulldir, '^.*\.(nii|img)$'); - + + ph_1tp = nifti(ph_files(1,:)); + ph = zeros([size(ph_1tp.dat) size(TEs,2)]); + mag = zeros([size(ph_1tp.dat) size(TEs,2)]); for t = 1:size(TEs,2) ph_1tp = nifti(ph_files(t,:)); ph(:,:,:,t) = ph_1tp.dat(:,:,:) ; - - if size(TEs,2) >= 6 % for mtw acquisition we cannot perform complex fit so we don't need magnitude data for each TE - mag_1tp = nifti(mag_files(t,:)); - mag(:,:,:,t) = mag_1tp.dat(:,:,:) ; - - end + mag_1tp = nifti(mag_files(t,:)); + mag(:,:,:,t) = mag_1tp.dat(:,:,:) ; end + ph(~isfinite(ph))=0; + mag(~isfinite(mag))=0; - % rescaling phase into [0,2*pi] phase range - ph = 2*pi*single(ph - min(vector(ph)))/single(max(vector(ph))-min(vector(ph))) ; - - clear ph_V % otherwise if there is already this structure in memory it does not like it - ph_V(1) = spm_vol((ph_files(1,:))) ; - ph_V(1).fname = 'ph.nii' ; - ph_V(1).dt = [16 0] ; - ph_V(1).descrip = ''; - ph_V(2) = ph_V(1) ; - ph_V(2).n = [2 1] ; - - clear mag_V % otherwise if there is already this structure in memory it does not like it - mag_V(1) = spm_vol(mag_files(end,:)); - mag_V(1).fname = sprintf('mag_TE%i.nii',size(TEs,2)); - mag_V(1).descrip = ''; - mag_V(2) = mag_V(1); - mag_V(2).n = [2 1]; - mag_ref = nifti(mag_files(end,:)); - mag_ref = repmat(mag_ref.dat(:,:,:),[1 1 1 2]); - for read_dir = 1:2 - - if size(TEs,2) < 6 % complex fit is only possible if at least 3 echoes per readout direction available - disp('calculating phase difference') - FM(:,:,:,read_dir) = angle(exp(1i*(ph(:,:,:,read_dir+2)-ph(:,:,:,read_dir)))) ; - else - disp('complex fitting phase') - compl = single(mag).*exp(-1i*ph); - [FM_1, ~, ~, ~] = Fit_ppm_complex_TE(compl(:,:,:,read_dir:2:end),TEs(read_dir:2:end)); - FM(:,:,:,read_dir) = FM_1 ; - end - - spm_write_vol(ph_V(read_dir), FM(:,:,:,read_dir)) ; - spm_write_vol(mag_V(read_dir), mag_ref(:,:,:,read_dir)) ; - - end - clear mag ph FM mag_ref - - + ph_merged_file = fullfile(output_fulldir, 'ph.nii') ; + mag_merged_file = fullfile(output_fulldir, 'mag.nii') ; + + createNifti(ph, ph_merged_file, ph_1tp.mat) + createNifti(mag, mag_merged_file, mag_1tp.mat) + + disp('phase unwrapping with ROMEO and:') disp('...removing global mean value') disp('......field map calculation') disp('.........saving quality map for masking') - TE = (TEs(3)-TEs(1))*10^3; % effective echo time difference after phase complex fitting in seconds - [~, FM_name,~] = fileparts(ph_V(1).fname) ; - FM_romeo_file = sprintf('%s_romeo.nii',FM_name) ; - if isunix - status = unix(sprintf('%s %s -m %s -o %s -t [%i,%i] -k nomask -g -q -B', para.romeo_command, ph_V(1).fname, mag_V(1).fname, FM_romeo_file, TE, TE)) ; - elseif ispc - status = system(sprintf('%s %s -m %s -o %s -t [%i,%i] -k nomask -g -q -B', para.romeo_command, ph_V(1).fname, mag_V(1).fname, FM_romeo_file, TE, TE)) ; - end + ph_romeo_file = fullfile(output_fulldir, 'ph_romeo.nii') ; + + status = system(sprintf('%s -p %s -m %s -o %s -t [%s] -q -B --weights 111111 --phase-offset-correction bipolar', para.romeo_command, ph_merged_file, mag_merged_file, ph_romeo_file, num2str(TEs))) ; + if status == 1 error('ROMEO did not run properly - check your installation path') end - - - %% field map rotation to scanner space + %% field map rotation % defining affine matrix in scanner space for data rotation to scanner % space with mantaining the same image origin (i.e. no translation) - - data_dim = size(ph_1tp.dat) ; - Z = spm_imatrix(ph_1tp.mat) ; + FM_V = spm_vol('B0.nii') ; + FM = nifti('B0.nii'); + disp('field map rotation to scanner space') + data_dim = size(FM.dat) ; + Z = spm_imatrix(FM.mat) ; pixdim = Z(7:9); - Maff_image = ph_1tp.mat ; - O = Maff_image\[0 0 0 1]' ; + mat_image = FM.mat ; + O = mat_image\[0 0 0 1]' ; O = O(1:3)' ; - Maff_scanner(1,:) = [0 0 pixdim(3) -pixdim(3)*O(3)] ; - Maff_scanner(2,:) = [pixdim(1) 0 0 -pixdim(1)*O(1)] ; - Maff_scanner(3,:) = [0 pixdim(2) 0 -pixdim(2)*O(2)] ; - Maff_scanner(4,:) = [0 0 0 1] ; + mat_scanner(1,:) = [0 0 pixdim(3) -pixdim(3)*O(3)] ; + mat_scanner(2,:) = [pixdim(1) 0 0 -pixdim(1)*O(1)] ; + mat_scanner(3,:) = [0 pixdim(2) 0 -pixdim(2)*O(2)] ; + mat_scanner(4,:) = [0 0 0 1] ; - FM = nifti('B0.nii') ; - fm_data = FM.dat(:,:,:) ; - fm_data(isnan(fm_data)) = 0 ; - FM_V = spm_vol('B0.nii'); - spm_write_vol(FM_V, fm_data) ; - clear FM - img2scanner_mat = Maff_image\Maff_scanner ; + img2scanner_mat = mat_image\mat_scanner ; FMrot = zeros(data_dim) ; data_dim_xy = data_dim(1:2); @@ -194,41 +151,34 @@ FMrot(:,:,slice) = spm_slice_vol(FM_V, img2scanner_mat*spm_matrix([0 0 slice]), data_dim_xy, -7) ; end - FMrot(isnan(FMrot)) = 0 ; - FMrot_V = FM_V ; - FMrot_V.mat = Maff_scanner ; - FMrot_V.fname = 'B0_rot.nii'; - spm_write_vol(FMrot_V, FMrot); + FMrot(~isfinite(FMrot)) = 0 ; + createNifti(FMrot, 'B0_rot.nii', mat_scanner) %% creating mask for QSM calculation disp('quality masking') qmask = nifti('quality.nii') ; - mask_V = spm_vol('quality.nii') ; qmask = qmask.dat(:,:,:) ; - qmask(isnan(qmask)) = 0 ; - qmask(qmask>0.3) = 1 ; - qmask(qmask<=0.3) = 0 ; + qmask(~isfinite(qmask)) = 0; + qmask(qmask>para.mask_thr) = 1 ; + qmask(qmask<=para.mask_thr) = 0 ; + % filling holes in the mask qmask = imfill(qmask,6,'holes') ; qmask = smooth3(qmask, 'gaussian') ; qmask(qmask>0.6) = 1 ; qmask(qmask<=0.6) = 0 ; + qmask = int16(qmask) ; - mask_V.dt = [spm_type('int16') 0] ; - mask_V.fname = 'mask.nii' ; - spm_write_vol(mask_V, qmask) ; + createNifti(qmask, 'mask.nii', mat_image) clear qmask + + mask_V = spm_vol('mask.nii') ; qmask_rot = zeros(data_dim) ; - for slice = 1 : data_dim(3) qmask_rot(:,:,slice) = spm_slice_vol(mask_V, img2scanner_mat*spm_matrix([0 0 slice]), data_dim_xy, -7) ; end - qmask_rot(isnan(qmask_rot)) = 0 ; - mask_V.mat = Maff_scanner ; - mask_V.fname = 'mask_rot.nii'; - spm_write_vol(mask_V, qmask_rot) ; - - + qmask_rot(~isfinite(qmask_rot)) = 0 ; + createNifti(int16(qmask_rot), 'mask_rot.nii', mat_scanner) %% SEPIA - background field removal and dipole inversion yielding final QSM @@ -236,7 +186,8 @@ B0 = para.B0 ; CF = B0*42.58*1e6; % imaging frequency, in Hz (B0*gyromagnetic_ratio*1e6) delta_TE = 1; % echo spacing, in second - we have already combined data, in such situation set to 1 - B0_dir = [0;1;0]; % main magnetic field direction, it's always [0,1,0] because the images are resliced so that 2nd dimension is aligned with B0 + TE = 1 ; + B0_dir = para.B0_dir; % main magnetic field direction, it's always [0,1,0] because the images are resliced so that 2nd dimension is aligned with B0 matrixSize = data_dim ; % image matrix size voxelSize = pixdim ; % spatial resolution of the data, in mm header_fullfile = fullfile(output_fulldir, 'header_sepia.mat') ; @@ -255,7 +206,7 @@ input(4).name = header_fullfile ; algorParam.bfr.refine = 0 ; - algorParam.bfr.erode_radius = 0 ; + algorParam.bfr.erode_radius = 1 ; algorParam.bfr.method = 'pdf' ; algorParam.bfr.tol = 0.1 ; algorParam.bfr.iteration = 50 ; @@ -273,7 +224,7 @@ elseif strcmp(algorParam.qsm.method , 'Star-QSM') - algorParam.qsm.padsize = ones(1,3)*12 ; + algorParam.qsm.padsize = ones(1,3)*6 ; end @@ -291,7 +242,7 @@ QSM = nifti('sepia_QSM.nii') ; QSM = QSM.dat(:,:,:) ; - scanner2img_mat = Maff_scanner\Maff_image ; + scanner2img_mat = mat_scanner\mat_image ; QSMinvrot = zeros(data_dim) ; QSM_V = spm_vol('sepia_QSM.nii'); for slice = 1 : data_dim(3) @@ -302,9 +253,13 @@ QSMinvrot_V.fname = 'sepia_QSM_invrot.nii'; spm_write_vol(QSMinvrot_V, QSMinvrot); - delete(sprintf('mag_TE%i.nii',size(TEs,2))); - delete('sepia_mask-qsm.nii.gz') - delete('sepia_QSM.nii.gz') + delete sepia_mask-qsm.nii.gz sepia_QSM.nii.gz + if strcmp(para.data_cleanup,'small') || strcmp(para.data_cleanup,'big') + delete mag.nii corrected_phase.nii mask.nii mask_rot.nii ph.nii ph_romeo.nii quality.nii sepia_local-field.nii.gz + end + if strcmp(para.data_cleanup,'big') + delete B0.nii B0_rot.nii + end sprintf('finished after %s' , secs2hms(toc(tstart))) end diff --git a/MPM_QSM_caller.m b/MPM_QSM_caller.m index e768d1a..7ea257c 100644 --- a/MPM_QSM_caller.m +++ b/MPM_QSM_caller.m @@ -35,9 +35,6 @@ % (also known as iterative Tikhonov), % may give more contrast than Star-QSM but is less robust to noise % 'Star-QSM' - is very robust to noise and quick -% calc_mean_qsm : 'yes' or 'no', if 'yes' mean QSM from the three -% MPM acquisitions will be calculated, ATTENTION: currently nor -% coregistration between the scans is implemented %%%% Inputs - directories, parameters and files specific to given contrast % ATTENTION: ensure only niftis you want to use are in that folder, with increasing echo numbering: @@ -72,73 +69,50 @@ % script created by Barbara Dymerska % @ UCL FIL Physics -% last modifications 09/09/2021 totstart = tic ; %%%%% USER PARAMETERS %%%%% para.romeo_command = '/your_path/romeo_linux_3.2.0/bin/romeo' ; para.in_root_dir = '/your/root/path' ; -para.out_root_dir = '/your/output/path'; +para.out_root_dir = '/your/output/path'; -para.B0 = 7; +para.B0 = 3; +para.B0_dir = [0;1;0]; % main magnetic field direction after reslicing the data para.dipole_inv = 'Star-QSM' ; -calc_mean_qsm = 'yes' ; +para.data_cleanup = 'big' ; % 'small' leaves B0 maps and QSMs, 'big' leaves only QSMs -% directories, parameters and files specific to given contrast -% ensure they are in the right order (PDw, T1w, MTw) -% otherwise mean PDw+T1w QSM will be something different when calc_mean_qsm = 'yes': for run = 1:3 switch run - case 1 %PDw - para.mag_dir = 'pdw_mfc_3dflash_v1k_0025' ; % folder with magnitude niftis - para.ph_dir = 'pdw_mfc_3dflash_v1k_0026' ; % folder with phase inftis - para.TEs = [2.2 4.58 6.96 9.34 11.72 14.1] ; % echo time in ms - para.output_dir = 'pdw_25_26' ; % output QSM directory for a specific MPM contrast + case 1 % PDw + para.mag_dir = 'pdw/mag/' ; % folder with magnitude niftis + para.ph_dir = 'pdw/ph/' ; % folder with phase inftis + para.TEs = 2.3*[1 2 3 4 5 6 7 8] ; % echo time in ms + para.output_dir = 'QSM_pdw' ; % output directory for a specific submeasurement from MPM + para.mask_thr = 0.15 ; % larger threshold smaller mask case 2 % T1w - para.mag_dir = 't1w_mfc_3dflash_v1k_0022' ; - para.ph_dir = 't1w_mfc_3dflash_v1k_0023' ; - para.TEs = [2.3 4.68 7.06 9.44 11.82 14.2] ; - para.output_dir = 't1w_22_23' ; + para.mag_dir = 't1w/mag/' ; + para.ph_dir = 't1w/ph/' ; + para.TEs = 2.3*[1 2 3 4 5 6 7 8] ; + para.output_dir = 'QSM_t1w' ; + para.mask_thr = 0.15 ; case 3 % MTw - para.mag_dir = 'mtw_mfc_3dflash_v1k_180deg_0031' ; - para.ph_dir = 'mtw_mfc_3dflash_v1k_180deg_0032' ; - para.TEs = [2.2 4.58 6.96 9.34] ; - para.output_dir = 'mtw_31_32' ; + para.mag_dir = 'mtw/mag/' ; + para.ph_dir = 'mtw/ph/' ; + para.TEs = 2.3*[1 2 3 4 5 6] ; + para.output_dir = 'QSM_mtw' ; + para.mask_thr = 0.1 ; end %%%%% END OF USER PARAMETERS %%%%% +MPM_QSM(para) ; - [QSM_V, QSM_all(:,:,:,run), QSMinvrot_V, QSM_all_invrot(:,:,:,run)] = MPM_QSM(para) ; - -end - - -if strcmp(calc_mean_qsm, 'yes') - - disp('calculating mean QSM across acquisitions') - QSM_pdw_t1w_mean = mean(QSM_all(:,:,:,1:2), 4) ; - QSM_all_mean = mean(QSM_all, 4) ; - - QSM_pdw_t1w_invrot_mean = mean(QSM_all_invrot(:,:,:,1:2), 4) ; - QSM_all_invrot_mean = mean(QSM_all_invrot, 4) ; - - - QSM_V.fname = fullfile(para.out_root_dir,'QSM_pdw_t1w_mean.nii') ; - spm_write_vol(QSM_V, QSM_pdw_t1w_mean) ; - QSM_V.fname = fullfile(para.out_root_dir,'QSM_all_mean.nii') ; - spm_write_vol(QSM_V, QSM_all_mean) ; - - - QSMinvrot_V.fname = fullfile(para.out_root_dir,'QSM_pdw_t1w_invrot_mean.nii') ; - spm_write_vol(QSMinvrot_V, QSM_pdw_t1w_invrot_mean) ; - QSMinvrot_V.fname = fullfile(para.out_root_dir,'QSM_all_invrot_mean.nii') ; - spm_write_vol(QSMinvrot_V, QSM_all_invrot_mean) ; end -sprintf('total processing finished after %s' , secs2hms(toc(totstart))) \ No newline at end of file +sprintf('total processing finished after %s' , secs2hms(toc(totstart))) +clear \ No newline at end of file diff --git a/README.md b/README.md index 969e60b..710bc01 100644 --- a/README.md +++ b/README.md @@ -19,7 +19,7 @@ QSM pipeline for Multi Parametric Mapping acquisitions 1. Download Zip with all the files from: - https://github.com/baskadym/MPM_QSM.git + https://github.com/fil-physics/MPM_QSM or clone it to your GitHub repository 2. Download compiled version of ROMEO either for windows or linux: diff --git a/other_dependencies/createNifti.m b/other_dependencies/createNifti.m new file mode 100644 index 0000000..eee8e34 --- /dev/null +++ b/other_dependencies/createNifti.m @@ -0,0 +1,88 @@ +%************************************************************************** +% +% Writes out data in nifti format. If the data is complex then it is +% written out as two files with the suffix _Magnitude & _Phase. +% +% MFC 14.12.2012 +% +% MFC 01.05.2013 Oddly I had dat(:,:,:) = data so changed now to +% N.dat(:,:,:) = data +% +% MFC 21.05.2013 Extra input parameter of spmMatrix so as to have +% orientational information. +% +%************************************************************************** + +function createNifti(data, fn, spmMatrix) + +% If a filename hasn't been passed in, where should I put it? +if nargin < 2 + [fn pn] = uiputfile('*.nii', 'Save As...'); + [~, fn] = fileparts(fn); +else + [pn fn] = fileparts(fn); + if isempty(pn) + pn = pwd; + end + +end + +% Is the data real or complex? +if isreal(data) + + dat = file_array; + dat.fname = [pn filesep fn '.nii']; + dat.dim = size(data); + dat.offset = ceil(348/8)*8; + + if isfloat(data) + dat.dtype = 'FLOAT32'; + N.descrip = 'FLOAT data'; + else + dat.dtype = 'UINT8'; + N.descrip = 'Non-float data'; + end + + % Create the NIFTI structure + N = nifti; + N.dat = dat; + if exist('spmMatrix', 'var') + N.mat = spmMatrix; + end + + + create(N); % Writes hdr + N.dat(:,:,:,:) = data; % Writes data + +else + + % General file_array fields: + dat = file_array; + dat.dim = size(data); + dat.dtype = 'FLOAT32'; + dat.offset = ceil(348/8)*8; + + % Create the NIFTI structure for the magnitude: + dat.fname = [pn filesep fn '_Magnitude.nii']; + N = nifti; + N.dat = dat; + N.descrip = 'Magnitude (from complex)'; + if exist('spmMatrix', 'var') + N.mat = spmMatrix; + end + create(N); % Writes hdr info + N.dat(:,:,:) = abs(data); % writes data to disk + + + % Create the NIFTI structure for the phase: + N = nifti; + dat.fname = [pn filesep fn '_Phase.nii']; + N.dat = dat; + N.descrip = 'Phase (from complex)'; + create(N); + N.dat(:,:,:) = angle(data); + if exist('spmMatrix', 'var') + N.mat = spmMatrix; + end + +end \ No newline at end of file