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/data/test_GenerateData.py

@@ -0,0 +1,71 @@
+import numpy as np
+import numpy.testing as npt
+import pytest
+import torch
+
+from utils.data_simulation.GenerateData import GenerateData
+
+#run using python -m pytest from the root folder
+
+test_monoexponential_data = [
+    pytest.param(0, np.linspace(0, 1000, 11), id='0'),
+    pytest.param(0.1, np.linspace(0, 1000, 11), id='0.1'),
+    pytest.param(0.2, np.linspace(0, 1000, 11), id='0.2'),
+    pytest.param(0.3, np.linspace(0, 1000, 11), id='0.3'),
+    pytest.param(0.4, np.linspace(0, 1000, 11), id='0.4'),
+    pytest.param(0.5, np.linspace(0, 1000, 11), id='0.5'),
+    pytest.param(0.8, np.linspace(0, 1000, 11), id='0.8'),
+    pytest.param(1, np.linspace(0, 1000, 11), id='1'),
+]
+@pytest.mark.parametrize("D, bvals", test_monoexponential_data)
+def test_monoexponential(D, bvals):
+    gd = GenerateData()
+    gd_signal = gd.exponential_signal(D, bvals)
+    testing_signal = np.exp(-D * np.asarray(bvals, dtype='float64'))
+    npt.assert_allclose(gd_signal, testing_signal)
+    assert(gd_signal[0] >= testing_signal[0])
+
+test_ivim_data = [
+    pytest.param(0.01, 0.0, 1, 1, np.linspace(0, 1000, 11), None),
+    pytest.param(0.01, 0.1, 0.05, 1, np.linspace(0, 1000, 11), None),
+    pytest.param(0.05, 0.2, 0.1, 1, np.linspace(0, 800, 11), None),
+    pytest.param(0.04, 0.15, 0.25, 1.5, np.linspace(0, 1000, 2), 10),
+    pytest.param(0.1, 0.5, 0.5, 0.5, np.linspace(0, 1500, 5), 100),
+    pytest.param(0.01, 0.2, 0.1, 1, np.linspace(10, 1500, 11), 100),
+    pytest.param(0.1, 0.15, 0.05, 1, np.linspace(10, 1000, 8), 5)
+]
+@pytest.mark.parametrize('D, Dp, f, S0, bvals, snr', test_ivim_data)
+def test_ivim(D, Dp, f, S0, bvals, snr):
+    gd = GenerateData()
+    gd_signal = gd.ivim_signal(D, Dp, f, S0, bvals, snr)
+    testing_signal = S0 * ((1 - f) * np.exp(-D * bvals) + f * np.exp(-Dp * bvals))
+    atol = 0.0
+    if snr is not None:
+        atol = 4 / snr
+    npt.assert_allclose(gd_signal, testing_signal, atol=atol)
+
+
+test_linear_data = [
+    pytest.param(0, np.linspace(0, 1000, 11), 0, id='0'),
+    pytest.param(0.1, np.linspace(0, 1000, 11), 10, id='0.1'),
+    pytest.param(0.2, np.linspace(0, 1000, 11), -10, id='0.2'),
+    pytest.param(0.3, np.linspace(0, 1000, 11), 0, id='0.3'),
+    pytest.param(0.4, np.linspace(0, 1000, 11), 0, id='0.4'),
+    pytest.param(0.5, np.linspace(0, 1000, 11), 0, id='0.5'),
+    pytest.param(0.8, np.linspace(0, 1000, 11), 0, id='0.8'),
+    pytest.param(1, np.linspace(0, 1000, 11), 0, id='1'),
+]
+@pytest.mark.parametrize("D, bvals, offset", test_linear_data)
+def test_linear(D, bvals, offset):
+    gd = GenerateData()
+    gd_signal = gd.linear_signal(D, bvals, offset)
+    testing_signal = -D * np.asarray(bvals, dtype='float64')
+    testing_signal += offset
+    npt.assert_allclose(gd_signal, testing_signal)
+    assert(gd_signal[0] >= testing_signal[0])
+
+    gd_exponential = gd.exponential_signal(D, bvals)
+    gd_log_exponential = np.log(gd_exponential) + offset
+    real_mask = np.isfinite(gd_log_exponential)
+
+    npt.assert_allclose(gd_log_exponential[real_mask], gd_signal[real_mask])

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)