From a2d510b2a1728ed5d81e3245d0458dc92664e8be Mon Sep 17 00:00:00 2001 From: Paulien Voorter Date: Sat, 8 Nov 2025 21:31:36 +0100 Subject: [PATCH] fix wrong initial guess --- src/original/PV_MUMC/two_step_IVIM_fit.py | 36 +++++++++++------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/src/original/PV_MUMC/two_step_IVIM_fit.py b/src/original/PV_MUMC/two_step_IVIM_fit.py index 7007a606..ac413223 100644 --- a/src/original/PV_MUMC/two_step_IVIM_fit.py +++ b/src/original/PV_MUMC/two_step_IVIM_fit.py @@ -1,6 +1,6 @@ """ January 2022 by Paulien Voorter -p.voorter@maastrichtuniversity.nl +p.voorter@maastrichtuniversity.nl https://www.github.com/paulienvoorter requirements: @@ -20,11 +20,11 @@ def two_exp_noS0(bvalues, Dpar, Fmv, Dmv): """ bi-exponential IVIM function, and S0 set to 1""" return Fmv * np.exp(-bvalues * Dmv) + (1 - Fmv ) * np.exp(-bvalues * Dpar) - + def two_exp(bvalues, S0, Dpar, Fmv, Dmv): """ bi-exponential IVIM function""" return S0 * (Fmv * np.exp(-bvalues * Dmv) + (1 - Fmv ) * np.exp(-bvalues * Dpar)) - + def fit_least_squares_array(bvalues, dw_data, fitS0=False, bounds=([0.9, 0.0001, 0.0, 0.0025], [1.1, 0.0025, 0.5, 0.2]), cutoff=200): @@ -34,8 +34,8 @@ def fit_least_squares_array(bvalues, dw_data, fitS0=False, bounds=([0.9, 0.0001, is done on an array. :param bvalues: 1D Array with the b-values :param dw_data: 2D Array with diffusion-weighted signal in different voxels at different b-values - :param bounds: Array with fit bounds ([S0min, Dparmin, Fmvmin, Dmvmin],[S0max, Dparmax, Fmvmax, Dmvmax]). - :param cutoff: cutoff b-value used in step 1 + :param bounds: Array with fit bounds ([S0min, Dparmin, Fmvmin, Dmvmin],[S0max, Dparmax, Fmvmax, Dmvmax]). + :param cutoff: cutoff b-value used in step 1 :return Dpar: 1D Array with Dpar in each voxel :return Fmv: 1D Array with Fmv in each voxel :return Dmv: 1D Array with Dmv in each voxel @@ -60,48 +60,48 @@ def fit_least_squares(bvalues, dw_data, S0_output=False, fitS0=False, is done on an array. It fits a single curve :param bvalues: 1D Array with the b-values :param dw_data: 1D Array with diffusion-weighted signal in different voxels at different b-values - :param S0_output: Boolean determining whether to output (often a dummy) variable S0; + :param S0_output: Boolean determining whether to output (often a dummy) variable S0; :param fitS0: Boolean determining whether to fix S0 to 1; - :param bounds: Array with fit bounds ([S0min, Dparmin, Fmvmin, Dmvmin],[S0max, Dparmax, Fmvmax, Dmvmax]). - :param cutoff: cutoff b-value used in step 1 + :param bounds: Array with fit bounds ([S0min, Dparmin, Fmvmin, Dmvmin],[S0max, Dparmax, Fmvmax, Dmvmax]). + :param cutoff: cutoff b-value used in step 1 :return S0: optional 1D Array with S0 in each voxel :return Dpar: scalar with Dpar of the specific voxel :return Fmv: scalar with Fmv of the specific voxel :return Dmv: scalar with Dmv of the specific voxel """ - + #try: def monofit(bvalues, Dpar, Fmv): return (1-Fmv)*np.exp(-bvalues * Dpar) - + high_b = bvalues[bvalues >= cutoff] high_dw_data = dw_data[bvalues >= cutoff] boundsmonoexp = ([bounds[0][1], bounds[0][2]], [bounds[1][1], bounds[1][2]]) - params, _ = curve_fit(monofit, high_b, high_dw_data, p0=[(bounds[1][1]-bounds[0][1])/2,(bounds[1][2]-bounds[0][2])/2], bounds=boundsmonoexp) + params, _ = curve_fit(monofit, high_b, high_dw_data, p0=[(bounds[1][1]+bounds[0][1])/2,(bounds[1][2]+bounds[0][2])/2], bounds=boundsmonoexp) Dpar1 = params[0] if not fitS0: boundsupdated=([bounds[0][2] , bounds[0][3] ], - [bounds[1][2] , bounds[1][3] ]) + [bounds[1][2] , bounds[1][3] ]) params, _ = curve_fit(lambda b, Fmv, Dmv: two_exp_noS0(b, Dpar1, Fmv, Dmv), bvalues, dw_data, p0=[(bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated) Fmv, Dmv = params[0], params[1] #when the fraction of a compartment equals zero (or very very small), the corresponding diffusivity is non-existing (=NaN) #if Fmv < 1e-4: # Dmv = float("NaN") - - else: + + else: boundsupdated = ([bounds[0][0] , bounds[0][2] , bounds[0][3] ], - [bounds[1][0] , bounds[1][2] , bounds[1][3] ]) + [bounds[1][0] , bounds[1][2] , bounds[1][3] ]) params, _ = curve_fit(lambda b, S0, Fmv, Dmv: two_exp(b, S0, Dpar1, Fmv, Dmv), bvalues, dw_data, p0=[1, (bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated) S0 = params[0] Fmv, Dmv = params[1] , params[2] #when the fraction of a compartment equals zero (or very very small), the corresponding diffusivity is non-existing (=NaN) #if Fmv < 1e-4: - # Dmv = float("NaN") - + # Dmv = float("NaN") + if S0_output: return Dpar1, Fmv, Dmv, S0 else: return Dpar1, Fmv, Dmv - +