Skip to content

fix wrong initial guess #128

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
+18 −18
Open

fix wrong initial guess #128

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
36 changes: 18 additions & 18 deletions src/original/PV_MUMC/two_step_IVIM_fit.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
"""
January 2022 by Paulien Voorter
p.voorter@maastrichtuniversity.nl
p.voorter@maastrichtuniversity.nl
https://www.github.com/paulienvoorter

requirements:
Expand All @@ -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):
Expand All @@ -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
Expand All @@ -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



Loading