Merge pull request #79 from OSIPI/Unit_Testing_Tolerances

Unit testing tolerances + more consistent definition of SNR

Author: etpeterson

WrapImage/nifti_wrapper.py

@@ -111,7 +111,7 @@ def loop_over_first_n_minus_1_dimensions(arr):
         for idx, view in tqdm(loop_over_first_n_minus_1_dimensions(data), desc=f"{args.algorithm} is fitting", dynamic_ncols=True, total=total_iteration):
             fit_result = fit.osipi_fit(view, bvals)
             f_image.append(fit_result["f"])
-            Dp_image.append(fit_result["D*"])
+            Dp_image.append(fit_result["Dp"])
             D_image.append(fit_result["D"])
 
         # Convert lists to NumPy arrays

doc/Introduction_to_TF24_IVIM-MRI_CodeCollection_github_and_IVIM_Analysis_using_Python.ipynb

@@ -532,7 +532,7 @@
     {
      "data": {
       "text/plain": [
-       "{'f': array(0.04609779), 'D*': array(0.01136011), 'D': array(0.00071134)}"
+       "{'f': array(0.04609779), 'Dp': array(0.01136011), 'D': array(0.00071134)}"
       ]
      },
      "execution_count": 15,
@@ -569,7 +569,7 @@
     {
      "data": {
       "text/plain": [
-       "{'f': array(0.04611801), 'D*': array(0.0113541), 'D': array(0.0007113)}"
+       "{'f': array(0.04611801), 'Dp': array(0.0113541), 'D': array(0.0007113)}"
       ]
      },
      "execution_count": 16,
@@ -636,10 +636,10 @@
     "#plot the results of algorithm 1\n",
     "plt.subplot(121)\n",
     "plt.plot(np.unique(bval),signal_1dir,'x')\n",
-    "plt.plot(np.unique(bval),fit['f']*np.exp(-np.unique(bval)*fit['D*'])+(1-fit['f'])*np.exp(-np.unique(bval)*fit['D']))\n",
-    "plt.plot(np.unique(bval),fit['f']*np.exp(-np.unique(bval)*fit['D*']))\n",
+    "plt.plot(np.unique(bval),fit['f']*np.exp(-np.unique(bval)*fit['Dp'])+(1-fit['f'])*np.exp(-np.unique(bval)*fit['D']))\n",
+    "plt.plot(np.unique(bval),fit['f']*np.exp(-np.unique(bval)*fit['Dp']))\n",
     "plt.plot(np.unique(bval),(1-fit['f'])*np.exp(-np.unique(bval)*fit['D']))\n",
-    "plt.legend(['measured data','model fit','D*','D'])\n",
+    "plt.legend(['measured data','model fit','Dp','D'])\n",
     "plt.ylabel('S/S0')\n",
     "plt.xlabel('b-value [s/mm^2]')\n",
     "plt.title('algorithm 1')\n",
@@ -650,10 +650,10 @@
     "#plot the results of algorithm 2\n",
     "plt.subplot(122)\n",
     "plt.plot(np.unique(bval),signal_1dir,'x')\n",
-    "plt.plot(np.unique(bval),fit['f']*np.exp(-np.unique(bval)*fit['D*'])+(1-fit['f'])*np.exp(-np.unique(bval)*fit['D']))\n",
-    "plt.plot(np.unique(bval),fit['f']*np.exp(-np.unique(bval)*fit['D*']))\n",
+    "plt.plot(np.unique(bval),fit['f']*np.exp(-np.unique(bval)*fit['Dp'])+(1-fit['f'])*np.exp(-np.unique(bval)*fit['D']))\n",
+    "plt.plot(np.unique(bval),fit['f']*np.exp(-np.unique(bval)*fit['Dp']))\n",
     "plt.plot(np.unique(bval),(1-fit['f'])*np.exp(-np.unique(bval)*fit['D']))\n",
-    "plt.legend(['measured data','model fit','D*','D'])\n",
+    "plt.legend(['measured data','model fit','Dp','D'])\n",
     "plt.ylabel('S/S0')\n",
     "plt.xlabel('b-value [s/mm^2]')\n",
     "plt.title('algorithm 2')\n"
@@ -818,7 +818,7 @@
      "data": {
       "text/plain": [
        "{'f': array([0., 0., 0., ..., 0., 0., 0.]),\n",
-       " 'D*': array([0., 0., 0., ..., 0., 0., 0.]),\n",
+       " 'Dp': array([0., 0., 0., ..., 0., 0., 0.]),\n",
        " 'D': array([0., 0., 0., ..., 0., 0., 0.])}"
       ]
      },

phantoms/MR_XCAT_qMRI/sim_ivim_sig.py

@@ -459,9 +459,9 @@ def parse_bvalues_file(file_path):
             signals = np.squeeze(voxels[int(voxels.shape[0] * voxel_selector_fraction)]).tolist()
             generic_data[name] = {
                 'noise': noise,
-                'D': np.mean(Dim[selector], axis=0),
-                'f': np.mean(fim[selector], axis=0),
-                'Dp': np.mean(Dpim[selector], axis=0),
+                'D': np.median(Dim[selector], axis=0),
+                'f': np.median(fim[selector], axis=0),
+                'Dp': np.median(Dpim[selector], axis=0),
                 'data': signals
             }
         generic_data['config'] = {

src/original/OGC_AmsterdamUMC/LSQ_fitting.py

@@ -124,29 +124,32 @@ 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
+        dw_data=dw_data/np.mean(dw_data[bvalues==0])
         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., 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=(p0[0], p0[3]-p0[1]/10),
+        assert bounds[0][1] < p0[1]
+        assert bounds[1][1] > p0[1]
+        bounds1 = ([bounds[0][0], 0], [bounds[1][0], 10000000000])
+        params, _ = curve_fit(lambda b, Dt, int: int * np.exp(-b * Dt ), high_b, high_dw_data,
+                              p0=(p0[0], p0[3]-p0[1]),
                               bounds=bounds1)
-        Dt, Fp = params[0] / 1000, 1 - params[1]
+        Dt, Fp = params[0], 1 - params[1]
+        if Fp < bounds[0][1] : Fp = np.float64(bounds[0][1])
+        if Fp > bounds[1][1] : Fp = np.float64(bounds[1][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]*10, bounds[1][2]*10)
+        bounds2 = (bounds[0][2], bounds[1][2])
         # fit for D*
         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:
         # if fit fails, return zeros
-        # print('segnetned fit failed')
+        # print('segmented fit failed')
         return 0., 0., 0.
 
 
@@ -235,17 +238,17 @@ def fit_least_squares(bvalues, dw_data, S0_output=False, fitS0=True,
     try:
         if not fitS0:
             # 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],
+            bounds2 = ([bounds[0][0] * 1000, bounds[0][1] * 10, bounds[0][2] * 10],
                       [bounds[1][0] * 1000, bounds[1][1] * 10, bounds[1][2] * 10])
-            p0=[p0[0]*1000,p0[1]*10,p0[2]*10]
-            params, _ = curve_fit(ivimN_noS0, bvalues, dw_data, p0=p0, bounds=bounds)
+            p1=[p0[0]*1000,p0[1]*10,p0[2]*10]
+            params, _ = curve_fit(ivimN_noS0, bvalues, dw_data, p0=p1, bounds=bounds2)
             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]],
+            bounds2 = ([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]])
-            p0=[p0[0]*1000,p0[1]*10,p0[2]*10,p0[3]]
-            params, _ = curve_fit(ivimN, bvalues, dw_data, p0=p0, bounds=bounds)
+            p1=[p0[0]*1000,p0[1]*10,p0[2]*10,p0[3]]
+            params, _ = curve_fit(ivimN, bvalues, dw_data, p0=p1, bounds=bounds2)
             S0 = params[3]
         # correct for the rescaling of parameters
         Dt, Fp, Dp = params[0] / 1000, params[1] / 10, params[2] / 10
@@ -258,10 +261,10 @@ def fit_least_squares(bvalues, dw_data, S0_output=False, fitS0=True,
         # if fit fails, then do a segmented fit instead
         # print('lsq fit failed, trying segmented')
         if S0_output:
-            Dt, Fp, Dp = fit_segmented(bvalues, dw_data, bounds=bounds)
+            Dt, Fp, Dp = fit_segmented(bvalues, dw_data, bounds=bounds,p0=p0)
             return Dt, Fp, Dp, 1
         else:
-            return fit_segmented(bvalues, dw_data)
+            return fit_segmented(bvalues, dw_data, bounds=bounds,p0=p0)
 
 
 def fit_least_squares_array_tri_exp(bvalues, dw_data, S0_output=True, fitS0=True, njobs=4,
@@ -561,19 +564,19 @@ def neg_log_prior(p):
             Dt, Fp, Dp = p[0], p[1], p[2]
         # make D*<D very unlikely
         if (Dp < Dt):
-            return 1e3
+            return 1e10
         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
+                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
+                    return 1e10
             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 1e10
 
     return neg_log_prior
 
@@ -638,7 +641,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, 1], fitS0=True):
+def fit_bayesian(bvalues, dw_data, neg_log_prior, x0=[0.001, 0.2, 0.05, 1], fitS0=True, bounds=([0,0,0,0],[0.005,1.5,2,2.5])):
     '''
     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
@@ -655,7 +658,7 @@ def fit_bayesian(bvalues, dw_data, neg_log_prior, x0=[0.001, 0.2, 0.05, 1], fitS
     '''
     try:
         # define fit bounds
-        bounds = [(0, 0.005), (0, 1.5), (0, 2), (0, 2.5)]
+        bounds = [(bounds[0][0], bounds[1][0]), (bounds[0][1], bounds[1][1]), (bounds[0][2], bounds[1][2]), (bounds[0][3], bounds[1][3])]
         # 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)

src/standardized/ETP_SRI_LinearFitting.py

@@ -41,8 +41,13 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             Our OsipiBase object could contain functions that compare the inputs with
             the requirements.
         """
-        super(ETP_SRI_LinearFitting, self).__init__(bvalues=bvalues, thresholds=thresholds, bounds=bounds, initial_guess=initial_guess)
-        
+
+        super(ETP_SRI_LinearFitting, self).__init__(bvalues, thresholds, bounds, initial_guess)
+        if bounds is not None:
+            print('warning, bounds from wrapper are not (yet) used in this algorithm')
+        self.use_bounds = False
+        self.use_initial_guess = False
+
         # Could be a good idea to have all the submission-specfic variable be 
         # defined with initials?
         self.ETP_weighting = weighting
@@ -62,6 +67,7 @@ def ivim_fit(self, signals, bvalues=None, linear_fit_option=False, **kwargs):
         Returns:
             _type_: _description_
         """
+        signals[signals<0.0000001]=0.0000001
         if bvalues is None:
             bvalues = self.bvalues
 
@@ -75,14 +81,14 @@ def ivim_fit(self, signals, bvalues=None, linear_fit_option=False, **kwargs):
             f, Dstar = ETP_object.linear_fit(bvalues, signals, self.ETP_weighting, self.ETP_stats)
 
             results["f"] = f
-            results["D*"] = Dstar
+            results["Dp"] = Dstar
 
             return results
         else: 
             f, D, Dstar = ETP_object.ivim_fit(bvalues, signals)
 
             results["f"] = f
-            results["D*"] = Dstar
+            results["Dp"] = Dstar
             results["D"] = D
 
             return results

src/standardized/IAR_LU_biexp.py

@@ -42,7 +42,10 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             the requirements.
         """
         super(IAR_LU_biexp, self).__init__(bvalues, thresholds, bounds, initial_guess)
-        
+        if bounds is not None:
+            print('warning, bounds from wrapper are not (yet) used in this algorithm')
+        self.use_bounds = False
+        self.use_initial_guess = False
         # Check the inputs
 
         # Initialize the algorithm
@@ -83,7 +86,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
 
         results = {}
         results["f"] = fit_results.model_params[1]
-        results["D*"] = fit_results.model_params[2]
+        results["Dp"] = fit_results.model_params[2]
         results["D"] = fit_results.model_params[3]
 
         return results
@@ -115,7 +118,7 @@ def ivim_fit_full_volume(self, signals, bvalues, **kwargs):
 
         results = {}
         results["f"] = fit_results.model_params[..., 1]
-        results["D*"] = fit_results.model_params[..., 2]
+        results["Dp"] = fit_results.model_params[..., 2]
         results["D"] = fit_results.model_params[..., 3]
 
         return results

src/standardized/IAR_LU_modified_mix.py

@@ -42,7 +42,10 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             the requirements.
         """
         super(IAR_LU_modified_mix, self).__init__(bvalues, thresholds, bounds, initial_guess)
-        
+        if bounds is not None:
+            print('warning, bounds from wrapper are not (yet) used in this algorithm')
+        self.use_bounds = False
+        self.use_initial_guess = False
         # Check the inputs
 
         # Initialize the algorithm
@@ -86,7 +89,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         #D = fit_results.model_params[3]
         results = {}
         results["f"] = fit_results.model_params[1]
-        results["D*"] = fit_results.model_params[2]
+        results["Dp"] = fit_results.model_params[2]
         results["D"] = fit_results.model_params[3]
 
         return results

src/standardized/IAR_LU_modified_topopro.py

@@ -42,7 +42,10 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             the requirements.
         """
         super(IAR_LU_modified_topopro, self).__init__(bvalues, thresholds, bounds, initial_guess)
-        
+        if bounds is not None:
+            print('warning, bounds from wrapper are not (yet) used in this algorithm')
+        self.use_bounds = False
+        self.use_initial_guess = False
         # Check the inputs
 
         # Initialize the algorithm
@@ -88,7 +91,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         #return f, Dstar, D
         results = {}
         results["f"] = fit_results.model_params[1]
-        results["D*"] = fit_results.model_params[2]
+        results["Dp"] = fit_results.model_params[2]
         results["D"] = fit_results.model_params[3]
 
         return results

src/standardized/IAR_LU_segmented_2step.py

@@ -42,7 +42,10 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             the requirements.
         """
         super(IAR_LU_segmented_2step, self).__init__(bvalues, thresholds, bounds, initial_guess)
-        
+        if bounds is not None:
+            print('warning, bounds from wrapper are not (yet) used in this algorithm')
+        self.use_bounds = False
+        self.use_initial_guess = False
         # Check the inputs
 
         # Initialize the algorithm
@@ -92,7 +95,7 @@ def ivim_fit(self, signals, bvalues, thresholds=None, **kwargs):
         #return f, Dstar, D
         results = {}
         results["f"] = fit_results.model_params[1]
-        results["D*"] = fit_results.model_params[2]
+        results["Dp"] = fit_results.model_params[2]
         results["D"] = fit_results.model_params[3]
 
         return results

src/standardized/IAR_LU_segmented_3step.py

@@ -42,7 +42,10 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             the requirements.
         """
         super(IAR_LU_segmented_3step, self).__init__(bvalues, thresholds, bounds, initial_guess)
-        
+        if bounds is not None:
+            print('warning, bounds from wrapper are not (yet) used in this algorithm')
+        self.use_bounds = False
+        self.use_initial_guess = False
         # Check the inputs
 
         # Initialize the algorithm
@@ -88,7 +91,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         #return f, Dstar, D
         results = {}
         results["f"] = fit_results.model_params[1]
-        results["D*"] = fit_results.model_params[2]
+        results["Dp"] = fit_results.model_params[2]
         results["D"] = fit_results.model_params[3]
 
         return results

src/standardized/IAR_LU_subtracted.py

@@ -42,7 +42,10 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             the requirements.
         """
         super(IAR_LU_subtracted, self).__init__(bvalues, thresholds, bounds, initial_guess)
-        
+        if bounds is not None:
+            print('warning, bounds from wrapper are not (yet) used in this algorithm')
+        self.use_bounds = False
+        self.use_initial_guess = False
         # Check the inputs
 
         # Initialize the algorithm
@@ -88,7 +91,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         #return f, Dstar, D
         results = {}
         results["f"] = fit_results.model_params[1]
-        results["D*"] = fit_results.model_params[2]
+        results["Dp"] = fit_results.model_params[2]
         results["D"] = fit_results.model_params[3]
 
         return results