From c79ff0c95f052b8e30cbd450fd6a86b9bdab89fb Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Thu, 18 Jul 2024 15:20:56 +0300 Subject: [PATCH 01/12] Begin DRO basic functions --- poetry.lock | 16 +++++++++- pyproject.toml | 1 + src/osipi/DRO/DICOM_processing.py | 45 +++++++++++++++++++++++++++ src/osipi/DRO/Model.py | 50 ++++++++++++++++++++++++++++++ src/osipi/DRO/__init__.py | 0 src/osipi/DRO/filters_and_noise.py | 26 ++++++++++++++++ src/osipi/DRO/roi_selection.py | 48 ++++++++++++++++++++++++++++ 7 files changed, 185 insertions(+), 1 deletion(-) create mode 100644 src/osipi/DRO/DICOM_processing.py create mode 100644 src/osipi/DRO/Model.py create mode 100644 src/osipi/DRO/__init__.py create mode 100644 src/osipi/DRO/filters_and_noise.py create mode 100644 src/osipi/DRO/roi_selection.py diff --git a/poetry.lock b/poetry.lock index 99c34be..d7e5951 100644 --- a/poetry.lock +++ b/poetry.lock @@ -1754,6 +1754,20 @@ doc = ["ablog (>=0.11.8)", "colorama", "graphviz", "ipykernel", "ipyleaflet", "i i18n = ["Babel", "jinja2"] test = ["pytest", "pytest-cov", "pytest-regressions", "sphinx[test]"] +[[package]] +name = "pydicom" +version = "2.4.4" +description = "A pure Python package for reading and writing DICOM data" +optional = false +python-versions = ">=3.7" +files = [ + {file = "pydicom-2.4.4-py3-none-any.whl", hash = "sha256:f9f8e19b78525be57aa6384484298833e4d06ac1d6226c79459131ddb0bd7c42"}, + {file = "pydicom-2.4.4.tar.gz", hash = "sha256:90b4801d851ce65be3df520e16bbfa3d6c767cf2a3a3b1c18f6780e6b670b87a"}, +] + +[package.extras] +docs = ["matplotlib", "numpy", "numpydoc", "pillow", "sphinx", "sphinx-copybutton", "sphinx-gallery", "sphinx_rtd_theme", "sphinxcontrib-napoleon"] + [[package]] name = "pyflakes" version = "3.2.0" @@ -2599,4 +2613,4 @@ test = ["big-O", "importlib-resources", "jaraco.functools", "jaraco.itertools", [metadata] lock-version = "2.0" python-versions = "^3.9" -content-hash = "6b6056507de0b7072e84eff4a02da7c1ac369939f4f234c861fa44df3b3c18df" +content-hash = "0665b4c4d2ef51dc8d6b4144a4953cc11fe83d3c3ac3f8a0e748d9202ed0a4bc" diff --git a/pyproject.toml b/pyproject.toml index 7ce1c7b..5357a3e 100644 --- a/pyproject.toml +++ b/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" diff --git a/src/osipi/DRO/DICOM_processing.py b/src/osipi/DRO/DICOM_processing.py new file mode 100644 index 0000000..7cefc0a --- /dev/null +++ b/src/osipi/DRO/DICOM_processing.py @@ -0,0 +1,45 @@ +import os + +import numpy as np +import pydicom + + +def read_dicom_slices_as_signal(folder_path): + """ + Read a DICOM series from a folder path. + """ + 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, spatial_shape, 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 + + +def signal_enhanecment(signal, baseline=0): + """ + Calculate signal enhancement. + """ + s0 = np.average(signal[:, :, :, :baseline], axis=-1) + return (signal - s0) / s0, s0 diff --git a/src/osipi/DRO/Model.py b/src/osipi/DRO/Model.py new file mode 100644 index 0000000..ca625c2 --- /dev/null +++ b/src/osipi/DRO/Model.py @@ -0,0 +1,50 @@ +import numpy as np +from scipy.integrate import cumtrapz + + +def tofts(cp, c_tiss, dt, datashape): + """ + Tofts model for DCE-MRI DRO + + ref : https://github.com/JCardenasRdz/Gage-repeatability-DCE-MRI?tab=readme-ov-file + + """ + ktrans = np.zeros(c_tiss.shape[:-1]) + kep = np.zeros(c_tiss.shape[:-1]) + + for k in range(0, datashape[2]): + for j in range(0, datashape[0]): + for i in range(0, datashape[1]): + c = c_tiss[j, i, k, :] + cp_integrated = cumtrapz(cp, dx=dt, initial=0) + c_tiss_integrated = -cumtrapz(c_tiss[i, j, k, :], dx=dt, initial=0) + A = np.column_stack((cp_integrated, c_tiss_integrated)) + ktrans_voxel, kep_voxel = np.linalg.lstsq(A, c, rcond=None)[0] + ktrans[j, i, k] = ktrans_voxel + kep[j, i, k] = kep_voxel + + return ktrans, kep + + +def exntended_tofts(cp, c_tiss, dt, datashape): + """ + Extended Tofts model for DCE-MRI DRO + + """ + ktrans = np.zeros(c_tiss.shape[:-1]) + kep = np.zeros(c_tiss.shape[:-1]) + vp = np.zeros(c_tiss.shape[:-1]) + + for k in range(0, datashape[2]): + for j in range(0, datashape[0]): + for i in range(0, datashape[1]): + c = c_tiss[j, i, k, :] + cp_integrated = cumtrapz(cp, dx=dt, initial=0) + c_tiss_integrated = -cumtrapz(c_tiss[i, j, k, :], dx=dt, initial=0) + A = np.column_stack((cp_integrated, c_tiss_integrated, cp)) + ktrans_voxel, kep_voxel, vp_voxel = np.linalg.lstsq(A, c, rcond=None)[0] + ktrans[j, i, k] = ktrans_voxel + kep[j, i, k] = kep_voxel + vp[j, i, k] = vp_voxel + + return ktrans, kep, vp diff --git a/src/osipi/DRO/__init__.py b/src/osipi/DRO/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/osipi/DRO/filters_and_noise.py b/src/osipi/DRO/filters_and_noise.py new file mode 100644 index 0000000..21cee3b --- /dev/null +++ b/src/osipi/DRO/filters_and_noise.py @@ -0,0 +1,26 @@ +import numpy as np +from scipy import ndimage + + +def median_filter(signal): + """ + Apply a median filter to a signal. + """ + return ndimage.median_filter(signal, size=(3, 3, 1, 1)) + + +def add_gaussian_noise(signal, mean=0, std=1): + """ + Add Gaussian noise to the 4 MR data. + """ + return signal + np.random.normal(mean, std, signal.shape) + + +def add_rician_noise(signal, mean=0, std=1): + """ + Add Rician noise to the 4D MR data. + """ + noise_real = np.random.normal(0, std, signal.shape) + noise_imag = np.random.normal(0, std, signal.shape) + noisy_signal = np.sqrt(signal**2 + noise_real**2 + noise_imag**2) + return noisy_signal diff --git a/src/osipi/DRO/roi_selection.py b/src/osipi/DRO/roi_selection.py new file mode 100644 index 0000000..0a5a301 --- /dev/null +++ b/src/osipi/DRO/roi_selection.py @@ -0,0 +1,48 @@ +import numpy as np +from matplotlib import path +from matplotlib import pyplot as plt +from matplotlib.widgets import LassoSelector + + +def roi(signal, slice_num): + """ + Select a region of interest (ROI) in a slice of a 4D signal. + """ + + fig = plt.figure() + ax1 = fig.add_subplot(121) + ax1.set_title("Slice:") + ax1.imshow(signal[:, :, slice_num, 0], cmap="gray", interpolation="nearest") + + # Empty array to be filled with lasso selector + array = np.zeros_like(signal[:, :, slice_num, 0]) + ax2 = fig.add_subplot(122) + ax2.set_title("numpy array:") + mask = ax2.imshow(array, vmax=1, interpolation="nearest") + + plt.show() + # Pixel coordinates + pix = np.arange(signal.shape[1]) + xv, yv = np.meshgrid(pix, pix) + pix = np.vstack((xv.flatten(), yv.flatten())).T + + def updateArray(array, indices): + lin = np.arange(array.size) + newArray = array.flatten() + newArray[lin[indices]] = 1 + return newArray.reshape(array.shape) + + def onselect(verts): + array = mask._A._data + p = path.Path(verts) + ind = p.contains_points(pix, radius=1) + array = updateArray(array, ind) + mask.set_data(array) + fig.canvas.draw_idle() + + lineprops = {"color": "red", "linewidth": 4, "alpha": 0.8} + lasso = LassoSelector(ax1, onselect, lineprops, button=1) + mask_4d = np.zeros(signal.shape) + mask_4d[:, :, slice_num, :] = mask + + return lasso, mask_4d From d653e40c1bdb4ee80daa86a48f6743ac8c6bb52d Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Fri, 19 Jul 2024 21:45:08 +0300 Subject: [PATCH 02/12] Create signal to concentration methods --- src/osipi/DRO/DICOM_processing.py | 63 +++++++++++++++++++++++++++++-- 1 file changed, 59 insertions(+), 4 deletions(-) diff --git a/src/osipi/DRO/DICOM_processing.py b/src/osipi/DRO/DICOM_processing.py index 7cefc0a..070b5c4 100644 --- a/src/osipi/DRO/DICOM_processing.py +++ b/src/osipi/DRO/DICOM_processing.py @@ -7,6 +7,7 @@ def read_dicom_slices_as_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): @@ -37,9 +38,63 @@ def read_dicom_slices_as_signal(folder_path): return signal -def signal_enhanecment(signal, baseline=0): +def calculate_baseline(signal, baseline): """ - Calculate signal enhancement. + 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=-1) + 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. + """ + R1 = -1 / TR * np.log(signal / S0) + 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. """ - s0 = np.average(signal[:, :, :, :baseline], axis=-1) - return (signal - s0) / s0, s0 + E = (R10 / r1) * (signal - S0) / S0 + return E From eb9b39f6ae79fcff925250780eb4cfe973f7d01f Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Fri, 19 Jul 2024 21:47:02 +0300 Subject: [PATCH 03/12] Create forward tofts function --- src/osipi/DRO/Model.py | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/src/osipi/DRO/Model.py b/src/osipi/DRO/Model.py index ca625c2..17b718b 100644 --- a/src/osipi/DRO/Model.py +++ b/src/osipi/DRO/Model.py @@ -1,5 +1,5 @@ import numpy as np -from scipy.integrate import cumtrapz +from scipy.integrate import cumtrapz, trapz def tofts(cp, c_tiss, dt, datashape): @@ -48,3 +48,30 @@ def exntended_tofts(cp, c_tiss, dt, datashape): vp[j, i, k] = vp_voxel return ktrans, kep, vp + + +def forward_tofts(ktrans, kep, cp, vp, dt): + """ + Forward Tofts model for DCE-MRI DRO + + Parameters: + ktrans (numpy.ndarray): The transfer constant Ktrans. + kep (numpy.ndarray): The rate constant kep. + cp (numpy.ndarray): The plasma concentration C_p(t). + vp (numpy.ndarray): The plasma volume fraction v_p. + dt (float): The time step between measurements. + + Returns: + numpy.ndarray: The tissue concentration C_tiss(t). + """ + time_points = cp.shape[-1] + c_tiss = np.zeros(ktrans.shape) + for t in range(time_points): + if t == 0: + c_tiss[..., t] = vp * cp[..., t] + else: + exp = np.exp(-kep * np.arange(t + 1)[::-1] * dt) + integral = trapz(cp[..., : t + 1] * exp, dx=dt, axis=-1) + c_tiss[..., t] = vp * cp[..., t] + ktrans * integral + + return c_tiss From 24afe12b146e3441d73e00910ee4231428b278a8 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Sun, 21 Jul 2024 17:32:36 +0300 Subject: [PATCH 04/12] display slices and select specific region --- .gitignore | 3 +++ src/osipi/DRO/DICOM_processing.py | 9 ++++---- src/osipi/DRO/Display.py | 35 +++++++++++++++++++++++++++++++ src/osipi/DRO/main.py | 25 ++++++++++++++++++++++ src/osipi/DRO/roi_selection.py | 24 ++++++++++++--------- 5 files changed, 82 insertions(+), 14 deletions(-) create mode 100644 src/osipi/DRO/Display.py create mode 100644 src/osipi/DRO/main.py diff --git a/.gitignore b/.gitignore index 701f9ad..4f359b8 100644 --- a/.gitignore +++ b/.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/aifmask.npy diff --git a/src/osipi/DRO/DICOM_processing.py b/src/osipi/DRO/DICOM_processing.py index 070b5c4..6ca1e3f 100644 --- a/src/osipi/DRO/DICOM_processing.py +++ b/src/osipi/DRO/DICOM_processing.py @@ -25,7 +25,7 @@ def read_dicom_slices_as_signal(folder_path): spatial_shape = slices[slice_location][0][1].pixel_array.shape - data_shape = (spatial_shape, spatial_shape, len(slices), len(slices[slice_location])) + data_shape = (spatial_shape[0], spatial_shape[1], len(slices), len(slices[slice_location])) signal = np.zeros(data_shape) @@ -35,7 +35,7 @@ def read_dicom_slices_as_signal(folder_path): ): # Sort by acquisition time signal[:, :, z, t] = slice.pixel_array - return signal + return signal, slices[slice_location][0][1] def calculate_baseline(signal, baseline): @@ -49,7 +49,7 @@ def calculate_baseline(signal, baseline): Returns: numpy.ndarray: The baseline signal (S0). """ - S0 = np.average(signal[..., :baseline], axis=-1) + S0 = np.average(signal[:, :, :, :baseline], axis=3, keepdims=True) return S0 @@ -65,7 +65,8 @@ def signal_to_R1(signal, S0, TR): Returns: numpy.ndarray: The R1 values. """ - R1 = -1 / TR * np.log(signal / S0) + 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 diff --git a/src/osipi/DRO/Display.py b/src/osipi/DRO/Display.py new file mode 100644 index 0000000..be815af --- /dev/null +++ b/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 diff --git a/src/osipi/DRO/main.py b/src/osipi/DRO/main.py new file mode 100644 index 0000000..fb0fd23 --- /dev/null +++ b/src/osipi/DRO/main.py @@ -0,0 +1,25 @@ +import numpy as np +from DICOM_processing import ( + calc_concentration, + calculate_baseline, + read_dicom_slices_as_signal, + signal_to_R1, +) +from Display import animate_mri +from roi_selection import roi + +slices, dicom_ref = read_dicom_slices_as_signal("data/subject1/13.000000-perfusion-23726") +anim = animate_mri(slices, mode="time", slice_index=7, time_index=5) +TR = 1e-3 * dicom_ref.RepetitionTime +TE = 1e-3 * dicom_ref.EchoTime +flip_angle = dicom_ref.FlipAngle +baseline = 5 +R10 = 1.18 +r1 = 3.9 + +baseline_signal = calculate_baseline(slices, baseline) +R1 = signal_to_R1(slices, baseline_signal, TR) +c_tissue = calc_concentration(R1, R10, r1) + +Max = np.max(c_tissue, axis=-1, keepdims=True) +mask, roi_sum, lasso = roi(Max, 5, c_tissue.shape, save=True) diff --git a/src/osipi/DRO/roi_selection.py b/src/osipi/DRO/roi_selection.py index 0a5a301..3752929 100644 --- a/src/osipi/DRO/roi_selection.py +++ b/src/osipi/DRO/roi_selection.py @@ -4,7 +4,7 @@ from matplotlib.widgets import LassoSelector -def roi(signal, slice_num): +def roi(signal, slice_num, data_shape, save=False): """ Select a region of interest (ROI) in a slice of a 4D signal. """ @@ -12,15 +12,14 @@ def roi(signal, slice_num): fig = plt.figure() ax1 = fig.add_subplot(121) ax1.set_title("Slice:") - ax1.imshow(signal[:, :, slice_num, 0], cmap="gray", interpolation="nearest") + ax1.imshow(signal[:, :, slice_num], cmap="gray", interpolation="nearest") # Empty array to be filled with lasso selector - array = np.zeros_like(signal[:, :, slice_num, 0]) + array = np.zeros_like(signal[:, :, slice_num]) ax2 = fig.add_subplot(122) ax2.set_title("numpy array:") mask = ax2.imshow(array, vmax=1, interpolation="nearest") - plt.show() # Pixel coordinates pix = np.arange(signal.shape[1]) xv, yv = np.meshgrid(pix, pix) @@ -40,9 +39,14 @@ def onselect(verts): mask.set_data(array) fig.canvas.draw_idle() - lineprops = {"color": "red", "linewidth": 4, "alpha": 0.8} - lasso = LassoSelector(ax1, onselect, lineprops, button=1) - mask_4d = np.zeros(signal.shape) - mask_4d[:, :, slice_num, :] = mask - - return lasso, mask_4d + lasso = LassoSelector(ax1, onselect, button=1) + plt.show() + mask2d = mask._A._data.copy() + roivox = np.sum(mask2d) + timemask = np.tile(mask2d[:, :, np.newaxis], (1, 1, signal.shape[-1])) + mask4d = np.zeros(data_shape) + mask4d[:, :, slice_num, :] = timemask + if save: + np.save("mask/roi_mask.npy", mask4d) + np.save("mask/roi_voxels.npy", roivox) + return mask4d, roivox, lasso From 327c892906d20203c3c1bef8b6046ec18f2d1f69 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Thu, 25 Jul 2024 21:33:14 +0300 Subject: [PATCH 05/12] All pipeline of DRO is finished --- .gitignore | 2 +- .pre-commit-config.yaml | 3 +- src/osipi/Conc2DROSignal.py | 48 ++++++++ src/osipi/DRO/DICOM_processing.py | 17 ++- src/osipi/DRO/Display.py | 2 +- src/osipi/DRO/Model.py | 96 ++++++++++++++++ src/osipi/DRO/filters_and_noise.py | 6 +- src/osipi/DRO/main.py | 173 +++++++++++++++++++++++++---- src/osipi/DRO/roi_selection.py | 12 +- 9 files changed, 330 insertions(+), 29 deletions(-) create mode 100644 src/osipi/Conc2DROSignal.py diff --git a/.gitignore b/.gitignore index 4f359b8..86b9bcf 100644 --- a/.gitignore +++ b/.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/ diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index c2429e4..2dd8c9b 100644 --- a/.pre-commit-config.yaml +++ b/.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 diff --git a/src/osipi/Conc2DROSignal.py b/src/osipi/Conc2DROSignal.py new file mode 100644 index 0000000..ad17259 --- /dev/null +++ b/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 diff --git a/src/osipi/DRO/DICOM_processing.py b/src/osipi/DRO/DICOM_processing.py index 6ca1e3f..86d3fea 100644 --- a/src/osipi/DRO/DICOM_processing.py +++ b/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): diff --git a/src/osipi/DRO/Display.py b/src/osipi/DRO/Display.py index be815af..7e16e18 100644 --- a/src/osipi/DRO/Display.py +++ b/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 diff --git a/src/osipi/DRO/Model.py b/src/osipi/DRO/Model.py index 17b718b..d32e384 100644 --- a/src/osipi/DRO/Model.py +++ b/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 diff --git a/src/osipi/DRO/filters_and_noise.py b/src/osipi/DRO/filters_and_noise.py index 21cee3b..1111cc0 100644 --- a/src/osipi/DRO/filters_and_noise.py +++ b/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): diff --git a/src/osipi/DRO/main.py b/src/osipi/DRO/main.py index fb0fd23..c66ad36 100644 --- a/src/osipi/DRO/main.py +++ b/src/osipi/DRO/main.py @@ -1,25 +1,158 @@ import numpy as np from DICOM_processing import ( - calc_concentration, - calculate_baseline, + SignalEnhancementExtract, read_dicom_slices_as_signal, - signal_to_R1, ) from Display import animate_mri -from roi_selection import roi - -slices, dicom_ref = read_dicom_slices_as_signal("data/subject1/13.000000-perfusion-23726") -anim = animate_mri(slices, mode="time", slice_index=7, time_index=5) -TR = 1e-3 * dicom_ref.RepetitionTime -TE = 1e-3 * dicom_ref.EchoTime -flip_angle = dicom_ref.FlipAngle -baseline = 5 -R10 = 1.18 -r1 = 3.9 - -baseline_signal = calculate_baseline(slices, baseline) -R1 = signal_to_R1(slices, baseline_signal, TR) -c_tissue = calc_concentration(R1, R10, r1) - -Max = np.max(c_tissue, axis=-1, keepdims=True) -mask, roi_sum, lasso = roi(Max, 5, c_tissue.shape, save=True) +from roi_selection import ICfromROI + +from osipi.Conc2DROSignal import Conc2Sig_Linear, STDmap, addnoise, calcR1_R2, createR10_withref +from osipi.DRO.filters_and_noise import median_filter +from osipi.DRO.Model import ForwardsModTofts, modifiedToftsMurase, modifiedToftsMurase1Vox + +signal, slices, dicom_ref = read_dicom_slices_as_signal("data/subject2/12.000000-perfusion-17557") +# anim = animate_mri(slices, mode="time", slice_index=7, time_index=5) +data_shape = signal.shape + +E, S0, S = SignalEnhancementExtract(signal, data_shape, 5) + +Max = np.max(E, axis=-1, keepdims=True) + +dt = 4.8 / 60 # mins from the RIDER website DCE description +t = np.linspace(0, dt * S.shape[-1], S.shape[-1]) # time series points + +aif_dir = "ROI_saved/" + +aifmask = np.load("{}aifmask.npy".format(aif_dir)) +sagmask = np.load("{}sagmask.npy".format(aif_dir)) +roivoxa = np.load("{}aifmaskvox.npy".format(aif_dir)) +roivoxv = np.load("{}sagmaskvox.npy".format(aif_dir)) + +z = 5 + +Ea = ICfromROI(E, aifmask, roivoxa, aifmask.ndim - 1) +Ev = ICfromROI(E, sagmask, roivoxv, sagmask.ndim - 1) +S0ref = ICfromROI(S0[:, :, z], sagmask[:, :, z, 0], roivoxv, sagmask.ndim - 2) + +Hct = 0.45 + +K1, k2, Vp = modifiedToftsMurase(Ea / (1 - Hct), E, dt, data_shape) + +pvc_K1, pvc_k2, pvc_Vp = modifiedToftsMurase1Vox(Ea / (1 - Hct), Ev, dt, data_shape) +pvc = abs((1 - Hct) / pvc_Vp) # sagittal sinus Vp should be 0.55, so calc correction factor if not + +# Apply correction factor to fitted parameters +cor_K1 = K1 * pvc +cor_k2 = k2 * pvc +cor_Vp = Vp * pvc +# Apply Median Filter to parameters all with footprint (3,3) + +mf_K1 = median_filter(cor_K1) +mf_k2 = median_filter(cor_k2) +mf_Vp = median_filter(cor_Vp) + +mf_Vp[mf_Vp <= 0] = 0 +mf_Vp[mf_Vp > 1] = 1.0 +mf_K1[mf_K1 <= 0] = 0 +mf_k2[mf_k2 <= 0] = 0 + +# evolve forwards model +aif_cor = np.max((Ea / (1 - Hct)) * pvc) / 6 # caluclate enhancement to concentration correction +Cp = ( + (Ea / (1 - Hct)) * pvc +) / aif_cor # calculate Cp using signal enhancement aif and concentration conversion factor + +c_tissue = ForwardsModTofts(mf_K1, mf_k2, mf_Vp, Cp, dt) + +r1 = 3.9 # longitudinal relaxivity Gd-DTPA (Hz/mM) source: (Pintaske,2006) +r2st = 10 # transverse relaxivity Gd-DTPA (Hz/mM) +# roughly estimated using (Pintaske,2006) and (Siemonsen, 2008) + +fa = np.deg2rad(float(dicom_ref.FlipAngle)) # flip angle (rads) + +Te = 1e-3 * float(dicom_ref.EchoTime) # Echo Time (s) +Tr = 1e-3 * float(dicom_ref.RepetitionTime) # Repetition Time (s) +T1b = 1.48 # T1 for blood measured in sagittal sinus @ 1.5T (s) (Zhang et al., 2013) + +R10_value = 1.18 # precontrast T1 relaxation rate (Hz) brain avg (radiopedia) +R20st_value = ( + 17.24 # precontrast T1 relaxation rate (Hz) brain avg using T2* from (Siemonsen, 2008) +) +R10 = createR10_withref( + S0ref, S0, Tr, fa, T1b, data_shape +) # precontrast R1 map (normalised to sagittal sinus) + +R1, R2st = calcR1_R2(R10, R20st_value, r1, r2st, c_tissue) # returns R10 and R2st maps +dro_S, M = Conc2Sig_Linear(Tr, Te, fa, R1, R2st, S, 1, 0) + +stdS = STDmap(signal, t0=5) # caluclate Standard deviation for original data +dro_Snoise = addnoise(1, stdS, dro_S, Nt=data_shape[-1]) + +trans_K1 = mf_K1.copy() +trans_k2 = mf_k2.copy() +trans_Vp = mf_Vp.copy() + +vmax_K1, vmax_k2, vmax_Vp = 1, 1, 0.2 +vmin_K1, vmin_k2, vmin_Vp = 0.2, 0.2, 0.01 +lb_K1, lb_k2, lb_Vp = 0.54, 0.52, 0.49 +ub_K1, ub_k2, ub_Vp = 1.52, 1.5, 1.43 +lim_K1, lim_k2, lim_Vp = vmax_K1 + 0.5, vmax_k2 + 0.1, vmax_Vp + 0.5 +ub_lim = 1.01 + +trans_K1[trans_K1 <= vmin_K1] = trans_K1[trans_K1 <= vmin_K1] * lb_K1 +trans_K1[trans_K1 >= lim_K1] = trans_K1[trans_K1 >= lim_K1] * ub_lim +trans_K1[(trans_K1 >= vmax_K1) & (trans_K1 < lim_K1)] = ( + trans_K1[(trans_K1 >= vmax_K1) & (trans_K1 < lim_K1)] * ub_K1 +) +trans_K1[(trans_K1 > vmin_K1) & (trans_K1 < vmax_K1)] = trans_K1[ + (trans_K1 > vmin_K1) & (trans_K1 < vmax_K1) +] * ( + lb_K1 + + ( + ((trans_K1[(trans_K1 > vmin_K1) & (trans_K1 < vmax_K1)] - vmin_K1) / (vmax_K1 - vmin_K1)) + * (ub_K1 - lb_K1) + ) +) + +trans_k2[trans_k2 <= vmin_k2] = trans_k2[trans_k2 <= vmin_k2] * lb_k2 +trans_k2[trans_k2 >= lim_k2] = trans_k2[trans_k2 >= lim_k2] * ub_lim +trans_k2[(trans_k2 >= vmax_k2) & (trans_k2 < lim_k2)] = ( + trans_k2[(trans_k2 >= vmax_k2) & (trans_k2 < lim_k2)] * ub_k2 +) +trans_k2[(trans_k2 > vmin_k2) & (trans_k2 < vmax_k2)] = trans_k2[ + (trans_k2 > vmin_k2) & (trans_k2 < vmax_k2) +] * ( + lb_k2 + + ( + ((trans_k2[(trans_k2 > vmin_k2) & (trans_k2 < vmax_k2)] - vmin_k2) / (vmax_k2 - vmin_k2)) + * (ub_k2 - lb_k2) + ) +) + +trans_Vp[trans_Vp <= vmin_Vp] = trans_Vp[trans_Vp <= vmin_Vp] * lb_Vp +trans_Vp[(trans_Vp >= lim_Vp)] = trans_Vp[trans_Vp >= lim_Vp] * ub_lim +trans_Vp[(trans_Vp >= vmax_Vp) & (trans_Vp < lim_Vp)] = ( + trans_Vp[(trans_Vp >= vmax_Vp) & (trans_Vp < lim_Vp)] * ub_Vp +) +trans_Vp[(trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp)] = trans_Vp[ + (trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp) +] * ( + lb_Vp + + ( + ((trans_Vp[(trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp)] - vmin_Vp) / (vmax_Vp - vmin_Vp)) + * (ub_Vp - lb_Vp) + ) +) + +trans_Vp[trans_Vp > 1] = 1 + +Ctiss_tr = ForwardsModTofts(trans_K1, trans_k2, trans_Vp, Cp, dt) + +R1_tr, R2st_tr = calcR1_R2(R10, R20st_value, r1, r2st, Ctiss_tr) +dro_S_tr, M_tr = Conc2Sig_Linear(Tr, Te, fa, R1_tr, R2st_tr, signal, 1, M) + +dro_Snoise_tr = addnoise(1, stdS, dro_S_tr, Nt=data_shape[-1]) + +animate_mri(signal, mode="time", slice_index=7, time_index=5) +animate_mri(dro_Snoise_tr, mode="time", slice_index=7, time_index=5) +animate_mri(dro_Snoise, mode="time", slice_index=7, time_index=5) diff --git a/src/osipi/DRO/roi_selection.py b/src/osipi/DRO/roi_selection.py index 3752929..2709f7c 100644 --- a/src/osipi/DRO/roi_selection.py +++ b/src/osipi/DRO/roi_selection.py @@ -47,6 +47,14 @@ def onselect(verts): mask4d = np.zeros(data_shape) mask4d[:, :, slice_num, :] = timemask if save: - np.save("mask/roi_mask.npy", mask4d) - np.save("mask/roi_voxels.npy", roivox) + np.save("ROI_saved/aif_mask.npy", mask4d) + np.save("ROI_saved/roi_voxels.npy", roivox) return mask4d, roivox, lasso + + +def ICfromROI(E, mask, roivox, numaxis): + Eroi = ( + np.sum(mask * E, axis=tuple(range(0, numaxis))) + ) / roivox # calculates average roi signal enhancement + + return Eroi From 331b23a54b685e98517e925d4970697eeac5be94 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Sat, 27 Jul 2024 01:17:38 +0300 Subject: [PATCH 06/12] Create ETM nonlinear fitting --- src/osipi/DRO/DICOM_processing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/osipi/DRO/DICOM_processing.py b/src/osipi/DRO/DICOM_processing.py index 86d3fea..e2781d9 100644 --- a/src/osipi/DRO/DICOM_processing.py +++ b/src/osipi/DRO/DICOM_processing.py @@ -6,7 +6,7 @@ from osipi.DRO.filters_and_noise import median_filter -def read_dicom_slices_as_signal(folder_path): +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). From 41f55cd4fa26978e561b17b861e0fb280154c12c Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Mon, 29 Jul 2024 01:30:22 +0300 Subject: [PATCH 07/12] create ETM non_linear fitting Added fitting for 1 voxel and for the whole image --- src/osipi/{ => DRO}/Conc2DROSignal.py | 0 src/osipi/DRO/Model.py | 81 ++++++++++++++++----------- src/osipi/DRO/main.py | 56 ++++++++++++++++-- 3 files changed, 97 insertions(+), 40 deletions(-) rename src/osipi/{ => DRO}/Conc2DROSignal.py (100%) diff --git a/src/osipi/Conc2DROSignal.py b/src/osipi/DRO/Conc2DROSignal.py similarity index 100% rename from src/osipi/Conc2DROSignal.py rename to src/osipi/DRO/Conc2DROSignal.py diff --git a/src/osipi/DRO/Model.py b/src/osipi/DRO/Model.py index d32e384..8ffbd3d 100644 --- a/src/osipi/DRO/Model.py +++ b/src/osipi/DRO/Model.py @@ -1,5 +1,7 @@ import numpy as np +from scipy import integrate from scipy.integrate import cumtrapz, trapz +from scipy.optimize import curve_fit def modifiedToftsMurase(Cp, Ctiss, dt, datashape): @@ -95,55 +97,66 @@ def tofts(cp, c_tiss, dt, datashape): return ktrans, kep -def exntended_tofts(cp, c_tiss, dt, datashape): +# Curve fitting function for a single voxel's time series data + + +def Extended_Tofts_Integral(t, Cp, Kt=0.1, ve=0.1, vp=0.2, uniform_sampling=True): + nt = len(t) + Ct = np.zeros(nt) + for k in range(nt): + tmp = vp * Cp[: k + 1] + integrate.cumtrapz( + np.exp(-Kt * (t[k] - t[: k + 1]) / ve) * Cp[: k + 1], t[: k + 1], initial=0.0 + ) + Ct[k] = tmp[-1] + return Ct + + +def FIT_single_voxel(ct, ca, time): + def fit_func(t, kt, ve, vp): + return Extended_Tofts_Integral(t, ca, Kt=kt, ve=ve, vp=vp) + + ini = [0.1, 0.1, 0.2] # Initial guess for [Kt, ve, vp] + popt, pcov = curve_fit(fit_func, time, ct, p0=ini) + return popt, pcov + + +def extended_tofts_model(ca, c_tiss, t, datashape): """ Extended Tofts model for DCE-MRI DRO - + ca -- arterial input function + c_tiss -- 4D array of tissue concentration data (x, y, z, time) + dt -- time interval between samples + datashape -- shape of the data """ ktrans = np.zeros(c_tiss.shape[:-1]) - kep = np.zeros(c_tiss.shape[:-1]) + ve = np.zeros(c_tiss.shape[:-1]) vp = np.zeros(c_tiss.shape[:-1]) for k in range(0, datashape[2]): + print(f"Processing slice {k+1}/{datashape[2]}") for j in range(0, datashape[0]): + print(f"Processing row {j+1}/{datashape[0]}") for i in range(0, datashape[1]): - c = c_tiss[j, i, k, :] - cp_integrated = cumtrapz(cp, dx=dt, initial=0) - c_tiss_integrated = -cumtrapz(c_tiss[i, j, k, :], dx=dt, initial=0) - A = np.column_stack((cp_integrated, c_tiss_integrated, cp)) - ktrans_voxel, kep_voxel, vp_voxel = np.linalg.lstsq(A, c, rcond=None)[0] - ktrans[j, i, k] = ktrans_voxel - kep[j, i, k] = kep_voxel - vp[j, i, k] = vp_voxel + ct = c_tiss[j, i, k, :] + popt, _ = FIT_single_voxel(ct, ca, t) + ktrans[j, i, k], ve[j, i, k], vp[j, i, k] = popt - return ktrans, kep, vp + return ktrans, ve, vp -def forward_tofts(ktrans, kep, cp, vp, dt): +def extended_tofts_model_1vox(ca, c_tiss, t): + """ + Extended Tofts model for DCE-MRI DRO + ca -- arterial input function + c_tiss -- 1D array of tissue concentration data (time) + dt -- time interval between samples """ - Forward Tofts model for DCE-MRI DRO - Parameters: - ktrans (numpy.ndarray): The transfer constant Ktrans. - kep (numpy.ndarray): The rate constant kep. - cp (numpy.ndarray): The plasma concentration C_p(t). - vp (numpy.ndarray): The plasma volume fraction v_p. - dt (float): The time step between measurements. + ct = c_tiss[:] + popt, _ = FIT_single_voxel(ct, ca, t) + ktrans, ve, vp = popt - Returns: - numpy.ndarray: The tissue concentration C_tiss(t). - """ - time_points = cp.shape[-1] - c_tiss = np.zeros(ktrans.shape) - for t in range(time_points): - if t == 0: - c_tiss[..., t] = vp * cp[..., t] - else: - exp = np.exp(-kep * np.arange(t + 1)[::-1] * dt) - integral = trapz(cp[..., : t + 1] * exp, dx=dt, axis=-1) - c_tiss[..., t] = vp * cp[..., t] + ktrans * integral - - return c_tiss + return ktrans, ve, vp def ForwardsModTofts(K1, k2, Vp, Cp, dt): diff --git a/src/osipi/DRO/main.py b/src/osipi/DRO/main.py index c66ad36..76e8463 100644 --- a/src/osipi/DRO/main.py +++ b/src/osipi/DRO/main.py @@ -1,16 +1,27 @@ +import time + import numpy as np from DICOM_processing import ( SignalEnhancementExtract, - read_dicom_slices_as_signal, + read_dicom_slices_as_4d_signal, ) from Display import animate_mri +from matplotlib import pyplot as plt from roi_selection import ICfromROI -from osipi.Conc2DROSignal import Conc2Sig_Linear, STDmap, addnoise, calcR1_R2, createR10_withref +import osipi +from osipi.DRO.Conc2DROSignal import Conc2Sig_Linear, STDmap, addnoise, calcR1_R2, createR10_withref from osipi.DRO.filters_and_noise import median_filter -from osipi.DRO.Model import ForwardsModTofts, modifiedToftsMurase, modifiedToftsMurase1Vox +from osipi.DRO.Model import ( + ForwardsModTofts, + extended_tofts_model_1vox, + modifiedToftsMurase, + modifiedToftsMurase1Vox, +) -signal, slices, dicom_ref = read_dicom_slices_as_signal("data/subject2/12.000000-perfusion-17557") +signal, slices, dicom_ref = read_dicom_slices_as_4d_signal( + "data/subject2/12.000000-perfusion-17557" +) # anim = animate_mri(slices, mode="time", slice_index=7, time_index=5) data_shape = signal.shape @@ -34,16 +45,49 @@ Ev = ICfromROI(E, sagmask, roivoxv, sagmask.ndim - 1) S0ref = ICfromROI(S0[:, :, z], sagmask[:, :, z, 0], roivoxv, sagmask.ndim - 2) +max_index = np.unravel_index(np.argmax(E * aifmask, axis=None), E.shape) +print(max_index) + Hct = 0.45 -K1, k2, Vp = modifiedToftsMurase(Ea / (1 - Hct), E, dt, data_shape) +E_vox = E[max_index[0], max_index[1], max_index[2], :] + +start_time = time.time() + +k1, ve, vp = extended_tofts_model_1vox(Ea, E_vox, t) + +end_time = time.time() + +execution_time = end_time - start_time + +print(f"The execution time of the line is: {execution_time} seconds") + +K1, K2, Vp = modifiedToftsMurase(Ea / (1 - Hct), E, dt, data_shape) + +k1_vox = K1[max_index[0], max_index[1], max_index[2]] +k2_vox = K2[max_index[0], max_index[1], max_index[2]] +Vp_vox = Vp[max_index[0], max_index[1], max_index[2]] + +ct = osipi.extended_tofts(t, Ea, k1_vox, k1_vox / k2_vox, Vp_vox) +ct_real = E[max_index[0], max_index[1], max_index[2], :] +ct_extended = osipi.extended_tofts(t, Ea, k1, ve, vp) + +plt.figure(figsize=(10, 6)) +plt.plot(t, ct, label="ct_Mudified_Tofts") +plt.plot(t, ct_real, label="ct_Raw") +plt.plot(t, ct_extended, label="ct_Extended_Tofts") +plt.xlabel("Time") +plt.ylabel("Concentration") +plt.title("ct_raw vs model") +plt.legend() +plt.show() pvc_K1, pvc_k2, pvc_Vp = modifiedToftsMurase1Vox(Ea / (1 - Hct), Ev, dt, data_shape) pvc = abs((1 - Hct) / pvc_Vp) # sagittal sinus Vp should be 0.55, so calc correction factor if not # Apply correction factor to fitted parameters cor_K1 = K1 * pvc -cor_k2 = k2 * pvc +cor_k2 = K2 * pvc cor_Vp = Vp * pvc # Apply Median Filter to parameters all with footprint (3,3) From 47131c98ed397aca88cac0542527dafff32be524 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Sat, 3 Aug 2024 01:33:59 +0300 Subject: [PATCH 08/12] finished full pipeline using different models --- src/osipi/DRO/Display.py | 2 +- src/osipi/DRO/Model.py | 165 ++++++++++++------ src/osipi/DRO/main.py | 364 ++++++++++++++++++++------------------- 3 files changed, 300 insertions(+), 231 deletions(-) diff --git a/src/osipi/DRO/Display.py b/src/osipi/DRO/Display.py index 7e16e18..be815af 100644 --- a/src/osipi/DRO/Display.py +++ b/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=10000, blit=False + fig=fig, func=animate, frames=frames, init_func=init, interval=100, blit=False ) plt.show() return anim diff --git a/src/osipi/DRO/Model.py b/src/osipi/DRO/Model.py index 8ffbd3d..e0d77d9 100644 --- a/src/osipi/DRO/Model.py +++ b/src/osipi/DRO/Model.py @@ -1,8 +1,11 @@ +import multiprocessing as mp + import numpy as np -from scipy import integrate -from scipy.integrate import cumtrapz, trapz +from scipy.integrate import cumtrapz, cumulative_trapezoid, trapz from scipy.optimize import curve_fit +import osipi + def modifiedToftsMurase(Cp, Ctiss, dt, datashape): # Fit Modified Tofts (Linear from Murase, 2004) @@ -73,73 +76,67 @@ def modifiedToftsMurase1Vox(Cp, Ctiss, dt, datashape): return K1, k2, Vp -def tofts(cp, c_tiss, dt, datashape): - """ - Tofts model for DCE-MRI DRO +def Extended_Tofts_Integral(t, Cp, Kt=0.1, ve=0.1, vp=0.2, uniform_sampling=True): + exp_term = np.exp(-Kt * (t[:, None] - t[None, :]) / ve) + exp_term = np.tril(exp_term, k=0) + integral_term = cumulative_trapezoid(exp_term * Cp[None, :], t, axis=1, initial=0) + Ct = vp * Cp + integral_term[:, -1] - ref : https://github.com/JCardenasRdz/Gage-repeatability-DCE-MRI?tab=readme-ov-file + return Ct - """ - ktrans = np.zeros(c_tiss.shape[:-1]) - kep = np.zeros(c_tiss.shape[:-1]) - for k in range(0, datashape[2]): - for j in range(0, datashape[0]): - for i in range(0, datashape[1]): - c = c_tiss[j, i, k, :] - cp_integrated = cumtrapz(cp, dx=dt, initial=0) - c_tiss_integrated = -cumtrapz(c_tiss[i, j, k, :], dx=dt, initial=0) - A = np.column_stack((cp_integrated, c_tiss_integrated)) - ktrans_voxel, kep_voxel = np.linalg.lstsq(A, c, rcond=None)[0] - ktrans[j, i, k] = ktrans_voxel - kep[j, i, k] = kep_voxel +def Tofts_Integral(t, Cp, Kt=0.1, ve=0.1): + exp_term = np.exp(-Kt * (t[:, None] - t[None, :]) / ve) + exp_term = np.tril(exp_term, k=0) + integral_term = cumulative_trapezoid(exp_term * Cp[None, :], t, axis=1, initial=0) + Ct = integral_term[:, -1] + return Ct - return ktrans, kep +def FIT_single_voxel(ct, ca, time): + def fit_func_ET(t, kt, ve, vp): + return Extended_Tofts_Integral(t, ca, Kt=kt, ve=ve, vp=vp) -# Curve fitting function for a single voxel's time series data + ini = [0.1, 0.1, 0.2] + popt, pcov = curve_fit(fit_func_ET, time, ct, p0=ini) + return popt -def Extended_Tofts_Integral(t, Cp, Kt=0.1, ve=0.1, vp=0.2, uniform_sampling=True): - nt = len(t) - Ct = np.zeros(nt) - for k in range(nt): - tmp = vp * Cp[: k + 1] + integrate.cumtrapz( - np.exp(-Kt * (t[k] - t[: k + 1]) / ve) * Cp[: k + 1], t[: k + 1], initial=0.0 - ) - Ct[k] = tmp[-1] - return Ct +def FIT_single_voxel_tofts(ct, ca, time): + def fit_func_T(t, kt, ve): + return Tofts_Integral(t, ca, Kt=kt, ve=ve) + ini = [0.1, 0.1] + popt, pcov = curve_fit(fit_func_T, time, ct, p0=ini) + return popt -def FIT_single_voxel(ct, ca, time): - def fit_func(t, kt, ve, vp): - return Extended_Tofts_Integral(t, ca, Kt=kt, ve=ve, vp=vp) - ini = [0.1, 0.1, 0.2] # Initial guess for [Kt, ve, vp] - popt, pcov = curve_fit(fit_func, time, ct, p0=ini) - return popt, pcov +def process_voxel(j, i, k, c_tiss, ca, t, type="ET"): + ct = c_tiss[j, i, k, :] + if type == "ET": + popt = FIT_single_voxel(ct, ca, t) + elif type == "T": + popt = FIT_single_voxel_tofts(ct, ca, t) + return j, i, k, popt -def extended_tofts_model(ca, c_tiss, t, datashape): - """ - Extended Tofts model for DCE-MRI DRO - ca -- arterial input function - c_tiss -- 4D array of tissue concentration data (x, y, z, time) - dt -- time interval between samples - datashape -- shape of the data - """ +def extended_tofts_model(ca, c_tiss, t): ktrans = np.zeros(c_tiss.shape[:-1]) ve = np.zeros(c_tiss.shape[:-1]) vp = np.zeros(c_tiss.shape[:-1]) - for k in range(0, datashape[2]): - print(f"Processing slice {k+1}/{datashape[2]}") - for j in range(0, datashape[0]): - print(f"Processing row {j+1}/{datashape[0]}") - for i in range(0, datashape[1]): - ct = c_tiss[j, i, k, :] - popt, _ = FIT_single_voxel(ct, ca, t) - ktrans[j, i, k], ve[j, i, k], vp[j, i, k] = popt + tasks = [ + (j, i, k, c_tiss, ca, t) + for k in range(c_tiss.shape[2]) + for j in range(c_tiss.shape[0]) + for i in range(c_tiss.shape[1]) + ] + + with mp.Pool(processes=mp.cpu_count()) as pool: + results = pool.starmap(process_voxel, tasks) + + for j, i, k, popt in results: + ktrans[j, i, k], ve[j, i, k], vp[j, i, k] = popt return ktrans, ve, vp @@ -154,9 +151,28 @@ def extended_tofts_model_1vox(ca, c_tiss, t): ct = c_tiss[:] popt, _ = FIT_single_voxel(ct, ca, t) - ktrans, ve, vp = popt - return ktrans, ve, vp + return popt + + +def tofts_model(ca, c_tiss, t): + ktrans = np.zeros(c_tiss.shape[:-1]) + ve = np.zeros(c_tiss.shape[:-1]) + + tasks = [ + (j, i, k, c_tiss, ca, t, "T") + for k in range(c_tiss.shape[2]) + for j in range(c_tiss.shape[0]) + for i in range(c_tiss.shape[1]) + ] + + with mp.Pool(processes=mp.cpu_count()) as pool: + results = pool.starmap(process_voxel, tasks) + + for j, i, k, popt in results: + ktrans[j, i, k], ve[j, i, k] = popt + + return ktrans, ve def ForwardsModTofts(K1, k2, Vp, Cp, dt): @@ -184,3 +200,44 @@ def ForwardsModTofts(K1, k2, Vp, Cp, dt): Ctiss[:, :, :, tk] = np.matmul(B, A).squeeze() return Ctiss + + +def ForwardsModTofts_1vox(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 + + t = Cp.shape[0] + + Ctiss = np.zeros(t) + + b1 = K1 + k2 * Vp # define combined parameter + B = np.zeros((1, 3)) + A = np.zeros((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 + + +def forward_extended_tofts(K1, Ve, Vp, Ca, time): + x, y, z = K1.shape + t = Ca.shape[0] + c_tiss = np.zeros((y, x, z, t)) + + for k in range(0, K1.shape[2]): + for j in range(0, K1.shape[0]): + for i in range(0, K1.shape[1]): + c_tiss[i, j, k, :] = osipi.extended_tofts( + time, Ca, K1[i, j, k], Ve[i, j, k], Vp[i, j, k] + ) + + return c_tiss diff --git a/src/osipi/DRO/main.py b/src/osipi/DRO/main.py index 76e8463..1eefd6a 100644 --- a/src/osipi/DRO/main.py +++ b/src/osipi/DRO/main.py @@ -6,197 +6,209 @@ read_dicom_slices_as_4d_signal, ) from Display import animate_mri -from matplotlib import pyplot as plt from roi_selection import ICfromROI -import osipi from osipi.DRO.Conc2DROSignal import Conc2Sig_Linear, STDmap, addnoise, calcR1_R2, createR10_withref from osipi.DRO.filters_and_noise import median_filter from osipi.DRO.Model import ( - ForwardsModTofts, - extended_tofts_model_1vox, - modifiedToftsMurase, + extended_tofts_model, + forward_extended_tofts, modifiedToftsMurase1Vox, ) -signal, slices, dicom_ref = read_dicom_slices_as_4d_signal( - "data/subject2/12.000000-perfusion-17557" -) -# anim = animate_mri(slices, mode="time", slice_index=7, time_index=5) -data_shape = signal.shape - -E, S0, S = SignalEnhancementExtract(signal, data_shape, 5) - -Max = np.max(E, axis=-1, keepdims=True) - -dt = 4.8 / 60 # mins from the RIDER website DCE description -t = np.linspace(0, dt * S.shape[-1], S.shape[-1]) # time series points - -aif_dir = "ROI_saved/" - -aifmask = np.load("{}aifmask.npy".format(aif_dir)) -sagmask = np.load("{}sagmask.npy".format(aif_dir)) -roivoxa = np.load("{}aifmaskvox.npy".format(aif_dir)) -roivoxv = np.load("{}sagmaskvox.npy".format(aif_dir)) - -z = 5 - -Ea = ICfromROI(E, aifmask, roivoxa, aifmask.ndim - 1) -Ev = ICfromROI(E, sagmask, roivoxv, sagmask.ndim - 1) -S0ref = ICfromROI(S0[:, :, z], sagmask[:, :, z, 0], roivoxv, sagmask.ndim - 2) - -max_index = np.unravel_index(np.argmax(E * aifmask, axis=None), E.shape) -print(max_index) - -Hct = 0.45 - -E_vox = E[max_index[0], max_index[1], max_index[2], :] - -start_time = time.time() - -k1, ve, vp = extended_tofts_model_1vox(Ea, E_vox, t) - -end_time = time.time() - -execution_time = end_time - start_time - -print(f"The execution time of the line is: {execution_time} seconds") - -K1, K2, Vp = modifiedToftsMurase(Ea / (1 - Hct), E, dt, data_shape) - -k1_vox = K1[max_index[0], max_index[1], max_index[2]] -k2_vox = K2[max_index[0], max_index[1], max_index[2]] -Vp_vox = Vp[max_index[0], max_index[1], max_index[2]] - -ct = osipi.extended_tofts(t, Ea, k1_vox, k1_vox / k2_vox, Vp_vox) -ct_real = E[max_index[0], max_index[1], max_index[2], :] -ct_extended = osipi.extended_tofts(t, Ea, k1, ve, vp) - -plt.figure(figsize=(10, 6)) -plt.plot(t, ct, label="ct_Mudified_Tofts") -plt.plot(t, ct_real, label="ct_Raw") -plt.plot(t, ct_extended, label="ct_Extended_Tofts") -plt.xlabel("Time") -plt.ylabel("Concentration") -plt.title("ct_raw vs model") -plt.legend() -plt.show() - -pvc_K1, pvc_k2, pvc_Vp = modifiedToftsMurase1Vox(Ea / (1 - Hct), Ev, dt, data_shape) -pvc = abs((1 - Hct) / pvc_Vp) # sagittal sinus Vp should be 0.55, so calc correction factor if not - -# Apply correction factor to fitted parameters -cor_K1 = K1 * pvc -cor_k2 = K2 * pvc -cor_Vp = Vp * pvc -# Apply Median Filter to parameters all with footprint (3,3) - -mf_K1 = median_filter(cor_K1) -mf_k2 = median_filter(cor_k2) -mf_Vp = median_filter(cor_Vp) - -mf_Vp[mf_Vp <= 0] = 0 -mf_Vp[mf_Vp > 1] = 1.0 -mf_K1[mf_K1 <= 0] = 0 -mf_k2[mf_k2 <= 0] = 0 - -# evolve forwards model -aif_cor = np.max((Ea / (1 - Hct)) * pvc) / 6 # caluclate enhancement to concentration correction -Cp = ( - (Ea / (1 - Hct)) * pvc -) / aif_cor # calculate Cp using signal enhancement aif and concentration conversion factor - -c_tissue = ForwardsModTofts(mf_K1, mf_k2, mf_Vp, Cp, dt) - -r1 = 3.9 # longitudinal relaxivity Gd-DTPA (Hz/mM) source: (Pintaske,2006) -r2st = 10 # transverse relaxivity Gd-DTPA (Hz/mM) -# roughly estimated using (Pintaske,2006) and (Siemonsen, 2008) - -fa = np.deg2rad(float(dicom_ref.FlipAngle)) # flip angle (rads) - -Te = 1e-3 * float(dicom_ref.EchoTime) # Echo Time (s) -Tr = 1e-3 * float(dicom_ref.RepetitionTime) # Repetition Time (s) -T1b = 1.48 # T1 for blood measured in sagittal sinus @ 1.5T (s) (Zhang et al., 2013) - -R10_value = 1.18 # precontrast T1 relaxation rate (Hz) brain avg (radiopedia) -R20st_value = ( - 17.24 # precontrast T1 relaxation rate (Hz) brain avg using T2* from (Siemonsen, 2008) -) -R10 = createR10_withref( - S0ref, S0, Tr, fa, T1b, data_shape -) # precontrast R1 map (normalised to sagittal sinus) - -R1, R2st = calcR1_R2(R10, R20st_value, r1, r2st, c_tissue) # returns R10 and R2st maps -dro_S, M = Conc2Sig_Linear(Tr, Te, fa, R1, R2st, S, 1, 0) - -stdS = STDmap(signal, t0=5) # caluclate Standard deviation for original data -dro_Snoise = addnoise(1, stdS, dro_S, Nt=data_shape[-1]) - -trans_K1 = mf_K1.copy() -trans_k2 = mf_k2.copy() -trans_Vp = mf_Vp.copy() - -vmax_K1, vmax_k2, vmax_Vp = 1, 1, 0.2 -vmin_K1, vmin_k2, vmin_Vp = 0.2, 0.2, 0.01 -lb_K1, lb_k2, lb_Vp = 0.54, 0.52, 0.49 -ub_K1, ub_k2, ub_Vp = 1.52, 1.5, 1.43 -lim_K1, lim_k2, lim_Vp = vmax_K1 + 0.5, vmax_k2 + 0.1, vmax_Vp + 0.5 -ub_lim = 1.01 - -trans_K1[trans_K1 <= vmin_K1] = trans_K1[trans_K1 <= vmin_K1] * lb_K1 -trans_K1[trans_K1 >= lim_K1] = trans_K1[trans_K1 >= lim_K1] * ub_lim -trans_K1[(trans_K1 >= vmax_K1) & (trans_K1 < lim_K1)] = ( - trans_K1[(trans_K1 >= vmax_K1) & (trans_K1 < lim_K1)] * ub_K1 -) -trans_K1[(trans_K1 > vmin_K1) & (trans_K1 < vmax_K1)] = trans_K1[ - (trans_K1 > vmin_K1) & (trans_K1 < vmax_K1) -] * ( - lb_K1 - + ( - ((trans_K1[(trans_K1 > vmin_K1) & (trans_K1 < vmax_K1)] - vmin_K1) / (vmax_K1 - vmin_K1)) - * (ub_K1 - lb_K1) +if __name__ == "__main__": + signal, slices, dicom_ref = read_dicom_slices_as_4d_signal( + "data/subject2/12.000000-perfusion-17557" + ) + # anim = animate_mri(slices, mode="time", slice_index=7, time_index=5) + data_shape = signal.shape + + E, S0, S = SignalEnhancementExtract(signal, data_shape, 5) + + Max = np.max(E, axis=-1, keepdims=True) + + dt = 4.8 / 60 # mins from the RIDER website DCE description + t = np.linspace(0, dt * S.shape[-1], S.shape[-1]) # time series points + + aif_dir = "ROI_saved/" + + aifmask = np.load("{}aifmask.npy".format(aif_dir)) + sagmask = np.load("{}sagmask.npy".format(aif_dir)) + roivoxa = np.load("{}aifmaskvox.npy".format(aif_dir)) + roivoxv = np.load("{}sagmaskvox.npy".format(aif_dir)) + + z = 5 + + Ea = ICfromROI(E, aifmask, roivoxa, aifmask.ndim - 1) + Ev = ICfromROI(E, sagmask, roivoxv, sagmask.ndim - 1) + S0ref = ICfromROI(S0[:, :, z], sagmask[:, :, z, 0], roivoxv, sagmask.ndim - 2) + max_index = np.unravel_index(np.argmax(E * aifmask, axis=None), E.shape) + print(max_index) + + Hct = 0.45 + + E_vox = E[max_index[0], max_index[1], max_index[2], :] + + start_time = time.time() + + k1, ve, vp = extended_tofts_model(Ea, E[:, :, :8, :], t) + k1_vox1 = k1[0, 0, 0] + ve_vox1 = ve[0, 0, 0] + vp_vox1 = vp[0, 0, 0] + end_time = time.time() + + execution_time = end_time - start_time + + print(f"The execution time of the line is: {execution_time} seconds") + + # K1, K2, Vp = modifiedToftsMurase(Ea / (1 - Hct), E, dt, data_shape) + # + # k1_vox = K1[max_index[0], max_index[1], max_index[2]] + # k2_vox = K2[max_index[0], max_index[1], max_index[2]] + # Vp_vox = Vp[max_index[0], max_index[1], max_index[2]] + # + # ct = ForwardsModTofts_1vox(k1_vox, k2_vox, Vp_vox, Ea/(1-Hct), dt) + # ct_real = E[max_index[0], max_index[1], max_index[2], :] + # ct_tofts = osipi.extended_tofts(t, Ea, k1_vox1, ve_vox1, vp_vox1) + # + # plt.figure(figsize=(10, 6)) + # plt.plot(t, ct, label="ct_Mudified_Tofts") + # plt.plot(t, ct_real, label="ct_Raw") + # plt.plot(t, ct_tofts, label="ct_ETofts") + # plt.xlabel("Time") + # plt.ylabel("Concentration") + # plt.title("ct_raw vs model") + # plt.legend() + # plt.show() + + pvc_K1, pvc_k2, pvc_Vp = modifiedToftsMurase1Vox(Ea / (1 - Hct), Ev, dt, data_shape) + pvc = abs( + (1 - Hct) / pvc_Vp + ) # sagittal sinus Vp should be 0.55, so calc correction factor if not + + # Apply correction factor to fitted parameters + cor_K1 = k1 * pvc + cor_k2 = ve * pvc + cor_Vp = vp * pvc + # Apply Median Filter to parameters all with footprint (3,3) + + mf_K1 = median_filter(cor_K1) + mf_k2 = median_filter(cor_k2) + mf_Vp = median_filter(cor_Vp) + + mf_Vp[mf_Vp <= 0] = 0 + mf_Vp[mf_Vp > 1] = 1.0 + mf_K1[mf_K1 <= 0] = 0 + mf_k2[mf_k2 <= 0] = 0 + + # evolve forwards model + aif_cor = ( + np.max((Ea / (1 - Hct)) * pvc) / 6 + ) # caluclate enhancement to concentration correction + Cp = ( + (Ea / (1 - Hct)) * pvc + ) / aif_cor # calculate Cp using signal enhancement aif and concentration conversion factor + + c_tissue = forward_extended_tofts(mf_K1, mf_k2, mf_Vp, Ea, t) + + r1 = 3.9 # longitudinal relaxivity Gd-DTPA (Hz/mM) source: (Pintaske,2006) + r2st = 10 # transverse relaxivity Gd-DTPA (Hz/mM) + # roughly estimated using (Pintaske,2006) and (Siemonsen, 2008) + + fa = np.deg2rad(float(dicom_ref.FlipAngle)) # flip angle (rads) + + Te = 1e-3 * float(dicom_ref.EchoTime) # Echo Time (s) + Tr = 1e-3 * float(dicom_ref.RepetitionTime) # Repetition Time (s) + T1b = 1.48 # T1 for blood measured in sagittal sinus @ 1.5T (s) (Zhang et al., 2013) + + R10_value = 1.18 # precontrast T1 relaxation rate (Hz) brain avg (radiopedia) + R20st_value = ( + 17.24 # precontrast T1 relaxation rate (Hz) brain avg using T2* from (Siemonsen, 2008) + ) + R10 = createR10_withref( + S0ref, S0[:, :, :8], Tr, fa, T1b, data_shape + ) # precontrast R1 map (normalised to sagittal sinus) + + R1, R2st = calcR1_R2(R10, R20st_value, r1, r2st, c_tissue) # returns R10 and R2st maps + dro_S, M = Conc2Sig_Linear(Tr, Te, fa, R1, R2st, S, 1, 0) + + stdS = STDmap(signal[:, :, :8, :], t0=5) # caluclate Standard deviation for original data + dro_Snoise = addnoise(1, stdS, dro_S, Nt=data_shape[-1]) + + trans_K1 = mf_K1.copy() + trans_k2 = mf_k2.copy() + trans_Vp = mf_Vp.copy() + + vmax_K1, vmax_k2, vmax_Vp = 1, 1, 0.2 + vmin_K1, vmin_k2, vmin_Vp = 0.2, 0.2, 0.01 + lb_K1, lb_k2, lb_Vp = 0.54, 0.52, 0.49 + ub_K1, ub_k2, ub_Vp = 1.52, 1.5, 1.43 + lim_K1, lim_k2, lim_Vp = vmax_K1 + 0.5, vmax_k2 + 0.1, vmax_Vp + 0.5 + ub_lim = 1.01 + + trans_K1[trans_K1 <= vmin_K1] = trans_K1[trans_K1 <= vmin_K1] * lb_K1 + trans_K1[trans_K1 >= lim_K1] = trans_K1[trans_K1 >= lim_K1] * ub_lim + trans_K1[(trans_K1 >= vmax_K1) & (trans_K1 < lim_K1)] = ( + trans_K1[(trans_K1 >= vmax_K1) & (trans_K1 < lim_K1)] * ub_K1 + ) + trans_K1[(trans_K1 > vmin_K1) & (trans_K1 < vmax_K1)] = trans_K1[ + (trans_K1 > vmin_K1) & (trans_K1 < vmax_K1) + ] * ( + lb_K1 + + ( + ( + (trans_K1[(trans_K1 > vmin_K1) & (trans_K1 < vmax_K1)] - vmin_K1) + / (vmax_K1 - vmin_K1) + ) + * (ub_K1 - lb_K1) + ) ) -) -trans_k2[trans_k2 <= vmin_k2] = trans_k2[trans_k2 <= vmin_k2] * lb_k2 -trans_k2[trans_k2 >= lim_k2] = trans_k2[trans_k2 >= lim_k2] * ub_lim -trans_k2[(trans_k2 >= vmax_k2) & (trans_k2 < lim_k2)] = ( - trans_k2[(trans_k2 >= vmax_k2) & (trans_k2 < lim_k2)] * ub_k2 -) -trans_k2[(trans_k2 > vmin_k2) & (trans_k2 < vmax_k2)] = trans_k2[ - (trans_k2 > vmin_k2) & (trans_k2 < vmax_k2) -] * ( - lb_k2 - + ( - ((trans_k2[(trans_k2 > vmin_k2) & (trans_k2 < vmax_k2)] - vmin_k2) / (vmax_k2 - vmin_k2)) - * (ub_k2 - lb_k2) + trans_k2[trans_k2 <= vmin_k2] = trans_k2[trans_k2 <= vmin_k2] * lb_k2 + trans_k2[trans_k2 >= lim_k2] = trans_k2[trans_k2 >= lim_k2] * ub_lim + trans_k2[(trans_k2 >= vmax_k2) & (trans_k2 < lim_k2)] = ( + trans_k2[(trans_k2 >= vmax_k2) & (trans_k2 < lim_k2)] * ub_k2 + ) + trans_k2[(trans_k2 > vmin_k2) & (trans_k2 < vmax_k2)] = trans_k2[ + (trans_k2 > vmin_k2) & (trans_k2 < vmax_k2) + ] * ( + lb_k2 + + ( + ( + (trans_k2[(trans_k2 > vmin_k2) & (trans_k2 < vmax_k2)] - vmin_k2) + / (vmax_k2 - vmin_k2) + ) + * (ub_k2 - lb_k2) + ) ) -) -trans_Vp[trans_Vp <= vmin_Vp] = trans_Vp[trans_Vp <= vmin_Vp] * lb_Vp -trans_Vp[(trans_Vp >= lim_Vp)] = trans_Vp[trans_Vp >= lim_Vp] * ub_lim -trans_Vp[(trans_Vp >= vmax_Vp) & (trans_Vp < lim_Vp)] = ( - trans_Vp[(trans_Vp >= vmax_Vp) & (trans_Vp < lim_Vp)] * ub_Vp -) -trans_Vp[(trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp)] = trans_Vp[ - (trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp) -] * ( - lb_Vp - + ( - ((trans_Vp[(trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp)] - vmin_Vp) / (vmax_Vp - vmin_Vp)) - * (ub_Vp - lb_Vp) + trans_Vp[trans_Vp <= vmin_Vp] = trans_Vp[trans_Vp <= vmin_Vp] * lb_Vp + trans_Vp[(trans_Vp >= lim_Vp)] = trans_Vp[trans_Vp >= lim_Vp] * ub_lim + trans_Vp[(trans_Vp >= vmax_Vp) & (trans_Vp < lim_Vp)] = ( + trans_Vp[(trans_Vp >= vmax_Vp) & (trans_Vp < lim_Vp)] * ub_Vp + ) + trans_Vp[(trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp)] = trans_Vp[ + (trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp) + ] * ( + lb_Vp + + ( + ( + (trans_Vp[(trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp)] - vmin_Vp) + / (vmax_Vp - vmin_Vp) + ) + * (ub_Vp - lb_Vp) + ) ) -) -trans_Vp[trans_Vp > 1] = 1 + trans_Vp[trans_Vp > 1] = 1 -Ctiss_tr = ForwardsModTofts(trans_K1, trans_k2, trans_Vp, Cp, dt) + Ctiss_tr = forward_extended_tofts(trans_K1, trans_k2, trans_Vp, Ea, t) -R1_tr, R2st_tr = calcR1_R2(R10, R20st_value, r1, r2st, Ctiss_tr) -dro_S_tr, M_tr = Conc2Sig_Linear(Tr, Te, fa, R1_tr, R2st_tr, signal, 1, M) + R1_tr, R2st_tr = calcR1_R2(R10, R20st_value, r1, r2st, Ctiss_tr) + dro_S_tr, M_tr = Conc2Sig_Linear(Tr, Te, fa, R1_tr, R2st_tr, signal[:, :, :8, :], 1, M) -dro_Snoise_tr = addnoise(1, stdS, dro_S_tr, Nt=data_shape[-1]) + dro_Snoise_tr = addnoise(1, stdS, dro_S_tr, Nt=data_shape[-1]) -animate_mri(signal, mode="time", slice_index=7, time_index=5) -animate_mri(dro_Snoise_tr, mode="time", slice_index=7, time_index=5) -animate_mri(dro_Snoise, mode="time", slice_index=7, time_index=5) + animate_mri(signal, mode="time", slice_index=7, time_index=5) + animate_mri(dro_Snoise_tr, mode="time", slice_index=7, time_index=5) + animate_mri(dro_Snoise, mode="time", slice_index=7, time_index=5) From d9c2f9921580a94e62db56c99c8b1787a3aaf04b Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Tue, 20 Aug 2024 17:54:57 +0300 Subject: [PATCH 09/12] Illustrate dro code steps, use osipi models for the fitting --- .pre-commit-config.yaml | 1 + src/osipi/DRO/Model.py | 67 ++++------ src/osipi/DRO/main.py | 281 ++++++++++++++++++++-------------------- 3 files changed, 170 insertions(+), 179 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 2dd8c9b..0f7ecb4 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -3,6 +3,7 @@ repos: rev: v0.4.8 # Use the latest version hooks: - id: ruff + args: ["--fix"] - id: ruff-format - repo: https://github.com/pycqa/flake8 rev: 7.0.0 diff --git a/src/osipi/DRO/Model.py b/src/osipi/DRO/Model.py index e0d77d9..96cdb8e 100644 --- a/src/osipi/DRO/Model.py +++ b/src/osipi/DRO/Model.py @@ -1,7 +1,7 @@ import multiprocessing as mp import numpy as np -from scipy.integrate import cumtrapz, cumulative_trapezoid, trapz +from scipy.integrate import cumtrapz, trapz from scipy.optimize import curve_fit import osipi @@ -76,37 +76,20 @@ def modifiedToftsMurase1Vox(Cp, Ctiss, dt, datashape): return K1, k2, Vp -def Extended_Tofts_Integral(t, Cp, Kt=0.1, ve=0.1, vp=0.2, uniform_sampling=True): - exp_term = np.exp(-Kt * (t[:, None] - t[None, :]) / ve) - exp_term = np.tril(exp_term, k=0) - integral_term = cumulative_trapezoid(exp_term * Cp[None, :], t, axis=1, initial=0) - Ct = vp * Cp + integral_term[:, -1] - - return Ct - - -def Tofts_Integral(t, Cp, Kt=0.1, ve=0.1): - exp_term = np.exp(-Kt * (t[:, None] - t[None, :]) / ve) - exp_term = np.tril(exp_term, k=0) - integral_term = cumulative_trapezoid(exp_term * Cp[None, :], t, axis=1, initial=0) - Ct = integral_term[:, -1] - return Ct - - -def FIT_single_voxel(ct, ca, time): +def fit_single_voxel_extended_tofts(ct, ca, time): def fit_func_ET(t, kt, ve, vp): - return Extended_Tofts_Integral(t, ca, Kt=kt, ve=ve, vp=vp) + return osipi.extended_tofts(t, ca, kt, ve, vp) - ini = [0.1, 0.1, 0.2] + ini = [0, 0, 0] popt, pcov = curve_fit(fit_func_ET, time, ct, p0=ini) return popt -def FIT_single_voxel_tofts(ct, ca, time): +def fit_single_voxel_tofts(ct, ca, time): def fit_func_T(t, kt, ve): - return Tofts_Integral(t, ca, Kt=kt, ve=ve) + return osipi.tofts(t, ca, kt, ve) - ini = [0.1, 0.1] + ini = [0, 0] popt, pcov = curve_fit(fit_func_T, time, ct, p0=ini) return popt @@ -114,9 +97,9 @@ def fit_func_T(t, kt, ve): def process_voxel(j, i, k, c_tiss, ca, t, type="ET"): ct = c_tiss[j, i, k, :] if type == "ET": - popt = FIT_single_voxel(ct, ca, t) + popt = fit_single_voxel_extended_tofts(ct, ca, t) elif type == "T": - popt = FIT_single_voxel_tofts(ct, ca, t) + popt = fit_single_voxel_tofts(ct, ca, t) return j, i, k, popt @@ -150,11 +133,26 @@ def extended_tofts_model_1vox(ca, c_tiss, t): """ ct = c_tiss[:] - popt, _ = FIT_single_voxel(ct, ca, t) + popt = fit_single_voxel_extended_tofts(ct, ca, t) return popt +def forward_extended_tofts(K1, Ve, Vp, Ca, time): + x, y, z = K1.shape + t = Ca.shape[0] + c_tiss = np.zeros((y, x, z, t)) + + for k in range(0, K1.shape[2]): + for j in range(0, K1.shape[0]): + for i in range(0, K1.shape[1]): + c_tiss[i, j, k, :] = osipi.extended_tofts( + time, Ca, K1[i, j, k], Ve[i, j, k], Vp[i, j, k] + ) + + return c_tiss + + def tofts_model(ca, c_tiss, t): ktrans = np.zeros(c_tiss.shape[:-1]) ve = np.zeros(c_tiss.shape[:-1]) @@ -226,18 +224,3 @@ def ForwardsModTofts_1vox(K1, k2, Vp, Cp, dt): Ctiss[tk] = np.matmul(B, A).squeeze() return Ctiss - - -def forward_extended_tofts(K1, Ve, Vp, Ca, time): - x, y, z = K1.shape - t = Ca.shape[0] - c_tiss = np.zeros((y, x, z, t)) - - for k in range(0, K1.shape[2]): - for j in range(0, K1.shape[0]): - for i in range(0, K1.shape[1]): - c_tiss[i, j, k, :] = osipi.extended_tofts( - time, Ca, K1[i, j, k], Ve[i, j, k], Vp[i, j, k] - ) - - return c_tiss diff --git a/src/osipi/DRO/main.py b/src/osipi/DRO/main.py index 1eefd6a..7cad3ba 100644 --- a/src/osipi/DRO/main.py +++ b/src/osipi/DRO/main.py @@ -1,34 +1,37 @@ -import time - import numpy as np from DICOM_processing import ( SignalEnhancementExtract, read_dicom_slices_as_4d_signal, ) from Display import animate_mri +from matplotlib import pyplot as plt from roi_selection import ICfromROI +import osipi from osipi.DRO.Conc2DROSignal import Conc2Sig_Linear, STDmap, addnoise, calcR1_R2, createR10_withref from osipi.DRO.filters_and_noise import median_filter from osipi.DRO.Model import ( extended_tofts_model, + extended_tofts_model_1vox, forward_extended_tofts, - modifiedToftsMurase1Vox, ) if __name__ == "__main__": signal, slices, dicom_ref = read_dicom_slices_as_4d_signal( "data/subject2/12.000000-perfusion-17557" ) - # anim = animate_mri(slices, mode="time", slice_index=7, time_index=5) + # shape of the mri data data_shape = signal.shape - E, S0, S = SignalEnhancementExtract(signal, data_shape, 5) - - Max = np.max(E, axis=-1, keepdims=True) + # extract the signal enhancement + # E is the signal enhancement, S0 is the baseline signal, signal is the original signal + E, S0, signal = SignalEnhancementExtract(signal, data_shape, 5) dt = 4.8 / 60 # mins from the RIDER website DCE description - t = np.linspace(0, dt * S.shape[-1], S.shape[-1]) # time series points + t = np.linspace(0, dt * signal.shape[-1], signal.shape[-1]) # time series points + + # load the ROI masks + # Will be added to create ROI masks from the GUI not just saved ones aif_dir = "ROI_saved/" @@ -42,74 +45,72 @@ Ea = ICfromROI(E, aifmask, roivoxa, aifmask.ndim - 1) Ev = ICfromROI(E, sagmask, roivoxv, sagmask.ndim - 1) S0ref = ICfromROI(S0[:, :, z], sagmask[:, :, z, 0], roivoxv, sagmask.ndim - 2) + + # choose a voxel with maximum enhancement to display the fitting process max_index = np.unravel_index(np.argmax(E * aifmask, axis=None), E.shape) print(max_index) - Hct = 0.45 - - E_vox = E[max_index[0], max_index[1], max_index[2], :] - - start_time = time.time() + # hematocrit value for partial volume correction + hct = 0.45 - k1, ve, vp = extended_tofts_model(Ea, E[:, :, :8, :], t) - k1_vox1 = k1[0, 0, 0] - ve_vox1 = ve[0, 0, 0] - vp_vox1 = vp[0, 0, 0] - end_time = time.time() + # fitting using extended tofts model for first 8 slices + # choose only first 8 slices due to memory constrains + kt, ve, vp = extended_tofts_model(Ea, E[:, :, :8, :], t) - execution_time = end_time - start_time - - print(f"The execution time of the line is: {execution_time} seconds") - - # K1, K2, Vp = modifiedToftsMurase(Ea / (1 - Hct), E, dt, data_shape) - # - # k1_vox = K1[max_index[0], max_index[1], max_index[2]] - # k2_vox = K2[max_index[0], max_index[1], max_index[2]] - # Vp_vox = Vp[max_index[0], max_index[1], max_index[2]] - # - # ct = ForwardsModTofts_1vox(k1_vox, k2_vox, Vp_vox, Ea/(1-Hct), dt) - # ct_real = E[max_index[0], max_index[1], max_index[2], :] - # ct_tofts = osipi.extended_tofts(t, Ea, k1_vox1, ve_vox1, vp_vox1) - # - # plt.figure(figsize=(10, 6)) - # plt.plot(t, ct, label="ct_Mudified_Tofts") - # plt.plot(t, ct_real, label="ct_Raw") - # plt.plot(t, ct_tofts, label="ct_ETofts") - # plt.xlabel("Time") - # plt.ylabel("Concentration") - # plt.title("ct_raw vs model") - # plt.legend() - # plt.show() - - pvc_K1, pvc_k2, pvc_Vp = modifiedToftsMurase1Vox(Ea / (1 - Hct), Ev, dt, data_shape) + # A partial volume correction was applied using the sagittal sinus signal. + pvc_K1, pvc_Ve, pvc_Vp = extended_tofts_model_1vox(Ea, Ev, t) pvc = abs( - (1 - Hct) / pvc_Vp + (1 - hct) / pvc_Vp ) # sagittal sinus Vp should be 0.55, so calc correction factor if not # Apply correction factor to fitted parameters - cor_K1 = k1 * pvc - cor_k2 = ve * pvc + cor_Kt = kt * pvc + cor_Ve = ve * pvc cor_Vp = vp * pvc - # Apply Median Filter to parameters all with footprint (3,3) - mf_K1 = median_filter(cor_K1) - mf_k2 = median_filter(cor_k2) + # Apply Median Filter to parameters all with footprint (3,3) + mf_Kt = median_filter(cor_Kt) + mf_Ve = median_filter(cor_Ve) mf_Vp = median_filter(cor_Vp) + # Apply thresholds to remove negative values + # limit volume fraction to max value of 1 mf_Vp[mf_Vp <= 0] = 0 mf_Vp[mf_Vp > 1] = 1.0 - mf_K1[mf_K1 <= 0] = 0 - mf_k2[mf_k2 <= 0] = 0 + mf_Ve[mf_Ve > 1] = 1.0 + mf_Kt[mf_Kt <= 0] = 0 + mf_Ve[mf_Ve <= 0] = 0 + + # the aif signal is corrected with hematocrit, scaled to 6mM as a realistic value + + aif_cor = np.max((Ea / (1 - hct)) * pvc) / 6 + + # caluclate enhancement to concentration correction + # Cp = ( + # (Ea / (1 - hct)) * pvc + # ) / aif_cor # calculate Cp using signal enhancement aif and concentration conversion factor + + # calculate the concentration in the tissue using the fitted parameters + c_tissue = forward_extended_tofts(mf_Kt, mf_Ve, mf_Vp, (Ea / aif_cor), t) + + # Choosing a specific voxel to plot Concentration curves and the fitting process + kt_vox1 = mf_Kt[96, 118, 5] + ve_vox1 = mf_Ve[96, 118, 5] + vp_vox1 = mf_Vp[96, 118, 5] - # evolve forwards model - aif_cor = ( - np.max((Ea / (1 - Hct)) * pvc) / 6 - ) # caluclate enhancement to concentration correction - Cp = ( - (Ea / (1 - Hct)) * pvc - ) / aif_cor # calculate Cp using signal enhancement aif and concentration conversion factor + # calculate the fitted concentration to compare with the real one + c_tissue_tofts = osipi.extended_tofts(t, Ea, kt_vox1, ve_vox1, vp_vox1) - c_tissue = forward_extended_tofts(mf_K1, mf_k2, mf_Vp, Ea, t) + plt.figure(figsize=(10, 6)) + plt.plot(t, c_tissue_tofts, label="ct_Tofts") + plt.plot(t, E[96, 118, 5, :], label="Ct_raw") + plt.xlabel("Time") + plt.ylabel("Concentration") + plt.title("ct_raw vs model") + plt.legend() + plt.show() + + # calculate the relaxation rates R1 and R2* for the tissue r1 = 3.9 # longitudinal relaxivity Gd-DTPA (Hz/mM) source: (Pintaske,2006) r2st = 10 # transverse relaxivity Gd-DTPA (Hz/mM) @@ -130,85 +131,91 @@ ) # precontrast R1 map (normalised to sagittal sinus) R1, R2st = calcR1_R2(R10, R20st_value, r1, r2st, c_tissue) # returns R10 and R2st maps - dro_S, M = Conc2Sig_Linear(Tr, Te, fa, R1, R2st, S, 1, 0) + + # calculate the signal using the concentration and relaxation rates + # spoiled gradient echo sequence + dro_S, M = Conc2Sig_Linear(Tr, Te, fa, R1, R2st, signal[:, :, :8, :], 1, 0) stdS = STDmap(signal[:, :, :8, :], t0=5) # caluclate Standard deviation for original data + # add noise to the signal dro_Snoise = addnoise(1, stdS, dro_S, Nt=data_shape[-1]) - - trans_K1 = mf_K1.copy() - trans_k2 = mf_k2.copy() - trans_Vp = mf_Vp.copy() - - vmax_K1, vmax_k2, vmax_Vp = 1, 1, 0.2 - vmin_K1, vmin_k2, vmin_Vp = 0.2, 0.2, 0.01 - lb_K1, lb_k2, lb_Vp = 0.54, 0.52, 0.49 - ub_K1, ub_k2, ub_Vp = 1.52, 1.5, 1.43 - lim_K1, lim_k2, lim_Vp = vmax_K1 + 0.5, vmax_k2 + 0.1, vmax_Vp + 0.5 - ub_lim = 1.01 - - trans_K1[trans_K1 <= vmin_K1] = trans_K1[trans_K1 <= vmin_K1] * lb_K1 - trans_K1[trans_K1 >= lim_K1] = trans_K1[trans_K1 >= lim_K1] * ub_lim - trans_K1[(trans_K1 >= vmax_K1) & (trans_K1 < lim_K1)] = ( - trans_K1[(trans_K1 >= vmax_K1) & (trans_K1 < lim_K1)] * ub_K1 - ) - trans_K1[(trans_K1 > vmin_K1) & (trans_K1 < vmax_K1)] = trans_K1[ - (trans_K1 > vmin_K1) & (trans_K1 < vmax_K1) - ] * ( - lb_K1 - + ( - ( - (trans_K1[(trans_K1 > vmin_K1) & (trans_K1 < vmax_K1)] - vmin_K1) - / (vmax_K1 - vmin_K1) - ) - * (ub_K1 - lb_K1) - ) - ) - - trans_k2[trans_k2 <= vmin_k2] = trans_k2[trans_k2 <= vmin_k2] * lb_k2 - trans_k2[trans_k2 >= lim_k2] = trans_k2[trans_k2 >= lim_k2] * ub_lim - trans_k2[(trans_k2 >= vmax_k2) & (trans_k2 < lim_k2)] = ( - trans_k2[(trans_k2 >= vmax_k2) & (trans_k2 < lim_k2)] * ub_k2 - ) - trans_k2[(trans_k2 > vmin_k2) & (trans_k2 < vmax_k2)] = trans_k2[ - (trans_k2 > vmin_k2) & (trans_k2 < vmax_k2) - ] * ( - lb_k2 - + ( - ( - (trans_k2[(trans_k2 > vmin_k2) & (trans_k2 < vmax_k2)] - vmin_k2) - / (vmax_k2 - vmin_k2) - ) - * (ub_k2 - lb_k2) - ) - ) - - trans_Vp[trans_Vp <= vmin_Vp] = trans_Vp[trans_Vp <= vmin_Vp] * lb_Vp - trans_Vp[(trans_Vp >= lim_Vp)] = trans_Vp[trans_Vp >= lim_Vp] * ub_lim - trans_Vp[(trans_Vp >= vmax_Vp) & (trans_Vp < lim_Vp)] = ( - trans_Vp[(trans_Vp >= vmax_Vp) & (trans_Vp < lim_Vp)] * ub_Vp - ) - trans_Vp[(trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp)] = trans_Vp[ - (trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp) - ] * ( - lb_Vp - + ( - ( - (trans_Vp[(trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp)] - vmin_Vp) - / (vmax_Vp - vmin_Vp) - ) - * (ub_Vp - lb_Vp) - ) - ) - - trans_Vp[trans_Vp > 1] = 1 - - Ctiss_tr = forward_extended_tofts(trans_K1, trans_k2, trans_Vp, Ea, t) - - R1_tr, R2st_tr = calcR1_R2(R10, R20st_value, r1, r2st, Ctiss_tr) - dro_S_tr, M_tr = Conc2Sig_Linear(Tr, Te, fa, R1_tr, R2st_tr, signal[:, :, :8, :], 1, M) - - dro_Snoise_tr = addnoise(1, stdS, dro_S_tr, Nt=data_shape[-1]) - - animate_mri(signal, mode="time", slice_index=7, time_index=5) - animate_mri(dro_Snoise_tr, mode="time", slice_index=7, time_index=5) animate_mri(dro_Snoise, mode="time", slice_index=7, time_index=5) + + # + # trans_K1 = mf_K1.copy() + # trans_k2 = mf_k2.copy() + # trans_Vp = mf_Vp.copy() + # + # vmax_K1, vmax_k2, vmax_Vp = 1, 1, 0.2 + # vmin_K1, vmin_k2, vmin_Vp = 0.2, 0.2, 0.01 + # lb_K1, lb_k2, lb_Vp = 0.54, 0.52, 0.49 + # ub_K1, ub_k2, ub_Vp = 1.52, 1.5, 1.43 + # lim_K1, lim_k2, lim_Vp = vmax_K1 + 0.5, vmax_k2 + 0.1, vmax_Vp + 0.5 + # ub_lim = 1.01 + # + # trans_K1[trans_K1 <= vmin_K1] = trans_K1[trans_K1 <= vmin_K1] * lb_K1 + # trans_K1[trans_K1 >= lim_K1] = trans_K1[trans_K1 >= lim_K1] * ub_lim + # trans_K1[(trans_K1 >= vmax_K1) & (trans_K1 < lim_K1)] = ( + # trans_K1[(trans_K1 >= vmax_K1) & (trans_K1 < lim_K1)] * ub_K1 + # ) + # trans_K1[(trans_K1 > vmin_K1) & (trans_K1 < vmax_K1)] = trans_K1[ + # (trans_K1 > vmin_K1) & (trans_K1 < vmax_K1) + # ] * ( + # lb_K1 + # + ( + # ( + # (trans_K1[(trans_K1 > vmin_K1) & (trans_K1 < vmax_K1)] - vmin_K1) + # / (vmax_K1 - vmin_K1) + # ) + # * (ub_K1 - lb_K1) + # ) + # ) + # + # trans_k2[trans_k2 <= vmin_k2] = trans_k2[trans_k2 <= vmin_k2] * lb_k2 + # trans_k2[trans_k2 >= lim_k2] = trans_k2[trans_k2 >= lim_k2] * ub_lim + # trans_k2[(trans_k2 >= vmax_k2) & (trans_k2 < lim_k2)] = ( + # trans_k2[(trans_k2 >= vmax_k2) & (trans_k2 < lim_k2)] * ub_k2 + # ) + # trans_k2[(trans_k2 > vmin_k2) & (trans_k2 < vmax_k2)] = trans_k2[ + # (trans_k2 > vmin_k2) & (trans_k2 < vmax_k2) + # ] * ( + # lb_k2 + # + ( + # ( + # (trans_k2[(trans_k2 > vmin_k2) & (trans_k2 < vmax_k2)] - vmin_k2) + # / (vmax_k2 - vmin_k2) + # ) + # * (ub_k2 - lb_k2) + # ) + # ) + # + # trans_Vp[trans_Vp <= vmin_Vp] = trans_Vp[trans_Vp <= vmin_Vp] * lb_Vp + # trans_Vp[(trans_Vp >= lim_Vp)] = trans_Vp[trans_Vp >= lim_Vp] * ub_lim + # trans_Vp[(trans_Vp >= vmax_Vp) & (trans_Vp < lim_Vp)] = ( + # trans_Vp[(trans_Vp >= vmax_Vp) & (trans_Vp < lim_Vp)] * ub_Vp + # ) + # trans_Vp[(trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp)] = trans_Vp[ + # (trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp) + # ] * ( + # lb_Vp + # + ( + # ( + # (trans_Vp[(trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp)] - vmin_Vp) + # / (vmax_Vp - vmin_Vp) + # ) + # * (ub_Vp - lb_Vp) + # ) + # ) + # + # trans_Vp[trans_Vp > 1] = 1 + # + # Ctiss_tr = forward_extended_tofts(trans_K1, trans_k2, trans_Vp, Ea, t) + # + # R1_tr, R2st_tr = calcR1_R2(R10, R20st_value, r1, r2st, Ctiss_tr) + # dro_S_tr, M_tr = Conc2Sig_Linear(Tr, Te, fa, R1_tr, R2st_tr, signal[:, :, :8, :], 1, M) + # + # dro_Snoise_tr = addnoise(1, stdS, dro_S_tr, Nt=data_shape[-1]) + # + # animate_mri(signal, mode="time", slice_index=7, time_index=5) + # animate_mri(dro_Snoise_tr, mode="time", slice_index=7, time_index=5) + # animate_mri(dro_Snoise, mode="time", slice_index=7, time_index=5) From b05437956158845d867032ba87a6b7c80446d6e7 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Thu, 5 Sep 2024 01:54:33 +0300 Subject: [PATCH 10/12] name refactoring, add imports - add imports as osipi dro - refactor function name to be snake case - added save function in utils --- .gitignore | 1 + src/osipi/DRO/DICOM_processing.py | 116 --------- src/osipi/DRO/Display.py | 35 --- src/osipi/DRO/Model.py | 226 ------------------ src/osipi/DRO/__init__.py | 18 ++ ...onc2DROSignal.py => _conc_2_dro_signal.py} | 2 - src/osipi/DRO/_dicom_processing.py | 55 +++++ ...ers_and_noise.py => _filters_and_noise.py} | 0 src/osipi/DRO/_model.py | 103 ++++++++ .../{roi_selection.py => _roi_selection.py} | 8 +- src/osipi/DRO/main.py | 213 +++++++---------- 11 files changed, 269 insertions(+), 508 deletions(-) delete mode 100644 src/osipi/DRO/DICOM_processing.py delete mode 100644 src/osipi/DRO/Display.py delete mode 100644 src/osipi/DRO/Model.py rename src/osipi/DRO/{Conc2DROSignal.py => _conc_2_dro_signal.py} (96%) create mode 100644 src/osipi/DRO/_dicom_processing.py rename src/osipi/DRO/{filters_and_noise.py => _filters_and_noise.py} (100%) create mode 100644 src/osipi/DRO/_model.py rename src/osipi/DRO/{roi_selection.py => _roi_selection.py} (92%) diff --git a/.gitignore b/.gitignore index 86b9bcf..b250319 100644 --- a/.gitignore +++ b/.gitignore @@ -14,3 +14,4 @@ docs-mk/docs/generated src/osipi/DRO/data src/osipi/DRO/__pycache__/ src/osipi/DRO/ROI_saved/ +output diff --git a/src/osipi/DRO/DICOM_processing.py b/src/osipi/DRO/DICOM_processing.py deleted file mode 100644 index e2781d9..0000000 --- a/src/osipi/DRO/DICOM_processing.py +++ /dev/null @@ -1,116 +0,0 @@ -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 diff --git a/src/osipi/DRO/Display.py b/src/osipi/DRO/Display.py deleted file mode 100644 index be815af..0000000 --- a/src/osipi/DRO/Display.py +++ /dev/null @@ -1,35 +0,0 @@ -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 diff --git a/src/osipi/DRO/Model.py b/src/osipi/DRO/Model.py deleted file mode 100644 index 96cdb8e..0000000 --- a/src/osipi/DRO/Model.py +++ /dev/null @@ -1,226 +0,0 @@ -import multiprocessing as mp - -import numpy as np -from scipy.integrate import cumtrapz, trapz -from scipy.optimize import curve_fit - -import osipi - - -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 fit_single_voxel_extended_tofts(ct, ca, time): - def fit_func_ET(t, kt, ve, vp): - return osipi.extended_tofts(t, ca, kt, ve, vp) - - ini = [0, 0, 0] - popt, pcov = curve_fit(fit_func_ET, time, ct, p0=ini) - return popt - - -def fit_single_voxel_tofts(ct, ca, time): - def fit_func_T(t, kt, ve): - return osipi.tofts(t, ca, kt, ve) - - ini = [0, 0] - popt, pcov = curve_fit(fit_func_T, time, ct, p0=ini) - return popt - - -def process_voxel(j, i, k, c_tiss, ca, t, type="ET"): - ct = c_tiss[j, i, k, :] - if type == "ET": - popt = fit_single_voxel_extended_tofts(ct, ca, t) - elif type == "T": - popt = fit_single_voxel_tofts(ct, ca, t) - return j, i, k, popt - - -def extended_tofts_model(ca, c_tiss, t): - ktrans = np.zeros(c_tiss.shape[:-1]) - ve = np.zeros(c_tiss.shape[:-1]) - vp = np.zeros(c_tiss.shape[:-1]) - - tasks = [ - (j, i, k, c_tiss, ca, t) - for k in range(c_tiss.shape[2]) - for j in range(c_tiss.shape[0]) - for i in range(c_tiss.shape[1]) - ] - - with mp.Pool(processes=mp.cpu_count()) as pool: - results = pool.starmap(process_voxel, tasks) - - for j, i, k, popt in results: - ktrans[j, i, k], ve[j, i, k], vp[j, i, k] = popt - - return ktrans, ve, vp - - -def extended_tofts_model_1vox(ca, c_tiss, t): - """ - Extended Tofts model for DCE-MRI DRO - ca -- arterial input function - c_tiss -- 1D array of tissue concentration data (time) - dt -- time interval between samples - """ - - ct = c_tiss[:] - popt = fit_single_voxel_extended_tofts(ct, ca, t) - - return popt - - -def forward_extended_tofts(K1, Ve, Vp, Ca, time): - x, y, z = K1.shape - t = Ca.shape[0] - c_tiss = np.zeros((y, x, z, t)) - - for k in range(0, K1.shape[2]): - for j in range(0, K1.shape[0]): - for i in range(0, K1.shape[1]): - c_tiss[i, j, k, :] = osipi.extended_tofts( - time, Ca, K1[i, j, k], Ve[i, j, k], Vp[i, j, k] - ) - - return c_tiss - - -def tofts_model(ca, c_tiss, t): - ktrans = np.zeros(c_tiss.shape[:-1]) - ve = np.zeros(c_tiss.shape[:-1]) - - tasks = [ - (j, i, k, c_tiss, ca, t, "T") - for k in range(c_tiss.shape[2]) - for j in range(c_tiss.shape[0]) - for i in range(c_tiss.shape[1]) - ] - - with mp.Pool(processes=mp.cpu_count()) as pool: - results = pool.starmap(process_voxel, tasks) - - for j, i, k, popt in results: - ktrans[j, i, k], ve[j, i, k] = popt - - return ktrans, ve - - -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 - - -def ForwardsModTofts_1vox(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 - - t = Cp.shape[0] - - Ctiss = np.zeros(t) - - b1 = K1 + k2 * Vp # define combined parameter - B = np.zeros((1, 3)) - A = np.zeros((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 diff --git a/src/osipi/DRO/__init__.py b/src/osipi/DRO/__init__.py index e69de29..59bfd88 100644 --- a/src/osipi/DRO/__init__.py +++ b/src/osipi/DRO/__init__.py @@ -0,0 +1,18 @@ +from ._model import (extended_tofts_model, + tofts_model, + extended_tofts_model_1vox, + forward_extended_tofts) + +from _dicom_processing import (read_dicom_slices_as_4d_signal, + signal_enhancement_extract) + +from _roi_selection import (ic_from_roi, roi) + +from _utils import (animate_mri, save_dicoms) +from _filters_and_noise import median_filter + +from _conc_2_dro_signal import (Conc2Sig_Linear, + calcR1_R2, + STDmap, + createR10_withref, + addnoise) diff --git a/src/osipi/DRO/Conc2DROSignal.py b/src/osipi/DRO/_conc_2_dro_signal.py similarity index 96% rename from src/osipi/DRO/Conc2DROSignal.py rename to src/osipi/DRO/_conc_2_dro_signal.py index ad17259..5a9cdd1 100644 --- a/src/osipi/DRO/Conc2DROSignal.py +++ b/src/osipi/DRO/_conc_2_dro_signal.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ @author: eveshalom """ diff --git a/src/osipi/DRO/_dicom_processing.py b/src/osipi/DRO/_dicom_processing.py new file mode 100644 index 0000000..8ba2b10 --- /dev/null +++ b/src/osipi/DRO/_dicom_processing.py @@ -0,0 +1,55 @@ +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 + slices[z] = slices.pop(slice_location) + for t, (_, slice) in enumerate( + sorted(slices[z], key=lambda x: x[0]) + ): # Sort by acquisition time + signal[:, :, z, t] = slice.pixel_array + slices[z][t] = slice + + return signal, slices, slices[0][0] + + +def signal_enhancement_extract(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 diff --git a/src/osipi/DRO/filters_and_noise.py b/src/osipi/DRO/_filters_and_noise.py similarity index 100% rename from src/osipi/DRO/filters_and_noise.py rename to src/osipi/DRO/_filters_and_noise.py diff --git a/src/osipi/DRO/_model.py b/src/osipi/DRO/_model.py new file mode 100644 index 0000000..d6781dc --- /dev/null +++ b/src/osipi/DRO/_model.py @@ -0,0 +1,103 @@ +import multiprocessing as mp + +import numpy as np +from scipy.optimize import curve_fit + +import osipi + + +def fit_single_voxel_extended_tofts(ct, ca, time): + def fit_func_ET(t, kt, ve, vp): + return osipi.extended_tofts(t, ca, kt, ve, vp) + + ini = [0, 0, 0] + popt, pcov = curve_fit(fit_func_ET, time, ct, p0=ini) + return popt + + +def fit_single_voxel_tofts(ct, ca, time): + def fit_func_T(t, kt, ve): + return osipi.tofts(t, ca, kt, ve) + + ini = [0, 0] + popt, pcov = curve_fit(fit_func_T, time, ct, p0=ini) + return popt + + +def process_voxel(j, i, k, c_tiss, ca, t, type="ET"): + ct = c_tiss[j, i, k, :] + if type == "ET": + popt = fit_single_voxel_extended_tofts(ct, ca, t) + elif type == "T": + popt = fit_single_voxel_tofts(ct, ca, t) + return j, i, k, popt + + +def extended_tofts_model(ca, c_tiss, t): + ktrans = np.zeros(c_tiss.shape[:-1]) + ve = np.zeros(c_tiss.shape[:-1]) + vp = np.zeros(c_tiss.shape[:-1]) + + tasks = [ + (j, i, k, c_tiss, ca, t) + for k in range(c_tiss.shape[2]) + for j in range(c_tiss.shape[0]) + for i in range(c_tiss.shape[1]) + ] + + with mp.Pool(processes=mp.cpu_count()) as pool: + results = pool.starmap(process_voxel, tasks) + + for j, i, k, popt in results: + ktrans[j, i, k], ve[j, i, k], vp[j, i, k] = popt + + return ktrans, ve, vp + + +def extended_tofts_model_1vox(ca, c_tiss, t): + """ + Extended Tofts model for DCE-MRI DRO + ca -- arterial input function + c_tiss -- 1D array of tissue concentration data (time) + dt -- time interval between samples + """ + + ct = c_tiss[:] + popt = fit_single_voxel_extended_tofts(ct, ca, t) + + return popt + + +def forward_extended_tofts(K1, Ve, Vp, Ca, time): + x, y, z = K1.shape + t = Ca.shape[0] + c_tiss = np.zeros((y, x, z, t)) + + for k in range(0, K1.shape[2]): + for j in range(0, K1.shape[0]): + for i in range(0, K1.shape[1]): + c_tiss[i, j, k, :] = osipi.extended_tofts( + time, Ca, K1[i, j, k], Ve[i, j, k], Vp[i, j, k] + ) + + return c_tiss + + +def tofts_model(ca, c_tiss, t): + ktrans = np.zeros(c_tiss.shape[:-1]) + ve = np.zeros(c_tiss.shape[:-1]) + + tasks = [ + (j, i, k, c_tiss, ca, t, "T") + for k in range(c_tiss.shape[2]) + for j in range(c_tiss.shape[0]) + for i in range(c_tiss.shape[1]) + ] + + with mp.Pool(processes=mp.cpu_count()) as pool: + results = pool.starmap(process_voxel, tasks) + + for j, i, k, popt in results: + ktrans[j, i, k], ve[j, i, k] = popt + + return ktrans, ve diff --git a/src/osipi/DRO/roi_selection.py b/src/osipi/DRO/_roi_selection.py similarity index 92% rename from src/osipi/DRO/roi_selection.py rename to src/osipi/DRO/_roi_selection.py index 2709f7c..32018fd 100644 --- a/src/osipi/DRO/roi_selection.py +++ b/src/osipi/DRO/_roi_selection.py @@ -52,9 +52,9 @@ def onselect(verts): return mask4d, roivox, lasso -def ICfromROI(E, mask, roivox, numaxis): - Eroi = ( +def ic_from_roi(E, mask, roi_vox, numaxis): + E_roi = ( np.sum(mask * E, axis=tuple(range(0, numaxis))) - ) / roivox # calculates average roi signal enhancement + ) / roi_vox # calculates average roi signal enhancement - return Eroi + return E_roi diff --git a/src/osipi/DRO/main.py b/src/osipi/DRO/main.py index 7cad3ba..3d84ee0 100644 --- a/src/osipi/DRO/main.py +++ b/src/osipi/DRO/main.py @@ -1,23 +1,14 @@ import numpy as np -from DICOM_processing import ( - SignalEnhancementExtract, - read_dicom_slices_as_4d_signal, -) -from Display import animate_mri -from matplotlib import pyplot as plt -from roi_selection import ICfromROI - -import osipi -from osipi.DRO.Conc2DROSignal import Conc2Sig_Linear, STDmap, addnoise, calcR1_R2, createR10_withref -from osipi.DRO.filters_and_noise import median_filter -from osipi.DRO.Model import ( - extended_tofts_model, - extended_tofts_model_1vox, - forward_extended_tofts, -) +import pydicom + +import osipi.DRO as dro if __name__ == "__main__": - signal, slices, dicom_ref = read_dicom_slices_as_4d_signal( + uidprefix = "1.3.6.1.4.1.9328.50.16." + study_instance_uid = pydicom.uid.generate_uid(prefix=uidprefix) + dro_IDnum = "9215224289" + + signal, slices, dicom_ref = dro.read_dicom_slices_as_4d_signal( "data/subject2/12.000000-perfusion-17557" ) # shape of the mri data @@ -25,7 +16,7 @@ # extract the signal enhancement # E is the signal enhancement, S0 is the baseline signal, signal is the original signal - E, S0, signal = SignalEnhancementExtract(signal, data_shape, 5) + E, S0, signal = dro.signal_enhancement_extract(signal, data_shape, 5) dt = 4.8 / 60 # mins from the RIDER website DCE description t = np.linspace(0, dt * signal.shape[-1], signal.shape[-1]) # time series points @@ -42,23 +33,22 @@ z = 5 - Ea = ICfromROI(E, aifmask, roivoxa, aifmask.ndim - 1) - Ev = ICfromROI(E, sagmask, roivoxv, sagmask.ndim - 1) - S0ref = ICfromROI(S0[:, :, z], sagmask[:, :, z, 0], roivoxv, sagmask.ndim - 2) + Ea = dro.ic_from_roi(E, aifmask, roivoxa, aifmask.ndim - 1) + Ev = dro.ic_from_roi(E, sagmask, roivoxv, sagmask.ndim - 1) + S0ref = dro.ic_from_roi(S0[:, :, z], sagmask[:, :, z, 0], roivoxv, sagmask.ndim - 2) # choose a voxel with maximum enhancement to display the fitting process max_index = np.unravel_index(np.argmax(E * aifmask, axis=None), E.shape) - print(max_index) # hematocrit value for partial volume correction hct = 0.45 # fitting using extended tofts model for first 8 slices # choose only first 8 slices due to memory constrains - kt, ve, vp = extended_tofts_model(Ea, E[:, :, :8, :], t) + kt, ve, vp = dro.extended_tofts_model(Ea, E[:, :, :8, :], t) # A partial volume correction was applied using the sagittal sinus signal. - pvc_K1, pvc_Ve, pvc_Vp = extended_tofts_model_1vox(Ea, Ev, t) + pvc_K1, pvc_Ve, pvc_Vp = dro.extended_tofts_model_1vox(Ea, Ev, t) pvc = abs( (1 - hct) / pvc_Vp ) # sagittal sinus Vp should be 0.55, so calc correction factor if not @@ -69,9 +59,9 @@ cor_Vp = vp * pvc # Apply Median Filter to parameters all with footprint (3,3) - mf_Kt = median_filter(cor_Kt) - mf_Ve = median_filter(cor_Ve) - mf_Vp = median_filter(cor_Vp) + mf_Kt = dro.median_filter(cor_Kt) + mf_Ve = dro.median_filter(cor_Ve) + mf_Vp = dro.median_filter(cor_Vp) # Apply thresholds to remove negative values # limit volume fraction to max value of 1 @@ -91,25 +81,13 @@ # ) / aif_cor # calculate Cp using signal enhancement aif and concentration conversion factor # calculate the concentration in the tissue using the fitted parameters - c_tissue = forward_extended_tofts(mf_Kt, mf_Ve, mf_Vp, (Ea / aif_cor), t) + c_tissue = dro.forward_extended_tofts(mf_Kt, mf_Ve, mf_Vp, (Ea / aif_cor), t) # Choosing a specific voxel to plot Concentration curves and the fitting process kt_vox1 = mf_Kt[96, 118, 5] ve_vox1 = mf_Ve[96, 118, 5] vp_vox1 = mf_Vp[96, 118, 5] - # calculate the fitted concentration to compare with the real one - c_tissue_tofts = osipi.extended_tofts(t, Ea, kt_vox1, ve_vox1, vp_vox1) - - plt.figure(figsize=(10, 6)) - plt.plot(t, c_tissue_tofts, label="ct_Tofts") - plt.plot(t, E[96, 118, 5, :], label="Ct_raw") - plt.xlabel("Time") - plt.ylabel("Concentration") - plt.title("ct_raw vs model") - plt.legend() - plt.show() - # calculate the relaxation rates R1 and R2* for the tissue r1 = 3.9 # longitudinal relaxivity Gd-DTPA (Hz/mM) source: (Pintaske,2006) @@ -126,96 +104,81 @@ R20st_value = ( 17.24 # precontrast T1 relaxation rate (Hz) brain avg using T2* from (Siemonsen, 2008) ) - R10 = createR10_withref( + R10 = dro.createR10_withref( S0ref, S0[:, :, :8], Tr, fa, T1b, data_shape ) # precontrast R1 map (normalised to sagittal sinus) - R1, R2st = calcR1_R2(R10, R20st_value, r1, r2st, c_tissue) # returns R10 and R2st maps + R1, R2st = dro.calcR1_R2(R10, R20st_value, r1, r2st, c_tissue) # returns R10 and R2st maps # calculate the signal using the concentration and relaxation rates # spoiled gradient echo sequence - dro_S, M = Conc2Sig_Linear(Tr, Te, fa, R1, R2st, signal[:, :, :8, :], 1, 0) + dro_s, M = dro.Conc2Sig_Linear(Tr, Te, fa, R1, R2st, signal[:, :, :8, :], 1, 0) - stdS = STDmap(signal[:, :, :8, :], t0=5) # caluclate Standard deviation for original data + stdS = dro.STDmap(signal[:, :, :8, :], t0=5) # caluclate Standard deviation for original data # add noise to the signal - dro_Snoise = addnoise(1, stdS, dro_S, Nt=data_shape[-1]) - animate_mri(dro_Snoise, mode="time", slice_index=7, time_index=5) - - # - # trans_K1 = mf_K1.copy() - # trans_k2 = mf_k2.copy() - # trans_Vp = mf_Vp.copy() - # - # vmax_K1, vmax_k2, vmax_Vp = 1, 1, 0.2 - # vmin_K1, vmin_k2, vmin_Vp = 0.2, 0.2, 0.01 - # lb_K1, lb_k2, lb_Vp = 0.54, 0.52, 0.49 - # ub_K1, ub_k2, ub_Vp = 1.52, 1.5, 1.43 - # lim_K1, lim_k2, lim_Vp = vmax_K1 + 0.5, vmax_k2 + 0.1, vmax_Vp + 0.5 - # ub_lim = 1.01 - # - # trans_K1[trans_K1 <= vmin_K1] = trans_K1[trans_K1 <= vmin_K1] * lb_K1 - # trans_K1[trans_K1 >= lim_K1] = trans_K1[trans_K1 >= lim_K1] * ub_lim - # trans_K1[(trans_K1 >= vmax_K1) & (trans_K1 < lim_K1)] = ( - # trans_K1[(trans_K1 >= vmax_K1) & (trans_K1 < lim_K1)] * ub_K1 - # ) - # trans_K1[(trans_K1 > vmin_K1) & (trans_K1 < vmax_K1)] = trans_K1[ - # (trans_K1 > vmin_K1) & (trans_K1 < vmax_K1) - # ] * ( - # lb_K1 - # + ( - # ( - # (trans_K1[(trans_K1 > vmin_K1) & (trans_K1 < vmax_K1)] - vmin_K1) - # / (vmax_K1 - vmin_K1) - # ) - # * (ub_K1 - lb_K1) - # ) - # ) - # - # trans_k2[trans_k2 <= vmin_k2] = trans_k2[trans_k2 <= vmin_k2] * lb_k2 - # trans_k2[trans_k2 >= lim_k2] = trans_k2[trans_k2 >= lim_k2] * ub_lim - # trans_k2[(trans_k2 >= vmax_k2) & (trans_k2 < lim_k2)] = ( - # trans_k2[(trans_k2 >= vmax_k2) & (trans_k2 < lim_k2)] * ub_k2 - # ) - # trans_k2[(trans_k2 > vmin_k2) & (trans_k2 < vmax_k2)] = trans_k2[ - # (trans_k2 > vmin_k2) & (trans_k2 < vmax_k2) - # ] * ( - # lb_k2 - # + ( - # ( - # (trans_k2[(trans_k2 > vmin_k2) & (trans_k2 < vmax_k2)] - vmin_k2) - # / (vmax_k2 - vmin_k2) - # ) - # * (ub_k2 - lb_k2) - # ) - # ) - # - # trans_Vp[trans_Vp <= vmin_Vp] = trans_Vp[trans_Vp <= vmin_Vp] * lb_Vp - # trans_Vp[(trans_Vp >= lim_Vp)] = trans_Vp[trans_Vp >= lim_Vp] * ub_lim - # trans_Vp[(trans_Vp >= vmax_Vp) & (trans_Vp < lim_Vp)] = ( - # trans_Vp[(trans_Vp >= vmax_Vp) & (trans_Vp < lim_Vp)] * ub_Vp - # ) - # trans_Vp[(trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp)] = trans_Vp[ - # (trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp) - # ] * ( - # lb_Vp - # + ( - # ( - # (trans_Vp[(trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp)] - vmin_Vp) - # / (vmax_Vp - vmin_Vp) - # ) - # * (ub_Vp - lb_Vp) - # ) - # ) - # - # trans_Vp[trans_Vp > 1] = 1 - # - # Ctiss_tr = forward_extended_tofts(trans_K1, trans_k2, trans_Vp, Ea, t) - # - # R1_tr, R2st_tr = calcR1_R2(R10, R20st_value, r1, r2st, Ctiss_tr) - # dro_S_tr, M_tr = Conc2Sig_Linear(Tr, Te, fa, R1_tr, R2st_tr, signal[:, :, :8, :], 1, M) - # - # dro_Snoise_tr = addnoise(1, stdS, dro_S_tr, Nt=data_shape[-1]) - # - # animate_mri(signal, mode="time", slice_index=7, time_index=5) - # animate_mri(dro_Snoise_tr, mode="time", slice_index=7, time_index=5) - # animate_mri(dro_Snoise, mode="time", slice_index=7, time_index=5) + dro_s_noise = dro.addnoise(1, stdS, dro_s, Nt=data_shape[-1]) + + trans_Kt = mf_Kt.copy() + trans_Ve = mf_Ve.copy() + trans_Vp = mf_Vp.copy() + + vmax_Kt, vmax_ke, vmax_Vp = 1, 1, 0.2 + vmin_Kt, vmin_ke, vmin_Vp = 0.2, 0.2, 0.01 + lb_Kt, lb_ke, lb_Vp = 0.54, 0.52, 0.49 + ub_Kt, ub_ke, ub_Vp = 1.52, 1.5, 1.43 + lim_Kt, lim_ke, lim_Vp = vmax_Kt + 0.5, vmax_ke + 0.1, vmax_Vp + 0.5 + ub_lim = 1.01 + + trans_Kt[trans_Kt <= vmin_Kt] = trans_Kt[trans_Kt <= vmin_Kt] * lb_Kt + trans_Kt[trans_Kt >= lim_Kt] = trans_Kt[trans_Kt >= lim_Kt] * ub_lim + trans_Kt[(trans_Kt >= vmax_Kt) & (trans_Kt < lim_Kt)] = ( + trans_Kt[(trans_Kt >= vmax_Kt) & (trans_Kt < lim_Kt)] * ub_Kt + ) + trans_Kt[(trans_Kt > vmin_Kt) & (trans_Kt < vmax_Kt)] = trans_Kt[ + (trans_Kt > vmin_Kt) & (trans_Kt < vmax_Kt) + ] * ( + lb_Kt + + ( + ( + (trans_Kt[(trans_Kt > vmin_Kt) & (trans_Kt < vmax_Kt)] - vmin_Kt) + / (vmax_Kt - vmin_Kt) + ) + * (ub_Kt - lb_Kt) + ) + ) + + trans_Vp[trans_Vp > 1] = 1 + + trans_Vp[trans_Vp <= vmin_Vp] = trans_Vp[trans_Vp <= vmin_Vp] * lb_Vp + trans_Vp[(trans_Vp >= lim_Vp)] = trans_Vp[trans_Vp >= lim_Vp] * ub_lim + trans_Vp[(trans_Vp >= vmax_Vp) & (trans_Vp < lim_Vp)] = ( + trans_Vp[(trans_Vp >= vmax_Vp) & (trans_Vp < lim_Vp)] * ub_Vp + ) + trans_Vp[(trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp)] = trans_Vp[ + (trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp) + ] * ( + lb_Vp + + ( + ( + (trans_Vp[(trans_Vp > vmin_Vp) & (trans_Vp < vmax_Vp)] - vmin_Vp) + / (vmax_Vp - vmin_Vp) + ) + * (ub_Vp - lb_Vp) + ) + ) + + c_tissue_tr = dro.forward_extended_tofts(trans_Kt, trans_Ve, trans_Vp, (Ea / aif_cor), t) + + R1_tr, R2st_tr = dro.calcR1_R2(R10, R20st_value, r1, r2st, c_tissue_tr) + dro_s_tr, M_tr = dro.Conc2Sig_Linear(Tr, Te, fa, R1_tr, R2st_tr, signal[:, :, :8, :], 1, M) + + dro_s_noise_tr = dro.addnoise(1, stdS, dro_s_tr, Nt=data_shape[-1]) + + dro.animate_mri(signal, mode="time", slice_index=7, time_index=5) + dro.animate_mri(dro_s_noise_tr, mode="time", slice_index=7, time_index=5) + dro.animate_mri(dro_s_noise, mode="time", slice_index=7, time_index=5) + + save = True + + if save: + dro.save_dicoms("output", slices, dro_s_noise, dro_IDnum, study_instance_uid) From 4df8bcd4bb73daa3f66bdf2e289f5d63ddf808f4 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Thu, 5 Sep 2024 01:58:58 +0300 Subject: [PATCH 11/12] Add save dicoms function --- src/osipi/DRO/_utils.py | 83 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 83 insertions(+) create mode 100644 src/osipi/DRO/_utils.py diff --git a/src/osipi/DRO/_utils.py b/src/osipi/DRO/_utils.py new file mode 100644 index 0000000..048d85a --- /dev/null +++ b/src/osipi/DRO/_utils.py @@ -0,0 +1,83 @@ +import os + +import numpy as np +import pydicom +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 + + +def save_dicoms(outdir, original_dicom, signal, patient_id_num, study_instance_uid): + new_dicoms = original_dicom.copy() + patient_id_num = "RIDER Neuro MRI-{}".format(patient_id_num) + + uidprefix = "1.3.6.1.4.1.9328.50.16." + + SeriesInstanceUID = pydicom.uid.generate_uid(prefix=uidprefix) + StorageMediaFileSetUID = pydicom.uid.generate_uid(prefix=uidprefix) + FrameOfReferenceUID = pydicom.uid.generate_uid(prefix=uidprefix) + SeriesID = SeriesInstanceUID[-5:] + Seriesdir = "{}.000000-perfusion-{}".format(original_dicom[0][0].SeriesNumber, SeriesID) + StudyDate = "19040323" + DateID = study_instance_uid[-5:] + Datedir = "23-03-1904-BRAINRESEARCH-{}".format(DateID) + + if not os.path.exists(outdir): + os.makedirs("{}/{}/{}/{}".format(outdir, patient_id_num, Datedir, Seriesdir)) + + z = 0 + for entries in new_dicoms: + t = 0 + for snap in new_dicoms["{}".format(entries)]: + SOPInstanceUID = pydicom.uid.generate_uid(prefix=uidprefix) + snap.PatientID = patient_id_num + snap.SOPInstanceUID = SOPInstanceUID + snap.StudyInstanceUID = study_instance_uid + snap.SeriesInstanceUID = SeriesInstanceUID + snap.StorageMediaFileSetUID = StorageMediaFileSetUID + snap.FrameOfReferenceUID = FrameOfReferenceUID + snap.StudyDate = StudyDate + snap.ContentDate = StudyDate + Sout = abs(signal[:, :, z, t]) + Sout = Sout.astype(np.uint16) + snap.PixelData = Sout.tobytes() + fname = str(t + 1).zfill(2) + "-" + str(z + 1).zfill(4) + snap.save_as( + "{}/{}/{}/{}/{}.dcm".format(outdir, patient_id_num, Datedir, Seriesdir, fname), + write_like_original=False, + ) + t += 1 + z += 1 + return From 9b7a275666fd6b4d5a5c910e5a255cc60d7a8c02 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Thu, 5 Sep 2024 02:15:20 +0300 Subject: [PATCH 12/12] solving save function bug --- src/osipi/DRO/_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/osipi/DRO/_utils.py b/src/osipi/DRO/_utils.py index 048d85a..5b88e4e 100644 --- a/src/osipi/DRO/_utils.py +++ b/src/osipi/DRO/_utils.py @@ -60,7 +60,7 @@ def save_dicoms(outdir, original_dicom, signal, patient_id_num, study_instance_u z = 0 for entries in new_dicoms: t = 0 - for snap in new_dicoms["{}".format(entries)]: + for snap in new_dicoms[z]: SOPInstanceUID = pydicom.uid.generate_uid(prefix=uidprefix) snap.PatientID = patient_id_num snap.SOPInstanceUID = SOPInstanceUID