Parallel

.github/workflows/unit_test.yml

@@ -33,4 +33,5 @@ jobs:
       - name: Test with pytest
         run: |
           pip install pytest pytest-cov
-          python -m pytest --doctest-modules --junitxml=junit/test-results.xml --cov=. --cov-report=xml --cov-report=html
+          python -m pytest --doctest-modules --junitxml=junit/test-results.xml --cov=. --cov-report=xml --cov-report=html -r w
+

.gitignore

@@ -5,6 +5,7 @@
 __pycache__/
 *.nii.gz
 *.nii
+*.npy
 *.dcm
 *.npy
 *.mat

conftest.py

@@ -3,7 +3,8 @@
 import json
 import csv
 # import datetime
-
+import numpy as np
+from phantoms.MR_XCAT_qMRI.sim_ivim_sig import phantom
 
 def pytest_addoption(parser):
     parser.addoption(
@@ -263,7 +264,8 @@ def algorithmlist(algorithms):
     for algorithm in algorithms["algorithms"]:
         algorithm_dict = algorithms.get(algorithm, {})
         requires_matlab = algorithm_dict.get("requires_matlab", False)
-        yield algorithm, requires_matlab, algorithm_dict.get('deep_learning', False)
+        if not algorithm_dict.get('deep_learning', False):
+            yield algorithm, requires_matlab
 
 
 def bound_input(datafile, algorithms):
@@ -302,3 +304,16 @@ def deep_learning_algorithms(datafile, algorithms):
             requires_matlab = algorithm_dict.get("requires_matlab", False)
             tolerances = algorithm_dict.get("tolerances", {"atol":{"f": 2e-1, "D": 8e-4, "Dp": 8e-2},"rtol":{"f": 0.2, "D": 0.3, "Dp": 0.4}})
             yield algorithm, all_data, bvals, kwargs, requires_matlab, tolerances
+
+
+@pytest.fixture(scope="session")
+def threeddata(request):
+    current_folder = pathlib.Path.cwd()
+    datafile = request.config.getoption("dataFile")
+    generic = current_folder / datafile
+    with generic.open() as f:
+        all_data = json.load(f)
+    bvals = all_data.pop('config')
+    bvals = np.array(bvals['bvalues'])
+    sig, _, Dim, fim, Dpim, _=phantom(bvals, 1/1000, TR=3000, TE=40, motion=False, rician=False, interleaved=False, T1T2=True)
+    return sig[::16,::8,::6,:], Dim[::16,::8,::6], fim[::16,::8,::6], Dpim[::16,::8,::6], bvals

phantoms/MR_XCAT_qMRI/sim_ivim_sig.py

@@ -5,28 +5,33 @@
 import argparse
 import os
 from utilities.data_simulation.Download_data import download_data
+import pathlib
 
 ##########
 # code written by Oliver J Gurney-Champion
 # code adapted from MAtlab code by Eric Schrauben: https://github.com/schrau24/XCAT-ERIC
 # This code generates a 4D IVIM phantom as nifti file
 
 def phantom(bvalue, noise, TR=3000, TE=40, motion=False, rician=False, interleaved=False,T1T2=True):
+    download_data()
     np.random.seed(42)
     if motion:
         states = range(1,21)
     else:
         states = [1]
     for state in states:
         # Load the .mat file
-        mat_data = loadmat('../../download/Phantoms/XCAT_MAT_RESP/XCAT5D_RP_' + str(state) + '_CP_1.mat')
+        project_root = pathlib.Path(__file__).resolve().parent.parent
+        filename = f'XCAT5D_RP_{state}_CP_1.mat'
+        mat_path = project_root / '..' /'download' / 'Phantoms' / 'XCAT_MAT_RESP' / filename
+        mat_data = loadmat(mat_path)
 
         # Access the variables in the loaded .mat file
         XCAT = mat_data['IMG']
         XCAT = XCAT[-1:0:-2,-1:0:-2,10:160:4]
 
         D, f, Ds = contrast_curve_calc()
-        S, Dim, fim, Dpim, legend = XCAT_to_MR_DCE(XCAT, TR, TE, bvalue, D, f, Ds,T1T2=T1T2)
+        S, Dim, fim, Dpim, legend = XCAT_to_MR_IVIM(XCAT, TR, TE, bvalue, D, f, Ds,T1T2=T1T2)
         if state == 1:
             Dim_out = Dim
             fim_out = fim
@@ -187,7 +192,7 @@ def contrast_curve_calc():
     return D, f, Ds
 
 
-def XCAT_to_MR_DCE(XCAT, TR, TE, bvalue, D, f, Ds, b0=3, ivim_cont = True, T1T2=True):
+def XCAT_to_MR_IVIM(XCAT, TR, TE, bvalue, D, f, Ds, b0=3, ivim_cont = True, T1T2=True):
     ###########################################################################################
     # This script converts XCAT tissue values to MR contrast based on the SSFP signal equation.
     # Christopher W. Roy 2018-12-04 # [email protected]
@@ -436,7 +441,6 @@ def parse_bvalues_file(file_path):
     motion = args.motion
     interleaved = args.interleaved
     T1T2 = args.T1T2
-    download_data()
     for key, bvalue in bvalues.items():
         bvalue = np.array(bvalue)
         sig, XCAT, Dim, fim, Dpim, legend = phantom(bvalue, noise, motion=motion, interleaved=interleaved,T1T2=T1T2)

pytest.ini

@@ -3,4 +3,7 @@ markers =
     slow: marks tests as slow (deselect with '-m "not slow"')
 addopts =  
     -m 'not slow'
-testpaths = tests
+testpaths = tests
+filterwarnings =
+    ignore::Warning
+    always::tests.IVIMmodels.unit_tests.test_ivim_fit.PerformanceWarning

src/original/ETP_SRI/LinearFitting.py

@@ -1,5 +1,6 @@
 import numpy as np
 import numpy.polynomial.polynomial as poly
+import warnings
 
 from utilities.data_simulation.GenerateData import GenerateData
 
@@ -82,11 +83,17 @@ def ivim_fit(self, bvalues, signal):
                 # print(Dp_prime)
 
             if np.any(np.asarray(Dp_prime) < 0) or not np.all(np.isfinite(Dp_prime)):
-                print('Perfusion fit failed')
+                warnings.warn('Perfusion fit failed',
+                    category=UserWarning,
+                    stacklevel=2  # Ensures correct file/line info in the warning
+                              )
                 Dp_prime = [0, 0]
             f = signal[0] - D[0]
         else:
-            print("This doesn't seem to be an IVIM set of b-values")
+            warnings.warn('This doesn\'t seem to be an IVIM set of b-values',
+                          category=UserWarning,
+                          stacklevel=2  # Ensures correct file/line info in the warning
+                          )
             f = 1
             Dp_prime = [0, 0]
         D = D[1]

src/original/PvH_KB_NKI/DWI_functions_standalone.py

@@ -10,6 +10,8 @@
 """
 
 import numpy as np
+import warnings
+
 
 def generate_ADC_standalone(DWIdata, bvalues, bmin: int = 150, bmax: int = 1000, specificBvals: list = []):
     """
@@ -116,9 +118,11 @@ def generate_IVIMmaps_standalone(DWIdata, bvalues, bminADC=150, bmaxADC=1000, bm
             raise Exception('b=0 was not acquired')
 
         if len(bvalues) > 2:
-            print('More than two b-values were detected for D* calculation. ' +
-                  'Note that only the first two will be used for D*.')
-
+            warnings.warn(
+                'More than two b-values were detected for D* calculation. ' + 'Note that only the first two will be used for D*.',
+                category=UserWarning,
+                stacklevel=2  # Ensures correct file/line info in the warning
+            )
     except Exception as e:
         print("Could not calculate D* due to an error: " + str(e))
         return

src/original/TF_reference/segmented_IVIMfit.py

@@ -24,6 +24,7 @@ def segmented_IVIM_fit(bvalues, dw_data, b_cutoff = 200, bounds=([0.0001, 0.0, 0
         """
     bvalues_D = bvalues[bvalues >= b_cutoff]
     dw_data_D = dw_data[bvalues >= b_cutoff]
+    dw_data_D = np.clip(dw_data_D, 1e-6, None) # ensure we do not get 0 or negative values that cause nans/infs
     log_data_D = np.log(dw_data_D)
 
     D, b0_intercept = d_fit_iterative_wls(bvalues_D, log_data_D)

src/standardized/ETP_SRI_LinearFitting.py

@@ -1,6 +1,8 @@
 import numpy as np
 from src.wrappers.OsipiBase import OsipiBase
 from src.original.ETP_SRI.LinearFitting import LinearFit
+import warnings
+warnings.simplefilter('once', UserWarning)
 
 
 class ETP_SRI_LinearFitting(OsipiBase):

src/standardized/IAR_LU_biexp.py

@@ -114,8 +114,9 @@ def ivim_fit_full_volume(self, signals, bvalues, **kwargs):
             gtab = gradient_table(bvalues, bvec, b0_threshold=0)
 
             self.IAR_algorithm = IvimModelBiExp(gtab, bounds=self.bounds, initial_guess=self.initial_guess)
-
-        fit_results = self.IAR_algorithm.fit(signals)
+        b0_index = np.where(bvalues == 0)[0][0]
+        mask = signals[...,b0_index]>0
+        fit_results = self.IAR_algorithm.fit(signals, mask=mask)
 
         results = {}
         results["f"] = fit_results.model_params[..., 1]

src/standardized/IAR_LU_modified_mix.py

@@ -34,6 +34,7 @@ class IAR_LU_modified_mix(OsipiBase):
     supported_dimensions = 1
     supported_priors = False
 
+
     def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, weighting=None, stats=False):
         """
             Everything this algorithm requires should be implemented here.
@@ -47,6 +48,10 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             print('warning, bounds from wrapper are not (yet) used in this algorithm')
         self.use_bounds = False
         self.use_initial_guess = False
+
+        # Additional options
+        self.stochastic = True
+
         # Check the inputs
 
         # Initialize the algorithm

src/standardized/IVIM_NEToptim.py

@@ -56,6 +56,9 @@ def initialize(self, bounds, initial_guess, fitS0, traindata, SNR):
         self.fitS0=fitS0
         self.deep_learning = True
         self.supervised = False
+        # Additional options
+        self.stochastic = True
+
         if traindata is None:
             warnings.warn('no training data provided (traindata = None). Training data will be simulated')
             if SNR is None:

src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py

@@ -80,7 +80,7 @@ def initialize(self, bounds=None, initial_guess=None, fitS0=True, prior_in=None,
             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]])
         else:
             print('warning, bounds are not used, as a prior is used instead')
-            if len(prior_in) is 4:
+            if len(prior_in) == 4:
                 self.neg_log_prior = empirical_neg_log_prior(prior_in[0], prior_in[1], prior_in[2],prior_in[3])
             else:
                 self.neg_log_prior = empirical_neg_log_prior(prior_in[0], prior_in[1], prior_in[2])
@@ -113,37 +113,61 @@ def ivim_fit(self, signals, bvalues, initial_guess=None, **kwargs):
 
         return results
 
-    def ivim_fit_full_volume(self, signals, bvalues, initial_guess=None, **kwargs):
+    def ivim_fit_full_volume(self, signals, bvalues, njobs=4, **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.reshape_to_voxelwise(signals)
+        # normalize signals
+        # Get index of b=0
+        shape=np.shape(signals)
+
+        b0_index = np.where(bvalues == 0)[0][0]
+        # Mask of voxels where signal at b=0 >= 0.5
+        valid_mask = signals[..., b0_index] >= 0
+        # Select only valid voxels for fitting
+        signals = signals[valid_mask]
+
+        minimum_bvalue = np.min(bvalues) # We normalize the signal to the minimum bvalue. Should be 0 or very close to 0.
+        b0_indices = np.where(bvalues == minimum_bvalue)[0]
+        normalization_factor = np.mean(signals[..., b0_indices],axis=-1)
+        signals = signals / np.repeat(normalization_factor[...,np.newaxis],np.shape(signals)[-1],-1)
 
         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)
+        fit_results = np.array(fit_segmented_array(bvalues, signals, bounds=self.bounds, cutoff=self.thresholds, p0=self.initial_guess))
+        #fit_results=np.array(fit_results+(1,))
+        # Loop over parameters (rows)
 
+        for i in range(4):
+            if i == 3:
+                fit_results[i] = np.random.normal(1,0.2,np.shape(fit_results[i]))
+            else:
+                below = fit_results[i] < self.bounds[0][i]
+                above = fit_results[i] > self.bounds[1][i]
+
+                fit_results[i, below] = self.bounds[0][i] + epsilon
+                fit_results[i, above] = self.bounds[1][i] - epsilon
+        self.jobs=njobs
+        fit_results = self.OGC_algorithm_array(bvalues, signals,fit_results, self)
+
+        D=np.zeros(shape[0:-1])
+        D[valid_mask]=fit_results[0]
+        f=np.zeros(shape[0:-1])
+        f[valid_mask]=fit_results[1]
+        Dp=np.zeros(shape[0:-1])
+        Dp[valid_mask]=fit_results[2]
         results = {}
-        results["D"] = fit_results[0]
-        results["f"] = fit_results[1]
-        results["Dp"] = fit_results[2]
-
+        results["D"] = D
+        results["f"] = f
+        results["Dp"] = Dp
         return results
 
 
-
     def reshape_to_voxelwise(self, data):
         """
         reshapes multi-D input (spatial dims, bvvalue) data to 2D voxel-wise array
@@ -154,4 +178,4 @@ def reshape_to_voxelwise(self, data):
         """
         B = data.shape[-1]
         voxels = int(np.prod(data.shape[:-1]))  # e.g., X*Y*Z
-        return data.reshape(voxels, B)
+        return data.reshape(voxels, B)

src/standardized/OGC_AmsterdamUMC_biexp.py

@@ -84,25 +84,4 @@ 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

@@ -86,4 +86,4 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         results["f"] = fit_results[1]
         results["Dp"] = fit_results[2]
 
-        return results
+        return results

src/standardized/PvH_KB_NKI_IVIMfit.py

@@ -1,6 +1,8 @@
 from src.wrappers.OsipiBase import OsipiBase
 from src.original.PvH_KB_NKI.DWI_functions_standalone import generate_IVIMmaps_standalone, generate_ADC_standalone
 import numpy as np
+import warnings
+warnings.simplefilter('once', UserWarning)
 
 class PvH_KB_NKI_IVIMfit(OsipiBase):
     """

src/standardized/Super_IVIM_DC.py

@@ -62,6 +62,10 @@ def initialize(self, bounds, initial_guess, fitS0, SNR, working_dir=os.getcwd(),
         self.use_bounds = False
         self.deep_learning = True
         self.supervised = True
+
+        # Additional options
+        self.stochastic = True
+
         modeldir = Path(working_dir)  # Ensure it's a Path object
         modeldir = modeldir / "models"
         modeldir.mkdir(parents=True, exist_ok=True)