From 4bf3cc4845f296d124972ee16a815a542288cd7f Mon Sep 17 00:00:00 2001 From: s-rauh Date: Tue, 21 Jul 2026 10:52:27 +0200 Subject: [PATCH 1/5] Add Matlab IVIM and IVIM-DTI fitting pipeline SR_LUMC to original/fitting code contributions. --- src/original/fitting/SR_LUMC/__init__.py | 0 .../diffusion_mri_tools-main/README.md | 24 ++ .../analyze_data/calc_dti_parameters.m | 89 +++++ .../diffusion_mri_tools-main/fit_dti.m | 162 +++++++++ .../diffusion_mri_tools-main/fit_ivim.m | 230 +++++++++++++ .../diffusion_mri_tools-main/fit_ivimdti.m | 312 ++++++++++++++++++ .../helper_fun/bval_scaling.m | 21 ++ .../helper_fun/calc_bmat.m | 28 ++ .../helper_fun/delta_fun.m | 17 + .../helper_fun/lower_triangular2tensor.m | 41 +++ .../helper_fun/mask_diffdata.m | 34 ++ .../helper_fun/norm_diffdata.m | 43 +++ .../helper_fun/remove_zeros.m | 25 ++ .../helper_fun/reshape_diffdata_for_fit.m | 49 +++ 14 files changed, 1075 insertions(+) create mode 100644 src/original/fitting/SR_LUMC/__init__.py create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/README.md create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/analyze_data/calc_dti_parameters.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_dti.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_ivim.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_ivimdti.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/bval_scaling.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/calc_bmat.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/delta_fun.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/lower_triangular2tensor.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/mask_diffdata.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/norm_diffdata.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/remove_zeros.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/reshape_diffdata_for_fit.m diff --git a/src/original/fitting/SR_LUMC/__init__.py b/src/original/fitting/SR_LUMC/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/README.md b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/README.md new file mode 100644 index 00000000..f3be8839 --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/README.md @@ -0,0 +1,24 @@ +# Diffusion MRI Tools + +Diffusion MRI Tools contains Matlab code to analyze and process diffusion MRI data. + +Author: Susi Rauh (s.s.rauh@lumc.nl) + +Three types of fitting are currently implemented: +* IVIM fit +* DTI fit +* IVIM-DTI fit + +The fit is performed voxel-wise. Only voxels with sufficient data are considered, i.e. zeros in the data are excluded. If a voxel contains to many zeros (i.e. not enough b-values or diffusion directions), the voxel is excluded from the fit. This can occur due to registration of data at the edges of FOV for specific b-values/directions. + +### IVIM fit +The IVIM model estimates the tissue diffusion coefficient D, the perfusion fraction f and the pseudo-diffusion coefficient D*. Three fit methods are supported: free, two-step approach or segmented/simplified IVIM. + +### DTI fit +The diffusion tensor is estimated from the provided data. A contrained and an unconstrained fit are implemented. The constrained fit uses a Cholesky decomposition to ensure the diffusion tensor is positive semi-definite. + +### IVIM-DTI fit +A combined IVIM-DTI fit is performed. The following methods are currently implemented: Free fit, two-step fit, segmented/simplified IVIM-DTI fit, IVIM-corrected DTI fit. + +## Analysis tools +The function calc_dti_parameters calculated the MD, FA, eigenvalues and eigenvectors from a diffusion tensor. diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/analyze_data/calc_dti_parameters.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/analyze_data/calc_dti_parameters.m new file mode 100644 index 00000000..6f526a8f --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/analyze_data/calc_dti_parameters.m @@ -0,0 +1,89 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% This function calculates MD, FA and Eigenvalues of a diffusion tensor. +% The tensor elements should be given in the following form: +% dimensions data (= tensor): (x,y,z,tensor-elements) +% tensor-elements(1) = Dxx +% tensor-elements(2) = Dyy +% tensor-elements(3) = Dzz +% tensor-elements(4) = Dxy +% tensor-elements(5) = Dxz +% tensor-elements(6) = Dyz +% +% Returns a structure dti_params with the following fields: +% - MD: Mean diffusivity +% - FA: Fractional anisotropy +% - eigenval: The 3 eigenvalues of the diffusion tensor +% - eigenvec: The corresponding 3 eigenvectors of the diffusion tensor +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function dti_params = calc_dti_parameters(tensor, options) +arguments + tensor + options.clip {mustBeNumericOrLogical} = 0 +end + +%reshape tensor +if size(tensor, ndims(tensor)) ~= 6 + error('Tensor elements are expected in last dimension.') +end + +tensor = permute(tensor, [length(size(tensor)), 1:length(size(tensor))-1]); +tensor_sz = size(tensor); +tensor_dim = ndims(tensor); + +tensor = reshape(tensor, 6, []); + +MD = zeros(size(tensor,2),1); +FA = zeros(size(tensor,2),1); +eval = zeros([size(tensor,2),3], 'single'); +evec = zeros([size(tensor,2),3,3], 'single'); + +fprintf('Calculate pyqmri DTI parameters for each voxel... \n'); +tic; +for x = 1:size(tensor,2) + temp = squeeze(tensor(:,x)); + %put tensor elements together + dti = [temp(1) temp(4) temp(5); temp(4) temp(2) temp(6); temp(5) temp(6) temp(3)]; + %calculate eigenvectors, eigenvalues, FA and MD for each voxel + [Eigenvectors, D] = eig(dti); + EigenValues = diag(D); + [~, index] = sort(EigenValues, 'descend'); + EigenValues = EigenValues(index); Eigenvectors = Eigenvectors(:,index); + + MDv = (EigenValues(1)+EigenValues(2)+EigenValues(3))/3; + FAv = sqrt(1.5)*(sqrt((EigenValues(1)- MDv).^2 + (EigenValues(2) - MDv).^2 + (EigenValues(3) - MDv).^2) ... + ./ sqrt(EigenValues(1).^2+EigenValues(2).^2+EigenValues(3).^2)); + + % store MD, FA and Eigenvalues + MD(x) = MDv; + FA(x) = FAv; + eval(x,:) = EigenValues; + evec(x,:,:) = Eigenvectors; +end + +%Assign output arguments. +%Reshape to matrix size. Not necessary for 1-voxel analysis +if tensor_sz(2) >1 + dti_params.MD = reshape(MD, tensor_sz(2:tensor_dim)); + dti_params.FA = reshape(FA, tensor_sz(2:tensor_dim)); + dti_params.eigenval = reshape(eval, [tensor_sz(2:tensor_dim), 3]); + dti_params.eigenvec = reshape(evec, [tensor_sz(2:tensor_dim), 3, 3]); +else + dti_params.MD = squeeze(MD); + dti_params.FA = squeeze(FA); + dti_params.eigenval = squeeze(eval); + dti_params.eigenvec = squeeze(evec); +end + +if options.clip + dti_params.MD(dti_params.MD<0) = 0; + dti_params.MD(dti_params.MD>5) = 0; + dti_params.eigenval(dti_params.eigenval>10) = 0; + dti_params.eigenval(dti_params.eigenval<0) = 0; + dti_params.FA(dti_params.FA>1) = 0; +end + +fprintf('Tensor calculation finished in %.2f seconds. \n', toc); +end diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_dti.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_dti.m new file mode 100644 index 00000000..6fd3690a --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_dti.m @@ -0,0 +1,162 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to perform an DTI fit to data. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input values: +% - diffusion data: Array containing the data. Can be single voxel, a +% slice or volumentric data. The last dimension needs +% to match the number of b-value-diffusion direction +% combinations. +% +% - b-values: Vector containing the b-values +% +% - diffusion +% directions: nx3 vector containing the n diffusion directions. +% +% Options: +% - initialguess: Initial guess for D +% If not provided, default guess will be used (1). +% For the diffusion tensor, the initial guess +% is set to D on the diagonal and 0 on the +% off-diagonal elements. +% The initial guess for S0 is either 1 (if data is +% normalized) or the mean b-0 signal (if data is not +% normalized). +% +% - mask: Default 1 +% Background is masked using a threshold cut-off +% value +% +% - constrained: Default 1 +% Constrained: fit the lower triangular matrix of a +% Cholesky decomposition. The fit boundaries differ +% from the unconstrained fit. After the tensor fit, +% the diffusion tensor needs to be calculated from the +% lower triangular matrix. +% Unconstrained: fit the 6 diffusion tensor elements +% directly. +% +% - normalized: Default 1 +% Normalize data to min(bval)-signal +% +% Output: +% Structure dti_fit containing the following fields: +% - S0 signal (fitted). In case of normalization S0 will be close to 1 +% for all voxels. +% - diffusion tensor D +% vector containing the 6 elements of the diffusion tensor. +% Order of the tensor elements is: +% Dxx, Dyy, Dzz, Dxy, Dxz, Dyz. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function dti_fit = fit_dti(data, bval, diffdir, options) + +arguments + data + bval + diffdir + options.initialguess (1,1) = 1 + options.mask {mustBeNumericOrLogical} = 1 + options.constrained {mustBeNumericOrLogical} = 1 + options.normalize {mustBeNumericOrLogical} = 1 +end + +%rearrange and reshape data to a bval x n array +[data, sz] = reshape_diffdata_for_fit(data, bval); + +if size(bval,1)==1 + bval = bval.'; +end +%% optional: mask background voxels +%define datafit +if options.mask + [datafit, sels] = mask_diffdata(data, bval); +else + sels = true(1,size(data,2)); + datafit = data; +end + +%% Set initial guess for fit, optional: normalize data +if options.normalize + datafit = norm_diffdata(datafit, bval); + %set initial guess + x0(1) = 1; + x0(2:4) = options.initialguess; + x0(5:7) = 0; +else + x0 = mean(datafit(bval==min(bval),:), 'all'); + x0(2:4) = options.initialguess; + x0(5:7) = 0; +end + +%% calculate b-matrix +bval = bval_scaling(bval); +bmat = calc_bmat(bval, diffdir); + +%% set fit options +if options.constrained + %Fit Cholesky lower triangular matrix of diffusion tensor + lb = [0 -9 -9 -9 -9 -9 -9]; + ub = [inf 9 9 9 9 9 9]; +else + lb = [0 0 0 0 0 0 0]; + ub = [inf 3 3 3 3 3 3]; +end + +if options.normalize + ub(1) = 10; +end +fitoptions = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt', ... + 'FunctionTolerance', 1e-10, 'MaxIterations', 1000, 'OptimalityTolerance', 1e-10, 'StepTolerance', 1e-10, ... + 'Display', 'off'); + +%initialize diffusion tensor and S0 +tensor = zeros(6, size(datafit,2)); +S0_fit = zeros(1, size(datafit,2)); + +%% perform fit +fprintf('Perform voxel-wise DTI fit...\n'); +tic; +fail = 0; +for v = 1:size(datafit,2) + tmpdat = datafit(:,v); + tmpbmat = bmat; + [tmpdat, tmpbmat] = remove_zeros(tmpdat, tmpbmat); + [~,~,udiffdir] = unique(tmpbmat, 'rows', 'stable'); + %only perform fit if at least 7 unique diffusion directions + % (scans) are used + if max(udiffdir) > 6 + if options.constrained + t = lsqcurvefit(@dtifun_constr, x0, tmpbmat, tmpdat, lb, ub, fitoptions); + % Calculate tensor elements from Cholesky lower triangular matrix + t(2:7) = lower_triangular2tensor(t(2:7)); + else + t = lsqcurvefit(@dtifun, x0, tmpbmat, tmpdat, lb, ub, fitoptions); + end + + tensor(:,v) = t(2:7); + S0_fit(v) = t(1); + else + fail = fail + 1; + t(1:7) = 0; + end +end +fprintf('%d pixels were not fitted due to too much missing data. \n', fail); + +fprintf('Fit performed in %.2f seconds. \n', toc); + +%% rearrange output +sz_fit = sz(2:end); +if length(sz_fit)==1 + sz_fit = [1 sz_fit]; +end + +dti_fit.S0 = zeros(sz_fit); +dti_fit.S0(sels) = S0_fit; + +dti_fit.D = squeeze(zeros([6, sz_fit])); +dti_fit.D(:,sels) = tensor; +dti_fit.D = permute(dti_fit.D, [2:length(size(dti_fit.D)) 1]); +end diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_ivim.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_ivim.m new file mode 100644 index 00000000..e766dcb7 --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_ivim.m @@ -0,0 +1,230 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to perform an IVIM fit to data. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input values: +% - Diffusion data: Array containing the data. Can be single voxel, a +% slice or volumetric data. The last dimension needs +% to match the number of b-values. +% +% - b-values: Vector containing the b-values +% +% Options: +% - initialguess: D, f, Ds +% If not provided, default guess will be used [1, +% 0.2, 50]. +% The initial guess for S0 is either 1 (if data is +% normalized) or the mean b-0 signal (if data is not +% normalized). +% +% - mask: Default 1 +% Background is masked by using a threshold cut-off +% value. +% +% - normalize: Default 1 +% Normalize data to the min(bval)-signal +% +% - fit_method: Default 'two_step' +% 'free', 'two_step', 'segmented' +% 'free': The full data with all b-values is fitted +% to the IVIM equation. +% 'two_step': A two-step fit will be performed, +% estimating D and f from high b-values (b >= bcut) +% only and performing a bi-exponential fit with fixed +% D to obtain f and D*. +% 'segmented': A segmented IVIM fit is performed, +% using b = 0 and b >= bcut to estimate D and f. +% D* is not fitted with this method. +% +% - bcut: Default 200 +% Cut-off b-value for the two_step and segmented fit. +% +% - seg_data: Default 0 +% Select only data with bval >= bcut and b = 0 for +% fitting. Only used for segmented IVIM fit +% (fit_method = 'segmented') +% +% Output: +% Structure ivim_fit containing the following fields: +% - S0 signal (fitted). In case of normalization S0 will be close to 1 +% for all voxels. +% - diffusion coefficient D +% - perfusion fraction f +% - pseudo-diffusion coefficient Ds +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function ivim_fit = fit_ivim(data, bval, options) + +arguments + data + bval + options.initialguess (1,3) = [1, 0.2, 50]; + options.mask {mustBeNumericOrLogical} = 1 + options.normalize {mustBeNumericOrLogical} = 1 + options.fit_method {mustBeMember(options.fit_method, {'free', 'two_step', 'segmented'})} = 'two_step' + options.bcut {mustBeNumeric} = 200 + options.seg_data {mustBeNumericOrLogical} = 0 + options.Ds_fix {mustBeNumeric} = 0 + +end + +%% rearrange and reshape data to a bval x n array +[data, sz] = reshape_diffdata_for_fit(data, bval); + +%% optional: mask background voxels +%define datafit +if options.mask + [datafit, sels] = mask_diffdata(data, bval); +else + sels = true(1,size(data,2)); + datafit = data; +end + +%% Set initial guess for fit, optional: normalize data +if options.normalize + datafit = norm_diffdata(datafit, bval); + x0 = [1, options.initialguess]; +else + x0 = [mean(datafit(bval==min(bval),:), 'all'), options.initialguess]; +end + +%% scale b-value +%scale b-value +if size(bval,1)==1 + bval = bval.'; +end +bval = bval_scaling(bval); +options.bcut = bval_scaling(options.bcut); +%% set fit options +%boundaries +lb = [0, 0, 0, 5+1e-5]; +ub = [inf, 5, 1, 300]; + +if options.normalize + ub(1) = 10; +end + +fitoptions = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt', ... + 'FunctionTolerance', 1e-10, 'MaxIterations', 1000, 'OptimalityTolerance', 1e-10, 'StepTolerance', 1e-10, ... + 'Display', 'off'); + +%initialize fit parameters +S0_fit = zeros(1, size(datafit,2)); +D_fit = zeros(1, size(datafit,2)); +f_fit = zeros(1, size(datafit,2)); +Ds_fit = zeros(1, size(datafit,2)); +%% fit data voxelwise +fprintf('Perform voxel-wise IVIM fit...\n') +tic; +fail = 0; +for v = 1:size(datafit, 2) + tmpdat = datafit(:,v); + tmpb = bval; + [tmpdat, tmpb] = remove_zeros(tmpdat, tmpb); + + switch options.fit_method + %%%%%%%%%%%%%%%%%%% + % free fit + %%%%%%%%%%%%%%%%%%% + case 'free' + %only perform fit if at least 4 b-values are used + if numel(unique(tmpb)) > 3 + x = lsqcurvefit(@ivimfun, x0, tmpb, tmpdat, lb, ub, fitoptions); + S0_fit(v) = x(1); + D_fit(v) = x(2); + f_fit(v) = x(3); + Ds_fit(v) = x(4); + else + fail = fail+1; + end + %%%%%%%%%%%%%%%%%%% + % two_step fit + %%%%%%%%%%%%%%%%%%% + case 'two_step' + if v == 1 + %adjust boundaries and initial guess + lb_twostep = lb([1 3 4]); ub_twostep = ub([1 3 4]); x0_twostep = x0([1 3 4]); + end + %only perform fit if at least 4 b-values are used + if numel(unique(tmpb)) > 3 + %First, perform mono-exponential fit to b >= bcut to get an estimate for D + [param(1:2)] = lsqcurvefit(@(x, xdat) x(1)*exp(-xdat*x(2)), x0([1,2]), tmpb(tmpb>=options.bcut), tmpdat(tmpb>=options.bcut), ... + lb([1,2]), ub([1,2]), fitoptions); + + x(2) = param(2); %D + + %Estimate f and use it as initial guess for IVIM fit + S0_tmp = mean(tmpdat(tmpb==min(tmpb))); + f_guess = 1 - (param(1) / S0_tmp); + if f_guess < 0 || isnan(f_guess) + f_guess = 0; + end + x0_twostep(2) = f_guess; + + %perform fit. D is fixed for the bi-exponential IVIM fit. + input.bval = tmpb; + input.D_fix = x(2); + if options.Ds_fix==0 %fit D* + x([1,3:4]) = lsqcurvefit(@ivimfun, x0_twostep, input, tmpdat, lb_twostep, ub_twostep, fitoptions); + else %fixed D* given as input + x(4) = options.Ds_fix; + input.Ds_fix = x(4); + x([1,3]) = lsqcurvefit(@ivimfun, x0_twostep([1 2]), input, tmpdat, lb_twostep([1 2]), ub_twostep([1 2]), fitoptions); + end + S0_fit(v) = x(1); + D_fit(v) = x(2); + f_fit(v) = x(3); + Ds_fit(v) = x(4); + else + fail = fail+1; + end + %%%%%%%%%%%%%%%%%%% + % segmented fit + %%%%%%%%%%%%%%%%%%% + case 'segmented' + %optional: select desired b-value range (if full range was + %acquired) + if options.seg_data && v == 1 + tmpdat = [tmpdat(tmpb==0); tmpdat(tmpb>=options.bcut)]; + tmpb = [tmpb(tmpb==0); tmpb(tmpb>=options.bcut)]; + end + %only perform fit if at least 2 b-values are used + if numel(unique(tmpb)) > 1 + %perform fit + x = lsqcurvefit(@ivimfun_seg, x0(1:3), tmpb, tmpdat, lb(1:3), ub(1:3), fitoptions); + S0_fit(v) = x(1); + D_fit(v) = x(2); + f_fit(v) = x(3); + + else + fail = fail+1; + x(1:4) = 0; + end + end +end +fprintf('%d pixels were not fitted due to too much missing data. \n', fail); + +fprintf('Fit performed in %.2f seconds. \n', toc); + +%% rearrange output +sz_fit = sz(2:end); +if length(sz_fit)==1 + sz_fit = [1 sz_fit]; +end +ivim_fit.S0 = zeros(sz_fit); +ivim_fit.S0(sels) = S0_fit; + +ivim_fit.D = zeros(sz_fit); +ivim_fit.D(sels) = D_fit; + +ivim_fit.f = zeros(sz_fit); +ivim_fit.f(sels) = f_fit; + +if ~isequal(options.fit_method, 'segmented') + ivim_fit.Ds = zeros(sz_fit); + ivim_fit.Ds(sels) = Ds_fit; +end + +end diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_ivimdti.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_ivimdti.m new file mode 100644 index 00000000..04834cea --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_ivimdti.m @@ -0,0 +1,312 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to perform a combined IVIM-DTI fit to data. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input values: +% - Diffusion data: Array containing the data. Can be single voxel, a +% slice or volumentric data. The last dimension needs +% to match the number of b-values/diffusion direction +% combinations. +% +% - b-values: Vector containing the b-values +% +% - diffusion +% directions: nx3 vector containing the diffusion directions +% +% Options: +% - initialguess: Initial guess for D, f and Ds +% If not provided, default guess will be used [1, +% 0.2, 50] +% For the diffusion tensor, the initial guess +% is set to D on the diagonal and 0 on the +% off-diagonal elements. +% The initial guess for S0 is either 1 (if data is +% normalized) or the mean b-0 signal (if data is not +% normalized). +% +% - mask: Default 1 +% Background is masked by using a threshold cut-off +% value based on the minimum b-value data. +% +% - normalize: Default 1 +% Normalize data to the min (bval)-signal +% +% - fit_method Default 'free' +% 'free', 'two_step', 'segmented', 'IVIM_correct' +% 'free': The data is fitted to the full IVIM-DTI +% equation. +% 'two_step': First, the diffusion tensor and f are +% estimated using the high b-values only (b >= bcut). +% Then, Dtensor is kept constant and a bi-exponential +% fit is performed to estimate f and Ds. The +% previously calculated value for f is used as +% initial guess. +% 'segmented': A segmented IVIM-DTI fit is performed, +% using b = 0 and b >= bcut to estimate the diffusion +% tensor and f. D* ist not fitted with this method. +% 'IVIM_correct': First an IVIM fit is performed to +% obtain f and Ds. The data is then IVIM-corrected by +% substracting the IVIM component. A pure DTI fit is +% performed on the remaining data. +% +% - dti_constrained Default 1 +% Constrained fit for the diffusion tensor elements. +% The lower triangular matrix of a Cholesky +% decomposition of the diffusion tensor is fitted. +% After the fit, the diffusion tensor can be +% calculated. The fit boundaries differ from the +% unconstrained fit. +% +% - bcut: Default 200 +% Cut-off b-value for the two_step fit. The cut-off +% b-value is included for analysis. +% +% Output: +% Structure ivimdti_fit containing the following fields: +% - S0 signal (fitted). In case of normalization S0 will be close to 1 +% for all voxels. +% - tensor: Diffusion tensor containing the 6 elements in the order +% Dxx, Dyy, Dzz, Dxy, Dxz, Dyz +% - perfusion fraction f +% - pseudo-diffusion coefficient Ds +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function ivimdti_fit = fit_ivimdti(data, bval, diffdir, options) + +arguments + data + bval + diffdir + options.initialguess (3,1) = [1, 0.2, 50] + options.mask {mustBeNumericOrLogical} = 1 + options.normalize {mustBeNumericOrLogical} = 1 + options.fit_method {mustBeMember(options.fit_method, {'free', 'two_step', 'segmented', 'IVIM_correct'})} = 'free' + options.dti_constrained = 1; + options.bcut {mustBeNumeric} = 200 +end + +%rearrange and reshape data to a bval x n array +[data, sz] = reshape_diffdata_for_fit(data, bval); + +%% optional: mask background voxels +%define datafit +if options.mask + [datafit, sels] = mask_diffdata(data, bval); +else + sels = true(1, size(data,2)); + datafit = data; +end + +%% Set initial guess for fit, optional: normalize data +% x0 = [S0, Dxx, Dyy, Dzz, Dxy, Dxz, Dyz, f, D*] +if options.normalize + datafit = norm_diffdata(datafit, bval); + %set initial guess + x0(1) = 1; + x0(2:4) = options.initialguess(1); + x0(5:7) = 0; + x0(8:9) = options.initialguess(2:3); +else + x0(1) = mean(datafit(bval==min(bval),:), 'all'); + x0(2:4) = options.initialguess(1); + x0(5:7) = 0; + x0(8:9) = options.initialguess(2:3); +end + +%% scale b-value +if size(bval,1)==1 + bval = bval.'; +end +bval = bval_scaling(bval); +options.bcut = bval_scaling(options.bcut); +%% set fit options +diffparams.bval = bval; +diffparams.diffdir = diffdir; + +if options.dti_constrained + %Fit Cholesky lower triangular matrix coefficients + lb = [0 -9 -9 -9 -9 -9 -9 0 5+1e-5]; + ub = [inf 9 9 9 9 9 9 1 300]; +else + lb = [0 0 0 0 0 0 0 0 5]; + ub = [inf 3 3 3 3 3 3 1 300]; +end +if options.normalize + ub(1) = 10; +end + +fitoptions = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt', ... + 'FunctionTolerance', 1e-10, 'maxIterations', 1000, 'OptimalityTolerance', 1e-10, ... + 'Display', 'off'); + +%initialize fit parameters +S0_fit = zeros(1, size(datafit, 2)); +tensor_fit = zeros(6, size(datafit, 2)); +f_fit = zeros(1, size(datafit, 2)); +Ds_fit = zeros(1, size(datafit, 2)); + +%% perform fit +fprintf('Perform voxel-wise IVIM-DTI fit...\n'); +tic; +fail = 0; +%Select fit method +switch options.fit_method + + %%%%%%%%%%%%%%%%%%% + % free fit + %%%%%%%%%%%%%%%%%% + case 'free' + %loop over voxels + for v = 1:size(datafit,2) + tmpdat = datafit(:,v); + tmpdiff = diffparams; + [tmpdat, tmpdiff] = remove_zeros(tmpdat, tmpdiff); + [~,~,udiffdir] = unique(tmpdiff.diffdir(tmpdiff.bval>0,:), 'rows', 'stable'); + %only perform fit if at least 4 b-values and 6 unique diffusion + %directions (without b0) are available + if numel(unique(tmpdiff.bval)) > 3 && max(udiffdir) >= 6 + if options.dti_constrained + x = lsqcurvefit(@ivimdtifun_constr, x0, tmpdiff, tmpdat, lb, ub, fitoptions); + x(2:7) = lower_triangular2tensor(x(2:7)); + else + x = lsqcurvefit(@ivimdtifun, x0, tmpdiff, tmpdat, lb, ub, fitoptions); + end + S0_fit(:,v) = x(1); + tensor_fit(:,v) = x(2:7); + f_fit(:,v) = x(8); + Ds_fit(:,v) = x(9); + else + fail = fail+1; + x(1:9) = 0; + end + end + + %%%%%%%%%%%%%%%%%%% + % two-step fit + %%%%%%%%%%%%%%%%%% + case 'two_step' + %loop over voxels + for v = 1:size(datafit,2) + tmpdat = datafit(:,v); + tmpdiff = diffparams; + [tmpdat, tmpdiff] = remove_zeros(tmpdat, tmpdiff); + [~,~,udiffdir] = unique(tmpdiff.diffdir(tmpdiff.bval>0,:), 'rows', 'stable'); + if v == 1 + lb_seg = lb([1 8 9]); ub_seg = ub([1 8 9]); x0_seg = x0([1 8 9]); + end + %only perform fit if at least 2 b-values and 6 unique diffusion + %directions (without b0) are available + if numel(unique(tmpdiff.bval)) > 1 && max(udiffdir) >= 6 + %First, fit diffusion tensor to b >= bcut + if options.dti_constrained + [param(1:7)] = lsqcurvefit(@dtifun_constr, x0(1:7), calc_bmat(tmpdiff.bval(tmpdiff.bval>=options.bcut), tmpdiff.diffdir(tmpdiff.bval>=options.bcut,:)), ... + tmpdat(tmpdiff.bval>=options.bcut), lb(1:7), ub(1:7), fitoptions); + x(2:7) = lower_triangular2tensor(param(2:7)); + else + [param(1:7)] = lsqcurvefit(@dtifun, x0(1:7), calc_bmat(tmpdiff.bval(tmpdiff.bval>=options.bcut), tmpdiff.diffdir(tmpdiff.bval>=options.bcut,:)), ... + tmpdat(tmpdiff.bval>=options.bcut), lb(1:7), ub(1:7), fitoptions); + x(2:7) = param(2:7); + end + %Estimate f and use it as initial guess for IVIM fit + S0_tmp = mean(tmpdat(tmpdiff.bval==min(tmpdiff.bval))); + f_guess = 1 - (param(1) / S0_tmp); + if f_guess < 0 || isnan(f_guess) + f_guess = 0; + end + x0_seg(2) = f_guess; + + %perform fit with fixed diffusion tensor + tmpdiff.tensor = x(2:7); + x([1,8:9]) = lsqcurvefit(@ivimdtifun, x0_seg, tmpdiff, tmpdat, lb_seg, ub_seg, fitoptions); + + + S0_fit(:,v) = x(1); + tensor_fit(:,v) = x(2:7); + f_fit(:,v) = x(8); + Ds_fit(:,v) = x(9); + else + fail = fail+1; + x(1:9) = 0; + end + end + + %%%%%%%%%%%%%%%%%%% + % segmented fit + %%%%%%%%%%%%%%%%%%% + case 'segmented' + %loop over voxels + for v = 1:size(datafit,2) + tmpdat = datafit(:,v); + tmpdiff = diffparams; + [tmpdat, tmpdiff] = remove_zeros(tmpdat, tmpdiff); + [~,~,udiffdir] = unique(tmpdiff.diffdir(tmpdiff.bval>0,:),'rows', 'stable'); + %only perform fit if at least 2 b-values and 6 unique diffusion + %directions (without b0) are available + if numel(unique(tmpdiff.bval)) > 1 && max(udiffdir) >= 6 + if options.dti_constrained + x = lsqcurvefit(@ivimdtifun_seg_constr, x0(1:8), tmpdiff, tmpdat, lb(1:8), ub(1:8), fitoptions); + x(2:7) = lower_triangular2tensor(x(2:7)); + else + x = lsqcurvefit(@ivimdtifun_seg, x0(1:8), tmpdiff, tmpdat, lb(1:8), ub(1:8), fitoptions); + end + S0_fit(:,v) = x(1); + tensor_fit(:,v) = x(2:7); + f_fit(:,v) = x(8); + else + fail = fail+1; + x(1:8) = 0; + end + end + + %%%%%%%%%%%%%%%%%%% + % IVIM-correct DTI data + %%%%%%%%%%%%%%%%%%% + % Based on QMRITools in Mathematica by Martijn Froeling + case 'IVIM_correct' + %first, fit IVIM to mean signal to obtain D* to fix + [mdat, ~] = mask_diffdata(data, bval); + ivim_mean = fit_ivim(mean(mdat,2), bval, 'mask', 0, 'normalize', 1); + fprintf('Fixed D* value obtained from mean fit: D* = %d \n', ivim_mean.Ds); + %then, fit IVIM (voxel-wise) to the signal with a fixed D* value (ivimfun) + ivim_pars = fit_ivim(datafit.', bval, 'Ds_fix', ivim_mean.Ds, ... + 'mask', 0, 'normalize', 0); + %IVIM-correct DTI signal by substracting IVIM component + datafit_cor = datafit - ivim_pars.S0.*ivim_pars.f.*exp(-bval.*ivim_pars.Ds); + %fit DTI to remaining DTI signal (dtifun) + dti_pars = fit_dti(datafit_cor.', bval, diffdir, ... + 'constrained', options.dti_constrained, 'mask', 0, 'normalize', 0); + + S0_fit = ivim_pars.S0; + f_fit = ivim_pars.f; + Ds_fit = ivim_pars.Ds; + tensor_fit = dti_pars.D.'; +end +fprintf('%d pixels were not fitted due to too much missing data. \n', fail); + +fprintf('Fit performed in %.2f seconds. \n', toc); + +%% rearrange output +sz_fit = sz(2:end); +if length(sz_fit)==1 + sz_fit = [1 sz_fit]; +end + +ivimdti_fit.S0 = zeros(sz_fit); +ivimdti_fit.S0(sels) = S0_fit; + +ivimdti_fit.tensor = squeeze(zeros([6, sz_fit])); +ivimdti_fit.tensor(:,sels) = tensor_fit; +ivimdti_fit.tensor = permute(ivimdti_fit.tensor, ... + [2:length(size(ivimdti_fit.tensor)) 1]); + +ivimdti_fit.f = zeros(sz_fit); +ivimdti_fit.f(sels) = f_fit; + +if ~isequal(options.fit_method, 'segmented') + ivimdti_fit.Ds = zeros(sz_fit); + ivimdti_fit.Ds(sels) = Ds_fit; +end + +end diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/bval_scaling.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/bval_scaling.m new file mode 100644 index 00000000..b4613d4b --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/bval_scaling.m @@ -0,0 +1,21 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to scale the b-value to the unit 10^-3 s/mm² to match with the +% IVIM/DTI parameter units. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input: b-value (array/list) +% Output: scaled b-value +% +% If the b-value is scaled correctly already, the original b-value is +% returned. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function bval_scaled = bval_scaling(bval) + +if max(bval) > 100 + bval_scaled = bval./1000; +else + bval_scaled = bval; +end \ No newline at end of file diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/calc_bmat.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/calc_bmat.m new file mode 100644 index 00000000..1505294d --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/calc_bmat.m @@ -0,0 +1,28 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to calculate the b-matrix from the b-values and the diffusion +% gradient directions. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input: b-values (1 x n vector), diffusion directions (n x 3 vector) +% Output: +% bmatrix = bval.* [gx.^2, gy.^2, gz.^2, 2*gx*gy, 2*gx*gz, 2*gy*gz]; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function bmat = calc_bmat(bval, diffdir) + if size(bval,1)==1 + bval = bval.'; + end + if size(diffdir,2) ~= 3 + diffdir = diffdir.'; + end + %scale b-value + bval = bval_scaling(bval); + + dir = [diffdir(:,1).^2, diffdir(:,2).^2, diffdir(:,3).^2, ... + 2*diffdir(:,1).*diffdir(:,2), ... + 2*diffdir(:,1).*diffdir(:,3), ... + 2*diffdir(:,2).*diffdir(:,3)]; + bmat = bval.*dir; +end \ No newline at end of file diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/delta_fun.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/delta_fun.m new file mode 100644 index 00000000..76947870 --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/delta_fun.m @@ -0,0 +1,17 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Dirac Delta-function +% +% Returns 1 if the input value is 0 and returns 0 otherwise. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function val = delta_fun(x) + +val = zeros(size(x)); + +val(x==0) = 1; +val(x~=0) = 0; + +end \ No newline at end of file diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/lower_triangular2tensor.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/lower_triangular2tensor.m new file mode 100644 index 00000000..ed561699 --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/lower_triangular2tensor.m @@ -0,0 +1,41 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% This function calculates the (diffusion) tensor from the given Cholesky +% lower triangular matrix (obtained by Cholesky decomposition). +% Input is a n x 6 vector containing the matrix elements. +% +% Output are the tensor elements in the order Dxx, Dyy, Dzz, Dxy, Dxz, Dyz +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function cd = lower_triangular2tensor(r) + % Cholesky decomposition + % cd = (r1 0 0 ) (r1 r4 r5 ) + % (r4 r2 0 ) * (0 r2 r6 ) + % (r5 r6 r3 ) (0 0 r3 ) + % + % + % = ( r1^2 r1*r4 r1*r5 ) + % ( r4^2+r2^2 r4*r5 + r2*r6 ) + % ( r5^2+r6^2+r3^2 ) + +%check dimension +if ndims(r) ~= 2 %#ok + error('wrong matrix size') +end +if isempty(find(size(r)==6, 1)) + error('no valid tensor provided') +end + +if min(size(r))==1 + %single voxel + cd = [r(1).^2, r(4).^2+r(2).^2, r(5).^2+r(6).^2+r(3).^2, ... + r(1).*r(4), r(1).*r(5), r(4).*r(5)+r(2).*r(6)]; +else + %multiple voxels + if size(r,1) ~= 6 + r = r.'; + end + cd = [r(1,:).^2; r(4,:).^2+r(2,:).^2; r(5,:).^2+r(6,:).^2+r(3,:).^2; ... + r(1,:).*r(4,:); r(1,:).*r(5,:); r(4,:).*r(5,:)+r(2,:).*r(6,:)]; +end \ No newline at end of file diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/mask_diffdata.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/mask_diffdata.m new file mode 100644 index 00000000..24d7af3c --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/mask_diffdata.m @@ -0,0 +1,34 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to mask diffusion data. The masking is performed using a +% threshold on the mean (bval = min(bval)-data). The threshold can be set +% manually using options.thr +% +% The masked data is reshaped to a 2-D array and has the size +% (numel(bval), sels). +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function [masked, sels] = mask_diffdata(data, bval, options) + +arguments + data + bval + options.thr {mustBeNumeric} = 0.01 +end + +%check data size +try + bdim = find(size(data)==length(bval)); +catch + error('Datasize and number of b-values not matching!') +end + +data = permute(data, [bdim, 1:bdim-1 bdim+1:length(size(data))]); + +%print b-value used for masking +fprintf('Smallest b-value found and used for masking: b = %.2f \n', min(bval)); + +sels = mean(data(bval==min(bval),:),1) > options.thr*mean(data(bval==min(bval),:)); +masked = data(:, sels); + diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/norm_diffdata.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/norm_diffdata.m new file mode 100644 index 00000000..59e33a02 --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/norm_diffdata.m @@ -0,0 +1,43 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to normalize diffusion data. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input: diffusion data, b-values +% Output: normalized diffusion data +% +% The diffusion data is normalized to the mean b = 0 data. If no b = 0 data +% is provided the data is normalized to the smallest b-value data. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function normalized = norm_diffdata(data, bval) + +% check data +if ndims(data) > 2 %#ok + [data, sz, bdim_org, perm] = reshape_diffdata_for_fit(data, bval); +end + +% check if b0 data is provided +if min(bval) >= 1 + fprintf('No b = 0 data provided. Normalize to minimum b-value data: %.2f s/mm². \n',... + min(bval)) + mS0 = mean(data((bval == min(bval)),:), 1); +else + fprintf('Normalize data to b = %.2f data. \n', min(bval)) + mS0 = mean(data(bval<1, :), 1); +end + +normalized = data ./ mS0; +normalized(isnan(normalized)) = 0; + +if exist('sz', 'var') + %reshape data to origianl size + normalized = reshape(normalized, sz); + bdim = find(size(normalized)==length(bval)); + if bdim ~= bdim_org + %reverse permutation + rp(perm) = 1:length(perm); + normalized = permute(normalized, rp); + end +end \ No newline at end of file diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/remove_zeros.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/remove_zeros.m new file mode 100644 index 00000000..1e7a0bf4 --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/remove_zeros.m @@ -0,0 +1,25 @@ +function [data, diffparams, varargout] = remove_zeros(data, diffparams, options) +arguments + data + diffparams + options.bval = [] +end + +if isstruct(diffparams) + %IVIM-DTI + diffparams.bval(data==0) = []; + diffparams.diffdir(data==0,:) = []; +elseif isvector(diffparams) + %IVIM, diffparams = bval + diffparams(data==0) = []; +else + %DTI, diffparams = b-matrix + diffparams(data==0,:) = []; +end +data(data==0) = []; + +if ~isempty(options.bval) + options.bval(data==0) = []; + varargout{1} = options.bval; +end + diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/reshape_diffdata_for_fit.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/reshape_diffdata_for_fit.m new file mode 100644 index 00000000..2cbba605 --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/reshape_diffdata_for_fit.m @@ -0,0 +1,49 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to reshape diffusion data to a bval x n array for data fitting. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input: +% - data +% - b-values +% +% Output: +% - diffdata +% Reshaped diffusion data +% The first dimension is the b-value and the second dimension +% contains the data. +% - data_sz +% The original array size (data_sz) is returned to allow reshaping +% of the data. +% - Optional: +% The b-value dimension can be returned (varargout{1}). +% Permutation vector can be returned (varargout{2}). +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function [diffdata, data_sz, varargout] = reshape_diffdata_for_fit(data, bval) + +%only for multi-voxel fit +if ~(size(data,2) == 1) + if isempty(find(size(data, ndims(data))==length(bval), 1)) + error('Datasize and number of b-values are not matching!') + end + bdim = ndims(data); + diffdata = permute(data, [bdim, 1:bdim-1 bdim+1:length(size(data))]); + data_sz = size(diffdata); + + diffdata = reshape(diffdata, length(bval), []); + diffdata = double(diffdata); + + if nargout > 2 + varargout{1} = bdim; + end + + if nargout > 3 + %permutation vector + varargout{2} = [bdim, 1:bdim-1 bdim+1:length(size(data))]; + end +else %single voxel case (e.g. simulations) + diffdata = data; + data_sz = size(data); +end \ No newline at end of file From 364900da0ee530331819cdf0c48d147337ac4441 Mon Sep 17 00:00:00 2001 From: s-rauh Date: Wed, 5 Aug 2026 08:38:17 +0200 Subject: [PATCH 2/5] Add IVIM and DTI models to SR_LUMC/diffusion_mri_tools fiting code --- .../models_ivim_dti/dtifun.m | 39 ++++++++ .../models_ivim_dti/dtifun_constr.m | 42 +++++++++ .../models_ivim_dti/ivimdtifun.m | 94 +++++++++++++++++++ .../models_ivim_dti/ivimdtifun_constr.m | 49 ++++++++++ .../models_ivim_dti/ivimdtifun_seg.m | 66 +++++++++++++ .../models_ivim_dti/ivimdtifun_seg_constr.m | 59 ++++++++++++ .../models_ivim_dti/ivimfun.m | 72 ++++++++++++++ .../models_ivim_dti/ivimfun_seg.m | 37 ++++++++ 8 files changed, 458 insertions(+) create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/dtifun.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/dtifun_constr.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun_constr.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun_seg.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun_seg_constr.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimfun.m create mode 100644 src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimfun_seg.m diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/dtifun.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/dtifun.m new file mode 100644 index 00000000..a04b7abf --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/dtifun.m @@ -0,0 +1,39 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to calculate the signal from a DTI measurement given the +% diffusion tensor. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input values: +% - diffusion tensor: given as n x 7 vector +% The first entry of the diffusion vector is S0. +% Order of the tensor elements tensor(2:7): +% Dxx, Dyy, Dzz, Dxy, Dxz, Dyz +% +% - b-matrix: given as (n x 6 vector) +% +% The signal is calculated according to S = S0*exp(-bmat*diff_tensor) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function signal = dtifun(tensor, bmat) + +if ndims(tensor) ~= 2 %#ok + error('wrong matrix size') +end +if isempty(find(size(tensor)==7, 1)) + error('no valid diffusion tensor provided') +end + +%check sizes and calculate signal +if min(size(tensor))==1 + %single voxel + signal = tensor(1)*exp(-bmat*tensor(2:7).'); +else + if size(tensor,1) ~= 7 + tensor = tensor.'; + end + %multiple voxels + signal = tensor(1,:).*exp(-bmat*tensor(2:7,:)); +end + diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/dtifun_constr.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/dtifun_constr.m new file mode 100644 index 00000000..d795b756 --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/dtifun_constr.m @@ -0,0 +1,42 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to calculate the diffusion signal given the lower triangular +% matrix (from the Cholesky decomposition) of the diffusion tensor. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input values: +% - lower triangular matrix of the diffusion tensor +% - b-matrix +% The diffusion tensor and b-matrix are given as n x 6 element vectors. The +% first entry of the diffusion tensor is the S0 signal. +% +% First, the diffusion tensor is calculated from the lower triangular +% matrix using the function lower_triangular2tensor. +% Then, the signal is calculated according to S = S0*exp(-bmat*diff_tensor) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function signal = dtifun_constr(tensor, bmat) + + +if ndims(tensor) ~= 2 %#ok + error('wrong matrix size') +end +if isempty(find(size(tensor)==7, 1)) + error('no valid diffusion tensor provided') +end + + +%Calculate tensor from Cholesky decomposition +if min(size(tensor))==1 + %single voxel + tensor(2:7) = lower_triangular2tensor(tensor(2:7)); +else + if size(tensor,1) ~= 7 + tensor = tensor.'; + end + tensor(2:7,:) = lower_triangular2tensor(tensor(2:7,:)); +end + +%Calculate signal +signal = dtifun(tensor, bmat); \ No newline at end of file diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun.m new file mode 100644 index 00000000..bf219dab --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun.m @@ -0,0 +1,94 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to calculate the IVIM-DTI signal from a measurement with many +% b-bvalues and diffusion directions. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input values: +% - params: model parameters +% params(1) = S0 signal +% params(2:7) = diffusion tensor, given as +% n x 6 vector. The order of the tensor elemts +% is: Dxx, Dyy, Dzz, Dxy, Dxz, Dyz +% params(8) = perfusion fraction f +% params(9) = pseudo-diffusion coefficient D* +% +% - diffparams: structure containing the fields 'bval' with the +% diffusion b-values and 'diffdir', containing the +% gradient directinos. +% In case of two-step fitting, the diffusion +% tensor is fixed and diffparams contains an +% additional field 'tensor'. +% +% The b-matrix is calculated from the b-values and gradient directions. +% The signal is calculated according to +% S = S0*(f*exp(-b*Ds) + (1-f)*exp(-bmat*diff_tensor)) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function signal = ivimdtifun(params, diffparams) + +%check input parameters +if ndims(params) ~= 2 %#ok + error('Wrong matrix size of input parameters') +end + +if isempty(find(size(params)==9,1)) && isempty(find(size(params)==3,1)) + error('Wrong number of input parameters') +end + +if ~isstruct(diffparams) + error('Diffparams must contain b-values and diffusion gradient directions') +end +if (~isfield(diffparams, 'bval') || ~isfield(diffparams, 'diffdir')) + error('Diffparams must contain b-values and diffusion gradient directions') +end +if isempty(find(size(params)==9,1)) && ~isfield(diffparams, 'tensor') + error('For two-step fitting diffparams must contain the diffusion tensor.') +end + +%calculate b-matrix +bval = bval_scaling(diffparams.bval); +diffdir = diffparams.diffdir; + +bmat = calc_bmat(bval, diffdir); + +%calculate signal +if min(size(params))==1 + %single voxel + S0 = params(1); + if length(params) ~= 9 + %tensor fixed + tensor = diffparams.tensor; + else + %tensor fitted + tensor = params(2:7); + end + + f = params(end-1); + Ds = params(end); + +else + %multiple voxels + if isempty(find(size(params)==9,1)) + %tensor fixed + if size(params,1) ~=3 + params = params.'; + end + tensor = diffparams.tensor; + else + %tensor fitted + if size(params,1) ~=9 + params = params.'; + end + tensor = params(2:7,:); + end + + S0 = params(1,:); + + f = params(end-1,:); + Ds = params(end,:); +end + +signal = S0*(f*exp(-bval*Ds) + (1-f)*exp(-bmat*tensor.')); + diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun_constr.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun_constr.m new file mode 100644 index 00000000..d54b5a77 --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun_constr.m @@ -0,0 +1,49 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to calculate the IVIM-DTI signal given the lower triangular +% matrix (from the Cholesky decomposition) of the diffusion tensor, +% perfusion fraction f and pseudo-diffusion coefficient D_star. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input values: +% - params: model parameters +% params(1) = S0 signal +% params(2:7) = lower triangular matrix of the +% diffusion tensor +% params(8) = perfusion fraction f +% params(9) = pseudo-diffusion coefficient D* +% +% - diffparams: structure containing the fields 'bval' with the +% diffusion b-values and 'diffdir', containing the +% gradient directinos. +% +% First, the diffusion tensor is calculated from the lower triangular +% matrix using the function lower_triangular2tensor. +% Then, the signal is calculated using the function ivimdtifun: +% S = S0*(f*exp(-b*Ds) + (1-f)*exp(-bmat*diff_tensor)) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function signal = ivimdtifun_constr(params, diffparams) + +if ndims(params) ~= 2 %#ok + error('Wrong matrix size of input parameters') +end +if isempty(find(size(params)==9,1)) + error('Wrong number of input parameters') +end + +%Calculate tensor from Cholesky decomposition +if min(size(params))==1 + %single voxel + params(2:7) = lower_triangular2tensor(params(2:7)); +else + if size(params,1) ~= 9 + params = params.'; + end + params(2:7,:) = lower_triangular2tensor(params(2:7,:)); +end + +%Calculate signal +signal = ivimdtifun(params, diffparams); + diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun_seg.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun_seg.m new file mode 100644 index 00000000..e5db8ef8 --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun_seg.m @@ -0,0 +1,66 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to calculate the IVIM-DTI signal from a measurement with many +% b-bvalues and diffusion directions using the segmented IVIM-DTI equation. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input values: +% - params: model parameters +% params(1) = S0 signal +% params(2:7) = diffusion tensor, given as +% n x 7 vector. The order of the tensor elemts +% is: Dxx, Dyy, Dzz, Dxy, Dxz, Dyz +% params(8) = perfusion fraction f +% +% - diffparams: structure containing the fields 'bval' with the +% diffusion b-values and 'diffdir', containing the +% gradient directinos. +% +% The b-matrix is calculated from the b-values and gradient directions. +% The signal is calculated according to +% S = S0*(f*delta_fun(bval) + (1-f)*exp(-bmat*diff_tensor)) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function signal = ivimdtifun_seg(params, diffparams) + +%check input parameters +if ndims(params) ~= 2 %#ok + error('Wrong matrix size of input parameters') +end + +if length(params)>8 + warning('Too many input parameters provided, additional parameters are ignored') +end + +if ~isstruct(diffparams) + error('Diffparams must be a structure containing b-values and diffusion gradient directions') +end +if (~isfield(diffparams, 'bval') || ~isfield(diffparams, 'diffdir')) + error('Diffparams must contain b-values and diffusion gradient directions') +end + +%calculate b-matrix +bval = bval_scaling(diffparams.bval); +diffdir = diffparams.diffdir; + +bmat = calc_bmat(bval, diffdir); + +%calculate signal +if min(size(params))==1 + %single voxel + S0 = params(1); + tensor = params(2:7); + f = params(8); +else + %multiple voxels + if size(params,1) ~= 8 + params = params.'; + end + tensor = params(2:7,:); + S0 = params(1,:); + f = params(8,:); +end + +signal = S0*(f*delta_fun(bval) + (1-f)*exp(-bmat*tensor.')); + diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun_seg_constr.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun_seg_constr.m new file mode 100644 index 00000000..71098dc6 --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimdtifun_seg_constr.m @@ -0,0 +1,59 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to calculate the IVIM-DTI signal given the lower triangular +% matrix (from the Cholesky decomposition) of the diffusion tensor and the +% perfusion fraction f. The segmented IVIM-DTI equation is used. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input values: +% - params: model parameters +% params(1) = S0 signal +% params(2:7) = lower triangular matrix of the +% diffusion tensor +% params(8) = perfusion fraction f +% +% - diffparams: structure containing the fields 'bval' with the +% diffusion b-values and 'diffdir', containing the +% gradient directinos. +% +% First, the diffusion tensor is calculated from the lower triangular +% matrix using the function lower_triangular2tensor. +% Then, the signal is calculated using the function ivimdtifun_seg: +% S = S0*(f*delta_fun(bval) + (1-f)*exp(-bmat*diff_tensor)) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function signal = ivimdtifun_seg_constr(params, diffparams) + +%check input parameters +if ndims(params) ~= 2 %#ok + error('Wrong matrix size of input parameters') +end + +if length(params)>8 + warning('Too many input parameters provided, additional parameters are ignored') +end + +if ~isstruct(diffparams) + error('Diffparams must be a structure containing b-values and diffusion gradient directions') +end +if (~isfield(diffparams, 'bval') || ~isfield(diffparams, 'diffdir')) + error('Diffparams must contain b-values and diffusion gradient directions') +end + +%Calculate tensor from Cholesky decomposition +if min(size(params)) == 1 + %single voxel + params(2:7) = lower_triangular2tensor(params(2:7)); +else + if size(params,1) ~= 8 + params = params.'; + end + params(2:7,:) = lower_triangular2tensor(params(2:7,:)); +end + +%Calculate signal +signal = ivimdtifun_seg(params, diffparams); + + + diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimfun.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimfun.m new file mode 100644 index 00000000..f78d2f92 --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimfun.m @@ -0,0 +1,72 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to calculate the IVIM signal given the IVIM parameters and +% b-values. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input values: +% - params: vector containing the IVIM parameters in the order +% [S0, D, f, Ds] +% +% - bval: b-values +% In case of standard IVIM signal calculation, bval +% needs to be a vector containing the b-values. +% In case of two-step IVIM fitting, bval is a +% structure with fields bval.D_fix (fixed value for +% the diffusion coefficient) and bval.bval (diffusion +% b-values). +% +% The signal is calculated according to +% S = S0*(f*exp(-b*Ds) + (1-f)*exp(-b*D)) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function signal = ivimfun(params, bval) + +switch length(params) + case 4 + %standard IVIM signal calculation + S0 = params(1); + D = params(2); + f = params(3); + Ds = params(4); + + if ~isnumeric(bval) + error('No valid b-values provided') + end + case 3 + %two-step approach for IVIM fitting. D is fixed during the fit and is + %provided as second argument in 'bval'. + S0 = params(1); + f = params(2); + Ds = params(3); + + if isstruct(bval) + D = bval.D_fix; + bval = bval.bval; + else + error('Incorrect input. No value for D is provided.') + end + + case 2 + %two-step approach for IVIM fitting with fixed D*. + %D is also fixed during the fit (pre-calculated) and is provided as + %second argument in 'bval'. + %D* is also stored in the structure 'bval'. + S0 = params(1); + f = params(2); + if isstruct(bval) + D = bval.D_fix; + Ds = bval.Ds_fix; + bval = bval.bval; + else + error('Incorrect input. No value for D or D* is provided.') + end + otherwise + error('Wrong number of input parameters for IVIM equation') +end + +%scale b-value +bval = bval_scaling(bval); + +signal = S0*(f*exp(-bval*Ds) + (1-f)*exp(-bval*D)); \ No newline at end of file diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimfun_seg.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimfun_seg.m new file mode 100644 index 00000000..9261f618 --- /dev/null +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/models_ivim_dti/ivimfun_seg.m @@ -0,0 +1,37 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to calculate the segmented IVIM signal given the IVIM parameters +% and b-values using the segmented IVIM equation. +% +% The pseudo-diffusion coefficient D* is not calculated with this method. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input values: +% - params: vector containing the segmented IVIM parameters in +% the order [S0, D, f] +% +% - bval: b-values +% bval is a vector containing the b-values. +% +% The signal is calculated according to +% S = S0*(f*delta_fun(bval) + (1-f)*exp(-b*D)) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function signal = ivimfun_seg(params, bval) + +if length(params) < 3 + error('Wrong number of input parameters for segmented IVIM equation') +end +if length(params) > 3 + warning('Too many input parameters provided, additional parameters are ignored') +end + +S0 = params(1); +D = params(2); +f = params(3); + +%scale b-value +bval = bval_scaling(bval); + +signal = S0*(f*delta_fun(bval) + (1-f)*exp(-bval*D)); \ No newline at end of file From 621d2f7d0758c354b4d9655fea947a9ecc8fbc99 Mon Sep 17 00:00:00 2001 From: s-rauh Date: Wed, 5 Aug 2026 08:41:36 +0200 Subject: [PATCH 3/5] Bug fix in original/fitting/SR_LUMC in remove_zeros.m --- .../helper_fun/remove_zeros.m | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/remove_zeros.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/remove_zeros.m index 1e7a0bf4..19f9c152 100644 --- a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/remove_zeros.m +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/remove_zeros.m @@ -5,21 +5,23 @@ options.bval = [] end +idx = data == 0; + if isstruct(diffparams) %IVIM-DTI - diffparams.bval(data==0) = []; - diffparams.diffdir(data==0,:) = []; + diffparams.bval(idx) = []; + diffparams.diffdir(idx,:) = []; elseif isvector(diffparams) %IVIM, diffparams = bval - diffparams(data==0) = []; + diffparams(idx) = []; else %DTI, diffparams = b-matrix - diffparams(data==0,:) = []; + diffparams(idx,:) = []; end -data(data==0) = []; +data(idx) = []; if ~isempty(options.bval) - options.bval(data==0) = []; + options.bval(idx) = []; varargout{1} = options.bval; end From 35dd735367d5a9f916c109f491c87b89f65e11ec Mon Sep 17 00:00:00 2001 From: s-rauh Date: Wed, 5 Aug 2026 08:43:52 +0200 Subject: [PATCH 4/5] Bugfix in original/fitting/SR_LUMC fit_ivim.m: Correct b-value selection in loop for -segmented fit option --- .../fitting/SR_LUMC/diffusion_mri_tools-main/fit_ivim.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_ivim.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_ivim.m index e766dcb7..89f927d0 100644 --- a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_ivim.m +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/fit_ivim.m @@ -186,7 +186,7 @@ case 'segmented' %optional: select desired b-value range (if full range was %acquired) - if options.seg_data && v == 1 + if options.seg_data tmpdat = [tmpdat(tmpb==0); tmpdat(tmpb>=options.bcut)]; tmpb = [tmpb(tmpb==0); tmpb(tmpb>=options.bcut)]; end From 727a7349c80b8ac9f410da8eba742b216210a06c Mon Sep 17 00:00:00 2001 From: s-rauh Date: Wed, 5 Aug 2026 08:57:57 +0200 Subject: [PATCH 5/5] Fix potential bug in norm_diffdata in original/fitting/SR_LUMC --- .../diffusion_mri_tools-main/helper_fun/norm_diffdata.m | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/norm_diffdata.m b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/norm_diffdata.m index 59e33a02..ff5f7e43 100644 --- a/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/norm_diffdata.m +++ b/src/original/fitting/SR_LUMC/diffusion_mri_tools-main/helper_fun/norm_diffdata.m @@ -19,15 +19,15 @@ end % check if b0 data is provided -if min(bval) >= 1 +if min(bval) > 0 fprintf('No b = 0 data provided. Normalize to minimum b-value data: %.2f s/mm². \n',... min(bval)) - mS0 = mean(data((bval == min(bval)),:), 1); else - fprintf('Normalize data to b = %.2f data. \n', min(bval)) - mS0 = mean(data(bval<1, :), 1); + fprintf('Normalize data to b = %.2f data. \n', min(bval)) end +mS0 = mean(data((bval == min(bval)),:), 1); + normalized = data ./ mS0; normalized(isnan(normalized)) = 0;