From 3fd75a736cc84b0a8805c2df22f12286cc04090c Mon Sep 17 00:00:00 2001 From: Eric Peterson Date: Tue, 23 May 2023 22:36:52 -0700 Subject: [PATCH 1/6] Add an example that uses linear fitting --- .github/workflows/unit_test.yml | 2 +- src/original/ETP_SRI/LinearFitting.py | 86 +++++++++++ src/original/ETP_SRI/__init__.py | 0 .../unit_tests/test_ivim_fit_linear.py | 47 ++++++ utils/data_simulation/GenerateData.py | 140 ++++++++++++++++++ 5 files changed, 274 insertions(+), 1 deletion(-) create mode 100644 src/original/ETP_SRI/LinearFitting.py create mode 100644 src/original/ETP_SRI/__init__.py create mode 100644 tests/IVIMmodels/unit_tests/test_ivim_fit_linear.py create mode 100644 utils/data_simulation/GenerateData.py diff --git a/.github/workflows/unit_test.yml b/.github/workflows/unit_test.yml index 6331cc4..441bdc9 100644 --- a/.github/workflows/unit_test.yml +++ b/.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 diff --git a/src/original/ETP_SRI/LinearFitting.py b/src/original/ETP_SRI/LinearFitting.py new file mode 100644 index 0000000..a89589d --- /dev/null +++ b/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] diff --git a/src/original/ETP_SRI/__init__.py b/src/original/ETP_SRI/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/IVIMmodels/unit_tests/test_ivim_fit_linear.py b/tests/IVIMmodels/unit_tests/test_ivim_fit_linear.py new file mode 100644 index 0000000..bcb9251 --- /dev/null +++ b/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) diff --git a/utils/data_simulation/GenerateData.py b/utils/data_simulation/GenerateData.py new file mode 100644 index 0000000..a41fe77 --- /dev/null +++ b/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 \ No newline at end of file From f84e432533dccacd80d35c319c324cb2af66aa72 Mon Sep 17 00:00:00 2001 From: Eric Peterson Date: Tue, 23 May 2023 22:56:05 -0700 Subject: [PATCH 2/6] Try to not fail on job errors --- .github/workflows/unit_test.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/unit_test.yml b/.github/workflows/unit_test.yml index 441bdc9..842e207 100644 --- a/.github/workflows/unit_test.yml +++ b/.github/workflows/unit_test.yml @@ -6,8 +6,9 @@ jobs: build: runs-on: ${{ matrix.os }} - + continue-on-error: true strategy: + fail-fast: false matrix: os: [ubuntu-latest, macos-latest, windows-latest] python-version: ["3.7", "3.8", "3.9", "3.10", "3.11"] From 9f48d4208b2cd2ae1ba9ac0b9a8d7015cd4ad0c4 Mon Sep 17 00:00:00 2001 From: Eric Peterson Date: Thu, 25 May 2023 21:38:40 -0700 Subject: [PATCH 3/6] Added requirements and syntax fix to allow testing --- requirements.txt | 4 +++- src/wrappers/ivim_fit.py | 4 ++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/requirements.txt b/requirements.txt index f13a38a..94d89bc 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,4 +2,6 @@ numpy scipy torchio torch -logging \ No newline at end of file +logging +joblib +dipy \ No newline at end of file diff --git a/src/wrappers/ivim_fit.py b/src/wrappers/ivim_fit.py index 13e3f82..1ca99c2 100644 --- a/src/wrappers/ivim_fit.py +++ b/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 = [] From a58f53159b6af7b484f1581e720ad7dd0740eef3 Mon Sep 17 00:00:00 2001 From: Eric Peterson Date: Thu, 25 May 2023 21:50:57 -0700 Subject: [PATCH 4/6] Added more requirements and trying no continue on error and file fixes --- .github/workflows/unit_test.yml | 2 +- requirements.txt | 3 ++- src/original/OGC_AmsterdamUMC/LSQ_fitting.py | 2 +- utils/data_simulation/ivim_simulation.py | 2 +- 4 files changed, 5 insertions(+), 4 deletions(-) diff --git a/.github/workflows/unit_test.yml b/.github/workflows/unit_test.yml index 842e207..01a0812 100644 --- a/.github/workflows/unit_test.yml +++ b/.github/workflows/unit_test.yml @@ -6,7 +6,7 @@ jobs: build: runs-on: ${{ matrix.os }} - continue-on-error: true + continue-on-error: false strategy: fail-fast: false matrix: diff --git a/requirements.txt b/requirements.txt index 94d89bc..3b1ef29 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,4 +4,5 @@ torchio torch logging joblib -dipy \ No newline at end of file +dipy +matplotlib \ No newline at end of file diff --git a/src/original/OGC_AmsterdamUMC/LSQ_fitting.py b/src/original/OGC_AmsterdamUMC/LSQ_fitting.py index 08bd042..c7e8d29 100644 --- a/src/original/OGC_AmsterdamUMC/LSQ_fitting.py +++ b/src/original/OGC_AmsterdamUMC/LSQ_fitting.py @@ -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): """ diff --git a/utils/data_simulation/ivim_simulation.py b/utils/data_simulation/ivim_simulation.py index c4d2848..c43841e 100644 --- a/utils/data_simulation/ivim_simulation.py +++ b/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): From c4ab12b9235f46b24223647ff31202b43507f002 Mon Sep 17 00:00:00 2001 From: Eric Peterson Date: Thu, 25 May 2023 22:02:05 -0700 Subject: [PATCH 5/6] Python version 3.7 seemed too old, add more requirements --- .github/workflows/unit_test.yml | 2 +- requirements.txt | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/.github/workflows/unit_test.yml b/.github/workflows/unit_test.yml index 01a0812..65f1929 100644 --- a/.github/workflows/unit_test.yml +++ b/.github/workflows/unit_test.yml @@ -11,7 +11,7 @@ jobs: 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" diff --git a/requirements.txt b/requirements.txt index 3b1ef29..1e52fbb 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,4 +5,5 @@ torch logging joblib dipy -matplotlib \ No newline at end of file +matplotlib +scienceplots \ No newline at end of file From ba55cd318a1540605ba0cac5cc725d7d7f04f3bb Mon Sep 17 00:00:00 2001 From: Eric Peterson Date: Thu, 25 May 2023 22:08:39 -0700 Subject: [PATCH 6/6] Added cvxpy requirement and code fixes --- requirements.txt | 3 ++- src/original/OGC_AmsterdamUMC/LSQ_fitting.py | 6 +++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/requirements.txt b/requirements.txt index 1e52fbb..acb0486 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,4 +6,5 @@ logging joblib dipy matplotlib -scienceplots \ No newline at end of file +scienceplots +cvxpy \ No newline at end of file diff --git a/src/original/OGC_AmsterdamUMC/LSQ_fitting.py b/src/original/OGC_AmsterdamUMC/LSQ_fitting.py index c7e8d29..fd3166b 100644 --- a/src/original/OGC_AmsterdamUMC/LSQ_fitting.py +++ b/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: