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