matlab wrapper implemented for IVIM_seg (OJ_GU)

.github/workflows/unit_test.yml

@@ -26,11 +26,17 @@ jobs:
       # You can test your matrix by printing the current Python version
       - name: Display Python version
         run: python -c "import sys; print(sys.version)"
+      # Sets up MATLAB on the GitHub Actions runner
+      - name: Setup MATLAB
+        uses: matlab-actions/setup-matlab@v1
+        with:
+          release: R2023a
+      # Installs other dependenvies
       - name: Install dependencies
         run: |
           python -m pip install --upgrade pip
           pip install -r requirements.txt
       - 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

requirements.txt

@@ -7,4 +7,5 @@ joblib
 dipy
 matplotlib
 scienceplots
-cvxpy
+cvxpy
+matlabengine

src/original/OJ_GU/IVIM_fitting.py

@@ -0,0 +1,18 @@
+import os
+import matlab.engine
+
+def IVIM_seg(Y, b, lim, blim, disp_prog):
+    eng = matlab.engine.start_matlab()
+    s = eng.genpath(os.path.dirname(__file__))
+    eng.addpath(s, nargout=0)
+    pars = eng.IVIM_seg(Y, b, lim, blim, disp_prog)
+    eng.quit()
+    return pars
+
+def IVIM_bayes(Y, f, D, Dstar, S0, b, lim, n, rician, prior, burns, meanonly):
+    eng = matlab.engine.start_matlab()
+    s = eng.genpath(os.path.dirname(__file__))
+    eng.addpath(s, nargout=0)
+    out = eng.IVIM_bayes(Y, f, D, Dstar, S0, b, lim, n, rician, prior, burns, meanonly)
+    eng.quit()
+    return out

src/original/OJ_GU/IVIM_seg.m

@@ -190,14 +190,18 @@
         maskDstar = maskf;
 
         % Optimization
-        [Dstar(maskf),maskDstar(maskf)] = optimizeD(Y(:,maskf)-repmat(A(maskf),length(b),1).*exp(-b*D(maskf)),b,optlim,lim(:,4),disp_prog);
-
+        if sum(maskf) > 0
+            [Dstar(maskf),maskDstar(maskf)] = optimizeD(Y(:,maskf)-repmat(A(maskf),length(b),1).*exp(-b*D(maskf)),b,optlim,lim(:,4),disp_prog);
+        end
+
         % Calculates f
         f = lim(1,1)*ones(1,n);
         f(f1 < lim(1,3)) = lim(1,3);
         f(f2 > lim(2,3)) = lim(2,3); 
-        f(maskDstar) = fcalc(Dstar(maskDstar),maskDstar);
-
+        if sum(maskDstar) > 0
+            f(maskDstar) = fcalc(Dstar(maskDstar),maskDstar);
+        end
+
         % Checks for f out of bounds
         maskf = maskf & (f > lim(1,3)) & (f < lim(2,3));
         Dstar(f < lim(1,3)) = lim(1,4);

src/original/OJ_GU/simple_test_of_fits.py

@@ -0,0 +1,65 @@
+import numpy as np
+from scipy.stats import norm
+from IVIM_fitting import IVIM_seg
+
+def ivim_signal(b, S0, f, D_star, D):
+    return S0*(f*np.exp(-b*D_star) + (1-f)*np.exp(-b*D))
+
+def diffusion_signal(b, S0, f, D):
+    return S0*(1-f)*np.exp(-b*D)
+
+def generate_noise(loc, sigma):
+    real_component = norm.rvs(loc=loc, scale=sigma/loc)
+    imaginary_component = norm.rvs(loc=loc, scale=sigma/loc)
+    return np.absolute(complex(real_component, imaginary_component))
+
+def add_rician_noise(signal, SNR):
+    sigma = signal[-1]/SNR
+    # Sample real and imaginary noise components from gaussian distributions
+    # Use the last b-value as the SNR baseline in order to avoid the noise floor
+    noise = np.array([generate_noise(signal_value, sigma) for signal_value in signal])
+
+    # Add the two components to the signal and take the magniutde of the result
+    noised_signal = signal + noise
+    noised_signal = np.absolute(noised_signal)
+
+    return noised_signal
+
+# Ground truth
+factor = 1
+S0 = 1
+f = 0.1
+D_star = 30e-3
+D = 1e-3
+rescale_units = False
+
+# Settings
+lower_bounds = (0, 5, 0)
+upper_bounds = (1, 100, 4)
+bounds_um = np.array((lower_bounds, upper_bounds))
+
+lower_bounds = (0, 0.005, 0)
+upper_bounds = (1, 0.1, 0.004)
+bounds_mm = (lower_bounds, upper_bounds)
+initial_guess_mm = (1, 0.2, 0.03, 0.001)
+
+# Create gtab containing b-value informations
+bvals = np.array([0, 20, 40, 60, 80, 100, 150, 200, 300, 400, 500, 600, 700, 800]).astype(float)
+
+# Signal
+signal = ivim_signal(bvals, S0, f, D_star, D)
+noised_signal = add_rician_noise(signal, 3)
+noised_signal /= noised_signal[0]
+
+noised_signal6 = add_rician_noise(signal, 6)
+noised_signal6 /= noised_signal6[0]
+
+
+
+blim = 200
+
+
+# IVIM_seg
+pars = IVIM_seg(signal, bvals, np.array([[0, 0, 0, 0],[3e-3, np.inf, 1, 1]]), blim, False)
+
+

src/standardized/OJ_GU_seg.py

@@ -0,0 +1,59 @@
+import numpy as np
+from src.wrappers.OsipiBase import OsipiBase
+from src.original.OJ_GU.IVIM_fitting import IVIM_seg
+
+
+class OJ_GU_seg(OsipiBase):
+    """
+    Segmented fitting algorithm by Oscar Jalnefjord, University of Gothenburg
+    """
+
+    # I'm thinking that we define default attributes for each submission like this
+    # And in __init__, we can call the OsipiBase control functions to check whether
+    # the user inputs fulfil the requirements
+
+    # Some basic stuff that identifies the algorithm
+    id_author = "Oscar Jalnefjord, GU"
+    id_algorithm_type = "Segmented fit"
+    id_return_parameters = "f, D*, D"
+    id_units = "seconds per milli metre squared"
+
+    # Algorithm requirements
+    required_bvalues = 4
+    required_thresholds = [0,0] # Interval from "at least" to "at most", in case submissions allow a custom number of thresholds
+    required_bounds = True
+    required_bounds_optional = True # Bounds may not be required but are optional
+    required_initial_guess = False
+    required_initial_guess_optional = False
+    accepted_dimensions = 1 # Not sure how to define this for the number of accepted dimensions. Perhaps like the thresholds, at least and at most?
+
+    def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None):
+        """
+            Everything this algorithm requires should be implemented here.
+            Number of segmentation thresholds, bounds, etc.
+            
+            Our OsipiBase object could contain functions that compare the inputs with
+            the requirements.
+        """
+        super(OJ_GU_seg, self).__init__(bvalues, thresholds, bounds)
+        if bounds is None:
+            self.bounds = np.array([[0, 0, 0, 0],[3e-3, np.inf, 1, 1]])
+
+    def ivim_fit(self, signals, bvalues=None, verbose=False):
+        """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_
+        """
+        if bvalues is None:
+            bvalues = self.bvalues
+        else:
+            bvalues = np.asarray(bvalues)
+
+        pars = IVIM_seg(signals, bvalues, self.bounds, self.thresholds[0], disp_prog=verbose)
+
+        return pars['f'], pars['Dstar'], pars['D']

tests/IVIMmodels/unit_tests/test_ivim_fit.py

@@ -97,6 +97,7 @@ def test_ivim_fit_saved(name, bvals, data, algorithm):
         #[f_fit, Dp_fit, D_fit] = fit.ivim_fit(signal, bvals)
         #npt.assert_allclose([data['f'], data['D']], [f_fit, D_fit], atol=tolerance)
         #npt.assert_allclose(data['Dp'], Dp_fit, atol=1e-1)  # go easy on the perfusion as it's a linear fake
+
     [f_fit, Dp_fit, D_fit] = fit.ivim_fit(signal, bvals)
     npt.assert_allclose([data['f'], data['D']], [f_fit, D_fit], atol=tolerance)
     npt.assert_allclose(data['Dp'], Dp_fit, atol=1e-1)  # go easy on the perfusion as it's a linear fake