From d26288b54b2e3a2a9c51fca12963ad38011b8b79 Mon Sep 17 00:00:00 2001 From: JoeBae Date: Wed, 11 Sep 2024 11:58:42 -0400 Subject: [PATCH 1/4] Add loganplot --- src/dynamicpet/kineticmodel/logan.py | 116 +++++++++++++++++++++++++++ tests/test_kineticmodel.py | 11 +++ 2 files changed, 127 insertions(+) create mode 100644 src/dynamicpet/kineticmodel/logan.py diff --git a/src/dynamicpet/kineticmodel/logan.py b/src/dynamicpet/kineticmodel/logan.py new file mode 100644 index 0000000..7f582ff --- /dev/null +++ b/src/dynamicpet/kineticmodel/logan.py @@ -0,0 +1,116 @@ +"""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 ..typing_utils import RealNumber +from .kineticmodel import KineticModel + + +class LoganPlot(KineticModel): + """Logan Plot. + + Reference: + Logan, J., Fowler, J. S., Volkow, N. D., Wang, G. J., + Ding, Y. S., & Alexoff, D. L. (1996). + Distribution volume ratios without blood sampling + from graphical analysis of PET data. + Journal of Cerebral Blood Flow & Metabolism, 16(5), 834-840. + """ + + @classmethod + def get_param_names(cls) -> list[str]: + """Get names of kinetic model parameters.""" + return [ + "dvr", + # "r1", + # "k2", + # "k2a", + # "r1_lrsc", + # "k2_lrsc", + # "k2a_lrsc", + "noise_var_eq_dvr", + # "noise_var_eq_r1", + ] + + 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", + fwhm: RealNumber | list[RealNumber] | 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. + fwhm: scalar or length 3 sequence, FWHM in mm over which to smooth + """ + # 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() + + 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) + + weights = tacs.get_weights(weight_by) + w = np.diag(weights) + + dvr = np.zeros((num_elements, 1)) + + for k in trange(num_elements): + # get TAC and its cumulative integral as 1-D vectors + tac = tacs_mat[k, :][:, np.newaxis] + + # special case when tac is the same as reftac + if np.allclose(tac, reftac): + dvr[k] = 1 + + continue + + int_tac = int_tacs_mat[k, :][:, np.newaxis] + + # ----- Get DVR ----- + # Set up the weighted linear regression model based on Eq. 7 in Logan et al. + + x = np.column_stack((np.ones(tac.shape), np.divide(int_reftac, tac))) + y = np.divide(int_tac, tac) + + b: NumpyNumberArray + try: + b = solve(x.T @ w @ x, x.T @ w @ y, assume_a="sym") + except LinAlgError: + b = np.ones((2, 1)) + + # distribution volume ratio + dvr[k] = b[1] + + self.set_parameter("dvr", dvr, mask) + + def fitted_tacs(self) -> TemporalMatrix | TemporalImage: + """Get fitted TACs based on estimated model parameters.""" + # there is no parametric model for SUVR, so we just return the tacs + return self.tacs diff --git a/tests/test_kineticmodel.py b/tests/test_kineticmodel.py index 6778997..974302a 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 LoganPlot from dynamicpet.kineticmodel.srtm import SRTMZhou2003 from dynamicpet.kineticmodel.suvr import SUVR from dynamicpet.temporalobject import TemporalImage @@ -103,3 +104,13 @@ 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_loganplot(reftac: TemporalMatrix, tacs_img: TemporalImage) -> None: + """Test Logan Plot using TemporalImage.""" + km = LoganPlot(reftac, tacs_img) + km.fit(integration_type="trapz") + # fitted_tacs = km.fitted_tacs() + + dvr_img: SpatialImage = km.get_parameter("dvr") # type: ignore + assert dvr_img.shape == (1, 1, 2) From 83a3fa0db090bfb7b90841915c6ecda8b2a02e88 Mon Sep 17 00:00:00 2001 From: Murat Bilgel Date: Thu, 12 Sep 2024 21:50:01 -0400 Subject: [PATCH 2/4] add k2 and tstar to logan, add sim test --- src/dynamicpet/kineticmodel/logan.py | 81 ++++++++++++++++------------ tests/test_kineticmodel.py | 7 ++- tests/test_kineticmodel_simulated.py | 26 +++++++++ 3 files changed, 77 insertions(+), 37 deletions(-) diff --git a/src/dynamicpet/kineticmodel/logan.py b/src/dynamicpet/kineticmodel/logan.py index 7f582ff..a631c92 100644 --- a/src/dynamicpet/kineticmodel/logan.py +++ b/src/dynamicpet/kineticmodel/logan.py @@ -10,42 +10,35 @@ from ..temporalobject.temporalobject import INTEGRATION_TYPE_OPTS from ..temporalobject.temporalobject import WEIGHT_OPTS from ..typing_utils import NumpyNumberArray -from ..typing_utils import RealNumber from .kineticmodel import KineticModel -class LoganPlot(KineticModel): - """Logan Plot. +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, J. S., Volkow, N. D., Wang, G. J., - Ding, Y. S., & Alexoff, D. L. (1996). - Distribution volume ratios without blood sampling - from graphical analysis of PET data. - Journal of Cerebral Blood Flow & Metabolism, 16(5), 834-840. + 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", - # "r1", - # "k2", - # "k2a", - # "r1_lrsc", - # "k2_lrsc", - # "k2a_lrsc", - "noise_var_eq_dvr", - # "noise_var_eq_r1", - ] + 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", - fwhm: RealNumber | list[RealNumber] | None = None, + tstar: float = 0, + k2: float | None = None, ) -> None: """Estimate model parameters. @@ -62,46 +55,68 @@ def fit( # noqa: max-complexity: 12 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. - fwhm: scalar or length 3 sequence, FWHM in mm over which to smooth + 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() + ).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 = np.diag(weights) + 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 = tacs_mat[k, :][:, np.newaxis] + tac_tstar = tacs_mat_tstar[k, :][:, np.newaxis] # special case when tac is the same as reftac - if np.allclose(tac, reftac): + if np.allclose(tac_tstar, reftac_tstar): dvr[k] = 1 continue - int_tac = int_tacs_mat[k, :][:, np.newaxis] + int_tac_tstar = int_tacs_mat_tstar[k, :][:, np.newaxis] # ----- Get DVR ----- - # Set up the weighted linear regression model based on Eq. 7 in Logan et al. - - x = np.column_stack((np.ones(tac.shape), np.divide(int_reftac, tac))) - y = np.divide(int_tac, tac) + # 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 @ x, x.T @ w @ y, assume_a="sym") + b = solve(x.T @ w_star @ x, x.T @ w_star @ y, assume_a="sym") except LinAlgError: b = np.ones((2, 1)) @@ -109,8 +124,8 @@ def fit( # noqa: max-complexity: 12 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.""" - # there is no parametric model for SUVR, so we just return the tacs - return self.tacs + raise NotImplementedError() diff --git a/tests/test_kineticmodel.py b/tests/test_kineticmodel.py index 974302a..1221f8d 100644 --- a/tests/test_kineticmodel.py +++ b/tests/test_kineticmodel.py @@ -6,7 +6,7 @@ from nibabel.spatialimages import SpatialImage from numpy.typing import NDArray -from dynamicpet.kineticmodel.logan import LoganPlot +from dynamicpet.kineticmodel.logan import LRTM from dynamicpet.kineticmodel.srtm import SRTMZhou2003 from dynamicpet.kineticmodel.suvr import SUVR from dynamicpet.temporalobject import TemporalImage @@ -106,11 +106,10 @@ def test_srtm_zhou2003_ti(reftac: TemporalMatrix, tacs_img: TemporalImage) -> No assert np.allclose(dvr_img.get_fdata(), bp_img.get_fdata() + 1) -def test_loganplot(reftac: TemporalMatrix, tacs_img: TemporalImage) -> None: +def test_logan_tm(reftac: TemporalMatrix, tacs_img: TemporalImage) -> None: """Test Logan Plot using TemporalImage.""" - km = LoganPlot(reftac, tacs_img) + km = LRTM(reftac, tacs_img) km.fit(integration_type="trapz") - # fitted_tacs = km.fitted_tacs() 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) From 109ed78edc4d4fc6b835ddcbdb6aa4ac889c546a Mon Sep 17 00:00:00 2001 From: Murat Bilgel Date: Thu, 12 Sep 2024 21:50:17 -0400 Subject: [PATCH 3/4] add logan reference tissue model --- README.md | 2 ++ 1 file changed, 2 insertions(+) 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 From 671a6f570bc83ba3e9116dde3ab499379b2b2faf Mon Sep 17 00:00:00 2001 From: Murat Bilgel Date: Thu, 12 Sep 2024 22:04:50 -0400 Subject: [PATCH 4/4] bump up dynamicpet version --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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"