Skip to content

Commit

Permalink
Merge pull request #2 from baskadym/main
Browse files Browse the repository at this point in the history 
Release for 3T qMAP-PD
  • Loading branch information
@mfcallaghan
mfcallaghan authored Jun 16, 2022
2 parents 54d024b + 5ec7dc0 commit cf7d68c
Show file tree
Hide file tree
Showing 5 changed files with 171 additions and 154 deletions.
2 changes: 1 addition & 1 deletion License
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,5 +13,5 @@
You should have received a copy of the GNU General Public License
along with this program. If not, see <https://www.gnu.org/licenses/>.

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.
161 changes: 58 additions & 103 deletions MPM_QSM.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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)

Expand All @@ -82,161 +78,116 @@
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);

for slice = 1 : data_dim(3)
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

disp('creating SEPIA header')
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') ;
Expand All @@ -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 ;
Expand All @@ -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

Expand All @@ -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)
Expand All @@ -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

Expand Down
72 changes: 23 additions & 49 deletions MPM_QSM_caller.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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)))
sprintf('total processing finished after %s' , secs2hms(toc(totstart)))
clear
2 changes: 1 addition & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
Loading

0 comments on commit cf7d68c

Please sign in to comment.