All pipeline of DRO is finished

.gitignore

@@ -13,4 +13,4 @@ docs-mk/site
 docs-mk/docs/generated
 src/osipi/DRO/data
 src/osipi/DRO/__pycache__/
-src/osipi/DRO/ROI_saved/aifmask.npy
+src/osipi/DRO/ROI_saved/

.pre-commit-config.yaml

@@ -3,7 +3,6 @@ repos:
     rev: v0.4.8  # Use the latest version
     hooks:
       - id: ruff
-        args: [--fix]  # Optional: to enable lint fixes
       - id: ruff-format
   - repo: https://github.com/pycqa/flake8
     rev: 7.0.0
@@ -12,7 +11,7 @@ repos:
         args:
           - --exclude=venv,.git,__pycache__,.eggs,.mypy_cache,.pytest_cache
           - --max-line-length=100
-          - --ignore=E203
+          - --ignore=E203, W503
           - --per-file-ignores=__init__.py:F401
   - repo: https://github.com/pre-commit/pre-commit-hooks
     rev: v4.6.0

src/osipi/Conc2DROSignal.py

@@ -0,0 +1,48 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+@author: eveshalom
+"""
+
+import numpy as np
+
+
+def createR10_withref(S0ref, S0, Tr, fa, T1ref, datashape):
+    R10_ref = 1 / T1ref
+    ref_frac = (1 - np.cos(fa) * np.exp(-Tr * R10_ref)) / (1 - np.exp(-Tr * R10_ref))
+    R10 = (-1 / Tr) * np.log((S0 - (ref_frac * S0ref)) / ((S0 * np.cos(fa)) - (ref_frac * S0ref)))
+    R10 = np.tile(R10[:, :, :, np.newaxis], (1, 1, 1, datashape[-1]))
+    return R10
+
+
+def calcR1_R2(R10, R20st, r1, r2st, Ctiss):
+    R1 = R10 + r1 * Ctiss
+    R2st = R20st + r2st * Ctiss
+    return R1, R2st
+
+
+def Conc2Sig_Linear(Tr, Te, fa, R1, R2st, S, scale, scalevisit1):
+    dro_S = ((1 - np.exp(-Tr * R1)) / (1 - (np.cos(fa) * np.exp(-Tr * R1)))) * (
+        np.sin(fa) * np.exp(-Te * R2st)
+    )
+
+    if scale == 1:
+        ScaleConst = np.percentile(S, 98) / np.percentile(dro_S, 98)
+    elif scale == 2:
+        ScaleConst = scalevisit1
+
+    dro_S = dro_S * ScaleConst
+    return dro_S, ScaleConst
+
+
+def STDmap(S, t0):
+    stdev = np.std(S[:, :, :, 0:t0], axis=3)
+    return stdev
+
+
+def addnoise(a, std, Sexact, Nt):
+    from numpy.random import normal as gaussian
+
+    std = np.tile(std[:, :, :, np.newaxis], (1, 1, 1, Nt))
+    Snoise = abs(Sexact + (a * std * gaussian(0, 1, Sexact.shape)))
+    return Snoise

src/osipi/DRO/DICOM_processing.py

@@ -3,6 +3,8 @@
 import numpy as np
 import pydicom
 
+from osipi.DRO.filters_and_noise import median_filter
+
 
 def read_dicom_slices_as_signal(folder_path):
     """
@@ -35,7 +37,20 @@ def read_dicom_slices_as_signal(folder_path):
         ):  # Sort by acquisition time
             signal[:, :, z, t] = slice.pixel_array
 
-    return signal, slices[slice_location][0][1]
+    return signal, slices, slices[slice_location][0][1]
+
+
+def SignalEnhancementExtract(S, datashape, baselinepoints):
+    # Take baseline average
+    S0 = np.average(S[:, :, :, 0:baselinepoints], axis=3)  # Take baseline signal
+    E = np.zeros_like(S)
+
+    # Calcualte siganl enhancement
+    for i in range(0, datashape[-1]):
+        E[:, :, :, i] = S[:, :, :, i] - S0
+        E[:, :, :, i] = median_filter(E[:, :, :, i])  # Median filter size (3,3)
+
+    return E, S0, S
 
 
 def calculate_baseline(signal, baseline):

src/osipi/DRO/Display.py

@@ -29,7 +29,7 @@ def animate(z):
             ax.set_title(f"Slice: {z}, Time: {time_index}")
 
     anim = FuncAnimation(
-        fig=fig, func=animate, frames=frames, init_func=init, interval=100, blit=False
+        fig=fig, func=animate, frames=frames, init_func=init, interval=10000, blit=False
     )
     plt.show()
     return anim

src/osipi/DRO/Model.py

@@ -2,6 +2,75 @@
 from scipy.integrate import cumtrapz, trapz
 
 
+def modifiedToftsMurase(Cp, Ctiss, dt, datashape):
+    # Fit Modified Tofts (Linear from Murase, 2004)
+    # Cp = Ea/0.45, Ctis=E/0.45
+    # Matrix equation C=AB (same notation as Murase)
+    # C: matrix of Ctis at distinct time steps
+    # A: 3 Coumns, rows of tk:
+    #   (1) Integral up to tk of Cp
+    #   (2) - Integral up to tk of Ctiss
+    #   (3) Cp at tk
+    # B: Array length 3 of parameters:
+    #   (1) K1 + k2 dot Vp
+    #   (2) k2
+    #   (3) Vp
+    # Use np.linalg.solve for equations form Zx=y aimed to find x
+    # np.linalg.solve(Z,y)=x  so need to use np.linalg.solve(A,C)
+    # solve only works for square matrices so use .lstsq for a least squares solve
+    #  Allocate parameter holding arrays
+
+    K1 = np.zeros(Ctiss.shape[:-1])  # only spatial maps
+    k2 = np.zeros(Ctiss.shape[:-1])  # only spatial maps
+    Vp = np.zeros(Ctiss.shape[:-1])  # only spatial maps
+
+    # Allocate matrices used from solver as defined above
+    C = np.zeros(datashape[-1])
+    A = np.zeros((datashape[-1], 3))
+
+    # iterate over slices
+    for k in range(0, datashape[2]):
+        # iterate over rows
+        for j in range(0, datashape[0]):
+            # iterate over columns
+            for i in range(0, datashape[1]):
+                # Build matrices for Modified Tofts for voxel
+                C = Ctiss[j, i, k, :]
+                A[:, 0] = cumtrapz(Cp, dx=dt, initial=0)
+                A[:, 1] = -cumtrapz(Ctiss[j, i, k, :], dx=dt, initial=0)
+                A[:, 2] = Cp
+                # Use least squares solver
+                sing_B1, sing_k2, sing_Vp = np.linalg.lstsq(A, C, rcond=None)[0]
+                sing_K1 = sing_B1 - (sing_k2 * sing_Vp)
+                # Assign Ouputs into parameter maps
+                K1[j, i, k] = sing_K1
+                k2[j, i, k] = sing_k2
+                Vp[j, i, k] = sing_Vp
+
+    return K1, k2, Vp
+
+
+def modifiedToftsMurase1Vox(Cp, Ctiss, dt, datashape):
+    K1 = np.zeros(Ctiss.shape[:-1])  # only spatial maps
+    k2 = np.zeros(Ctiss.shape[:-1])  # only spatial maps
+    Vp = np.zeros(Ctiss.shape[:-1])  # only spatial maps
+
+    # Allocate matrices used from solver as defined above
+    C = np.zeros(datashape[-1])
+    A = np.zeros((datashape[-1], 3))
+
+    # Build matrices for Modified Tofts for voxel
+    C = Ctiss
+    A[:, 0] = cumtrapz(Cp, dx=dt, initial=0)
+    A[:, 1] = -cumtrapz(Ctiss, dx=dt, initial=0)
+    A[:, 2] = Cp
+    # Use least squares solver
+    B1, k2, Vp = np.linalg.lstsq(A, C, rcond=None)[0]
+    K1 = B1 - (k2 * Vp)
+
+    return K1, k2, Vp
+
+
 def tofts(cp, c_tiss, dt, datashape):
     """
     Tofts model for DCE-MRI DRO
@@ -75,3 +144,30 @@ def forward_tofts(ktrans, kep, cp, vp, dt):
             c_tiss[..., t] = vp * cp[..., t] + ktrans * integral
 
     return c_tiss
+
+
+def ForwardsModTofts(K1, k2, Vp, Cp, dt):
+    # To be carried out as matmul C=BA
+    # Where C is the output Ctiss and B the parameters
+    # With A a matrix of cumulative integrals
+
+    x, y, z = K1.shape
+    t = Cp.shape[0]
+
+    Ctiss = np.zeros((y, x, z, t))
+
+    b1 = K1 + np.multiply(k2, Vp)  # define combined parameter
+    B = np.zeros((x, y, z, 1, 3))
+    A = np.zeros((x, y, z, 3, 1))
+
+    B[:, :, :, 0, 0] = b1
+    B[:, :, :, 0, 1] = -k2
+    B[:, :, :, 0, 2] = Vp
+
+    for tk in range(1, t):
+        A[:, :, :, 0, 0] = trapz(Cp[0 : tk + 1], dx=dt)
+        A[:, :, :, 1, 0] = trapz(Ctiss[:, :, :, 0 : tk + 1], dx=dt)
+        A[:, :, :, 2, 0] = Cp[tk]
+
+        Ctiss[:, :, :, tk] = np.matmul(B, A).squeeze()
+    return Ctiss

src/osipi/DRO/filters_and_noise.py

@@ -2,11 +2,13 @@
 from scipy import ndimage
 
 
-def median_filter(signal):
+def median_filter(param_map):
     """
     Apply a median filter to a signal.
     """
-    return ndimage.median_filter(signal, size=(3, 3, 1, 1))
+    for i in range(param_map.shape[-1]):
+        param_map[:, :, i] = ndimage.median_filter(param_map[:, :, i], size=(3, 3))
+    return param_map
 
 
 def add_gaussian_noise(signal, mean=0, std=1):