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)