update after review Ivan

- Now use ivim_fit_full_volume instead
- Skipped

Author: oliverchampion

src/standardized/IVIM_NEToptim.py

@@ -60,7 +60,7 @@ def initialize(self, bounds, initial_guess, fitS0, traindata, SNR):
             if SNR is None:
                 warnings.warn('No SNR indicated. Data simulated with SNR = (5-1000)')
                 SNR = (5, 1000)
-            self.training_data(self.bvalues,n=1000000,SNR=SNR)
+            self.osipi_training_data(self.bvalues,n=1000000,SNR=SNR)
         self.arg=Arg()
         if bounds is not None:
             self.arg.net_pars.cons_min = bounds[0]  # Dt, Fp, Ds, S0
@@ -77,7 +77,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
+            bvalues (array-like): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -94,6 +94,29 @@ def ivim_fit(self, signals, bvalues, **kwargs):
 
         return results
 
+
+    def ivim_fit_full_volume(self, signals, bvalues, **kwargs):
+        """Perform the IVIM fit
+
+        Args:
+            signals (array-like)
+            bvalues (array-like): b-values for the signals. If None, self.bvalues will be used. Default is None.
+
+        Returns:
+            _type_: _description_
+        """
+        if not np.array_equal(bvalues, self.bvalues):
+            raise ValueError("bvalue list at fitting must be identical as the one at initiation, otherwise it will not run")
+        signals = self.osipi_reshape_to_voxelwise(signals)
+        paramsNN = deep.predict_IVIM(signals, self.bvalues, self.net, self.arg)
+
+        results = {}
+        results["D"] = paramsNN[0]
+        results["f"] = paramsNN[1]
+        results["Dp"] = paramsNN[2]
+
+        return results
+
 class NetArgs:
     def __init__(self):
         self.optim = 'adam'  # these are the optimisers implementd. Choices are: 'sgd'; 'sgdr'; 'adagrad' adam

src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py

@@ -1,5 +1,5 @@
 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
+from src.original.OGC_AmsterdamUMC.LSQ_fitting import flat_neg_log_prior, fit_bayesian, empirical_neg_log_prior, fit_segmented, fit_bayesian_array, fit_segmented_array
 import numpy as np
 
 class OGC_AmsterdamUMC_Bayesian_biexp(OsipiBase):
@@ -34,7 +34,7 @@ class OGC_AmsterdamUMC_Bayesian_biexp(OsipiBase):
     supported_dimensions = 1
     supported_priors = True
 
-    def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, fitS0=True, prior_in=None):
+    def __init__(self, bvalues=None, thresholds=150, bounds=None, initial_guess=None, fitS0=True, prior_in=None):
 
         """
             Everything this algorithm requires should be implemented here.
@@ -44,14 +44,21 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             the requirements.
 
             Args:
-                 datain is a 2D array with values of D, f, D* (and S0) that will form the prior.
+                 datain (Array): is a 2D array with values of D, f, D* (and S
+                 ) that will form the prior.
+                 thresholds (Bolean, optional): a bolean indicating what threshold is used
+                 prior_in (array, optional): 2D array of D, f, D* and (optionally) S0 values which form the prior
+
         """
         super(OGC_AmsterdamUMC_Bayesian_biexp, self).__init__(bvalues=bvalues, bounds=bounds, thresholds=thresholds, initial_guess=initial_guess) #, fitS0, prior_in)
         self.OGC_algorithm = fit_bayesian
-        self.initialize(bounds, initial_guess, fitS0, prior_in)
+        self.OGC_algorithm_array = fit_bayesian_array
+        self.initialize(bounds, initial_guess, fitS0, prior_in, thresholds)
         self.fit_segmented=fit_segmented
 
-    def initialize(self, bounds=None, initial_guess=None, fitS0=True, prior_in=None):
+    def initialize(self, bounds=None, initial_guess=None, fitS0=True, prior_in=None, thresholds=None):
+
+
         if bounds is None:
             print('warning, no bounds were defined, so default bounds are used of [0, 0, 0.005, 0.7],[0.005, 1.0, 0.2, 1.3]')
             self.bounds=([0, 0, 0.005, 0.7],[0.005, 1.0, 0.2, 1.3])
@@ -64,7 +71,7 @@ def initialize(self, bounds=None, initial_guess=None, fitS0=True, prior_in=None)
             self.initial_guess = initial_guess
         self.use_initial_guess = True
         self.use_bounds = True
-        self.thresholds = 150
+        self.thresholds = thresholds
         if prior_in is None:
             print('using a flat prior between bounds')
             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]])
@@ -101,4 +108,31 @@ def ivim_fit(self, signals, bvalues, initial_guess=None, **kwargs):
         results["f"] = fit_results[1]
         results["Dp"] = fit_results[2]
 
+        return results
+
+    def ivim_fit_full_volume(self, signals, bvalues, initial_guess=None, **kwargs):
+        """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)
+
+        epsilon = 0.000001
+        fit_results = fit_segmented_array(bvalues, signals, bounds=self.bounds, cutoff=self.thresholds, p0=self.initial_guess)
+        fit_results=np.array(fit_results+(1,))
+        for i in range(4):
+            if fit_results[i] < self.bounds[0][i] : fit_results[0] = self.bounds[0][i]+epsilon
+            if fit_results[i] > self.bounds[1][i] : fit_results[0] = self.bounds[1][i]-epsilon
+        fit_results = self.OGC_algorithm_array(bvalues, signals, self.neg_log_prior, x0=fit_results, fitS0=self.fitS0, bounds=self.bounds)
+
+        results = {}
+        results["D"] = fit_results[0]
+        results["f"] = fit_results[1]
+        results["Dp"] = fit_results[2]
+
         return results

src/standardized/OGC_AmsterdamUMC_biexp.py

@@ -1,5 +1,5 @@
 from src.wrappers.OsipiBase import OsipiBase
-from src.original.OGC_AmsterdamUMC.LSQ_fitting import fit_least_squares
+from src.original.OGC_AmsterdamUMC.LSQ_fitting import fit_least_squares, fit_least_squares_array
 import numpy as np
 
 class OGC_AmsterdamUMC_biexp(OsipiBase):
@@ -45,6 +45,7 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
         #super(OGC_AmsterdamUMC_biexp, self).__init__(bvalues, bounds, initial_guess, fitS0)
         super(OGC_AmsterdamUMC_biexp, self).__init__(bvalues=bvalues, bounds=bounds, initial_guess=initial_guess)
         self.OGC_algorithm = fit_least_squares
+        self.OGC_algorithm_array = fit_least_squares_array
         self.fitS0=fitS0
         self.initialize(bounds, initial_guess, fitS0)
 
@@ -83,4 +84,25 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         results["f"] = fit_results[1]
         results["Dp"] = fit_results[2]
 
+        return results
+
+    def ivim_fit_full_volume(self, signals, bvalues, **kwargs):
+        """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 = self.OGC_algorithm_array(bvalues, signals, p0=self.initial_guess, bounds=self.bounds, fitS0=self.fitS0)
+
+        results = {}
+        results["D"] = fit_results[0]
+        results["f"] = fit_results[1]
+        results["Dp"] = fit_results[2]
+
         return results

src/standardized/OGC_AmsterdamUMC_biexp_segmented.py

@@ -1,5 +1,5 @@
 from src.wrappers.OsipiBase import OsipiBase
-from src.original.OGC_AmsterdamUMC.LSQ_fitting import fit_segmented
+from src.original.OGC_AmsterdamUMC.LSQ_fitting import fit_segmented, fit_segmented_array
 import numpy as np
 
 class OGC_AmsterdamUMC_biexp_segmented(OsipiBase):
@@ -44,6 +44,7 @@ def __init__(self, bvalues=None, thresholds=150, bounds=None, initial_guess=None
         """
         super(OGC_AmsterdamUMC_biexp_segmented, self).__init__(bvalues, thresholds, bounds, initial_guess)
         self.OGC_algorithm = fit_segmented
+        self.OGC_algorithm_array = fit_segmented_array
         self.initialize(bounds, initial_guess, thresholds)
 
 
@@ -65,6 +66,7 @@ def initialize(self, bounds, initial_guess, thresholds):
             print('warning, no thresholds were defined, so default bounds are used of  150')
         else:
             self.thresholds = thresholds
+
     def ivim_fit(self, signals, bvalues, **kwargs):
         """Perform the IVIM fit
 
@@ -84,4 +86,27 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         results["f"] = fit_results[1]
         results["Dp"] = fit_results[2]
 
+        return results
+
+
+    def ivim_fit_full_volume(self, signals, bvalues, **kwargs):
+        """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_
+        """
+        signals = self.osipi_reshape_to_voxelwise(signals)
+
+        bvalues=np.array(bvalues)
+        fit_results = self.OGC_algorithm_array(bvalues, signals, bounds=self.bounds, cutoff=self.thresholds, p0=self.initial_guess)
+
+        results = {}
+        results["D"] = fit_results[0]
+        results["f"] = fit_results[1]
+        results["Dp"] = fit_results[2]
+
         return results

src/standardized/Super_IVIM_DC.py

@@ -108,3 +108,30 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         results["Dp"] = Dp
 
         return results
+
+
+    def test_deep_learning_algorithms(self, signals, bvalues, **kwargs):
+        """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:
+            results: a dictionary containing "d", "f", and "Dp".
+        """
+        if not np.array_equal(bvalues, self.bvalues):
+            raise ValueError("bvalue list at fitting must be identical as the one at initiation, otherwise it will not run")
+
+        Dp, Dt, f, S0_superivimdc = infer_from_signal(
+            signal=signals,
+            bvalues=self.bvalues,
+            model_path=f"{self.working_dir}/{self.super_ivim_dc_filename}.pt",
+        )
+
+        results = {}
+        results["D"] = Dt
+        results["f"] = f
+        results["Dp"] = Dp
+
+        return results

src/wrappers/OsipiBase.py

@@ -110,31 +110,26 @@ def osipi_fit(self, data, bvalues=None, **kwargs):
         minimum_bvalue = np.min(use_bvalues) # We normalize the signal to the minimum bvalue. Should be 0 or very close to 0.
         b0_indices = np.where(use_bvalues == minimum_bvalue)[0]
 
-        if not self.deep_learning:
-          for ijk in tqdm(np.ndindex(data.shape[:-1]), total=np.prod(data.shape[:-1])):
-              # Normalize array
-              single_voxel_data = data[ijk]
-              single_voxel_data_normalization_factor = np.mean(single_voxel_data[b0_indices])
-              single_voxel_data_normalized = single_voxel_data/single_voxel_data_normalization_factor
-
-              args = [single_voxel_data_normalized, use_bvalues]
-              fit = self.ivim_fit(*args, **kwargs) # For single voxel fits, we assume this is a dict with a float value per key.
-              for key in list(fit.keys()):
-                  results[key][ijk] = fit[key]
-        else:
-          # Note that I am probably high-jacking the wrong part of the code as DL works best for 2D arrays (b-values vs signal decays).
-          # normalizing of signal still needs implementing for DL. 
-          args = [data, use_bvalues]
-          results = self.ivim_fit(*args,**kwargs)  # For single voxel fits, we assume this is a dict with a float value per key.
+        for ijk in tqdm(np.ndindex(data.shape[:-1]), total=np.prod(data.shape[:-1])):
+            # Normalize array
+            single_voxel_data = data[ijk]
+            single_voxel_data_normalization_factor = np.mean(single_voxel_data[b0_indices])
+            single_voxel_data_normalized = single_voxel_data/single_voxel_data_normalization_factor
+
+            args = [single_voxel_data_normalized, use_bvalues]
+            fit = self.ivim_fit(*args, **kwargs) # For single voxel fits, we assume this is a dict with a float value per key.
+            for key in list(fit.keys()):
+                results[key][ijk] = fit[key]
 
         #self.parameter_estimates = self.ivim_fit(data, bvalues)
         return results
-    
+
+
     def osipi_fit_full_volume(self, data, bvalues=None, **kwargs):
         """Sends a full volume in one go to the fitting algorithm. The osipi_fit method only sends one voxel at a time.
 
         Args:
-            data (array): 3D (single slice) or 4D (multi slice) DWI data.
+            data (array): 2D (data x b-values), 3D (single slice) or 4D (multi slice) DWI data, with last dimension the b-value dimension.
             bvalues (array, optional): The b-values of the DWI data. Defaults to None.
 
         Returns:
@@ -179,7 +174,20 @@ def osipi_fit_full_volume(self, data, bvalues=None, **kwargs):
                 print("Full volume fitting not supported for this algorithm")
 
             return False
-    
+
+
+    def osipi_reshape_to_voxelwise(self, data):
+        """
+        reshapes multi-D input (spatial dims, bvvalue) data to 2D voxel-wise array
+        Args:
+            data (array): mulit-D array (data x b-values)
+        Returns:
+            out (array): 2D array (voxel x b-value)
+        """
+        B = data.shape[-1]
+        voxels = int(np.prod(data.shape[:-1]))  # e.g., X*Y*Z
+        return data.reshape(voxels, B)
+
     def osipi_print_requirements(self):
         """
         Prints the requirements of the algorithm.
@@ -285,11 +293,11 @@ def osipi_check_required_initial_guess(self):
         return True
 
 
-    def osipi_check_required_bvalues():
+    def osipi_check_required_bvalues(self):
         """Minimum number of b-values required"""
         pass
 
-    def osipi_author():
+    def osipi_author(self):
         """Author identification"""
         return ''
 
@@ -323,7 +331,7 @@ def osipi_simple_bias_and_RMSE_test(self, SNR, bvalues, f, Dstar, D, noise_reali
         print(f"Dstar bias:\t{Dstar_bias}\nDstar RMSE:\t{Dstar_RMSE}")
         print(f"D bias:\t{D_bias}\nD RMSE:\t{D_RMSE}")
 
-    def training_data(self, bvalues, data=None, SNR=(5,1000), n=1000000,Drange=(0.0005,0.0034),frange=(0,1),Dprange=(0.005,0.1),rician_noise=False):
+    def osipi_training_data(self, bvalues, data=None, SNR=(5,1000), n=1000000,Drange=(0.0005,0.0034),frange=(0,1),Dprange=(0.005,0.1),rician_noise=False):
         rng = np.random.RandomState(42)
         if data is None:
             gen = GenerateData(rng=rng)

tests/IVIMmodels/unit_tests/algorithms.json

@@ -14,15 +14,11 @@
         "PV_MUMC_biexp",
         "PvH_KB_NKI_IVIMfit",
         "OJ_GU_seg",
-        "IVIM_NEToptim",
-        "Super_IVIM_DC"
+        "IVIM_NEToptim"
     ],
     "IVIM_NEToptim": {
         "deep_learning": true
     },
-    "Super_IVIM_DC": {
-        "deep_learning": true
-    },
     "ASD_MemorialSloanKettering_QAMPER_IVIM": {
         "requires_matlab": true
     },

tests/IVIMmodels/unit_tests/test_ivim_fit.py

@@ -163,7 +163,7 @@ def test_deep_learning_algorithms(deep_learning_algorithms, record_property):
 
     array_2d = np.array([dat["data"] for _, dat in data.items()])
     start_time = time.time()
-    fit_result = fit.osipi_fit(array_2d, bvals)
+    fit_result = fit.osipi_fit_full_volume(array_2d, bvals)
     elapsed_time = time.time() - start_time
 
     errors = []  # Collect all assertion errors