Add NESMA; update deprecated numpy dependencies

src/dynamicpet/__main__.py

@@ -11,6 +11,7 @@
 from nibabel.spatialimages import SpatialImage
 
 from dynamicpet.denoise import hypr
+from dynamicpet.denoise import nesma
 from dynamicpet.kineticmodel.kineticmodel import KineticModel
 from dynamicpet.kineticmodel.srtm import SRTMLammertsma1996
 from dynamicpet.kineticmodel.srtm import SRTMZhou2003
@@ -41,13 +42,44 @@
 
 @click.command()
 @click.argument("pet", type=str)
-@click.argument("fwhm", type=float)
 @click.option(
     "--method",
-    type=click.Choice(["HYPRLR"], case_sensitive=False),
+    type=click.Choice(["HYPRLR", "NESMA"], case_sensitive=False),
     required=True,
     help="Name of denoising method",
 )
+@click.option(
+    "--fwhm",
+    default=None,
+    type=float,
+    help="Full width at half max in mm for smoothing (only for HYPR-LR)",
+)
+@click.option(
+    "--mask",
+    default=None,
+    type=str,
+    help=(
+        "Binary mask specifying voxels where denoising should be performed. "
+        "(only for NESMA)"
+    ),
+)
+@click.option(
+    "--window_half_size",
+    default=None,
+    nargs=3,
+    type=int,
+    help=(
+        "The size of search window centered around a voxel will be "
+        "2 * window_half_size + 1 (any part of the window that extends beyond "
+        "the image will be truncated). (only for NESMA)"
+    ),
+)
+@click.option(
+    "--thresh",
+    default=0.05,
+    type=click.FloatRange(0, 1),
+    help="threshold (only for NESMA)",
+)
 @click.option(
     "--outputdir",
     default=None,
@@ -59,21 +91,48 @@
 )
 @click.option("--json", default=None, type=str, help="PET-BIDS json file")
 def denoise(
-    pet: str, fwhm: float, method: str, outputdir: str | None, json: str | None
+    pet: str,
+    method: str,
+    fwhm: float | None,
+    mask: str | None,
+    window_half_size: tuple[int, int, int] | None,
+    thresh: float | None,
+    outputdir: str | None,
+    json: str | None,
 ) -> None:
     """Perform dynamic PET denoising.
 
     Outputs will have a '_<method>' suffix.
 
     PET: 4-D PET image
-
-    FWHM: full width at half max, in mm, for smoothing filter
     """
     # load PET
     pet_img = petbidsimage_load(pet, json)
 
     if method == "HYPRLR":
-        res = hypr.hypr_lr(pet_img, fwhm)
+        if fwhm:
+            res = hypr.hypr_lr(pet_img, fwhm)
+        else:
+            raise ValueError("fwhm must be specified for HYPR-LR")
+    elif method == "NESMA":
+        if mask:
+            mask_img: SpatialImage = nib_load(mask)  # type: ignore
+            # check that mask is in the same space as pet
+            if not np.all(pet_img.img.affine == mask_img.affine):
+                raise ValueError("PET and mask are not in the same space")
+            mask_img_mat: NumpyRealNumberArray = mask_img.get_fdata().astype("bool")
+
+            if window_half_size:
+                if thresh:
+                    res, _ = nesma.nesma_semiadaptive(
+                        pet_img, mask_img_mat, window_half_size, thresh
+                    )
+                else:
+                    raise ValueError("thresh must be specified for NESMA")
+            else:
+                raise ValueError("window_half_size must be specified for NESMA")
+        else:
+            raise ValueError("mask must be specified for NESMA")
     else:
         raise NotImplementedError(f"Denoising method {method} is not supported")
 
@@ -82,7 +141,7 @@ def denoise(
         os.makedirs(outputdir, exist_ok=True)
         bname = os.path.basename(froot)
         froot = os.path.join(outputdir, bname)
-    output = froot + "_hyprlr" + ext + addext
+    output = froot + "_" + method.lower() + ext + addext
     res.to_filename(output)
 
 

src/dynamicpet/denoise/nesma.py

@@ -0,0 +1,96 @@
+"""NESMA denoising."""
+
+import numpy as np
+from tqdm import tqdm
+
+from ..petbids.petbidsimage import PETBIDSImage
+from ..temporalobject.temporalimage import image_maker
+from ..typing_utils import NumpyRealNumberArray
+
+
+def nesma_semiadaptive(
+    ti: PETBIDSImage,
+    mask: NumpyRealNumberArray,
+    window_half_size: tuple[int, int, int],
+    thresh: float = 0.05,
+) -> tuple[PETBIDSImage, NumpyRealNumberArray]:
+    """NESMA denoising.
+
+    NESMA is short for Nonlocal EStimation of multispectral MAgnitudes.
+
+    Reference:
+    Bouhrara M, Reiter DA, Maring MC, Bonny JM, Spencer RG. Use of the NESMA
+    Filter to Improve Myelin Water Fraction Mapping with Brain MRI.
+    J Neuroimaging. 2018;28(6):640–9. https://doi.org/10.1111/jon.12537
+
+    Args:
+        ti: 4-D image
+        mask: 3-D boolean mask indicating where smoothing should be performed
+        window_half_size: the size of search window centered around a voxel will
+                          be 2 * window_half_size + 1 (any part of the window
+                          that extends beyond the image will be truncated)
+        thresh: a value between 0 and 1, optional (default = 0.05)
+
+    Returns:
+        s_nesma: filtered output image, 4-D
+        vsm: voxel similarity map, 3-D indicates the number of voxels in the
+             search window that were below the threshold
+
+    Raises:
+        ValueError: incorrect input(s)
+    """
+    if thresh < 0 or thresh > 1:
+        raise ValueError("thresh must be between 0 and 1")
+    if ti.img.ndim != 4:
+        raise ValueError("input image must be 4-D")
+    if mask.ndim != 3:
+        raise ValueError("mask must be 3-D")
+    if not np.all(ti.img.shape[:-1] == mask.shape):
+        raise ValueError("input image and mask shapes must be compatible")
+
+    m, n, o, _ = ti.img.shape
+
+    if m <= window_half_size[0] or n <= window_half_size[1] or o <= window_half_size[2]:
+        raise ValueError("image size should be greater than window_half_size")
+
+    s_nesma = np.copy(ti.dataobj)
+    vsm = np.zeros(ti.shape[:-1])
+
+    indices = np.where(mask)
+
+    for i, j, k in tqdm(zip(*indices, strict=True), total=np.sum(mask)):  # type: ignore
+        tmin = max(k - window_half_size[2], 0)
+        tmax = min(k + window_half_size[2] + 1, o)
+
+        rmin = max(i - window_half_size[0], 0)
+        rmax = min(i + window_half_size[0] + 1, m)
+
+        # index voxel signal vector
+        x = ti.dataobj[i : i + 1, j : j + 1, k : k + 1, :]
+
+        smin = max(j - window_half_size[1], 0)
+        smax = min(j + window_half_size[1] + 1, n)
+
+        # get the neighborhood around voxel
+        nbhd = ti.dataobj[rmin:rmax, smin:smax, tmin:tmax, :]
+
+        # Relative Manhattan distance calculation
+        #
+        # In the NESMA paper, the denominator includes the L1 norm of x only.
+        # Here, we use the sum of the L1 norm of x and the L1 norm of the voxel
+        # x is being compared to (+ a small number to prevent division by 0).
+        #
+        # This choice of denominator is motivated by the inequality that
+        # |y - x| <= |x| + |y|
+        # This way, the threshold is bounded above by 1.
+        distance = np.sum(np.abs(nbhd - x), axis=-1) / (
+            np.sum(np.abs(x)) + np.sum(np.abs(nbhd), axis=-1) + np.finfo(float).eps
+        )
+
+        pos = distance < thresh
+        n_below_thresh = np.sum(pos)
+        if n_below_thresh > 0:
+            s_nesma[i, j, k, :] = np.sum(nbhd[pos, :], axis=0) / n_below_thresh
+            vsm[i, j, k] = n_below_thresh
+
+    return PETBIDSImage(image_maker(s_nesma, ti.img), ti.json_dict), vsm

src/dynamicpet/kineticmodel/kinfitr.py

@@ -64,7 +64,9 @@ def fit(
                         param_estimates.update(
                             {param_name: np.zeros((num_elements, 1))}
                         )
-                    param_estimates[param_name][i] = res["par"][param_name]
+                    param_estimates[param_name][i] = res["par"][  # noqa: B909
+                        param_name
+                    ]
 
         for param_name, param_estimate in param_estimates.items():
             self.set_parameter(param_name, param_estimate, mask)

src/dynamicpet/kineticmodel/srtm.py

@@ -301,7 +301,7 @@ def fit(  # noqa: max-complexity: 12
                 # based on Eq. 11 in Zhou et al.
                 x = np.column_stack((reftac, int_reftac, -int_tac))
                 h_d = np.diag(h[k, :])
-                b_sc = np.row_stack(
+                b_sc = np.vstack(
                     (smooth_r1_mat[k], smooth_k2_mat[k], smooth_k2a_mat[k])
                 )
                 try:

src/dynamicpet/temporalobject/temporalobject.py

@@ -10,6 +10,7 @@
 
 import numpy as np
 from scipy.integrate import cumulative_trapezoid  # type: ignore
+from scipy.integrate import trapezoid
 
 from ..typing_utils import NumpyRealNumberArray
 from ..typing_utils import RealNumber
@@ -329,7 +330,7 @@ def _dynamic_mean(
                     RuntimeWarning,
                     stacklevel=2,
                 )
-            dyn_mean = np.trapz(self.dataobj, self.frame_mid) / (
+            dyn_mean = trapezoid(self.dataobj, self.frame_mid) / (
                 self.frame_mid[-1] - self.frame_mid[0]
             )
         else:

tests/test_kinfitr.py

@@ -38,7 +38,7 @@ def simref0() -> TACPair:
     # this should not affect kinfitr functions as they will add the 0 back
     reftac_tm = TemporalMatrix(reftac[1:], frame_start[1:], frame_duration[1:])
     tac_tm = TemporalMatrix(
-        np.row_stack((tac1[1:], tac2[1:], tac3[1:])),
+        np.vstack((tac1[1:], tac2[1:], tac3[1:])),
         frame_start[1:],
         frame_duration[1:],
     )

tests/test_main.py

@@ -119,7 +119,7 @@ def images() -> dict[str, Path]:
     return fnames
 
 
-def test_denoise(images: dict[str, Path]) -> None:
+def test_denoise_hyprlr(images: dict[str, Path]) -> None:
     """Test denoise in __main__.py."""
     from dynamicpet.__main__ import denoise
 
@@ -129,7 +129,15 @@ def test_denoise(images: dict[str, Path]) -> None:
     runner = CliRunner()
     result = runner.invoke(
         denoise,
-        ["--method", "HYPRLR", str(pet_fname), "5.0", "--outputdir", str(outputdir)],
+        [
+            "--method",
+            "HYPRLR",
+            str(pet_fname),
+            "--fwhm",
+            "5.0",
+            "--outputdir",
+            str(outputdir),
+        ],
         catch_exceptions=False,
     )
 
@@ -138,6 +146,41 @@ def test_denoise(images: dict[str, Path]) -> None:
     assert os.path.isfile(outputdir / "pet_hyprlr.json")
 
 
+@pytest.mark.skip(reason="current NESMA implementation is too slow")
+def test_denoise_nesma(images: dict[str, Path]) -> None:
+    """Test denoise in __main__.py."""
+    from dynamicpet.__main__ import denoise
+
+    pet_fname = images["pet_fname"]
+    petmask_fname = images["petmask_fname"]
+    outputdir = pet_fname.parent / "test_output" / "nesma"
+
+    runner = CliRunner()
+    result = runner.invoke(
+        denoise,
+        [
+            "--method",
+            "NESMA",
+            str(pet_fname),
+            "--mask",
+            str(petmask_fname),
+            "--window_half_size",
+            "3",
+            "3",
+            "3",
+            "--thresh",
+            "0.05",
+            "--outputdir",
+            str(outputdir),
+        ],
+        catch_exceptions=False,
+    )
+
+    assert result.exit_code == 0
+    assert os.path.isfile(outputdir / "pet_nesma.nii")
+    assert os.path.isfile(outputdir / "pet_nesma.json")
+
+
 def test_kineticmodel_suvr(images: dict[str, Path]) -> None:
     """Test SUVR kineticmodel in __main__.py."""
     from dynamicpet.__main__ import kineticmodel