diff --git a/README.md b/README.md index 1b8a6a3..c58a9db 100644 --- a/README.md +++ b/README.md @@ -36,6 +36,7 @@ Methods implemented in the CLI include: - HighlY constrained backPRojection method constraining the backprojections to Local Regions of interest ([HYPR-LR]) - Reference tissue-based modeling - Standardized Uptake Value Ratio (SUVR) + - Logan Reference Tissue Model ([LRTM]) - Simplified Reference Tissue Model (SRTM) Several implementations of estimating SRTM parameters are available: @@ -44,6 +45,7 @@ Several implementations of estimating SRTM parameters are available: - Two-step SRTM (SRTM2), as described by [Wu and Carson, J Cereb Blood Flow Metab (2002)](https://doi.org/10.1097/01.WCB.0000033967.83623.34) - Linear Regression with Spatial Constraint (LRSC), as described by [Zhou et al., Neuroimage (2003)]() +[lrtm]: https://doi.org/10.1097/00004647-199609000-00008 [hypr-lr]: https://doi.org/10.2967/jnumed.109.073999 ## Requirements diff --git a/pyproject.toml b/pyproject.toml index 119f00d..ad41c44 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "dynamicpet" -version = "0.1.1" +version = "0.1.2" description = "Dynamic PET" authors = ["Murat Bilgel "] license = "MIT" diff --git a/src/dynamicpet/kineticmodel/logan.py b/src/dynamicpet/kineticmodel/logan.py new file mode 100644 index 0000000..a631c92 --- /dev/null +++ b/src/dynamicpet/kineticmodel/logan.py @@ -0,0 +1,131 @@ +"""Logan plot.""" + +import numpy as np +from numpy.linalg import LinAlgError +from scipy.linalg import solve # type: ignore +from tqdm import trange + +from ..temporalobject.temporalimage import TemporalImage +from ..temporalobject.temporalmatrix import TemporalMatrix +from ..temporalobject.temporalobject import INTEGRATION_TYPE_OPTS +from ..temporalobject.temporalobject import WEIGHT_OPTS +from ..typing_utils import NumpyNumberArray +from .kineticmodel import KineticModel + + +class LRTM(KineticModel): + """Logan reference tissue model. + + Also known as Logan plot, graphical Logan, Logan graphical analysis (or + some permutation thereof these words) with reference tissue. + "Non-invasive" can also be used instead of "with reference tissue" to convey + the same meaning. + + Reference: + Logan J, Fowler JS, Volkow ND, Wang GJ, Ding YS, Alexoff DL. Distribution + volume ratios without blood sampling from graphical analysis of PET data. + J Cereb Blood Flow Metab. 1996 Sep;16(5):834-40. + """ + + @classmethod + def get_param_names(cls) -> list[str]: + """Get names of kinetic model parameters.""" + return ["dvr"] + + def fit( # noqa: max-complexity: 12 + self, + mask: NumpyNumberArray | None = None, + integration_type: INTEGRATION_TYPE_OPTS = "trapz", + weight_by: WEIGHT_OPTS | NumpyNumberArray | None = "frame_duration", + tstar: float = 0, + k2: float | None = None, + ) -> None: + """Estimate model parameters. + + Args: + integration_type: If 'rect', rectangular integration is used for TACs. + If 'trapz', trapezoidal integration is used based + on middle timepoint of each frame. + weight_by: [optional] frame weights used in model fitting. + If weight_by == None, each frame is weighted equally. + If weight_by == 'frame_duration', each frame is weighted + proportionally to its duration (inverse variance weighting). + If weight_by is a 1-D array, then specified values are used. + mask: [optional] A 1-D (for TemporalMatrix TACs) or + 3-D (for TemporalImage TACs) binary mask that defines where + to fit the kinetic model. Elements outside the mask will + be set to to 0 in parametric estimate outputs. + tstar: time beyond which to assume linearity + k2: (avg.) effective tissue-to-plasma efflux constant, in unit of 1/min + """ + # get reference TAC as a 1-D vector + reftac: NumpyNumberArray = self.reftac.dataobj.flatten()[:, np.newaxis] + # numerical integration of reference TAC + int_reftac: NumpyNumberArray = self.reftac.cumulative_integral( + integration_type + ).flatten()[:, np.newaxis] + + tacs: TemporalMatrix = self.tacs.timeseries_in_mask(mask) + num_elements = tacs.num_elements + tacs_mat: NumpyNumberArray = tacs.dataobj + int_tacs_mat: NumpyNumberArray = tacs.cumulative_integral(integration_type) + + # time indexing should be done after integrating + t_idx = tacs.frame_start >= tstar + reftac_tstar = reftac[t_idx, :] + int_reftac_tstar = int_reftac[t_idx, :] + tacs_mat_tstar = tacs_mat[:, t_idx] + int_tacs_mat_tstar = int_tacs_mat[:, t_idx] + + weights = tacs.get_weights(weight_by) + w_star = np.diag(weights[t_idx]) + + dvr = np.zeros((num_elements, 1)) + + if not k2: + # TODO + # Check Eq. 7 assumption (i.e., that tac / reftac is reasonably + # constant) by calculating R2 etc. for each tac. + # Display warning if assumption is off. + pass + + for k in trange(num_elements): + # get TAC and its cumulative integral as 1-D vectors + tac_tstar = tacs_mat_tstar[k, :][:, np.newaxis] + + # special case when tac is the same as reftac + if np.allclose(tac_tstar, reftac_tstar): + dvr[k] = 1 + + continue + + int_tac_tstar = int_tacs_mat_tstar[k, :][:, np.newaxis] + + # ----- Get DVR ----- + # Set up the weighted linear regression model based on Logan et al.: + # - use Eq. 6 if k2 is provided + # - use Eq. 7 if k2 is not provided + + x = np.column_stack( + ( + np.ones_like(tac_tstar), + (int_reftac_tstar + (reftac_tstar / k2 if k2 else 0)) / tac_tstar, + ) + ) + y = int_tac_tstar / tac_tstar + + b: NumpyNumberArray + try: + b = solve(x.T @ w_star @ x, x.T @ w_star @ y, assume_a="sym") + except LinAlgError: + b = np.ones((2, 1)) + + # distribution volume ratio + dvr[k] = b[1] + + self.set_parameter("dvr", dvr, mask) + # should tstar (and k2?) also be stored? + + def fitted_tacs(self) -> TemporalMatrix | TemporalImage: + """Get fitted TACs based on estimated model parameters.""" + raise NotImplementedError() diff --git a/tests/test_kineticmodel.py b/tests/test_kineticmodel.py index 6778997..1221f8d 100644 --- a/tests/test_kineticmodel.py +++ b/tests/test_kineticmodel.py @@ -6,6 +6,7 @@ from nibabel.spatialimages import SpatialImage from numpy.typing import NDArray +from dynamicpet.kineticmodel.logan import LRTM from dynamicpet.kineticmodel.srtm import SRTMZhou2003 from dynamicpet.kineticmodel.suvr import SUVR from dynamicpet.temporalobject import TemporalImage @@ -103,3 +104,12 @@ def test_srtm_zhou2003_ti(reftac: TemporalMatrix, tacs_img: TemporalImage) -> No bp_img: SpatialImage = km.get_parameter("bp") # type: ignore assert np.allclose(dvr_img.get_fdata(), bp_img.get_fdata() + 1) + + +def test_logan_tm(reftac: TemporalMatrix, tacs_img: TemporalImage) -> None: + """Test Logan Plot using TemporalImage.""" + km = LRTM(reftac, tacs_img) + km.fit(integration_type="trapz") + + dvr_img: SpatialImage = km.get_parameter("dvr") # type: ignore + assert dvr_img.shape == (1, 1, 2) diff --git a/tests/test_kineticmodel_simulated.py b/tests/test_kineticmodel_simulated.py index 4e40a71..75cb1ef 100644 --- a/tests/test_kineticmodel_simulated.py +++ b/tests/test_kineticmodel_simulated.py @@ -9,6 +9,7 @@ from numpy.typing import NDArray from scipy.integrate import odeint # type: ignore +from dynamicpet.kineticmodel.logan import LRTM from dynamicpet.kineticmodel.srtm import SRTMLammertsma1996 from dynamicpet.kineticmodel.srtm import SRTMZhou2003 from dynamicpet.temporalobject import TemporalImage @@ -268,3 +269,28 @@ def test_srtmlammertsma1996_tm( assert np.allclose(bp, bp_true, rtol=relative_tol) assert np.allclose(r1, r1_true, rtol=relative_tol) assert np.allclose(fitted_tacs.dataobj, tac.dataobj, rtol=relative_tol) + + +def test_lrtm_tm( + frame_start: NDArray[np.double], + frame_duration: NDArray[np.double], +) -> None: + """Test Logan 1996 calculation using TemporalMatrix.""" + bp_true = 1.5 + r1_true = 1.2 + ct, cref = get_tacs_and_reftac_dataobj( + frame_start, frame_duration, bp_true, r1_true + ) + + tac = TemporalMatrix(ct, frame_start, frame_duration) + reftac = TemporalMatrix(cref, frame_start, frame_duration) + + km = LRTM(reftac, tac) + # because of the way data were simulated, we need to use trapz integration + # to fit to obtain the best possible recovery of true values + km.fit(integration_type="trapz") + + bp: NumpyRealNumberArray = km.get_parameter("bp") # type: ignore + + relative_tol = 0.02 # .02 means that 2% error is tolerated + assert np.allclose(bp, bp_true, rtol=relative_tol)