JointR1R2star

Joint R2-R2* single compartment relaxometry using variable flip angle and multi-echo GRE data

gpuJointR1R2star

With askadam solver

Usage

obj = gpuJointR1R2starMapping(te,tr,fa);
[out] = obj.estimate( data, mask, extraData, fitting);

Model parameters

% M0    : Proton density weighted signal
% R1    : (=1/T1) in s^-1
% R2star: R2* in s^-1
model_params    = {'M0';'R1';'R2star'};
ub              = [  2;  10;  200];
lb              = [  0; 0.1;  0.1];
startpoint      = [   1;   1;  30];

I/O overview

obj = gpuJointR1R2starMapping(te,tr,fa);

Input	Description
te	1xNte echo time [s]
tr	repetition time [s]
fa	1xNfa flip angle vector [degree]

[out] = obj.estimate( data, mask, extraData, fitting);

Input	Description
data	5D MRI data, [x,y,z,t,fa]
mask	3D mask, [x,y,z]
extraData	Structure array with additional data
extraData.b1	3D B1+ map [ratio], [x,y,z]
fitting	Structure array for model parameter estimation
fitting.optimiser	Algorithm for parameter update, ‘adam’ (default) | ‘sgdm’ | ‘rmsprop’
fitting.isdisplay	boolean, display optimisation process in graphic plot
fitting.convergenceValue	tolerance in loss gradient to stop the optimisation
fitting.convergenceWindow	# of elements in which ‘convergenceValue’ is computed
fitting.iteration	maximum # of optimisation iterations
fitting.initialLearnRate	initial learn rate of Adam optimiser
fitting.tol	tolerance in loss
fitting.lambda	regularisation parameter(s)
fitting.regmap	model parameter(s) in which regularisation is applied,
fitting.TVmode	Mode for total variation (TV) regularisation, ‘2D’ | ‘3D’
fitting.lossFunction	loss function, ‘L1’ | ‘L2’ | ‘huber’ | ‘mse’
fitting.isWeighted	is cost weighted, true|false, default = true
fitting.weightMethod	Weighting method, ‘1stecho’ (default) | ‘norm’
fitting.weightPower	power order of the weight, default = 2

Example

Example script for noise propagation:

addpath(genpath('../../gacelle/'));
clear;

%% Simulate data

% for reproducibility
seed        = 23439; rng(seed);
Nsample     = 1e3;  % #voxel
SNR         = 100;

Nt  = 15; 
TE1 = 1.5e-3;
tr  = 45e-3;
t   = linspace(TE1, tr-3e-3, Nt);
fa  = [5, 10, 20, 30, 40, 50, 70];
Nfa = numel(fa);

% Parameter raneg for forward simulation
M0_range        = [0.5, 2];
R1_range        = [0.2, 2.5];
R2star_range    = [5,   60];

% generate ground truth
M0_GT       = single(rand(1,Nsample) * diff(M0_range) + min(M0_range) );
R1_GT       = single(rand(1,Nsample) * diff(R1_range)  + min(R1_range));
R2star_GT   = single(rand(1,Nsample) * diff(R2star_range)  + min(R2star_range));

% forward signal simulation
pars.M0     = M0_GT;
pars.R1     = R1_GT;
pars.R2star = R2star_GT;

b1      = ones(size(M0_GT)); extraData.b1 = b1;
objGPU  = gpuJointR1R2starMapping(t,tr,fa);
s       = gather(extractdata(objGPU.FWD(pars, extraData)));
s       = permute(reshape(s,[Nt, Nfa, Nsample]),[3 4 5 1 2]);
mask    = ones(size(s,1:3),'logical');

% Let's assume Gaussian noise for simplicity
noiseLv = (M0_GT.')./SNR;
s       = s + randn(size(s)) .* noiseLv;

%% askadam estimation
fitting                     = [];
fitting.optimiser           = 'adam';
fitting.iteration           = 4000;
fitting.initialLearnRate    = 0.001;
fitting.convergenceValue    = 1e-8;
fitting.lossFunction        = 'l1';
fitting.tol                 = 1e-8;
fitting.isdisplay           = false;
fitting.start               = 'prior';   
extraData                   = [];
extraData.b1                = b1.';

objGPU  = gpuJointR1R2starMapping(t,tr,fa);
out     = objGPU.estimate(s, mask, extraData, fitting);

%% askadam estimation
fitting                     = [];
fitting.optimiser           = 'adam';
fitting.iteration           = 4000;
fitting.initialLearnRate    = 0.001;
fitting.convergenceValue    = 1e-8;
fitting.lossFunction        = 'l1';
fitting.tol                 = 1e-8;
fitting.isdisplay           = false;
fitting.start               = 'prior';   
fitting.isWeighted          = true;
fitting.weightMethod        = '1stecho';
fitting.weightPower         = 2;
extraData                   = [];
extraData.b1                = b1.';

objGPU          = gpuJointR1R2starMapping(t,tr,fa);
out_weighted    = objGPU.estimate(s, mask, extraData, fitting);

%% starting position
objGPU  = gpuJointR1R2starMapping(t,tr,fa);
pars0   = objGPU.estimate_prior(s, mask, extraData);

%% plot result
figure(99);
field = fieldnames(pars);
tiledlayout(2,numel(field));
for k = 1:numel(field)
    nexttile;
    scatter(pars.(field{k}),pars0.(field{k}),5,'filled','MarkerFaceAlpha',.4);hold on
    scatter(pars.(field{k}),out.final.(field{k}),5,'filled','MarkerFaceAlpha',.4);
    h = refline(1);
    h.Color = 'k';
    title(field{k});
    xlabel('GT');ylabel('Fitted');
end
for k = 1:numel(field)
    nexttile;
    scatter(pars.(field{k}),pars0.(field{k}),5,'filled','MarkerFaceAlpha',.4);hold on
    scatter(pars.(field{k}),out_weighted.final.(field{k}),5,'filled','MarkerFaceAlpha',.4);
    h = refline(1);
    h.Color = 'k';
    title([field{k} ' w-P2']);
    xlabel('GT');ylabel('Fitted');
end

legend('Start','Fitted','Location','northwest');

Example script for real data:

addpath(genpath('/autofs/space/linen_001/users/kwokshing/tools/gacelle'))
addpath('/autofs/space/linen_001/users/kwokshing/tools/sepia/sepia_master');
sepia_addpath;

%% Subject info and directories
subj_label = 'sub-022';

bids_dir            = '/autofs/cluster/connectome2/Bay8_C2/bids/';
derivatives_dir     = fullfile(bids_dir,'derivatives/');

jointR1R2star_dir       = fullfile(derivatives_dir,'jointR1R2star',subj_label); 
SEPIA_dir               = fullfile(derivatives_dir,'sepia',subj_label);
processed_vibe_dir      = fullfile(derivatives_dir,'processed_vibe',subj_label);

file_list = dir(fullfile(processed_vibe_dir,'*rec-offline*space-MEGRE.mat'));
flip_angle = zeros(1,numel(file_list));
for kfile = 1:numel(file_list)
    load(fullfile(file_list(kfile).folder,file_list(kfile).name),'FA');
    flip_angle(kfile) = FA;
end
flip_angle = sort(flip_angle,'ascend');

%% load data
counter         = 0;
img             = [];
sepia_header    = [];
unwrappedPhase  = [];
totalField      = [];
fa              = zeros(1,length(flip_angle));
for kfa = 1:length(flip_angle)
    counter = counter + 1;
    
     % general GRE basename
    acq_label   = strcat('acq-',['FA' num2str(flip_angle(kfa))]);
    prefix      = strcat(subj_label,'_',acq_label,'_rec-offline');

    % magnitude nifti image filename
    magn_fn         = dir(fullfile(processed_vibe_dir,strcat(prefix,'*part-mag*MEGRE_space-MEGRE.nii*')));
    sepia_header_fn = dir(fullfile(processed_vibe_dir,strcat(prefix,'*space-MEGRE.mat')));

    nii                 = load_untouch_nii(fullfile(magn_fn.folder, magn_fn.name));
    img                 = cat(5,img,nii.img);
    sepia_header{kfa}   = load(fullfile(sepia_header_fn.folder, sepia_header_fn.name));

    fa(kfa)  = sepia_header{kfa}.FA;
    tr      = sepia_header{kfa}.TR;

end
te = sepia_header{end}.TE;

% B1 info
true_flip_angle_fn      = dir(fullfile(processed_vibe_dir,'*acq-famp*TB1TFL*space-MEGRE.nii*'));
true_flip_angle_json    = dir(fullfile(processed_vibe_dir,'*acq-famp*TB1TFL*.json'));

true_flip_angle         = load_nii_img_only( fullfile( true_flip_angle_fn.folder, true_flip_angle_fn.name));
b1_header               = jsondecode( fileread( fullfile( true_flip_angle_json.folder,true_flip_angle_json.name)));

b1                      = true_flip_angle / 10 / b1_header.FlipAngle;

mask_fn         = dir(fullfile(processed_vibe_dir,strcat(subj_label,'*acq-FA10*rec-offline*mask_brain*space-MEGRE.nii*')));
mask_filename   = fullfile(mask_fn.folder, mask_fn.name);
mask            = load_nii_img_only(mask_filename);

clear true_flip_angle

%% Prepare data fot batch processing
% Magnitude fitting
% setup algorithm parameters
fitting                     = [];
% fitting option
fitting.iteration           = 4000;
fitting.convergenceWindow   = 25;     % at least 100 iterations and more robust convergance estimate
fitting.convergenceValue    = 1e-8;
fitting.initialLearnRate    = 0.001;
fitting.decayRate           = 0;
fitting.tol                 = 1e-3;
% residual option
fitting.isWeighted          = false;
extradata = [];
extradata.b1 = b1;

obj     = gpuJointR1R2starMapping(te,tr,fa);
[out]   = obj.estimate(img, mask, extradata, fitting);

%% Prepare data fot batch processing
% Magnitude fitting
% setup algorithm parameters
fitting                     = [];
% fitting option
fitting.iteration           = 4000;
fitting.convergenceWindow   = 25;     % at least 100 iterations and more robust convergance estimate
fitting.convergenceValue    = 1e-8;
fitting.initialLearnRate    = 0.001;
fitting.decayRate           = 0;
fitting.tol                 = 1e-4;
% residual option
fitting.isWeighted          = true;
fitting.weightPower         = 2;
fitting.weightMethod        = '1stecho';
extradata = [];
extradata.b1 = b1;

obj             = gpuJointR1R2starMapping(te,tr,fa);
[out_weighted]   = obj.estimate(img, mask, extradata, fitting);

gpuJointR1R2starmcmc

With MCMC solver

Usage

obj     = gpuJointR1R2starMappingmcmc(te,tr,fa);
[out]   = obj.estimate( data, mask, extraData, fitting);

Model parameters

% M0    : Proton density weighted signal
% R1    : (=1/T1) in s^-1
% R2star: R2* in s^-1
model_params    = {'M0';'R1';'R2star';'noise'};
ub              = [   2;  10;     200;    0.1];
lb              = [   0; 0.1;     0.1;  0.001];
startpoint      = [   1;   1;      30;   0.05];
step            = [0.01;0.01;       1;  0.005];

I/O overview

obj = gpuJointR1R2starMappingmcmc(te,tr,fa);

Input	Description
te	1xNte echo time [s]
tr	repetition time [s]
fa	1xNfa flip angle vector [degree]

[out] = obj.estimate( data, mask, extraData, fitting);

Input	Description
data	5D MRI data, [x,y,z,t,fa]
mask	3D mask, [x,y,z]
extraData	Structure array with additional data
extraData.b1	3D B1+ map [ratio], [x,y,z]
fitting	Structure array for model parameter estimation
fitting.algorithm	MCMC algorithm, ‘MH’ (Metropolis-Hastings)|’GW’ (Affline-invariant ensemble)
fitting.iteration	# MCMC iterations
fitting.repetition	# repetition of MCMC proposal
fitting.thinning	sampling interval between iterations
fitting.burnin	iterations to be discarded at the beginning, if >1, the exact number will be used; else iteration*burnin
fitting.xStepSize	step size of model parameter in MCMC proposal, same size and order as ‘model_params’ (‘MH’ only)
fitting.StepSize	step size for ‘GW’ in MCMC proposal (‘GW’ only)
fitting.Nwalker	# random walkers (‘GW’ only)
fitting.metric	cell variable, metric(s) derived from posterior distribution, ‘mean’|’std’|’median’|’iqr’ (can be multiple)
fitting.isWeighted	is cost weighted, true|false, default = true
fitting.weightMethod	Weighting method, ‘1stecho’ (default) | ‘norm’
fitting.weightPower	power order of the weight, default = 2

Output	Description
out	structure contains optimisation result
out.posterior	structure contains MCMC posterior samples
out.posterior.(model_params{k})	Model parameter MCMC posterior samples, masked and unshaped for memory preservation
out.{metric}.(model_params{k})	Posterior statistics chosen in fitting.metric

Example

Example script for noise propagation:

addpath(genpath('../../gacelle/'));
clear;

%% Simulate data

% for reproducibility
seed        = 23439; rng(seed);
Nsample     = 1e3;  % #voxel
SNR         = 100;

Nt  = 15; 
TE1 = 1.5e-3;
tr  = 45e-3;
t   = linspace(TE1, tr-3e-3, Nt);
fa  = [5, 10, 20, 30, 40, 50, 70];
Nfa = numel(fa);

% Parameter raneg for forward simulation
M0_range        = [0.5, 2];
R1_range        = [0.2, 4];
R2star_range    = [5, 100];

% generate ground truth
M0_GT       = single(rand(1,Nsample) * diff(M0_range) + min(M0_range) );
R1_GT       = single(rand(1,Nsample) * diff(R1_range)  + min(R1_range));
R2star_GT   = single(rand(1,Nsample) * diff(R2star_range)  + min(R2star_range));

% forward signal simulation
pars.M0     = M0_GT;
pars.R1     = R1_GT;
pars.R2star = R2star_GT;

b1      = ones(size(M0_GT)); extraData.b1 = b1;
objGPU  = gpuJointR1R2starMappingmcmc(t,tr,fa);
s       = gather(objGPU.FWD(pars, [], extraData));
s       = permute(reshape(s,[Nt, Nfa, Nsample]),[3 4 5 1 2]);
mask    = ones(size(s,1:3),'logical');

% Let's assume Gaussian noise for simplicity
noiseLv = (M0_GT.')./SNR;
s       = s + randn(size(s)) .* noiseLv;

%% askadam estimation
fitting             = [];
fitting.algorithm   = 'GW';
fitting.Nwalker     = 50;
fitting.StepSize    = 2;
fitting.iteration   = 1e4;
fitting.thinning    = 10;        % Sample every 10 iteration
fitting.metric      = {'mean','std'};
fitting.burnin      = 0.1;       % 10% burn-in
fitting.start       = 'default';
extraData           = [];
extraData.b1        = b1.';

objGPU  = gpuJointR1R2starMappingmcmc(t,tr,fa);
out     = objGPU.estimate(s, mask, extraData, fitting);

%% plot result
figure(2)
field = fieldnames(pars);
tiledlayout(1,numel(field));
for k = 1:numel(field)
    nexttile;
    scatter(pars.(field{k}),out.mean.(field{k}),5,'filled','MarkerFaceAlpha',.4);
    h = refline(1);
    h.Color = 'k';
    title(field{k});
    xlabel('GT');ylabel('Fitted');
end

Example script for real data:

addpath(genpath('/autofs/space/linen_001/users/kwokshing/tools/gacelle'))
addpath('/autofs/space/linen_001/users/kwokshing/tools/sepia/sepia_master');
sepia_addpath;

%% Subject info and directories
subj_label = 'sub-022';

bids_dir            = '/autofs/cluster/connectome2/Bay8_C2/bids/';
derivatives_dir     = fullfile(bids_dir,'derivatives/');

jointR1R2star_dir       = fullfile(derivatives_dir,'jointR1R2star',subj_label); 
SEPIA_dir               = fullfile(derivatives_dir,'sepia',subj_label);
processed_vibe_dir      = fullfile(derivatives_dir,'processed_vibe',subj_label);

file_list = dir(fullfile(processed_vibe_dir,'*rec-offline*space-MEGRE.mat'));
flip_angle = zeros(1,numel(file_list));
for kfile = 1:numel(file_list)
    load(fullfile(file_list(kfile).folder,file_list(kfile).name),'FA');
    flip_angle(kfile) = FA;
end
flip_angle = sort(flip_angle,'ascend');

%% load data
counter         = 0;
img             = [];
sepia_header    = [];
unwrappedPhase  = [];
totalField      = [];
fa              = zeros(1,length(flip_angle));
for kfa = 1:length(flip_angle)
    counter = counter + 1;
    
     % general GRE basename
    acq_label   = strcat('acq-',['FA' num2str(flip_angle(kfa))]);
    prefix      = strcat(subj_label,'_',acq_label,'_rec-offline');

    % magnitude nifti image filename
    magn_fn         = dir(fullfile(processed_vibe_dir,strcat(prefix,'*part-mag*MEGRE_space-MEGRE.nii*')));
    sepia_header_fn = dir(fullfile(processed_vibe_dir,strcat(prefix,'*space-MEGRE.mat')));

    nii                 = load_untouch_nii(fullfile(magn_fn.folder, magn_fn.name));
    img                 = cat(5,img,nii.img);
    sepia_header{kfa}   = load(fullfile(sepia_header_fn.folder, sepia_header_fn.name));

    fa(kfa)  = sepia_header{kfa}.FA;
    tr      = sepia_header{kfa}.TR;

end
te = sepia_header{end}.TE;

% B1 info
true_flip_angle_fn      = dir(fullfile(processed_vibe_dir,'*acq-famp*TB1TFL*space-MEGRE.nii*'));
true_flip_angle_json    = dir(fullfile(processed_vibe_dir,'*acq-famp*TB1TFL*.json'));

true_flip_angle         = load_nii_img_only( fullfile( true_flip_angle_fn.folder, true_flip_angle_fn.name));
b1_header               = jsondecode( fileread( fullfile( true_flip_angle_json.folder,true_flip_angle_json.name)));

b1                      = true_flip_angle / 10 / b1_header.FlipAngle;

mask_fn         = dir(fullfile(processed_vibe_dir,strcat(subj_label,'*acq-FA10*rec-offline*mask_brain*space-MEGRE.nii*')));
mask_filename   = fullfile(mask_fn.folder, mask_fn.name);
mask            = load_nii_img_only(mask_filename);

clear true_flip_angle

%% Prepare data fot batch processing
% Magnitude fitting
% setup algorithm parameters
fitting                 = [];
fitting.algorithm       = 'GW';
fitting.Nwalker         = 50;
fitting.StepSize        = 2;
fitting.iteration       = 1e4;
fitting.thinning        = 10;        % Sample every 20 iteration
fitting.metric          = {'median','iqr'};
fitting.burnin          = 0.1;       % 10% burn-in
% residual option
fitting.isWeighted      = true;
fitting.weightPower     = 2;
fitting.weightMethod    = '1stecho';
extradata = [];
extradata.b1 = b1;

obj     = gpuJointR1R2starMappingmcmc(te,tr,fa);
[out]   = obj.estimate(img, mask, extradata, fitting);