Add an example that uses linear fitting

Author: etpeterson

.github/workflows/unit_test.yml

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

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/ETP_SRI/__init__.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)

tests/IVIMmodels/unit_tests/test_ivim_fit_linear.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