diff --git a/docs/user-guide/gen_aif.md b/docs/user-guide/gen_aif.md index 8143615..599ed6e 100644 --- a/docs/user-guide/gen_aif.md +++ b/docs/user-guide/gen_aif.md @@ -9,6 +9,8 @@ import osipi t = np.arange(0, 6*60, 1) ca = osipi.aif_parker(t) plt.plot(t, ca) +plt.xlabel('Time (s)') +plt.ylabel('Indicator concentration (mM)') plt.show() ``` @@ -24,13 +26,55 @@ ca = osipi.aif_parker(t) Ktrans = 0.6 ve = 0.2 ct = osipi.tofts(t, ca, Ktrans=Ktrans/60, ve=ve) -plt.plot(t, ct) +fig, ax = plt.subplots(1, 2) +ax[0].plot(t, ca) +ax[0].set_xlabel('Time (s)') +ax[0].set_ylabel('Indicator concentration (mM)') +ax[0].set_title('AIF') +ax[1].plot(t, ct) +ax[1].set_xlabel('Time (s)') +ax[1].set_ylabel('Indicator concentration (mM)') +ax[1].set_title('Tissue') + +fig.tight_layout() plt.show() ``` ## Generating an MRI signal -!!! note "Coming Soon" - This section is under development and will be available soon. + +``` py +import numpy as np +import matplotlib.pyplot as plt +import osipi + +t = np.arange(0, 6*60, 1) +ca = osipi.aif_parker(t) +Ktrans = 0.6 +ve = 0.2 +ct = osipi.tofts(t, ca, Ktrans=Ktrans/60, ve=ve) + +R10 = 0.5 +r1 = 4.5 +R1t = osipi.C_to_R1_linear_relaxivity(ct, R10, r1) + +S0 = 1 +TR = 0.004 +a = 12 +St = osipi.signal_SPGR(R1t, S0, TR, a) + +fig, ax = plt.subplots(1, 2) +ax[0].plot(t, ct) +ax[0].set_xlabel('Time (s)') +ax[0].set_ylabel('Indicator concentration (mM)') +ax[0].set_title('Concentration') +ax[1].plot(t, St) +ax[1].set_xlabel('Time (s)') +ax[1].set_ylabel('MRI signal (a.u)') +ax[1].set_title('MRI signal') + +fig.tight_layout() +plt.show() +``` ## Adding measurement error !!! note "Coming Soon" diff --git a/src/osipi/__init__.py b/src/osipi/__init__.py index 454bb7e..68753a2 100644 --- a/src/osipi/__init__.py +++ b/src/osipi/__init__.py @@ -17,6 +17,10 @@ from ._signal_to_concentration import ( S_to_C_via_R1_SPGR, - S_to_R1_SPGR, - R1_to_C_linear_relaxivity + S_to_R1_SPGR +) + +from ._electromagnetic_property import ( + R1_to_C_linear_relaxivity, + C_to_R1_linear_relaxivity ) diff --git a/src/osipi/_electromagnetic_property.py b/src/osipi/_electromagnetic_property.py index a746c97..a14447b 100644 --- a/src/osipi/_electromagnetic_property.py +++ b/src/osipi/_electromagnetic_property.py @@ -38,3 +38,37 @@ def R1_to_C_linear_relaxivity( elif not (r1 >= 0): raise ValueError("r1 must be positive") return (R1 - R10) / r1 # C + + +def C_to_R1_linear_relaxivity( + C: NDArray[np.float64], R10: np.float64, r1: np.float64 +) -> NDArray[np.float64]: + """ + Electromagnetic property forward model: + - longitudinal relaxation rate, linear with relaxivity + + Converts tissue concentration to R1 + + Args: + C (1D array of np.float64): + Vector of indicator concentrations in units of mM. [OSIPI code Q.IC1.001] + R10 (np.float64): + Native longitudinal relaxation rate in units of /s. [OSIPI code Q.EL1.002] + r1 (np.float64): + Longitudinal relaxivity in units of /s/mM. [OSIPI code Q.EL1.015] + + Returns: + NDArray[np.float64]: + Vector of longitudinal relaxation rate in units of /s. [OSIPI code Q.EL1.001] + + References: + - Lexicon URL: https://osipi.github.io/OSIPI_CAPLEX/perfusionModels/# + - Lexicon code: M.EL1.003 + - Adapted from equation given in lexicon + """ + # Check C is a 1D array of floats + if not (isinstance(C, np.ndarray) and C.ndim == 1 and C.dtype == np.float64): + raise TypeError("C must be a 1D NumPy array of np.float64") + elif not (r1 >= 0): + raise ValueError("r1 must be positive") + return R10 + r1 * C # R1