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

[s-rauh]

Fitting code contribution in src/original/fitting. The code contains Matlab fitting routines for IVIM, IVIM-DTI, and DTI.

[Devguru-codes]

I went through the code and spotted a few things:

1. remove_zeros.m - zero mask applied too late
   The zeros get removed from data on line 13, but then the options.bval filtering on line 16 checks data==0 after the removal - so there are no zeros left and it's a no-op. The mask should be saved before modifying data
   code -
   % something like:
   idx = data==0;
   data(idx) = [];
   if ~isempty(options.bval)
   options.bval(idx) = [];
   varargout{1} = options.bval;
   end
   Doesn't affect the current fitting routines since none of them pass options.bval, but it'll bite anyone who tries to use that feature.

2. fit_ivim.m - segmented - v == 1 guard on data selection
   code -
   case 'segmented'
   if options.seg_data && v == 1
   tmpdat = [tmpdat(tmpb==0); tmpdat(tmpb>=options.bcut)];
   tmpb = [tmpb(tmpb==0); tmpb(tmpb>=options.bcut)];
   end

tmpdat and tmpb are re-initialized from scratch at the start of every loop iteration (tmpdat = datafit(:,v); tmpb = bval;), so the v == 1 condition means only the first voxel gets the b-value range selection. All other voxels fit the full range. Compare with how v == 1 is correctly used elsewhere in the same file - for one-time initialization of bounds that persist:
code -
case 'two_step'
if v == 1 % Correct: bounds persist across iterations
lb_twostep = lb([1 3 4]); ...
end

Looks like the v == 1 condition was carried over from the two_step case where it's used for one-time bounds init, but here the data selection needs to happen every iteration.

3. norm_diffdata.m - wrong normalization after b-value scaling
   if min(bval) >= 1
   mS0 = mean(data((bval == min(bval)),:), 1); % normalize to min bval
   else
   mS0 = mean(data(bval<1, :), 1); % normalize to "b=0"
   end

This is called after bval_scaling() divides b-values by 1000 (when max(bval) > 100). So typical clinical b-values [0, 50, 200, 800] become [0, 0.05, 0.2, 0.8]. Since min(bval) = 0 < 1, it uses bval < 1, which selects ALL b-values (they're all < 1 after scaling). The normalization factor becomes the mean of ALL signals across ALL b-values, not just b=0.

This effectively destroys the normalization - it should be bval == min(bval) or bval < 0.001 in the else branch.
Also, I noticed the fitting routines call @ivimfun, @dtifun, etc. but those model function files don't seem to be included - are they in a models/ directory that didn't get added?

[s-rauh]

Thanks for looking into the code. I adjusted a few things:

Models: The folder "models/" including the ivimfun, dtifun etc was indeed missing, thanks for spotting this. "models" appears in gitignore (which I missed). I renamed the folder and added it.

1. remove_zeros.m: I adjusted the code accordingly
2. fit_ivim.m option 'segmented': I corrected the data selection for this option
3. norm_diffdata.m: the normalization is performed before b-value scaling, so it should not be a problem. However, as the differentiation was redundant, I adjusted the code to always use normalization to the minimum b-value:
   if min(bval) > 0
   fprintf('No b = 0 data provided. Normalize to minimum b-value data: %.2f s/mm². \n',...
   min(bval))
   else
   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;

[Devguru-codes]

Good work. I just want to point to 2 things here -

1. MATLAB Search Path / Subfolder structure: Since helper_fun/ and models_ivim_dti/ are in subdirectories, running fit_ivim (or fit_dti/fit_ivimdti) will throw an Undefined function 'ivimfun' or 'reshape_diffdata_for_fit' error unless the user explicitly ran addpath(genpath(...)) beforehand.
   Adding something like addpath(genpath(fileparts(mfilename('fullpath')))) at the top of the entry points (or noting addpath(genpath(...)) in the README) would make sure MATLAB picks up the subfolders automatically.

2. Zero-division / NaN in calc_dti_parameters.m: In calc_dti_parameters.m (line 100), the FA formula divides by sqrt(sum(EigenValues.^2)). For zero/background voxels where all eigenvalues are 0, this evaluates to 0 / 0 = NaN.
   Since NaN > 1 returns false in MATLAB, the options.clip check on line 129 doesn't catch it, leaving NaNs in the FA map. A quick zero-guard before division would fix it.

Please look into these 2. Thank you.