Merge pull request #15 from OSIPI/OGC_AmsterdamUMC_standardized

Ogc amsterdam umc standardized

Author: oliverchampion

.gitignore

@@ -1,8 +1,9 @@
 .idea*
 .github*
 __pycache__/
-
-
+*.nii.gz
+bvals.txt
+*json
 # Unit test / coverage reports
 .tox/
 .coverage

phantoms/MR_XCAT_qMRI/sim_ivim_sig.py

@@ -8,7 +8,7 @@
 # code adapted from MAtlab code by Eric Schrauben: https://github.com/schrau24/XCAT-ERIC
 # This code generates a 4D IVIM phantom as nifti file
 
-def phantom(bvalue, noise, TR=3000, TE=40, motion=False, rician=False, interleaved=False):
+def phantom(bvalue, noise, TR=8000, TE=80, motion=False, rician=False, interleaved=False):
     np.random.seed(42)
     if motion:
         states = range(1,21)
@@ -367,7 +367,8 @@ def XCAT_to_MR_DCE(XCAT, TR, TE, bvalue, D, f, Ds, b0=3, ivim_cont = True):
     bvalue = np.array([0., 1, 2, 5, 10, 20, 30, 50, 75, 100, 150, 250, 350, 400, 550, 700, 850, 1000])
     noise = 0.0005
     motion = False
-    sig, XCAT, Dim, fim, Dpim, legend = phantom(bvalue, noise, motion=motion, interleaved=False)
+    interleaved = False
+    sig, XCAT, Dim, fim, Dpim, legend = phantom(bvalue, noise, motion=motion, interleaved=interleaved)
     # sig = np.flip(sig,axis=0)
     # sig = np.flip(sig,axis=1)
     res=np.eye(4)
@@ -402,12 +403,16 @@ def XCAT_to_MR_DCE(XCAT, TR, TE, bvalue, D, f, Ds, b0=3, ivim_cont = True):
     nifti_img = nib.Nifti1Image(sig, affine=res)  # Replace affine if necessary
     # Save the NIfTI image to a file
     nifti_img.header.set_data_dtype(np.float64)
-    if motion:
+    if not motion:
+        output_file = 'output.nii.gz'  # Replace with your desired output file name
+    elif interleaved:
         output_file = 'output_resp_int.nii.gz'  # Replace with your desired output file name
     else:
-        output_file = 'output.nii.gz'  # Replace with your desired output file name
+        output_file = 'output_resp.nii.gz'  # Replace with your desired output file name
+
     nib.save(nifti_img, output_file)
 
+
     nifti_img = nib.Nifti1Image(XCAT, affine=res)  # Replace affine if necessary
     # Save the NIfTI image to a file
     output_file = 'output_xcat.nii.gz'  # Replace with your desired output file name

src/original/OGC_AmsterdamUMC/LSQ_fitting.py

@@ -66,7 +66,7 @@ def order(Dt, Fp, Dp, S0=None):
         return Dt, Fp, Dp, S0
 
 
-def fit_segmented_array(bvalues, dw_data, njobs=4, bounds=([0, 0, 0.005],[0.005, 0.7, 0.2]), cutoff=75):
+def fit_segmented_array(bvalues, dw_data, njobs=4, bounds=([0, 0, 0.005],[0.005, 0.7, 0.2]), cutoff=75,p0=[0.001, 0.1, 0.01,1]):
     """
     This is an implementation of the segmented fit, in which we first estimate D using a curve fit to b-values>cutoff;
     then estimate f from the fitted S0 and the measured S0 and finally estimate D* while fixing D and f. This fit
@@ -90,7 +90,7 @@ def fit_segmented_array(bvalues, dw_data, njobs=4, bounds=([0, 0, 0.005],[0.005,
         try:
             # define the parallel function
             def parfun(i):
-                return fit_segmented(bvalues, dw_data[i, :], bounds=bounds, cutoff=cutoff)
+                return fit_segmented(bvalues, dw_data[i, :], bounds=bounds, cutoff=cutoff,p0=p0)
 
             output = Parallel(n_jobs=njobs)(delayed(parfun)(i) for i in tqdm(range(len(dw_data)), position=0, leave=True))
             Dt, Fp, Dp = np.transpose(output)
@@ -107,11 +107,11 @@ def parfun(i):
         Fp = np.zeros(len(dw_data))
         for i in tqdm(range(len(dw_data)), position=0, leave=True):
             # fill arrays with fit results on a per voxel base:
-            Dt[i], Fp[i], Dp[i] = fit_segmented(bvalues, dw_data[i, :], bounds=bounds, cutoff=cutoff)
+            Dt[i], Fp[i], Dp[i] = fit_segmented(bvalues, dw_data[i, :], bounds=bounds, cutoff=cutoff,p0=p0)
     return [Dt, Fp, Dp, S0]
 
 
-def fit_segmented(bvalues, dw_data, bounds=([0, 0, 0.005],[0.005, 0.7, 0.2]), cutoff=75):
+def fit_segmented(bvalues, dw_data, bounds=([0, 0, 0.005],[0.005, 0.7, 0.2]), cutoff=75,p0=[0.001, 0.1, 0.01,1]):
     """
     This is an implementation of the segmented fit, in which we first estimate D using a curve fit to b-values>cutoff;
     then estimate f from the fitted S0 and the measured S0 and finally estimate D* while fixing D and f.
@@ -124,23 +124,24 @@ def fit_segmented(bvalues, dw_data, bounds=([0, 0, 0.005],[0.005, 0.7, 0.2]), cu
     :return Dp: Fitted Dp
     :return S0: Fitted S0
     """
+    p0 = [p0[0] * 1000, p0[1] * 10, p0[2] * 10, p0[3]]
     try:
         # determine high b-values and data for D
         high_b = bvalues[bvalues >= cutoff]
         high_dw_data = dw_data[bvalues >= cutoff]
         # correct the bounds. Note that S0 bounds determine the max and min of f
-        bounds1 = ([bounds[0][0] * 1000., 1 - bounds[1][1]], [bounds[1][0] * 1000., 1. - bounds[0][
+        bounds1 = ([bounds[0][0] * 1000., 0.7 - bounds[1][1]], [bounds[1][0] * 1000., 1.3 - bounds[0][
             1]])  # By bounding S0 like this, we effectively insert the boundaries of f
         # fit for S0' and D
         params, _ = curve_fit(lambda b, Dt, int: int * np.exp(-b * Dt / 1000), high_b, high_dw_data,
-                              p0=(1, 1),
+                              p0=(p0[0], p0[3]-p0[1]/10),
                               bounds=bounds1)
         Dt, Fp = params[0] / 1000, 1 - params[1]
         # remove the diffusion part to only keep the pseudo-diffusion
         dw_data_remaining = dw_data - (1 - Fp) * np.exp(-bvalues * Dt)
-        bounds2 = (bounds[0][2], bounds[1][2])
+        bounds2 = (bounds[0][2]*10, bounds[1][2]*10)
         # fit for D*
-        params, _ = curve_fit(lambda b, Dp: Fp * np.exp(-b * Dp), bvalues, dw_data_remaining, p0=(0.1), bounds=bounds2)
+        params, _ = curve_fit(lambda b, Dp: Fp * np.exp(-b * Dp), bvalues, dw_data_remaining, p0=(p0[2]), bounds=bounds2)
         Dp = params[0]
         return Dt, Fp, Dp
     except:
@@ -150,7 +151,7 @@ def fit_segmented(bvalues, dw_data, bounds=([0, 0, 0.005],[0.005, 0.7, 0.2]), cu
 
 
 def fit_least_squares_array(bvalues, dw_data, S0_output=True, fitS0=True, njobs=4,
-                            bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3])):
+                            bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3]),p0=[0.001, 0.1, 0.01, 1]):
     """
     This is an implementation of the conventional IVIM fit. It is fitted in array form.
     :param bvalues: 1D Array with the b-values
@@ -175,7 +176,7 @@ def fit_least_squares_array(bvalues, dw_data, S0_output=True, fitS0=True, njobs=
             try:
                 # defining parallel function
                 def parfun(i):
-                    return fit_least_squares(bvalues, dw_data[i, :], S0_output=S0_output, fitS0=fitS0, bounds=bounds)
+                    return fit_least_squares(bvalues, dw_data[i, :], S0_output=S0_output, fitS0=fitS0, bounds=bounds,p0=p0)
 
                 output = Parallel(n_jobs=njobs)(delayed(parfun)(i) for i in tqdm(range(len(dw_data)), position=0, leave=True))
                 Dt, Fp, Dp, S0 = np.transpose(output)
@@ -192,14 +193,14 @@ def parfun(i):
             # running in a single loop and filling arrays
             for i in tqdm(range(len(dw_data)), position=0, leave=True):
                 Dt[i], Fp[i], Dp[i], S0[i] = fit_least_squares(bvalues, dw_data[i, :], S0_output=S0_output, fitS0=fitS0,
-                                                               bounds=bounds)
+                                                               bounds=bounds,p0=p0)
         return [Dt, Fp, Dp, S0]
     else:
         # if S0 is not exported
         if njobs > 1:
             try:
                 def parfun(i):
-                    return fit_least_squares(bvalues, dw_data[i, :], fitS0=fitS0, bounds=bounds)
+                    return fit_least_squares(bvalues, dw_data[i, :], fitS0=fitS0, bounds=bounds,p0=p0)
 
                 output = Parallel(n_jobs=njobs)(delayed(parfun)(i) for i in tqdm(range(len(dw_data)), position=0, leave=True))
                 Dt, Fp, Dp = np.transpose(output)
@@ -213,12 +214,12 @@ def parfun(i):
             Fp = np.zeros(len(dw_data))
             for i in range(len(dw_data)):
                 Dt[i], Fp[i], Dp[i] = fit_least_squares(bvalues, dw_data[i, :], S0_output=S0_output, fitS0=fitS0,
-                                                        bounds=bounds)
+                                                        bounds=bounds,p0=p0)
         return [Dt, Fp, Dp]
 
 
 def fit_least_squares(bvalues, dw_data, S0_output=False, fitS0=True,
-                      bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3])):
+                      bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3]), p0=[0.001, 0.1, 0.01, 1]):
     """
     This is an implementation of the conventional IVIM fit. It fits a single curve
     :param bvalues: Array with the b-values
@@ -236,13 +237,15 @@ def fit_least_squares(bvalues, dw_data, S0_output=False, fitS0=True,
             # bounds are rescaled such that each parameter changes at roughly the same rate to help fitting.
             bounds = ([bounds[0][0] * 1000, bounds[0][1] * 10, bounds[0][2] * 10],
                       [bounds[1][0] * 1000, bounds[1][1] * 10, bounds[1][2] * 10])
-            params, _ = curve_fit(ivimN_noS0, bvalues, dw_data, p0=[1, 1, 0.1], bounds=bounds)
+            p0=[p0[0]*1000,p0[1]*10,p0[2]*10]
+            params, _ = curve_fit(ivimN_noS0, bvalues, dw_data, p0=p0, bounds=bounds)
             S0 = 1
         else:
             # bounds are rescaled such that each parameter changes at roughly the same rate to help fitting.
             bounds = ([bounds[0][0] * 1000, bounds[0][1] * 10, bounds[0][2] * 10, bounds[0][3]],
                       [bounds[1][0] * 1000, bounds[1][1] * 10, bounds[1][2] * 10, bounds[1][3]])
-            params, _ = curve_fit(ivimN, bvalues, dw_data, p0=[1, 1, 0.1, 1], bounds=bounds)
+            p0=[p0[0]*1000,p0[1]*10,p0[2]*10,p0[3]]
+            params, _ = curve_fit(ivimN, bvalues, dw_data, p0=p0, bounds=bounds)
             S0 = params[3]
         # correct for the rescaling of parameters
         Dt, Fp, Dp = params[0] / 1000, params[1] / 10, params[2] / 10
@@ -463,6 +466,20 @@ def fit_segmented_tri_exp(bvalues, dw_data, bounds=([0, 0, 0, 0.005, 0, 0.06], [
         return 0., 0., 0., 0., 0., 0.
 
 
+def neg_log_likelihood(p, bvalues, dw_data):
+    """
+    This function determines the negative of the log of the likelihood of parameters p, given the data dw_data for the Bayesian fit
+    :param p: 1D Array with the estimates of D, f, D* and (optionally) S0
+    :param bvalues: 1D array with b-values
+    :param dw_data: 1D Array diffusion-weighted data
+    :returns: the log-likelihood of the parameters given the data
+    """
+    if len(p) == 4:
+        return 0.5 * (len(bvalues) + 1) * np.log(
+            np.sum((ivim(bvalues, p[0], p[1], p[2], p[3]) - dw_data) ** 2))  # 0.5*sum simplified
+    else:
+        return 0.5 * (len(bvalues) + 1) * np.log(
+            np.sum((ivim(bvalues, p[0], p[1], p[2], 1) - dw_data) ** 2))  # 0.5*sum simplified
 def empirical_neg_log_prior(Dt0, Fp0, Dp0, S00=None):
     """
     This function determines the negative of the log of the empirical prior probability of the IVIM parameters
@@ -516,23 +533,6 @@ def neg_log_prior(p):
 
     return neg_log_prior
 
-
-def neg_log_likelihood(p, bvalues, dw_data):
-    """
-    This function determines the negative of the log of the likelihood of parameters p, given the data dw_data for the Bayesian fit
-    :param p: 1D Array with the estimates of D, f, D* and (optionally) S0
-    :param bvalues: 1D array with b-values
-    :param dw_data: 1D Array diffusion-weighted data
-    :returns: the log-likelihood of the parameters given the data
-    """
-    if len(p) == 4:
-        return 0.5 * (len(bvalues) + 1) * np.log(
-            np.sum((ivim(bvalues, p[0], p[1], p[2], p[3]) - dw_data) ** 2))  # 0.5*sum simplified
-    else:
-        return 0.5 * (len(bvalues) + 1) * np.log(
-            np.sum((ivim(bvalues, p[0], p[1], p[2], 1) - dw_data) ** 2))  # 0.5*sum simplified
-
-
 def neg_log_posterior(p, bvalues, dw_data, neg_log_prior):
     """
     This function determines the negative of the log of the likelihood of parameters p, given the prior likelihood and the data
@@ -545,6 +545,39 @@ def neg_log_posterior(p, bvalues, dw_data, neg_log_prior):
     return neg_log_likelihood(p, bvalues, dw_data) + neg_log_prior(p)
 
 
+def flat_neg_log_prior(Dt_range, Fp_range, Dp_range, S0_range=None):
+    """
+    This function determines the negative of the log of the empirical prior probability of the IVIM parameters
+    :param Dt0: 1D Array with the initial D estimates
+    :param Dt0: 1D Array with the initial f estimates
+    :param Dt0: 1D Array with the initial D* estimates
+    :param Dt0: 1D Array with the initial S0 estimates (optional)
+    """
+    def neg_log_prior(p):
+        # depends on whether S0 is fitted or not
+        if len(p) == 4:
+            Dt, Fp, Dp, S0 = p[0], p[1], p[2], p[3]
+        else:
+            Dt, Fp, Dp = p[0], p[1], p[2]
+        # make D*<D very unlikely
+        if (Dp < Dt):
+            return 1e3
+        else:
+            # determine and return the prior for D, f and D* (and S0)
+            if len(p) == 4:
+                if Dt_range[0] < Dt < Dt_range[1] and Fp_range[0] < Fp < Fp_range[1] and Dp_range[0] < Dp < Dp_range[1]: # and S0_range[0] < S0 < S0_range[1]: << not sure whether this helps. Technically it should be here
+                    return 0
+                else:
+                    return 1e3
+            else:
+                if Dt_range[0] < Dt < Dt_range[1] and Fp_range[0] < Fp < Fp_range[1] and Dp_range[0] < Dp < Dp_range[1]:
+                    return 0
+                else:
+                    return 1e3
+
+    return neg_log_prior
+
+
 def fit_bayesian_array(bvalues, dw_data, paramslsq, arg):
     """
     This is an implementation of the Bayesian IVIM fit for arrays. The fit is taken from Barbieri et al. which was
@@ -563,7 +596,6 @@ def fit_bayesian_array(bvalues, dw_data, paramslsq, arg):
     :return Dp: Array with Dp in each voxel
     :return S0: Array with S0 in each voxel
     """
-    arg = checkarg_lsq(arg)
     # fill out missing args
     Dt0, Fp0, Dp0, S00 = paramslsq
     # determine prior
@@ -606,7 +638,7 @@ def parfun(i):
     return Dt_pred, Fp_pred, Dp_pred, S0_pred
 
 
-def fit_bayesian(bvalues, dw_data, neg_log_prior, x0=[0.001, 0.2, 0.05], fitS0=True):
+def fit_bayesian(bvalues, dw_data, neg_log_prior, x0=[0.001, 0.2, 0.05, 1], fitS0=True):
     '''
     This is an implementation of the Bayesian IVIM fit. It returns the Maximum a posterior probability.
     The fit is taken from Barbieri et al. which was initially introduced in http://arxiv.org/10.1002/mrm.25765 and
@@ -623,7 +655,7 @@ def fit_bayesian(bvalues, dw_data, neg_log_prior, x0=[0.001, 0.2, 0.05], fitS0=T
     '''
     try:
         # define fit bounds
-        bounds = [(0, 0.005), (0, 0.7), (0.005, 0.2), (0, 2.5)]
+        bounds = [(0, 0.005), (0, 1.5), (0, 2), (0, 2.5)]
         # Find the Maximum a posterior probability (MAP) by minimising the negative log of the posterior
         if fitS0:
             params = minimize(neg_log_posterior, x0=x0, args=(bvalues, dw_data, neg_log_prior), bounds=bounds)
@@ -704,7 +736,7 @@ def goodness_of_fit(bvalues, Dt, Fp, Dp, S0, dw_data, Fp2=None, Dp2=None):
         # plt.ion()
         # plt.show()
         # print(R2[vox])
-    return R2, adjust
+    return R2, adjusted_R2
 # ed_R2
 
 def MSE(bvalues, Dt, Fp, Dp, S0, dw_data):

src/original/PV_MUMC/two_step_IVIM_fit.py

@@ -50,7 +50,7 @@ def fit_least_squares_array(bvalues, dw_data, fitS0=True, bounds=([0.9, 0.0001,
     Fmv = np.zeros(len(dw_data))
     for i in tqdm.tqdm(range(len(dw_data)), position=0, leave=True):
         # fill arrays with fit results on a per voxel base:
-        Dpar[i], Fmv[i], Dmv[i], S0[i] = fit_least_squares(bvalues, dw_data[i, :], S0_output=True, fitS0=fitS0, bounds=bounds)
+        Dpar[i], Fmv[i], Dmv[i], S0[i] = fit_least_squares(bvalues, dw_data[i, :], S0_output=False, fitS0=fitS0, bounds=bounds)
     return [Dpar, Fmv, Dmv, S0]
 
 

src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py

@@ -0,0 +1,79 @@
+from src.wrappers.OsipiBase import OsipiBase
+from src.original.OGC_AmsterdamUMC.LSQ_fitting import flat_neg_log_prior, fit_bayesian, empirical_neg_log_prior, fit_segmented
+import numpy as np
+
+class OGC_AmsterdamUMC_Bayesian_biexp(OsipiBase):
+    """
+    Bayesian Bi-exponential fitting algorithm by Oliver Gurney-Champion, Amsterdam UMC
+    """
+
+    # I'm thinking that we define default attributes for each submission like this
+    # And in __init__, we can call the OsipiBase control functions to check whether
+    # the user inputs fulfil the requirements
+
+    # Some basic stuff that identifies the algorithm
+    id_author = "Oliver Gurney Champion, Amsterdam UMC"
+    id_algorithm_type = "Bi-exponential fit"
+    id_return_parameters = "f, D*, D, S0"
+    id_units = "seconds per milli metre squared or milliseconds per micro metre squared"
+
+    # Algorithm requirements
+    required_bvalues = 4
+    required_thresholds = [0,
+                           0]  # Interval from "at least" to "at most", in case submissions allow a custom number of thresholds
+    required_bounds = False
+    required_bounds_optional = True  # Bounds may not be required but are optional
+    required_initial_guess = False
+    required_initial_guess_optional = True
+    accepted_dimensions = 1  # Not sure how to define this for the number of accepted dimensions. Perhaps like the thresholds, at least and at most?
+
+    def __init__(self, bvalues=None, bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3]), initial_guess=None, fitS0=True, thresholds=None, prior_in=None):
+        """
+            Everything this algorithm requires should be implemented here.
+            Number of segmentation thresholds, bounds, etc.
+
+            Our OsipiBase object could contain functions that compare the inputs with
+            the requirements.
+
+            Args:
+                 datain is a 2D array with values of D, f, D* (and S0) that will form the prior.
+        """
+        super(OGC_AmsterdamUMC_Bayesian_biexp, self).__init__(bvalues, bounds,initial_guess,fitS0,prior_in)
+        if bounds is None:
+            self.bounds=([0, 0, 0.005, 0.7],[0.005, 1, 0.2, 1.3])
+        else:
+            self.bounds=bounds
+        if prior_in is None:
+            self.neg_log_prior=flat_neg_log_prior([self.bounds[0][0],self.bounds[1][0]],[self.bounds[0][1],self.bounds[1][1]],[self.bounds[0][2],self.bounds[1][2]],[self.bounds[0][3],self.bounds[1][3]])
+        else:
+            if len(prior_in) is 4:
+                self.neg_log_prior = empirical_neg_log_prior(prior_in[0], prior_in[1], prior_in[2],prior_in[3])
+            else:
+                self.neg_log_prior = empirical_neg_log_prior(prior_in[0], prior_in[1], prior_in[2])
+        self.OGC_algorithm = fit_bayesian
+        if initial_guess is None:
+            self.initial_guess = [0.001, 0.5, 0.1, 1]
+        else:
+            self.initial_guess = initial_guess
+        self.fitS0=fitS0
+
+    def ivim_fit(self, signals, bvalues=None):
+        """Perform the IVIM fit
+
+        Args:
+            signals (array-like)
+            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
+
+        Returns:
+            _type_: _description_
+        """
+        bvalues=np.array(bvalues)
+        fit_results = fit_segmented(bvalues, signals,bounds=self.bounds,cutoff=150,p0=self.initial_guess)
+        fit_results=fit_results+(1,)
+        fit_results = self.OGC_algorithm(bvalues, signals, self.neg_log_prior, x0=fit_results, fitS0=self.fitS0)
+
+        D = fit_results[0]
+        f = fit_results[1]
+        Dstar = fit_results[2]
+
+        return f, Dstar, D