Fitting code contribution: SR_LUMC Matlab code for IVIM, DTI, and IVIM-DTI fitting

src/original/fitting/SR_LUMC/diffusion_mri_tools-main/README.md

@@ -0,0 +1,26 @@
+# Diffusion MRI Tools
+
+Diffusion MRI Tools contains Matlab code to analyze and process diffusion MRI data. 
+
+Author: Susi Rauh (s.s.rauh@lumc.nl) 
+
+Run setup.m once before using any project functions to add the subfolders to the path. 
+
+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. 

src/original/fitting/SR_LUMC/diffusion_mri_tools-main/analyze_data/calc_dti_parameters.m

@@ -0,0 +1,90 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% 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;
+    dti_params.FA(~isfinite(dti_params.FA)) = 0;
+end
+
+fprintf('Tensor calculation finished in %.2f seconds. \n', toc);
+end

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