From 8e28306c947688e2580562193ec73eec77962412 Mon Sep 17 00:00:00 2001 From: Samuel Date: Mon, 7 Sep 2026 00:27:09 +0200 Subject: [PATCH 1/4] Add vTR-qMRI as a dependency for variable-TR MCR-MWI --- .gitmodules | 4 ++++ dependencies/vTR-qMRI | 1 + 2 files changed, 5 insertions(+) create mode 160000 dependencies/vTR-qMRI diff --git a/.gitmodules b/.gitmodules index f61abb90..a064eb80 100644 --- a/.gitmodules +++ b/.gitmodules @@ -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"] + path = dependencies/vTR-qMRI + url = https://github.com/samuelmelke/vTR-qMRI.git + branch = main diff --git a/dependencies/vTR-qMRI b/dependencies/vTR-qMRI new file mode 160000 index 00000000..22b05c8e --- /dev/null +++ b/dependencies/vTR-qMRI @@ -0,0 +1 @@ +Subproject commit 22b05c8e037d40866b49e0c4d7ed8e17c1a9f19a From ffca023dc1a04a69d8cf8bbb84eac79050eee23c Mon Sep 17 00:00:00 2001 From: Samuel Date: Mon, 7 Sep 2026 03:48:52 +0200 Subject: [PATCH 2/4] Add variable-TR protocol support to MCR_GPUWorker --- +qb/+workers/MCR_GPUWorker.m | 82 ++++++++++++++++++++++++++++-------- 1 file changed, 64 insertions(+), 18 deletions(-) diff --git a/+qb/+workers/MCR_GPUWorker.m b/+qb/+workers/MCR_GPUWorker.m index a93bcab9..8bc8baf4 100644 --- a/+qb/+workers/MCR_GPUWorker.m +++ b/+qb/+workers/MCR_GPUWorker.m @@ -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:" "----------------------" @@ -87,6 +88,13 @@ function get_work_done(obj, workitem) localfmask = obj.ask_team('localfmask'); % Multiple FA-images per run TB1map_GRE = obj.ask_team('TB1map_GRE'); % Single image per run + % >>> TEMPORARY for development testing + %testflips = [1 2 5]; % TR44/FA5, TR44/FA15, TR44/FA70 + testflips = [21 12 5]; % TR11/FA5, TR28/FA15, TR44/FA70 + keep = cellfun(@(f) ismember(str2double(bids.File(f).entities.flip), testflips), ME4Dmag); + ME4Dmag = ME4Dmag(keep); unwrapped = unwrapped(keep); + fieldmap = fieldmap(keep); localfmask = localfmask(keep); + % Check the number of items we got: TODO: FIXME: multi-run acquisitions if numel(unique([length(unwrapped), length(fieldmap)])) > 1 obj.logger.exception('%s received an ambiguous number of ME4Dmag, unwrapped or fieldmaps:%s', obj.name, ... @@ -102,31 +110,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 + TR(n) = bfile.metadata.RepetitionTime; %#ok + TE{n} = bfile.metadata.EchoTime(:); %#ok + 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 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; @@ -135,9 +162,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); From eef41b9b1dcd956f567af0abe6fcd77bd464a012 Mon Sep 17 00:00:00 2001 From: Samuel Date: Wed, 9 Sep 2026 20:16:54 +0200 Subject: [PATCH 3/4] VFAprepWorker: use despot1_VTR for variable-TR protocols --- +qb/+workers/VFAprepWorker.m | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/+qb/+workers/VFAprepWorker.m b/+qb/+workers/VFAprepWorker.m index 6aba6763..7d4bffb7 100644 --- a/+qb/+workers/VFAprepWorker.m +++ b/+qb/+workers/VFAprepWorker.m @@ -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}); @@ -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]}; @@ -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 From 4483315d9baa55c98f737cdca469a09b7d278166 Mon Sep 17 00:00:00 2001 From: samuelmelke Date: Wed, 9 Sep 2026 21:28:30 +0200 Subject: [PATCH 4/4] remove development only flip subset code --- +qb/+workers/MCR_GPUWorker.m | 7 ------- 1 file changed, 7 deletions(-) diff --git a/+qb/+workers/MCR_GPUWorker.m b/+qb/+workers/MCR_GPUWorker.m index 8bc8baf4..985925ea 100644 --- a/+qb/+workers/MCR_GPUWorker.m +++ b/+qb/+workers/MCR_GPUWorker.m @@ -88,13 +88,6 @@ function get_work_done(obj, workitem) localfmask = obj.ask_team('localfmask'); % Multiple FA-images per run TB1map_GRE = obj.ask_team('TB1map_GRE'); % Single image per run - % >>> TEMPORARY for development testing - %testflips = [1 2 5]; % TR44/FA5, TR44/FA15, TR44/FA70 - testflips = [21 12 5]; % TR11/FA5, TR28/FA15, TR44/FA70 - keep = cellfun(@(f) ismember(str2double(bids.File(f).entities.flip), testflips), ME4Dmag); - ME4Dmag = ME4Dmag(keep); unwrapped = unwrapped(keep); - fieldmap = fieldmap(keep); localfmask = localfmask(keep); - % Check the number of items we got: TODO: FIXME: multi-run acquisitions if numel(unique([length(unwrapped), length(fieldmap)])) > 1 obj.logger.exception('%s received an ambiguous number of ME4Dmag, unwrapped or fieldmaps:%s', obj.name, ...