Merge pull request #5 from OSIPI/peterson/linear-example

Peterson/linear example

.github/workflows/unit_test.yml

@@ -6,11 +6,12 @@ jobs:
   build:
 
     runs-on: ${{ matrix.os }}
-
+    continue-on-error: false
     strategy:
+      fail-fast: false
       matrix:
         os: [ubuntu-latest, macos-latest, windows-latest]
-        python-version: ["3.7", "3.8", "3.9", "3.10", "3.11"]
+        python-version: ["3.8", "3.9", "3.10", "3.11"]
         # exclude:
         #   - os: macos-latest
         #     python-version: "3.7"
@@ -32,4 +33,4 @@ jobs:
       - name: Test with pytest
         run: |
           pip install pytest pytest-cov
-          pytest tests.py --doctest-modules --junitxml=junit/test-results.xml --cov=com --cov-report=xml --cov-report=html
+          python -m pytest --doctest-modules --junitxml=junit/test-results.xml --cov=com --cov-report=xml --cov-report=html

requirements.txt

@@ -2,4 +2,9 @@ numpy
 scipy
 torchio
 torch
-logging
+logging
+joblib
+dipy
+matplotlib
+scienceplots
+cvxpy

src/original/ETP_SRI/LinearFitting.py

@@ -0,0 +1,86 @@
+import numpy as np
+import numpy.polynomial.polynomial as poly
+
+from utils.data_simulation.GenerateData import GenerateData
+
+
+
+class LinearFit:
+    """
+    Performs linear fits of exponential data
+    """
+    def __init__(self, linear_cutoff=500):
+        """
+        Parameters
+        ----------
+        linear_cutoff : float
+            The b-value after which it can be assumed that the perfusion value is negligible
+        """
+        self.linear_cutoff = linear_cutoff
+
+    def linear_fit(self, bvalues, signal, weighting=None, stats=False):
+        """
+        Fit a single line
+
+        Parameters
+        ----------
+        bvalues : list or array of float
+            The diffusion (b-values)
+        signal : list or array of float
+            The acquired signal to fit. It is assumed to be linearized at this point.
+        weighting : list or array fo float
+            Weights to pass into polyfit. None is no weighting.
+        stats : boolean
+            If true, return the polyfit statistics
+        """
+        assert bvalues.size == signal.size, "Signal and b-values don't have the same number of values"
+        if stats:
+            D, stats = poly.polyfit(np.asarray(bvalues), signal, 1, full=True, w=weighting)
+            return [np.exp(D[0]), *-D[1:]], stats
+        else:
+            D = poly.polyfit(np.asarray(bvalues), signal, 1, w=weighting)
+            return [np.exp(D[0]), *-D[1:]]
+
+    def ivim_fit(self, bvalues, signal):
+        """
+        Fit an IVIM curve
+        This fits a bi-exponential curve using linear fitting only
+
+
+        Parameters
+        ----------
+        bvalues : list or array of float
+            The diffusion (b-values)
+        signal : list or array of float
+            The acquired signal to fit. It is assumed to be exponential at this point
+        """
+        bvalues = np.asarray(bvalues)
+        assert bvalues.size > 1, 'Too few b-values'
+        signal = np.asarray(signal)
+        assert bvalues.size == signal.size, "Signal and b-values don't have the same number of values"
+        gd = GenerateData()
+        lt_cutoff = bvalues <= self.linear_cutoff
+        gt_cutoff = bvalues >= self.linear_cutoff
+        linear_signal = np.log(signal)
+        D = self.linear_fit(bvalues[gt_cutoff], linear_signal[gt_cutoff])
+
+        if lt_cutoff.sum() > 0:
+            signal_Dp = linear_signal[lt_cutoff] - gd.linear_signal(D[1], bvalues[lt_cutoff], np.log(D[0]))
+
+            Dp_prime = self.linear_fit(bvalues[lt_cutoff], np.log(signal_Dp))
+
+            if np.any(np.asarray(Dp_prime) < 0) or not np.all(np.isfinite(Dp_prime)):
+                print('Perfusion fit failed')
+                Dp_prime = [0, 0]
+            f = signal[0] - D[0]
+        else:
+            print("This doesn't seem to be an IVIM set of b-values")
+            f = 1
+            Dp_prime = [0, 0]
+        D = D[1]
+        Dp = D + Dp_prime[1]
+        if np.allclose(f, 0):
+            Dp = 0
+        elif np.allclose(f, 1):
+            D = 0
+        return [f, D, Dp]

src/original/OGC_AmsterdamUMC/LSQ_fitting.py

@@ -494,7 +494,7 @@ def empirical_neg_log_prior(Dt0, Fp0, Dp0, S00=None):
     # define the prior
     def neg_log_prior(p):
         # depends on whether S0 is fitted or not
-        if len(p) is 4:
+        if len(p) == 4:
             Dt, Fp, Dp, S0 = p[0], p[1], p[2], p[3]
         else:
             Dt, Fp, Dp = p[0], p[1], p[2]
@@ -507,7 +507,7 @@ def neg_log_prior(p):
             Dt_prior = stats.lognorm.pdf(Dt, Dt_shape, scale=Dt_scale)
             Fp_prior = stats.beta.pdf(Fp, Fp_a, Fp_b)
             # determine and return the prior for D, f and D* (and S0)
-            if len(p) is 4:
+            if len(p) == 4:
                 S0_prior = stats.beta.pdf(S0 / 2, S0_a, S0_b)
                 return -np.log(Dp_prior + eps) - np.log(Dt_prior + eps) - np.log(Fp_prior + eps) - np.log(
                     S0_prior + eps)
@@ -525,7 +525,7 @@ def neg_log_likelihood(p, bvalues, dw_data):
     :param dw_data: 1D Array diffusion-weighted data
     :returns: the log-likelihood of the parameters given the data
     """
-    if len(p) is 4:
+    if len(p) == 4:
         return 0.5 * (len(bvalues) + 1) * np.log(
             np.sum((ivim(bvalues, p[0], p[1], p[2], p[3]) - dw_data) ** 2))  # 0.5*sum simplified
     else:
@@ -705,7 +705,7 @@ def goodness_of_fit(bvalues, Dt, Fp, Dp, S0, dw_data, Fp2=None, Dp2=None):
         # plt.show()
         # print(R2[vox])
     return R2, adjust
-ed_R2
+# ed_R2
 
 def MSE(bvalues, Dt, Fp, Dp, S0, dw_data):
     """

src/wrappers/ivim_fit.py

@@ -24,8 +24,8 @@ def ivim_fit(author, signals=None, bvalues=None, data=None, initial_guess=None,
 
 
     # Some implementations can only fit a voxel at a time (i.e. all inputs must be 2-dimensional)
-    requires_2D = True if author in []
-    requires_4D = True if author in []
+    requires_2D = True if author in [] else False
+    requires_4D = True if author in [] else False
 
     # Bounds and initial guess for parameters
     initial_guess = []

tests/IVIMmodels/unit_tests/test_ivim_fit_linear.py

@@ -0,0 +1,47 @@
+import numpy as np
+import numpy.testing as npt
+import pytest
+import torch
+
+from utils.data_simulation.GenerateData import GenerateData
+from src.original.ETP_SRI.LinearFitting import LinearFit
+
+
+#run using python -m pytest from the root folder
+
+test_linear_data = [
+    pytest.param(0, np.linspace(0, 1000, 11), id='0'),
+    pytest.param(0.01, np.linspace(0, 1000, 11), id='0.1'),
+    pytest.param(0.02, np.linspace(0, 1000, 11), id='0.2'),
+    pytest.param(0.03, np.linspace(0, 1000, 11), id='0.3'),
+    pytest.param(0.04, np.linspace(0, 1000, 11), id='0.4'),
+    pytest.param(0.05, np.linspace(0, 1000, 11), id='0.5'),
+    pytest.param(0.08, np.linspace(0, 1000, 11), id='0.8'),
+    pytest.param(0.1, np.linspace(0, 1000, 11), id='1'),
+]
+@pytest.mark.parametrize("D, bvals", test_linear_data)
+def test_linear_fit(D, bvals):
+    gd = GenerateData()
+    gd_signal = gd.exponential_signal(D, bvals)
+    print(gd_signal)
+    fit = LinearFit()
+    D_fit = fit.linear_fit(bvals, np.log(gd_signal))
+    npt.assert_allclose([1, D], D_fit)
+
+test_ivim_data = [
+    pytest.param(0, 0.01, 0.05, np.linspace(0, 1000, 11), id='0'),
+    pytest.param(0.1, 0.01, 0.05, np.linspace(0, 1000, 11), id='0.1'),
+    pytest.param(0.2, 0.01, 0.05, np.linspace(0, 1000, 11), id='0.2'),
+    pytest.param(0.1, 0.05, 0.1, np.linspace(0, 1000, 11), id='0.3'),
+    pytest.param(0.4, 0.001, 0.05, np.linspace(0, 1000, 11), id='0.4'),
+    pytest.param(0.5, 0.001, 0.05, np.linspace(0, 1000, 11), id='0.5'),
+]
+@pytest.mark.parametrize("f, D, Dp, bvals", test_ivim_data)
+def test_ivim_fit(f, D, Dp, bvals):
+    gd = GenerateData()
+    gd_signal = gd.ivim_signal(D, Dp, f, 1, bvals)
+    fit = LinearFit()
+    [f_fit, D_fit, Dp_fit] = fit.ivim_fit(bvals, gd_signal)
+    npt.assert_allclose([f, D], [f_fit, D_fit], atol=1e-5)
+    if not np.allclose(f, 0):
+        npt.assert_allclose(Dp, Dp_fit, rtol=1e-2, atol=1e-3)

utils/data_simulation/GenerateData.py

@@ -0,0 +1,140 @@
+import numpy as np
+
+class GenerateData:
+    """
+    Generate exponential and linear data
+    """
+    def __init__(self, operator=None):
+        """
+        Parameters
+        ----------
+        operator : numpy or torch
+            Must provide mathematical operators
+        """
+        if operator is None:
+            self._op = np
+        else:
+            self._op = operator
+
+    def ivim_signal(self, D, Dp, f, S0, bvalues, snr=None):
+        """
+        Generates IVIM (biexponential) signal
+
+        Parameters
+        ----------
+        D : float
+            The tissue diffusion value
+        Dp : float
+            The pseudo perfusion value
+        f : float
+            The fraction of the signal from perfusion
+        S0 : float
+            The baseline signal (magnitude at no diffusion)
+        bvalues : list or array of float
+            The diffusion (b-values)
+        """
+        signal = self.multiexponential_signal([D, Dp], [1 - f, f], S0, self._op.asarray(bvalues, dtype='float64'))
+        return self.add_rician_noise(signal, snr)
+
+    def exponential_signal(self, D, bvalues):
+        """
+        Generates exponential signal
+
+        Parameters
+        ----------
+        D : float
+            The tissue diffusion value
+        bvalues : list or array of float
+            The diffusion (b-values)
+        """
+        assert D >= 0, 'D must be >= 0'
+        return self._op.exp(-self._op.asarray(bvalues, dtype='float64') * D)
+
+    def multiexponential_signal(self, D, F, S0, bvalues):
+        """
+        Generates multiexponential signal
+        The combination of exponential signals
+
+        Parameters
+        ----------
+        D : list or arrray of float
+            The tissue diffusion value
+        F : list or array of float
+            The fraction of the signal from perfusion
+        S0 : list or array of float
+            The baseline signal (magnitude at no diffusion)
+        bvalues : list or array of float
+            The diffusion (b-values)
+        """
+        assert len(D) == len(F), 'D and F must be the same length'
+        signal = self._op.zeros_like(bvalues)
+        for [d, f] in zip(D, F):
+            signal += f * self.exponential_signal(d, bvalues)
+        signal *= S0
+        return signal
+
+    def add_rician_noise(self, real_signal, snr=None, imag_signal=None):
+        """
+        Adds Rician noise to a real or complex signal
+
+        Parameters
+        ----------
+        real_signal : list or array of float
+            The real channel float
+        snr : float
+            The signal to noise ratio
+        imag_signal : list or array of float
+            The imaginary channel float
+        """
+        if imag_signal is None:
+            imag_signal = self._op.zeros_like(real_signal)
+        if snr is None:
+            noisy_data = self._op.sqrt(self._op.power(real_signal, 2) + self._op.power(imag_signal, 2))
+        else:
+            real_noise = self._op.random.normal(0, 1 / snr, real_signal.shape)
+            imag_noise = self._op.random.normal(0, 1 / snr, imag_signal.shape)
+            noisy_data = self._op.sqrt(self._op.power(real_signal + real_noise, 2) + self._op.power(imag_signal + imag_noise, 2))
+        return noisy_data
+
+    def linear_signal(self, D, bvalues, offset=0):
+        """
+        Generates linear signal
+
+        Parameters
+        ----------
+        D : float
+            The tissue diffusion value
+        bvalues : list or array of float
+            The diffusion (b-values)
+        offset : float
+            The signal offset
+        """
+        assert D >= 0, 'D must be >= 0'
+        data = -D * np.asarray(bvalues)
+        return data + offset
+
+    def multilinear_signal(self, D, F, S0, bvalues, offset=0):
+        """
+        Generates multilinear signal
+        The combination of multiple linear signals
+
+        Parameters
+        ----------
+        D : list or arrray of float
+            The tissue diffusion value
+        F : list or array of float
+            The fraction of the signal from perfusion
+        S0 : list or array of float
+            The baseline signal (magnitude at no diffusion)
+        bvalues : list or array of float
+            The diffusion (b-values)
+        offset : float
+            The signal offset
+        """
+        assert len(D) == len(F), 'D and F must be the same length'
+        signal = self._op.zeros_like(bvalues)
+        for [d, f] in zip(D, F):
+            signal += f * self.linear_signal(d, bvalues)
+        signal *= S0
+        signal += offset
+        return signal

utils/data_simulation/ivim_simulation.py

@@ -1,7 +1,7 @@
 import torch
 import numpy as np
 
-from utils.ivim.ivim_fit import ivim_parameters_to_signal
+from utils.ivim.forward_model import ivim_parameters_to_signal
 
 
 def simulate_ivim_signal(D, Dp, f, S0, bvalues, SNR_array, rg):