automatically use the base z stat image when the target image is the cluster corrected map

Author: jdkent

nimare/diagnostics.py

@@ -2,6 +2,7 @@
 
 import copy
 import logging
+import warnings
 from abc import abstractmethod
 
 import nibabel as nib
@@ -69,6 +70,124 @@ def _get_target_value_map(result):
     )
 
 
+def _resolve_target_threshold(target_threshold, voxel_thresh):
+    """Resolve diagnostics threshold aliases."""
+    if target_threshold is not None and voxel_thresh is not None:
+        raise ValueError(
+            "Only one of 'target_threshold' and deprecated 'voxel_thresh' may be provided."
+        )
+
+    if voxel_thresh is not None:
+        warnings.warn(
+            "'voxel_thresh' is deprecated for diagnostics and will be removed in a future "
+            "release. Use 'target_threshold' to threshold the selected target image before "
+            "diagnostics table/support generation.",
+            FutureWarning,
+            stacklevel=3,
+        )
+        return voxel_thresh
+
+    return target_threshold
+
+
+def _is_cluster_corrected_target(target_image):
+    """Determine whether a target image is a corrected cluster-level map."""
+    return "_level-cluster" in target_image and "_corr-" in target_image
+
+
+def _remove_cluster_stat_suffix(description):
+    """Remove size/mass cluster-stat suffixes from a map description."""
+    if description in {"size", "mass"}:
+        return None
+
+    for suffix in ("Size", "Mass", "size", "mass"):
+        if description.endswith(suffix):
+            description = description[: -len(suffix)]
+            break
+
+    return description or None
+
+
+def _candidate_peak_value_maps(target_image):
+    """Generate candidate original statistic maps for a corrected cluster target."""
+    uncorrected_target = target_image.split("_corr-", 1)[0]
+    candidate_maps = []
+
+    if "_desc-" in uncorrected_target:
+        description = uncorrected_target.split("_desc-", 1)[1].split("_", 1)[0]
+        description = _remove_cluster_stat_suffix(description)
+        if description is not None:
+            candidate_maps.extend([f"z_desc-{description}", f"stat_desc-{description}"])
+
+    candidate_maps.extend(
+        ["z", "stat", "est", "stat_desc-group1MinusGroup2", "z_desc-association"]
+    )
+    return candidate_maps
+
+
+def _get_peak_value_map_for_cluster_table(result, target_image):
+    """Select the original map used for peak statistics in corrected-cluster tables."""
+    for candidate_map in _candidate_peak_value_maps(target_image):
+        if "_level-cluster" not in candidate_map and candidate_map in result.maps:
+            return candidate_map
+
+    available_maps = ", ".join(sorted(result.maps.keys()))
+    raise ValueError(
+        "No supported original z/statistic map found for corrected-cluster table peaks. "
+        f"Target image was '{target_image}'. Available maps are: {available_maps}."
+    )
+
+
+def _get_cluster_support_data(label_maps, shape):
+    """Convert one or more cluster label maps to a binary support array."""
+    support = np.zeros(shape, dtype=bool)
+    for label_map in label_maps:
+        support |= np.asanyarray(label_map.dataobj) > 0
+    return support
+
+
+def _get_clusters_table_and_label_maps(
+    result,
+    target_img,
+    target_image,
+    threshold,
+    cluster_threshold,
+):
+    """Create diagnostics clusters from target image and peaks from original statistics."""
+    target_data = target_img.get_fdata(dtype=DEFAULT_FLOAT_DTYPE)
+    if hasattr(result.estimator, "two_sided"):
+        # Only present in Fisher's and Stouffer's estimators
+        two_sided = getattr(result.estimator, "two_sided")
+    else:
+        two_sided = (target_data < 0).any()
+
+    clusters_table, label_maps = get_clusters_table(
+        target_img,
+        0 if threshold is None else threshold,
+        cluster_threshold,
+        two_sided=two_sided,
+        return_label_maps=True,
+    )
+
+    if clusters_table.empty or not label_maps or not _is_cluster_corrected_target(target_image):
+        return clusters_table, label_maps
+
+    peak_value_map = _get_peak_value_map_for_cluster_table(result, target_image)
+    peak_img = result.get_map(peak_value_map, return_type="image")
+    peak_data = peak_img.get_fdata(dtype=DEFAULT_FLOAT_DTYPE)
+    support_data = _get_cluster_support_data(label_maps, peak_data.shape)
+    masked_peak_data = np.where(support_data, peak_data, 0).astype(DEFAULT_FLOAT_DTYPE, copy=False)
+    masked_peak_img = nib.Nifti1Image(masked_peak_data, peak_img.affine, peak_img.header)
+
+    peak_clusters_table = get_clusters_table(
+        masked_peak_img,
+        0,
+        0,
+        two_sided=two_sided or (masked_peak_data < 0).any(),
+    )
+    return peak_clusters_table, label_maps
+
+
 def _cluster_masker_kwargs():
     """Return standardized kwargs for label-based cluster summaries."""
     return _filter_kwargs(
@@ -210,17 +329,20 @@ class Diagnostics(NiMAREBase):
         The meta-analytic map for which clusters will be characterized.
         The default is z because log-p will not always have value of zero for non-cluster voxels.
     voxel_thresh : :obj:`float` or None, optional
-        An optional voxel-level threshold that may be applied to the ``target_image`` to define
-        clusters. This can be None if the ``target_image`` is already thresholded
-        (e.g., a cluster-level corrected map).
-        Default is None.
+        Deprecated alias for ``target_threshold``. Prefer ``target_threshold`` for new code.
     cluster_threshold : :obj:`int` or None, optional
         Cluster size threshold, in :term:`voxels<voxel>`.
         If None, then no cluster size threshold will be applied. Default=None.
     n_cores : :obj:`int`, optional
         Number of cores to use for parallelization.
         If <=0, defaults to using all available cores.
         Default is 1.
+    target_threshold : :obj:`float` or None, optional
+        Threshold applied to ``target_image`` before defining diagnostics clusters and tables.
+        For unthresholded Monte Carlo cluster-corrected maps, this should generally be the
+        corrected significance threshold in the target map's units. This is distinct from
+        :class:`~nimare.correct.FWECorrector` ``voxel_thresh``, which is the cluster-forming
+        threshold used during correction. Default=None.
 
     """
 
@@ -231,9 +353,11 @@ def __init__(
         cluster_threshold=None,
         display_second_group=False,
         n_cores=1,
+        target_threshold=None,
     ):
         self.target_image = target_image
         self.voxel_thresh = voxel_thresh
+        self.target_threshold = _resolve_target_threshold(target_threshold, voxel_thresh)
         self.cluster_threshold = cluster_threshold
         self.display_second_group = display_second_group
         self.n_cores = _check_ncores(n_cores)
@@ -312,21 +436,14 @@ def transform(self, result):
             )
 
         # Get clusters table and label maps
-        stat_threshold = self.voxel_thresh or 0
         cluster_threshold = 0 if self.cluster_threshold is None else self.cluster_threshold
 
-        if hasattr(result.estimator, "two_sided"):
-            # Only present in Fisher's and Stouffer's estimators
-            two_sided = getattr(result.estimator, "two_sided")
-        else:
-            two_sided = (target_img.get_fdata(dtype=DEFAULT_FLOAT_DTYPE) < 0).any()
-
-        clusters_table, label_maps = get_clusters_table(
+        clusters_table, label_maps = _get_clusters_table_and_label_maps(
+            result,
             target_img,
-            stat_threshold,
+            self.target_image,
+            self.target_threshold,
             cluster_threshold,
-            two_sided=two_sided,
-            return_label_maps=True,
         )
 
         n_clusters = clusters_table.shape[0]

nimare/tests/test_diagnostics.py

@@ -6,6 +6,7 @@
 
 import nibabel as nib
 import numpy as np
+import pandas as pd
 import pytest
 from nilearn.maskers import NiftiLabelsMasker
 
@@ -92,6 +93,100 @@ def test_get_target_value_map_raises_for_unsupported_maps():
         diagnostics._get_target_value_map(result)
 
 
+def test_diagnostics_voxel_thresh_deprecated_alias():
+    """voxel_thresh should remain a deprecated alias for target_threshold."""
+    with pytest.warns(FutureWarning, match="voxel_thresh"):
+        counter = diagnostics.FocusCounter(target_image="z", voxel_thresh=1.0)
+
+    assert counter.target_threshold == 1.0
+
+
+def test_diagnostics_target_and_voxel_threshold_error():
+    """Supplying both threshold names should fail explicitly."""
+    with pytest.raises(ValueError, match="target_threshold"):
+        diagnostics.FocusCounter(target_image="z", target_threshold=1.0, voxel_thresh=1.0)
+
+
+def test_corrected_cluster_table_uses_thresholded_support_and_original_z():
+    """Corrected cluster tables should report original-z peaks inside thresholded support."""
+    target_image = "z_desc-size_level-cluster_corr-FWE_method-montecarlo"
+    mask_data = np.ones((7, 7, 7), dtype=bool)
+    mask_img = nib.Nifti1Image(mask_data.astype(np.int8), affine=np.eye(4))
+
+    class DummyMasker:
+        def __init__(self, mask_img):
+            self.mask_img_ = mask_img
+
+        def transform(self, img):
+            return np.asanyarray(img.dataobj)[mask_data].reshape(1, -1)
+
+        def inverse_transform(self, data):
+            arr = np.asarray(data)
+            if arr.ndim > 1:
+                arr = np.squeeze(arr, axis=0)
+            out = np.zeros(mask_data.shape, dtype=arr.dtype)
+            out[mask_data] = arr
+            return nib.Nifti1Image(out, affine=mask_img.affine)
+
+    class DummyResult:
+        def __init__(self, maps, masker):
+            self.maps = maps
+            self.tables = {}
+            self.diagnostics = []
+            self.masker = masker
+            self.estimator = SimpleNamespace(
+                masker=masker,
+                inputs_={
+                    "id": ["study1"],
+                    "coordinates": pd.DataFrame(
+                        {"id": ["study1"], "x": [1.0], "y": [1.0], "z": [1.0]}
+                    ),
+                },
+            )
+
+        def get_map(self, name, return_type="image"):
+            values = self.maps[name]
+            if return_type == "array":
+                return values
+            return self.masker.inverse_transform(values)
+
+    corrected = np.zeros(mask_data.shape, dtype=float)
+    corrected[(1, 1, 1)] = 6.0
+    corrected[(1, 1, 2)] = 2.0
+    corrected[(1, 2, 1)] = 2.0
+    corrected[(5, 5, 5)] = 1.2
+    corrected[(5, 5, 4)] = 1.2
+
+    original_z = np.zeros(mask_data.shape, dtype=float)
+    original_z[(1, 1, 1)] = 3.0
+    original_z[(1, 1, 2)] = 9.0
+    original_z[(1, 2, 1)] = 2.0
+    original_z[(5, 5, 5)] = 8.0
+    original_z[(5, 5, 4)] = 8.0
+    original_z[(3, 3, 3)] = 99.0
+
+    masker = DummyMasker(mask_img)
+    result = DummyResult(
+        {
+            target_image: masker.transform(nib.Nifti1Image(corrected, affine=mask_img.affine))[0],
+            "z": masker.transform(nib.Nifti1Image(original_z, affine=mask_img.affine))[0],
+        },
+        masker,
+    )
+
+    counter = diagnostics.FocusCounter(target_image=target_image, target_threshold=1.64)
+    result = counter.transform(result)
+
+    clusters_table = result.tables[f"{target_image}_tab-clust"]
+    peak_row = clusters_table.loc[clusters_table["Peak Stat"].idxmax()]
+
+    assert clusters_table.shape[0] == 1
+    assert peak_row["Peak Stat"] == pytest.approx(9.0)
+    assert peak_row[["X", "Y", "Z"]].to_numpy().tolist() == [1.0, 1.0, 2.0]
+    assert not np.any(np.isclose(clusters_table["Peak Stat"], 8.0))
+    assert not np.any(np.isclose(clusters_table["Peak Stat"], 99.0))
+
+
 def test_is_voxelwise_masker_uses_round_trip_when_mask_count_mismatches():
     """Voxelwise detection should fall back to a round-trip feature-shape check."""
     mask_data = np.array([[[1], [0]], [[1], [1]]], dtype=np.int8)
@@ -150,7 +245,7 @@ def test_jackknife_smoke(
     testdata = testdata_ibma if meta_type == "ibma" else testdata_cbma_full
     res = meta.fit(dset1, dset2) if n_samples == "twosample" else meta.fit(testdata)
 
-    jackknife = diagnostics.Jackknife(target_image=target_image, voxel_thresh=voxel_thresh)
+    jackknife = diagnostics.Jackknife(target_image=target_image, target_threshold=voxel_thresh)
     results = jackknife.transform(res)
 
     image_name = "_".join(target_image.split("_")[1:])
@@ -177,7 +272,7 @@ def test_jackknife_with_zero_clusters(testdata_cbma_full):
     meta = cbma.ALE()
     res = meta.fit(testdata_cbma_full)
 
-    jackknife = diagnostics.Jackknife(target_image="z", voxel_thresh=10)
+    jackknife = diagnostics.Jackknife(target_image="z", target_threshold=10)
     results = jackknife.transform(res)
 
     contribution_table = results.tables["z_diag-Jackknife_tab-counts"]
@@ -200,14 +295,14 @@ def test_jackknife_with_custom_masker_smoke(testdata_ibma):
     meta = ibma.SampleSizeBasedLikelihood(mask=masker)
     res = meta.fit(testdata_ibma)
 
-    jackknife = diagnostics.Jackknife(target_image="z", voxel_thresh=0.5)
+    jackknife = diagnostics.Jackknife(target_image="z", target_threshold=0.5)
     results = jackknife.transform(res)
     contribution_table = results.tables["z_diag-Jackknife_tab-counts_tail-positive"]
     assert contribution_table.shape[0] == len(meta.inputs_["id"])
 
     # A Jackknife with a target_image that isn't present in the MetaResult raises a ValueError.
     with pytest.raises(ValueError):
-        jackknife = diagnostics.Jackknife(target_image="doggy", voxel_thresh=0.5)
+        jackknife = diagnostics.Jackknife(target_image="doggy", target_threshold=0.5)
         jackknife.transform(res)
 
 
@@ -227,7 +322,7 @@ def test_focuscounter_negative_tail_label_map_naming(testdata_cbma_full):
     neg_img = nib.Nifti1Image(neg_data, mask_img.affine)
     res.maps["z"] = np.squeeze(masker.transform(neg_img))
 
-    counter = diagnostics.FocusCounter(target_image="z", voxel_thresh=1.0)
+    counter = diagnostics.FocusCounter(target_image="z", target_threshold=1.0)
     results = counter.transform(res)
 
     assert "label_tail-negative" in results.maps
@@ -252,7 +347,7 @@ def test_focuscounter_positive_tail_label_map_naming(testdata_cbma_full):
     pos_img = nib.Nifti1Image(pos_data, mask_img.affine)
     res.maps["z"] = np.squeeze(masker.transform(pos_img))
 
-    counter = diagnostics.FocusCounter(target_image="z", voxel_thresh=1.0)
+    counter = diagnostics.FocusCounter(target_image="z", target_threshold=1.0)
     results = counter.transform(res)
 
     assert "label_tail-positive" in results.maps
@@ -283,7 +378,7 @@ def _fake_infer(_label_maps, _clusters_table, _n_clusters):
     monkeypatch.setattr(diagnostics, "_infer_label_map_tails", _fake_infer)
     caplog.set_level(logging.WARNING, logger="nimare.diagnostics")
 
-    counter = diagnostics.FocusCounter(target_image="z", voxel_thresh=1.0)
+    counter = diagnostics.FocusCounter(target_image="z", target_threshold=1.0)
     results = counter.transform(res)
 
     assert any("Mixed-sign clusters detected" in r.message for r in caplog.records)
@@ -308,7 +403,7 @@ def test_focuscounter_pairwise_negative_tail_uses_group2(testdata_cbma_full):
     neg_img = nib.Nifti1Image(neg_data, mask_img.affine)
     res.maps["z_desc-uniformity"] = np.squeeze(masker.transform(neg_img))
 
-    counter = diagnostics.FocusCounter(target_image="z_desc-uniformity", voxel_thresh=1.0)
+    counter = diagnostics.FocusCounter(target_image="z_desc-uniformity", target_threshold=1.0)
     results = counter.transform(res)
 
     table_key = "z_desc-uniformity_diag-FocusCounter_tab-counts_tail-negative"
@@ -342,7 +437,7 @@ def test_focuscounter_smoke(
     testdata = testdata_ibma if meta_type == "ibma" else testdata_cbma_full
     res = meta.fit(dset1, dset2) if n_samples == "twosample" else meta.fit(testdata)
 
-    counter = diagnostics.FocusCounter(target_image=target_image, voxel_thresh=1.65)
+    counter = diagnostics.FocusCounter(target_image=target_image, target_threshold=1.65)
     if meta_type == "ibma":
         with pytest.raises(ValueError):
             counter.transform(res)