DCE-dro

.gitignore

@@ -11,3 +11,6 @@ docs/source/generated
 docs/source/sg_execution_times.rst
 docs-mk/site
 docs-mk/docs/generated
+src/osipi/DRO/data
+src/osipi/DRO/__pycache__/
+src/osipi/DRO/ROI_saved/

.pre-commit-config.yaml

@@ -3,7 +3,7 @@ repos:
     rev: v0.4.8  # Use the latest version
     hooks:
       - id: ruff
-        args: [--fix]  # Optional: to enable lint fixes
+        args: ["--fix"]
       - id: ruff-format
   - repo: https://github.com/pycqa/flake8
     rev: 7.0.0
@@ -12,7 +12,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

pyproject.toml

@@ -114,6 +114,7 @@ numpy = "^1.21.2"
 scipy = "^1.7.3"
 matplotlib = "^3.4.3"
 requests = "^2.32.3"
+pydicom = "^2.4.4"
 
 [tool.poetry.group.dev.dependencies]
 flake8 = "^7.0.0"

src/osipi/DRO/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

@@ -0,0 +1,116 @@
+import os
+
+import numpy as np
+import pydicom
+
+from osipi.DRO.filters_and_noise import median_filter
+
+
+def read_dicom_slices_as_4d_signal(folder_path):
+    """
+    Read a DICOM series from a folder path.
+    Returns the signal data as a 4D numpy array (x, y, z, t).
+    """
+    slices = {}
+    for root, _, files in os.walk(folder_path):
+        for file in files:
+            if file.endswith(".dcm"):
+                dicom_file = os.path.join(root, file)
+                slice = pydicom.read_file(dicom_file)
+                if slice.SliceLocation not in slices:
+                    slices[slice.SliceLocation] = []
+                slices[slice.SliceLocation].append((slice.AcquisitionTime, slice))
+
+    # Sort each list of slices by the first element (AcquisitionTime)
+    for slice_location in slices:
+        slices[slice_location].sort(key=lambda x: x[0])
+
+    spatial_shape = slices[slice_location][0][1].pixel_array.shape
+
+    data_shape = (spatial_shape[0], spatial_shape[1], len(slices), len(slices[slice_location]))
+
+    signal = np.zeros(data_shape)
+
+    for z, slice_location in enumerate(sorted(slices.keys())):  # Sort by slice location
+        for t, (_, slice) in enumerate(
+            sorted(slices[slice_location], key=lambda x: x[0])
+        ):  # Sort by acquisition time
+            signal[:, :, z, t] = slice.pixel_array
+
+    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):
+    """
+    Calculate the baseline signal (S0) from pre-contrast time points.
+
+    Parameters:
+    signal (numpy.ndarray): The 4D signal data (x, y, z, t).
+    baseline (int): Number of time points before contrast injection.
+
+    Returns:
+    numpy.ndarray: The baseline signal (S0).
+    """
+    S0 = np.average(signal[:, :, :, :baseline], axis=3, keepdims=True)
+    return S0
+
+
+def signal_to_R1(signal, S0, TR):
+    """
+    Convert signal to R1 values using the Ernst equation.
+
+    Parameters:
+    signal (numpy.ndarray): The 4D signal data (x, y, z, t).
+    S0 (numpy.ndarray): The baseline signal (S0).
+    TR (float): Repetition time.
+
+    Returns:
+    numpy.ndarray: The R1 values.
+    """
+    epsilon = 1e-8  # Small constant to avoid division by zero and log of zero
+    R1 = -1 / TR * np.log((signal + epsilon) / (S0 + epsilon))
+    return R1
+
+
+def calc_concentration(R1, R10, r1):
+    """
+    Calculate the concentration of the contrast agent in tissue.
+
+    Parameters:
+    R1 (numpy.ndarray): The R1 values.
+    R10 (numpy.ndarray): The pre-contrast R1 values.
+    r1 (float): Relaxivity of the contrast agent.
+
+    Returns:
+    numpy.ndarray: The concentration of the contrast agent in the tissue.
+    """
+    Ctiss = (R1 - R10) / r1
+    return Ctiss
+
+
+def signal_enhancement(signal, S0, R10, r1):
+    """
+    Calculate the signal enhancement.
+
+    Parameters:
+    signal (numpy.ndarray): The 4D signal data (x, y, z, t).
+    other parameters same as previous function
+
+    Returns:
+    numpy.ndarray: The signal enhancement.
+    """
+    E = (R10 / r1) * (signal - S0) / S0
+    return E

src/osipi/DRO/Display.py

@@ -0,0 +1,35 @@
+from matplotlib import pyplot as plt
+from matplotlib.animation import FuncAnimation
+
+
+def animate_mri(slices, mode="time", slice_index=0, time_index=0):
+    fig, ax = plt.subplots()
+    if mode == "time":
+        frames = slices.shape[-1]
+
+        def init():
+            ax.imshow(slices[:, :, slice_index, 0], cmap="gray")
+            ax.set_title(f"Slice: {slice_index}, Time: 0")
+
+        def animate(t):
+            ax.clear()
+            ax.imshow(slices[:, :, slice_index, t], cmap="gray")
+            ax.set_title(f"Slice: {slice_index}, Time: {t}")
+
+    elif mode == "slice":
+        frames = slices.shape[2]
+
+        def init():
+            ax.imshow(slices[:, :, 0, time_index], cmap="gray")
+            ax.set_title(f"Slice: 0, Time: {time_index}")
+
+        def animate(z):
+            ax.clear()
+            ax.imshow(slices[:, :, z, time_index], cmap="gray")
+            ax.set_title(f"Slice: {z}, Time: {time_index}")
+
+    anim = FuncAnimation(
+        fig=fig, func=animate, frames=frames, init_func=init, interval=100, blit=False
+    )
+    plt.show()
+    return anim