Merge pull request #133 from bilgelm/pet-bids-derivatives

Align outputs with PET-BIDS derivatives

Author: bilgelm

README.md

@@ -1,4 +1,4 @@
-# Dynamic PET <img src="docs/logo.png" align="right" height="150" />
+# Dynamic PET <img src="docs/logo.png" align="right" height="150" alt="Dynamic PET logo"/>
 
 <!-- [![PyPI](https://img.shields.io/pypi/v/dynamicpet.svg)][pypi_]
 [![Status](https://img.shields.io/pypi/status/dynamicpet.svg)][status]
@@ -34,6 +34,7 @@ Methods implemented in the CLI include:
 
 - Denoising
   - HighlY constrained backPRojection method constraining the backprojections to Local Regions of interest ([HYPR-LR])
+  - Nonlocal EStimation of multispectral MAgnitudes ([NESMA])
 - Reference tissue-based modeling
   - Standardized Uptake Value Ratio (SUVR)
   - Logan Reference Tissue Model ([LRTM])
@@ -47,6 +48,7 @@ Several implementations of estimating SRTM parameters are available:
 
 [lrtm]: https://doi.org/10.1097/00004647-199609000-00008
 [hypr-lr]: https://doi.org/10.2967/jnumed.109.073999
+[nesma]: https://doi.org/10.1111/jon.12537
 
 ## Requirements
 
@@ -60,8 +62,8 @@ _Dynamic PET_ requires Python 3.11+ and the following modules:
 You can install _Dynamic PET_ via [pip] after cloning the repository:
 
 ```console
-$ git clone https://github.com/bilgelm/dynamicpet.git
-$ pip install -e dynamicpet
+git clone https://github.com/bilgelm/dynamicpet.git
+pip install -e dynamicpet
 ```
 
 ## Usage
@@ -78,13 +80,13 @@ You will then need to create a binary mask that is in the same space as the PET
 After installing _Dynamic PET_ as described above, execute:
 
 ```console
-$ kineticmodel PET --model SRTMZhou2003 --refmask <REFMASK> --outputdir <OUTPUTDIR> --fwhm 5
+kineticmodel PET --model SRTMZhou2003 --refmask <REFMASK> --outputdir <OUTPUTDIR> --fwhm 5
 ```
 
 where
 
 ```console
-$ PET=<OPENNEURODATA>/ds001705-download/sub-000101/ses-baseline/pet/sub-000101_ses-baseline_pet.nii
+PET=<OPENNEURODATA>/ds001705-download/sub-000101/ses-baseline/pet/sub-000101_ses-baseline_pet.nii
 ```
 
 Before running these commands, replace

pyproject.toml

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

src/dynamicpet/__init__.py

@@ -1 +1,6 @@
 """Dynamic PET."""
+
+import importlib.metadata
+
+
+__version__ = importlib.metadata.version("dynamicpet")

src/dynamicpet/__main__.py

@@ -2,14 +2,17 @@
 
 import csv
 import os
+import re
 import warnings
+from json import dump as json_dump
 
 import click
 import numpy as np
 from nibabel.filename_parser import splitext_addext
 from nibabel.loadsave import load as nib_load
 from nibabel.spatialimages import SpatialImage
 
+from dynamicpet import __version__
 from dynamicpet.denoise import hypr
 from dynamicpet.denoise import nesma
 from dynamicpet.kineticmodel.kineticmodel import KineticModel
@@ -100,19 +103,31 @@ def denoise(
 ) -> None:
     """Perform dynamic PET denoising.
 
-    Outputs will have a '_<method>' suffix.
+    Outputs will have a '_desc-<method>_pet' suffix.
 
     PET: 4-D PET image
     """
     # load PET
     pet_img = petbidsimage_load(pet, json)
 
     if method == "HYPRLR":
+        doc = hypr.hypr_lr.__doc__
+
+        # parameters not relevant to HYPR-LR
+        mask = None
+        window_half_size = None
+        thresh = None
+
         if fwhm:
             res = hypr.hypr_lr(pet_img, fwhm)
         else:
             raise ValueError("fwhm must be specified for HYPR-LR")
     elif method == "NESMA":
+        doc = nesma.nesma_semiadaptive.__doc__
+
+        # parameter not relevant to NESMA
+        fwhm = None
+
         if mask:
             mask_img: SpatialImage = nib_load(mask)  # type: ignore
             # check that mask is in the same space as pet
@@ -138,10 +153,37 @@ def denoise(
     if outputdir:
         os.makedirs(outputdir, exist_ok=True)
         bname = os.path.basename(froot)
-        froot = os.path.join(outputdir, bname)
-    output = froot + "_" + method.lower() + ext + addext
+        # if the input file name follows the PET-BIDS convention, it should end
+        # with "_pet". Need to move this to the end of the new file name to
+        # maintain compatibility with the PET-BIDS Derivatives convention.
+        froot = re.sub("_pet$", "", os.path.join(outputdir, bname))
+    output = froot + "_desc-" + method.lower() + "_pet" + ext + addext
+    output_json = froot + "_desc-" + method.lower() + "_pet.json"
     res.to_filename(output)
 
+    cmd = (
+        f"denoise --method {method} "
+        + (f"--fwhm {fwhm} " if fwhm else "")
+        + (f"--mask {mask} " if mask else "")
+        + (f"--window_half_size {window_half_size} " if window_half_size else "")
+        + (f"--thresh {thresh} " if thresh else "")
+        + (f"--outputdir {outputdir} " if outputdir else "")
+        + (f"--json {json} " if json else "")
+        + pet
+    )
+    derivative_json_dict = {
+        "Description": (
+            re.sub(r"\s+", " ", doc.split("Args:")[0]).strip() if doc else ""
+        ),
+        # "Sources": [pet],
+        "SoftwareName": "dynamicpet",
+        "SoftwareVersion": __version__,
+        "CommandLine": cmd,
+    }
+
+    with open(output_json, "w") as f:
+        json_dump(derivative_json_dict, f, indent=4)
+
 
 @click.command()
 @click.argument("pet", type=str)
@@ -200,9 +242,11 @@ def denoise(
     ),
 )
 @click.option(
-    "--start", default=None, type=float, help="Start of time window for model"
+    "--start", default=None, type=float, help="Start of time window for model in min"
+)
+@click.option(
+    "--end", default=None, type=float, help="End of time window for model in min"
 )
-@click.option("--end", default=None, type=float, help="End of time window for model")
 @click.option(
     "--fwhm",
     default=None,
@@ -226,7 +270,7 @@ def denoise(
         "'rect' is rectangular integration."
     ),
 )
-def kineticmodel(
+def kineticmodel(  # noqa: C901
     pet: str,
     model: str,
     refroi: str | None,
@@ -242,7 +286,7 @@ def kineticmodel(
 ) -> None:
     """Fit a reference tissue model to a dynamic PET image or TACs.
 
-    Outputs will have a '_km-<model>_kp-<parameter>' suffix.
+    Outputs will have a '_model-<model>_meas-<parameter>' suffix.
 
     PET: 4-D PET image (can be 3-D if model is SUVR) or a 2-D tabular TACs tsv file
     """
@@ -262,16 +306,27 @@ def kineticmodel(
         pet_img = pet_img.extract(start, end)
         reftac = reftac.extract(start, end)
 
+    if fwhm and model not in ["srtmzhou2003"]:
+        fwhm = None
+        warnings.warn(
+            "--fwhm argument is not relevant for this model, will be ignored",
+            RuntimeWarning,
+            stacklevel=2,
+        )
+
     # fit kinetic model
     km: KineticModel
     match model:
         case "suvr":
+            model_abbr = "SUVR"
             km = SUVR(reftac, pet_img)
             km.fit(mask=petmask_img_mat)
         case "srtmlammertsma1996":
+            model_abbr = "SRTM"
             km = SRTMLammertsma1996(reftac, pet_img)
             km.fit(mask=petmask_img_mat, weight_by=weight_by)
         case "srtmzhou2003":
+            model_abbr = "SRTM"
             km = SRTMZhou2003(reftac, pet_img)
             km.fit(
                 mask=petmask_img_mat,
@@ -288,7 +343,13 @@ def kineticmodel(
         froot = os.path.join(outputdir, bname)
 
     if isinstance(pet_img, PETBIDSMatrix):
-        output = froot + "_km-" + model.replace(".", "") + ext
+        # if the input file name follows the PET-BIDS Derivatives convention,
+        # it should end with "_tacs". Need to remove this to maintain
+        # compatibility with the PET-BIDS Derivatives convention.
+        froot = re.sub("_tacs$", "", froot)
+
+        output = froot + "_model-" + model_abbr + "_kinpar" + ext
+        output_json = froot + "_model-" + model_abbr + "_kinpar.json"
         data = np.empty((len(km.parameters), pet_img.num_elements))
         for i, param in enumerate(km.parameters.keys()):
             data[i] = km.get_parameter(param)
@@ -299,16 +360,73 @@ def kineticmodel(
             for i, elem in enumerate(pet_img.elem_names):
                 tsvwriter.writerow([elem] + datat[i].tolist())
     else:
+        # if the input file name follows the PET-BIDS convention, it should end
+        # with "_pet". Need to remove this to maintain compatibility with the
+        # PET-BIDS Derivatives convention.
+        froot = re.sub("_pet$", "", froot)
+
         # save estimated parameters as image
         for param in km.parameters.keys():
             res_img: SpatialImage = km.get_parameter(param)  # type: ignore
             output = (
-                froot + "_km-" + model.replace(".", "") + "_kp-" + param + ext + addext
+                froot
+                + "_model-"
+                + model_abbr
+                + "_meas-"
+                + param
+                + "_mimap"
+                + ext
+                + addext
             )
             res_img.to_filename(output)
-
-    # also need to save a json PET BIDS derivative file
-    # TODO
+        output_json = froot + "_model-" + model_abbr + "_mimap.json"
+
+    # save json PET BIDS derivative file
+    inputvalues = [start, end]
+    inputvalueslabels = [
+        "Start of time window for model",
+        "End of time window for model",
+    ]
+    inputvaluesunits = ["min", "min"]
+
+    if fwhm:
+        inputvalues += [fwhm]
+        inputvalueslabels += ["Full width at half max"]
+        inputvaluesunits += ["mm"]
+
+    cmd = (
+        f"kineticmodel --model {model} "
+        + (f"--refroi {refroi} " if refroi else f"--refmask {refmask} ")
+        + (f"--outputdir {outputdir} " if outputdir else "")
+        + (f"--json {json} " if json else "")
+        + (f"--petmask {petmask} " if petmask else "")
+        + f"--start {start} "
+        + f"--end {end} "
+        + (f"--fwhm {fwhm} " if fwhm else "")
+        + f"--weight_by {weight_by} "
+        + f"--integration_type {integration_type} "
+        + pet
+    )
+    doc = km.__class__.__doc__
+    derivative_json_dict = {
+        "Description": re.sub(r"\s+", " ", doc) if doc else "",
+        # "Sources": [pet],
+        "ModelName": model_abbr,
+        "ReferenceRegion": refroi if refroi else refmask,
+        "AdditionalModelDetails": (
+            f"Frame weighting by: {weight_by}. "
+            + f"Integration type: {integration_type}.",
+        ),
+        "InputValues": inputvalues,
+        "InputValuesLabels": inputvalueslabels,
+        "InputValuesUnits": inputvaluesunits,
+        "SoftwareName": "dynamicpet",
+        "SoftwareVersion": __version__,
+        "CommandLine": cmd,
+    }
+
+    with open(output_json, "w") as f:
+        json_dump(derivative_json_dict, f, indent=4)
 
 
 def parse_kineticmodel_inputs(

src/dynamicpet/kineticmodel/kineticmodel.py

@@ -35,6 +35,7 @@ class KineticModel(ABC):
     @abstractmethod
     def get_param_names(cls) -> list[str]:
         """Get names of kinetic model parameters."""
+        # parameter names should contain alphanumeric characters only
         raise NotImplementedError
 
     def __init__(
@@ -110,11 +111,11 @@ def get_parameter(self, param_name: str) -> SpatialImage | NumpyNumberArray:
             else:
                 param_vector: NumpyNumberArray = self.parameters[param_name]
                 return param_vector
-        elif param_name == "bp" and "dvr" in self.parameters:
-            self.parameters[param_name] = self.parameters["dvr"] - 1
+        elif param_name == "BPND" and "DVR" in self.parameters:
+            self.parameters[param_name] = self.parameters["DVR"] - 1
             return self.get_parameter(param_name)
-        elif param_name == "dvr" and "bp" in self.parameters:
-            self.parameters[param_name] = self.parameters["bp"] + 1
+        elif param_name == "DVR" and "BPND" in self.parameters:
+            self.parameters[param_name] = self.parameters["BPND"] + 1
             return self.get_parameter(param_name)
         else:
             raise AttributeError(

src/dynamicpet/kineticmodel/logan.py

@@ -30,15 +30,15 @@ class LRTM(KineticModel):
     @classmethod
     def get_param_names(cls) -> list[str]:
         """Get names of kinetic model parameters."""
-        return ["dvr"]
+        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,
+        k2prime: float | None = None,
     ) -> None:
         """Estimate model parameters.
 
@@ -56,7 +56,8 @@ def fit(  # noqa: max-complexity: 12
                   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
+            k2prime: (avg.) effective tissue-to-plasma efflux constant in the
+                     reference region, in unit of 1/min
         """
         # get reference TAC as a 1-D vector
         reftac: NumpyNumberArray = self.reftac.dataobj.flatten()[:, np.newaxis]
@@ -82,7 +83,7 @@ def fit(  # noqa: max-complexity: 12
 
         dvr = np.zeros((num_elements, 1))
 
-        if not k2:
+        if not k2prime:
             # TODO
             # Check Eq. 7 assumption (i.e., that tac / reftac is reasonably
             # constant) by calculating R2 etc. for each tac.
@@ -103,13 +104,14 @@ def fit(  # noqa: max-complexity: 12
 
             # ----- 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
+            # - use Eq. 6 if k2prime is provided
+            # - use Eq. 7 if k2prime is not provided
 
             x = np.column_stack(
                 (
                     np.ones_like(tac_tstar),
-                    (int_reftac_tstar + (reftac_tstar / k2 if k2 else 0)) / tac_tstar,
+                    (int_reftac_tstar + (reftac_tstar / k2prime if k2prime else 0))
+                    / tac_tstar,
                 )
             )
             y = int_tac_tstar / tac_tstar
@@ -123,8 +125,8 @@ def fit(  # noqa: max-complexity: 12
             # distribution volume ratio
             dvr[k] = b[1]
 
-        self.set_parameter("dvr", dvr, mask)
-        # should tstar (and k2?) also be stored?
+        self.set_parameter("DVR", dvr, mask)
+        # should tstar (and k2prime?) also be stored?
 
     def fitted_tacs(self) -> TemporalMatrix | TemporalImage:
         """Get fitted TACs based on estimated model parameters."""