Skip to content

Support variable-TR, variable-flip-angle protocols in MCR-MWI processing and fitting #42

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
samuelmelke wants to merge 4 commits into
base: dev
Choose a base branch
from
+77 −21
Open

Support variable-TR, variable-flip-angle protocols in MCR-MWI processing and fitting #42

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
8e28306
Add vTR-qMRI as a dependency for variable-TR MCR-MWI
samuelmelke Sep 6, 2026
ffca023
Add variable-TR protocol support to MCR_GPUWorker
samuelmelke Sep 7, 2026
eef41b9
VFAprepWorker: use despot1_VTR for variable-TR protocols
samuelmelke Sep 9, 2026
4483315
remove development only flip subset code
samuelmelke Sep 9, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
75 changes: 57 additions & 18 deletions +qb/+workers/MCR_GPUWorker.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,8 @@
""
"MCR_GPUWorker implements the MCR framework on GPU hardware, combining complex multi-echo GRE data (VFA or MPM)"
"with coregistered B1 transmit field maps to estimate myelin water fraction (MWF) and other quantitative"
"microstructural parameters."
"microstructural parameters. Protocols with a variable TR and/or a variable number of echoes across flip angles "
"are fitted with gpuMCRMWI_VFAVTR (https://github.com/samuelmelke/vTR-qMRI)."
""
"Theoretical Framework:"
"----------------------"
Expand Down Expand Up @@ -102,31 +103,50 @@ function get_work_done(obj, workitem)
obj.logger.exception('%s expected %d brainmasks but got:%s', obj.name, length(ME4Dmag), sprintf(' %s', localfmask{:}))
end

% Read the protocol of each acquisition. This is done before allocating because acquisitions
% may differ in TR and in the number of echoes
for n = 1:length(ME4Dmag)
bfile = bids.File(ME4Dmag{n}); % For reading metadata, parsing entities, etc
FA(n) = bfile.metadata.FlipAngle; %#ok<AGROW>
TR(n) = bfile.metadata.RepetitionTime; %#ok<AGROW>
TE{n} = bfile.metadata.EchoTime(:); %#ok<AGROW>
end
nTE = cellfun(@numel, TE);
te_indexrange = [cumsum(nTE(:)) - nTE(:) + 1, cumsum(nTE(:))]; % First/last index of each acquisition along the echo dimension
isVTR = ~all(abs(TR - TR(1)) < 1e-6*TR(1)) || ... % Different repetition times
~all(nTE == nTE(1)) || ... % Different numbers of echoes
~all(cellfun(@(t) isequal(t, TE{1}), TE)); % Different echo times

% Load the data + metadata
V = spm_vol(ME4Dmag{1}); % For reading the 3D image dimensions
dims = [V(1).dim length(V) length(ME4Dmag)]; % Dimensions: [x,y,z,TE,FA]
V = spm_vol(ME4Dmag{1}); % For reading the 3D image dimensions
dims = [V(1).dim sum(nTE)]; % Dimensions: [x,y,z,echo], all acquisitions concatenated
img = single(NaN(dims));
unwrappedPhase = single(NaN(dims));
totalField = single(NaN(dims([1:3 5]))); % Dimensions: [x,y,z,FA]
totalField = single(NaN([dims(1:3) length(ME4Dmag)])); % Dimensions: [x,y,z,FA]
mask = true;
for n = 1:dims(5)
bfile = bids.File(ME4Dmag{n}); % For reading metadata, parsing entities, etc
img(:,:,:,:,n) = spm_read_vols(spm_vol(ME4Dmag{n}));
unwrappedPhase(:,:,:,:,n) = spm_read_vols(spm_vol(unwrapped{n}));
for n = 1:length(ME4Dmag)
idx = te_indexrange(n,1):te_indexrange(n,2);
img(:,:,:,idx) = spm_read_vols(spm_vol(ME4Dmag{n}));
unwrappedPhase(:,:,:,idx) = spm_read_vols(spm_vol(unwrapped{n}));
totalField(:,:,:,n) = spm_read_vols(spm_vol(fieldmap{n}));
mask = spm_read_vols(spm_vol(localfmask{n})) & mask;
FA(n) = bfile.metadata.FlipAngle; %#ok<AGROW>
end
B1 = spm_read_vols(spm_vol(char(TB1map_GRE)));
TR = bfile.metadata.RepetitionTime;
TE = bfile.metadata.EchoTime;

% Obtain the initial estimation of the initial B1 phase

% Obtain the initial estimation of the initial B1 phase (NB: img is still 4D here, i.e. with all acquisitions concatenated)
img = img .* exp(1i*unwrappedPhase);
mask = mask & all(~isnan(img), [4 5]);
pini = squeeze(unwrappedPhase(:,:,:,1,:)) - 2*pi*totalField .* TE(1);
mask = mask & all(~isnan(img), 4);
TE1 = reshape(cellfun(@(te) te(1), TE), 1, 1, 1, []); % First echo time of each acquisition
pini = unwrappedPhase(:,:,:,te_indexrange(:,1)) - 2*pi*totalField .* TE1; % Dimensions: [x,y,z,FA]
pini = polyfit3D_NthOrder(double(mean(pini(:,:,:,1:end-1), 4)), mask, 6);

clear unwrappedPhase % not used after this line

if ~isVTR
img = reshape(img, [dims(1:3) nTE(1) length(ME4Dmag)]); % Back to [x,y,z,TE,FA] for gpuMCRMWI
dims = size(img);
end

% Construct the fixed parameters and extra data for the MCR model
fixed_params = obj.config.MCR_GPUWorker.fixed_params;
fixed_params.B0 = bfile.metadata.MagneticFieldStrength;
Expand All @@ -135,9 +155,28 @@ function get_work_done(obj, workitem)
extraData.pini = pini;
extraData.b1 = B1;

% Estimate the MCR model
objGPU = gpuMCRMWI(TE, TR, FA, fixed_params);
askadam_mcr = objGPU.estimate(img, mask, extraData, obj.config.MCR_GPUWorker.fitting);
% Variable-TR protocols need the echo index ranges and cannot use the EPG-X networks
% We need to override isEPG for VTR, and mutating obj.config would leak the change into other subjects processed in the same session
fitting = obj.config.MCR_GPUWorker.fitting;
if isVTR
extraData.te_indexrange = te_indexrange;
if fitting.isEPG
obj.logger.warning('%s: the EPG-X networks assume a single TR for all flip angles, using the Bloch-McConnell solution instead', obj.name)
fitting.isEPG = false;
end
end

% Estimate the MCR model (variable-TR protocols need the vTR-qMRI implementation)
if isVTR
if ~exist('gpuMCRMWI_VFAVTR', 'class')
obj.logger.exception('%s found a variable-TR protocol but gpuMCRMWI_VFAVTR is not on the MATLAB-path.\nPossible solution:\ngit submodule update --init dependencies/vTR-qMRI', obj.name)
end
obj.logger.info('%s detected a variable protocol (TR = [%s] ms), using gpuMCRMWI_VFAVTR', obj.name, num2str(TR*1e3, ' %.1f'))
objGPU = gpuMCRMWI_VFAVTR(TE, TR, FA, fixed_params);
else
objGPU = gpuMCRMWI(TE{1}, TR(1), FA, fixed_params);
end
askadam_mcr = objGPU.estimate(img, mask, extraData, fitting);

% Extract and save the output data
V(1).dim = dims(1:3);
Expand Down
18 changes: 15 additions & 3 deletions +qb/+workers/VFAprepWorker.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -198,7 +198,7 @@ function make_syntheticT1_M0(obj, bfilter)

% Compute T1 and M0 maps
obj.logger.info("--> Running despot1 to compute T1 and M0 maps from: " + VFA_e1{1})
flipangles = [];
flipangles = []; TRs = [];
VFAimg = NaN([Vref.dim length(VFA_e1)]);
for n = 1:length(VFA_e1)
VFAn = spm_vol(VFA_e1{n});
Expand All @@ -213,12 +213,23 @@ function make_syntheticT1_M0(obj, bfilter)
end
metadata = bids.File(VFA_e1{n}).metadata;
flipangles(n) = metadata.FlipAngle;
TRs(n) = metadata.RepetitionTime;
end

isVTR = ~all(abs(TRs - TRs(1)) < 1e-6*TRs(1));
if isVTR
if ~exist('despot1_VTR', 'class')
obj.logger.exception('%s found a variable-TR protocol but despot1_VTR is not on the MATLAB-path.\nPossible solution:\ngit submodule update --init dependencies/vTR-qMRI', obj.name)
end
obj.logger.info('%s detected a variable-TR protocol (TR = [%s] ms), using despot1_VTR', obj.name, num2str(TRs*1e3, ' %.2f'))
[T1, M0] = despot1_VTR(TRs, flipangles).estimate(VFAimg);
else
[T1, M0] = despot1_mapping(VFAimg, flipangles, TRs(1));
end
[T1, M0] = despot1_mapping(VFAimg, flipangles, metadata.RepetitionTime);

% Save T1w-like images in the work directory
for n = 1:length(VFA_e1)
T1w = M0 .* GRESignal(flipangles(n), metadata.RepetitionTime, T1);
T1w = M0 .* GRESignal(flipangles(n), TRs(n), T1);
T1w(~isfinite(T1w)) = 0;
bfile = obj.bfile_set(VFA_e1{n}, obj.bidsfilter.syntheticT1);
bfile.metadata.Sources = {['bids::' bfile.bids_path '/' bfile.filename]};
Expand All @@ -230,6 +241,7 @@ function make_syntheticT1_M0(obj, bfilter)
bfile = obj.bfile_set(Vref.fname, obj.bidsfilter.M0map_echo1);
bfile.metadata.Sources = strrep(VFA_e1, extractBefore(VFA_e1{1}, bfile.bids_path), 'bids::');
bfile.metadata.FlipAngle = flipangles;
bfile.metadata.RepetitionTime = TRs;
obj.logger.verbose("-> Saving M0 map " + fullfile(bfile.bids_path, bfile.filename))
write_vol(Vref, M0, bfile);
end
Expand Down
4 changes: 4 additions & 0 deletions .gitmodules
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,3 +42,7 @@
path = dependencies/Tensor-MP-PCA
url = https://github.com/Neurophysics-CFIN/Tensor-MP-PCA.git
branch = main
[submodule "dependencies/vTR-qMRI"]

@marcelzwiers marcelzwiers Sep 15, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Question, did you manually edit the .gitmodules file or did you use:

$ git submodule add -b main https://github.com/samuelmelke/vTR-qMRI.git dependencies/vTR-qMRI

I'm asking because I had created a lot of issues myself in the past when I edited the file manually...

path = dependencies/vTR-qMRI
url = https://github.com/samuelmelke/vTR-qMRI.git
branch = main
1 change: 1 addition & 0 deletions dependencies/vTR-qMRI
Open in desktop
Submodule vTR-qMRI added at 22b05c
Loading