Merge pull request #117 from baej3/logan

add logan reference tissue model [WIP]

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)](<https://doi.org/10.1016/S1053-8119(03)00017-X>)
 
+[lrtm]: https://doi.org/10.1097/00004647-199609000-00008
 [hypr-lr]: https://doi.org/10.2967/jnumed.109.073999
 
 ## Requirements

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "dynamicpet"
-version = "0.1.1"
+version = "0.1.2"
 description = "Dynamic PET"
 authors = ["Murat Bilgel <[email protected]>"]
 license = "MIT"

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()

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)

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)