Merge pull request #106 from OSIPI/wrapper_dev

Wrapper dev - attribute harmonization and signal normalization

Author: oliverchampion

doc/Introduction_to_TF24_IVIM-MRI_CodeCollection_github_and_IVIM_Analysis_using_Python.ipynb

@@ -0,0 +1 @@
+0 50 100 500 1000

ivim_simulation.bval

@@ -0,0 +1,3 @@
+1 0 0
+0 1 0
+0 0 1

ivim_simulation.bvec

@@ -25,13 +25,14 @@ class ETP_SRI_LinearFitting(OsipiBase):
     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 = False
     supported_initial_guess = False
     supported_thresholds = True
+    supported_dimensions = 1
+    supported_priors = False
 
     def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, weighting=None, stats=False):
         """

src/standardized/ETP_SRI_LinearFitting.py

@@ -26,12 +26,13 @@ class IAR_LU_biexp(OsipiBase):
     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?
 
     # Supported inputs in the standardized class
     supported_bounds = True
     supported_initial_guess = True
     supported_thresholds = False
+    supported_dimensions = 1
+    supported_priors = False
 
     def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, weighting=None, stats=False):
         """

src/standardized/IAR_LU_biexp.py

@@ -26,12 +26,13 @@ class IAR_LU_modified_mix(OsipiBase):
     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?
 
     # Supported inputs in the standardized class
     supported_bounds = True
     supported_initial_guess = False
     supported_thresholds = False
+    supported_dimensions = 1
+    supported_priors = False
 
     def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, weighting=None, stats=False):
         """

src/standardized/IAR_LU_modified_mix.py

@@ -26,12 +26,13 @@ class IAR_LU_modified_topopro(OsipiBase):
     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?
 
     # Supported inputs in the standardized class
     supported_bounds = True
     supported_initial_guess = False
     supported_thresholds = False
+    supported_dimensions = 1
+    supported_priors = False
 
     def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, weighting=None, stats=False):
         """

src/standardized/IAR_LU_modified_topopro.py

@@ -26,12 +26,13 @@ class IAR_LU_segmented_2step(OsipiBase):
     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?
 
     # Supported inputs in the standardized class
     supported_bounds = True
     supported_initial_guess = True
     supported_thresholds = True
+    supported_dimensions = 1
+    supported_priors = False
 
     def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None):
         """

src/standardized/IAR_LU_segmented_2step.py

@@ -26,12 +26,13 @@ class IAR_LU_segmented_3step(OsipiBase):
     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?
 
     # Supported inputs in the standardized class
     supported_bounds = True
     supported_initial_guess = True
     supported_thresholds = False
+    supported_dimensions = 1
+    supported_priors = False
 
     def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, weighting=None, stats=False):
         """

src/standardized/IAR_LU_segmented_3step.py

@@ -26,12 +26,13 @@ class IAR_LU_subtracted(OsipiBase):
     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?
 
     # Supported inputs in the standardized class
     supported_bounds = True
     supported_initial_guess = True
     supported_thresholds = False
+    supported_dimensions = 1
+    supported_priors = False
 
     def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None):
         """

src/standardized/IAR_LU_subtracted.py

@@ -25,14 +25,14 @@ class OGC_AmsterdamUMC_Bayesian_biexp(OsipiBase):
     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?
-    accepts_priors = True
 
 
     # Supported inputs in the standardized class
     supported_bounds = True
     supported_initial_guess = True
     supported_thresholds = True
+    supported_dimensions = 1
+    supported_priors = True
 
     def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, fitS0=True, prior_in=None):
 

src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py

@@ -25,13 +25,14 @@ class OGC_AmsterdamUMC_biexp(OsipiBase):
     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?
 
 
     # Supported inputs in the standardized class
     supported_bounds = True
     supported_initial_guess = True
     supported_thresholds = False
+    supported_dimensions = 1
+    supported_priors = False
 
     def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, fitS0=True):
         """

src/standardized/OGC_AmsterdamUMC_biexp.py

@@ -25,13 +25,14 @@ class OGC_AmsterdamUMC_biexp_segmented(OsipiBase):
     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?
 
 
     # Supported inputs in the standardized class
     supported_bounds = True
     supported_initial_guess = True
     supported_thresholds = True
+    supported_dimensions = 1
+    supported_priors = False
 
     def __init__(self, bvalues=None, thresholds=150, bounds=None, initial_guess=None):
         """

src/standardized/OGC_AmsterdamUMC_biexp_segmented.py

@@ -24,12 +24,13 @@ class OJ_GU_seg(OsipiBase):
     required_bounds_optional = False # 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 = False
     supported_initial_guess = False
     supported_thresholds = True
+    supported_dimensions = 1
+    supported_priors = False
 
     def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, weighting=None, stats=False):
         """

src/standardized/OJ_GU_seg.py

@@ -21,12 +21,13 @@ class PV_MUMC_biexp(OsipiBase):
     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?
 
     # Supported inputs in the standardized class
     supported_bounds = True
     supported_initial_guess = False
     supported_thresholds = True
+    supported_dimensions = 1
+    supported_priors = False
 
     def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, weighting=None, stats=False):
         """

src/standardized/PV_MUMC_biexp.py

@@ -25,12 +25,13 @@ class PvH_KB_NKI_IVIMfit(OsipiBase):
     required_bounds_optional = False  # 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 = False
     supported_initial_guess = False
     supported_thresholds = False
+    supported_dimensions = 1
+    supported_priors = False
 
     def __init__(self, bvalues=None, thresholds=None,bounds=None,initial_guess=None):
         """

src/standardized/PvH_KB_NKI_IVIMfit.py

@@ -103,9 +103,16 @@ def osipi_fit(self, data, bvalues=None, **kwargs):
             #args = [data[ijk], use_bvalues]
             #fit = list(self.ivim_fit(*args, **kwargs))
             #results[ijk] = fit
+        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]
 
         for ijk in tqdm(np.ndindex(data.shape[:-1]), total=np.prod(data.shape[:-1])):
-            args = [data[ijk], use_bvalues]
+            # 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]
@@ -141,7 +148,15 @@ def osipi_fit_full_volume(self, data, bvalues=None, **kwargs):
             for key in self.result_keys:
                 results[key] = np.empty(list(data.shape[:-1]))
 
-            args = [data, use_bvalues]
+            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]
+            b0_mean = np.mean(data[..., b0_indices], axis=-1)
+
+            normalization_factors = np.array([b0_mean for i in range(data.shape[-1])])
+            normalization_factors = np.moveaxis(normalization_factors, 0, -1)
+            data_normalized = data/normalization_factors
+
+            args = [data_normalized, use_bvalues]
             fit = self.ivim_fit_full_volume(*args, **kwargs) # Assume this is a dict with an array per key representing the parametric maps
             for key in list(fit.keys()):
                 results[key] = fit[key]

src/wrappers/OsipiBase.py

@@ -22,22 +22,22 @@ 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)
-    data = np.array([signals, signals, signals])
+    data = np.array([[signals, signals, signals], [signals, signals, signals]])
     #print(data)
     signals = data
 
     #model = ETP_SRI_LinearFitting(thresholds=[200])
     if kwargs:
-        results = model.osipi_fit(signals, bvalues, **kwargs)
+        results = model.osipi_fit_full_volume(signals, bvalues, **kwargs)
     else:
-        results = model.osipi_fit(signals, bvalues)
+        results = model.osipi_fit_full_volume(signals, bvalues)
     print(results)
     #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(bounds=([0,0,0,0], [1,1,1,1]))
+model2 = IAR_LU_biexp(bounds=([0,0,0,0], [1,1,1,1]))
 #model2 = IAR_LU_modified_mix()
-model2 = OGC_AmsterdamUMC_biexp()
+#model2 = OGC_AmsterdamUMC_biexp()
 
 #dev_test_run(model1)
 dev_test_run(model2)

tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py

@@ -25,7 +25,7 @@ def test_generated(algorithmlist, ivim_data, SNR, rtol, atol, fit_count, rician_
     Dp = data["Dp"]
     fit = OsipiBase(algorithm=ivim_algorithm)
     # here is a prior
-    if use_prior and hasattr(fit, "accepts_priors") and fit.accepts_priors:
+    if use_prior and hasattr(fit, "supported_priors") and fit.supported_priors:
         prior = [rng.normal(D, D/3, 10), rng.normal(f, f/3, 10), rng.normal(Dp, Dp/3, 10), rng.normal(1, 1/3, 10)]
         # prior = [np.repeat(D, 10)+np.random.normal(0,D/3,np.shape(np.repeat(D, 10))), np.repeat(f, 10)+np.random.normal(0,f/3,np.shape(np.repeat(D, 10))), np.repeat(Dp, 10)+np.random.normal(0,Dp/3,np.shape(np.repeat(D, 10))),np.repeat(np.ones_like(Dp), 10)+np.random.normal(0,1/3,np.shape(np.repeat(D, 10)))]  # D, f, D*
         fit.initialize(prior_in=prior)