Did some work on PV_MUMC. It now runs but the estimation errors are large. Needs investigation

Author: IvanARashid

src/original/PV_MUMC/two_step_IVIM_fit.py

@@ -75,45 +75,48 @@ def fit_least_squares(bvalues, dw_data, IR=False, S0_output=False, fitS0=True,
     :return Dmv: scalar with Dmv of the specific voxel
     """
 
-    try:
-        def monofit(bvalues, Dpar):
-             return np.exp(-bvalues * Dpar)
-        
-        high_b = bvalues[bvalues >= cutoff]
-        high_dw_data = dw_data[bvalues >= cutoff]
-        boundspar = ([bounds[0][1]], [bounds[1][1]])
-        params, _ = curve_fit(monofit, high_b, high_dw_data, p0=[(bounds[1][1]-bounds[0][1])/2], bounds=boundspar)
-        Dpar1 = params[0]
+    #try:
+    def monofit(bvalues, Dpar):
+            return np.exp(-bvalues * Dpar)
+    
+    high_b = bvalues[bvalues >= cutoff]
+    high_dw_data = dw_data[bvalues >= cutoff]
+    boundspar = ([bounds[0][1]], [bounds[1][1]])
+    params, _ = curve_fit(monofit, high_b, high_dw_data, p0=[(bounds[1][1]-bounds[0][1])/2], bounds=boundspar)
+    Dpar = params[0]
 
-        if not fitS0:
-            boundsupdated=([Dpar1 , bounds[0][2] , bounds[0][3] ],
-                      [Dpar1 , bounds[1][2] , bounds[1][3] ])    
-            params, _ = curve_fit(two_exp_noS0, bvalues, dw_data, p0=[Dpar1, (bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated)
-            Dpar, Fmv, Dmv = params[0], 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")
+    if not fitS0:
+        boundsupdated=([Dpar1 , bounds[0][2] , bounds[0][3] ],
+                    [Dpar1 , bounds[1][2] , bounds[1][3] ])    
+        params, _ = curve_fit(two_exp_noS0, bvalues, dw_data, p0=[Dpar1, (bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated)
+        Dpar1, Fmv, Dmv = params[0], 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")
+        
+    else: 
+        #boundsupdated = ([bounds[0][0] , Dpar1 , bounds[0][2] , bounds[0][3] ],
+        #            [bounds[1][0] , Dpar1, bounds[1][2] , bounds[1][3] ])   
+        #params, _ = curve_fit(two_exp, bvalues, dw_data, p0=[1, Dpar1, (bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated)
+        boundsupdated = ([bounds[0][0] , bounds[0][2] , bounds[0][3] ],
+                    [bounds[1][0] , bounds[1][2] , bounds[1][3] ])   
+        params, _ = curve_fit(lambda bvalues, S0, Fmv, Dmv: two_exp(bvalues, S0, Dpar, 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")     
 
-        else: 
-            boundsupdated = ([bounds[0][0] , Dpar1 , bounds[0][2] , bounds[0][3] ],
-                      [bounds[1][0] , Dpar1, bounds[1][2] , bounds[1][3] ])   
-            params, _ = curve_fit(two_exp, bvalues, dw_data, p0=[1, Dpar1, (bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated)
-            S0 = params[0]
-            Dpar, Fmv, Dmv = params[1] , params[2] , params[3]
-            #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")     
-                
-        if S0_output:
-            return Dpar, Fmv, Dmv, S0
-        else:
-            return Dpar, Fmv, Dmv
-    except:
+    if S0_output:
+        return Dpar, Fmv, Dmv, S0
+    else:
+        return Dpar, Fmv, Dmv
+    #except:
 
-        if S0_output:
-            return 0, 0, 0, 0, 0, 0
-        else:
-            return 0, 0, 0, 0, 0
+    if S0_output:
+        return 0, 0, 0, 0, 0, 0
+    else:
+        return 0, 0, 0, 0, 0
 
 
 

src/standardized/PV_MUMC_biexp.py

@@ -0,0 +1,56 @@
+import numpy as np
+from src.wrappers.OsipiBase import OsipiBase
+from src.original.PV_MUMC.two_step_IVIM_fit import fit_least_squares_array, fit_least_squares
+
+
+class PV_MUMC_biexp(OsipiBase):
+    """
+    Bi-exponential fitting algorithm by Paulien Voorter, Maastricht University
+    """
+    
+    # Some basic stuff that identifies the algorithm
+    id_author = "Paulien Voorter MUMC"
+    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,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, thresholds=None, bounds=None, initial_guess=None, weighting=None, stats=False):
+        """
+            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(PV_MUMC_biexp, self).__init__(bvalues=bvalues, thresholds=None, bounds=bounds, initial_guess=None)
+        self.PV_algorithm = fit_least_squares
+        
+    
+    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_
+        """
+        
+            
+        fit_results = self.PV_algorithm(bvalues, signals)
+        
+        f = fit_results[1]
+        Dstar = fit_results[2]
+        D = fit_results[0]
+        
+        return f, Dstar, D

tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py

@@ -2,14 +2,15 @@
 import os
 from pathlib import Path
 #from src.standardized.ETP_SRI_LinearFitting import ETP_SRI_LinearFitting
-from src.standardized.IAR_LU_biexp import IAR_LU_biexp
+#from src.standardized.IAR_LU_biexp import IAR_LU_biexp
+from src.standardized.PV_MUMC_biexp import PV_MUMC_biexp
 
 ## Simple test code... 
 # Used to just do a test run of an algorithm during development
 def dev_test_run(model, **kwargs):
-    bvalues = np.array([0, 200, 500, 800])
+    bvalues = np.array([0, 50, 100, 150, 200, 500, 800])
 
-    def ivim_model(b, S0=1, f=0.1, Dstar=0.03, D=0.001):
+    def ivim_model(b, S0=1, f=0.1, Dstar=0.01, D=0.001):
         return S0*(f*np.exp(-b*Dstar) + (1-f)*np.exp(-b*D))
 
     signals = ivim_model(bvalues)
@@ -23,7 +24,8 @@ def ivim_model(b, S0=1, f=0.1, Dstar=0.03, D=0.001):
     #test = model.osipi_simple_bias_and_RMSE_test(SNR=20, bvalues=bvalues, f=0.1, Dstar=0.03, D=0.001, noise_realizations=10)
 
 #model1 = ETP_SRI_LinearFitting(thresholds=[200])
-model2 = IAR_LU_biexp()
+#model2 = IAR_LU_biexp()
+model2 = PV_MUMC_biexp()
 
 #dev_test_run(model1, linear_fit_option=True)
 dev_test_run(model2)