Merge pull request #101 from OSIPI/reference_fitting_algorithm

Added a segmented-fitting algorithm based on the recommendation paper

Author: oliverchampion

phantoms/MR_XCAT_qMRI/sim_ivim_sig.py

@@ -429,7 +429,7 @@ def parse_bvalues_file(file_path):
     elif args.bvalue:
         bvalues = {"cmd": args.bvalue}
     else:
-        bvalues = parse_bvalues_file("b_values.json") 
+        bvalues = parse_bvalues_file("b_values.json")
 
 
     noise = args.noise

requirements.txt

@@ -16,4 +16,5 @@ pandas
 sphinx
 sphinx_rtd_theme
 pytest-json-report
+statsmodels
 ivimnet

src/original/TF_reference/segmented_IVIMfit.py

@@ -0,0 +1,78 @@
+"""
+Code to calculate IVIM maps
+
+Segmented fitting approach
+First fit D using a WLLS
+
+"""
+
+import numpy as np
+import statsmodels.api as sm
+import scipy
+
+def ivim_biexp(bvalues, D, f, Dp, S0=1):
+    return (S0 * (f * np.exp(-bvalues * Dp) + (1 - f) * np.exp(-bvalues * D)))
+
+
+def segmented_IVIM_fit(bvalues, dw_data, b_cutoff = 200, bounds=([0.0001, 0.0, 0.001], [0.004, 0.7, 0.01])):
+    """
+        A segmented fitting implementation for a bi-exponential model.
+        First D is fitted using a mono exponential model on all signal above the bvalue cutoff using an iterative WLLSVertrekpassage, Schiphol
+        Then f is fitted by using the b=0 intercept from the mono expontential fit and substracting this from the measured signal at b=0
+        Then D* is fitted using a bi-exponential model with fixed D and f
+
+        """
+    bvalues_D = bvalues[bvalues >= b_cutoff]
+    dw_data_D = dw_data[bvalues >= b_cutoff]
+    log_data_D = np.log(dw_data_D)
+
+    D, b0_intercept = d_fit_iterative_wls(bvalues_D, log_data_D)
+
+    D = np.clip(D, bounds[0][0], bounds[1][0])
+    S0 = dw_data[bvalues == 0].item()
+    f = (S0 - b0_intercept) / S0
+    f = np.clip(f, bounds[0][1], bounds[1][1])
+
+    Dp = scipy.optimize.least_squares(
+    fun=lambda Dp: ivim_biexp(bvalues, D, f, Dp[0]) - dw_data,
+    x0=np.clip(D * 10, bounds[0][2], bounds[1][2]), # Initial guess for D*
+    bounds=(bounds[0][2], bounds[1][2]))
+
+    Dp = Dp.x[0]
+
+    return D, f, Dp
+
+
+def d_fit_iterative_wls(bvalues_D, log_signal, max_iter=50):
+    """
+    Function to calculate D using an iterative wlls on the log(signal)
+    weights for the wlls are initialized from the predicted signal of a lls as described in
+    http://dx.doi.org/10.1016/j.neuroimage.2013.05.028 equation (7)
+
+    Parameters:
+    log_signal: the log() of the signal above the threshold for segmented fitting
+
+    bvalues_D: all bvalues above the threshold for fitting D
+
+    max_iter: the maximum number of iterations for WLS
+
+    """
+
+    bvals = sm.add_constant(-bvalues_D)
+    # First do a LLS to initialize the weights
+    beta_lls = np.linalg.lstsq(bvals, log_signal, rcond=None)[0]
+    # initialize the weights based on the predicted signals W=diag(exp(2*X*Beta_lls))
+    init_weights = np.exp(2*bvals@beta_lls)
+    weights = init_weights
+
+    for i in range(max_iter):
+        # Weighted linear regression
+        model = sm.WLS(log_signal, bvals, weights=weights)
+        results = model.fit()
+        # The weights for the next iteration are based on the predicted signal of this iteration
+        new_weights = np.exp(2*bvals@results.params)
+        weights = new_weights
+
+    ln_b0_intercept, D  = results.params
+    b0_intercept = np.exp(ln_b0_intercept)
+    return D, b0_intercept

src/standardized/TF_reference_IVIMfit.py

@@ -0,0 +1,88 @@
+from src.wrappers.OsipiBase import OsipiBase
+from src.original.TF_reference.segmented_IVIMfit import segmented_IVIM_fit
+import numpy as np
+
+class TF_reference_IVIMfit(OsipiBase):
+    """
+    Bi-exponential fitting algorithm by IVIM Task force
+    """
+
+    # 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 = "OSIPI IVIM TF"
+    id_algorithm_type = "Bi-exponential fit"
+    id_return_parameters = "f, D*, D"
+    id_units = "seconds per milli metre squared or milliseconds per micro metre squared"
+
+    # Algorithm requirements
+    required_bvalues = 4
+    required_thresholds = [0,
+                           1]  # 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 = False
+    accepted_dimensions = 1  # Not sure how to define this for the number of accepted dimensions. Perhaps like the thresholds, at least and at most?
+
+    # Supported inputs in the standardized class
+    supported_bounds = True
+    supported_initial_guess = False
+    supported_thresholds = True
+
+    def __init__(self, bvalues=None, thresholds=200, bounds=None, initial_guess=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.
+        """
+        super(TF_reference_IVIMfit, self).__init__(bvalues=bvalues, thresholds=thresholds,bounds=bounds,initial_guess=initial_guess)
+        self.TF_reference_algorithm = segmented_IVIM_fit
+        self.initialize(bounds, thresholds)
+        self.use_initial_guess = False
+
+
+    def initialize(self, bounds, thresholds):
+        if bounds is None:
+            print('warning, no bounds were defined, so default bounds are used of [0, 0, 0.005],[0.005, 1.0, 0.2]')
+            self.bounds=([0, 0, 0.005, 0.8],[0.005, 1.0, 0.2, 1.2])
+        else:
+            self.bounds=bounds
+        self.use_bounds = True
+        if thresholds is None:
+            self.thresholds = 200
+            print('warning, no thresholds were defined, so default threshold are used of  200')
+        else:
+            self.thresholds = thresholds
+
+
+    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)
+        bvalues = np.array(bvalues)
+
+        if np.any(signals < 0):
+            signals = np.clip(signals,0.01, None)
+            print('warning, negative values in signal: values clipped to 0.01')
+
+
+        fit_results = self.TF_reference_algorithm(bvalues,signals,b_cutoff=self.thresholds, bounds=self.bounds)
+
+        results = {}
+        results["D"] = fit_results[0]
+        results["f"] = fit_results[1]
+        results["Dp"] = fit_results[2]
+
+        return results

tests/IVIMmodels/unit_tests/algorithms.json

@@ -14,6 +14,7 @@
         "PV_MUMC_biexp",
         "PvH_KB_NKI_IVIMfit",
         "OJ_GU_seg",
+        "TF_reference_IVIMfit",
         "IVIM_NEToptim"
     ],
     "IVIM_NEToptim": {

tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py

@@ -12,6 +12,7 @@
 #from src.standardized.PvH_KB_NKI_IVIMfit import PvH_KB_NKI_IVIMfit
 #from src.standardized.PV_MUMC_biexp import PV_MUMC_biexp
 from src.standardized.OGC_AmsterdamUMC_biexp import OGC_AmsterdamUMC_biexp
+from src.standardized.TF_reference_IVIMfit import TF_reference_IVIMfit
 
 ## Simple test code... 
 # Used to just do a test run of an algorithm during development