Add -normalize-PAM50 argument to metric extraction script to allow interpolation to PAM50 space

spinalcordtoolbox/metrics_to_PAM50.py

@@ -33,10 +33,19 @@ def interpolate_metrics(metrics, fname_vert_levels_PAM50, fname_vert_levels):
     level_slices_PAM50 = [get_slices_from_vertebral_levels(im_seg_labeled_PAM50, level) for level in levels]
     level_slices_im = [get_slices_from_vertebral_levels(im_seg_labeled, level) for level in levels]
 
-    # Find the mean scaling between the image and PAM50 (excluding first and last levels)
-    scales = [len(slices_PAM50)/len(slices_im) for slices_PAM50, slices_im
-              in zip(level_slices_PAM50[1:-1], level_slices_im[1:-1])]
-    scale_mean = np.mean(scales)
+    # Find the mean scaling between the image and PAM50 (excluding first and last levels to avoid
+    # edge effects from potentially incomplete levels)
+    inner_scales = [len(slices_PAM50)/len(slices_im) for slices_PAM50, slices_im
+                    in zip(level_slices_PAM50[1:-1], level_slices_im[1:-1])
+                    if len(slices_im) > 0]
+    if inner_scales:
+        scale_mean = np.mean(inner_scales)
+    else:
+        # Fall back to all available levels when fewer than 3 levels are present
+        # This is done to pass test on mt/mtr.nii.gz and PAM50_levels.nii.gz that only has 2 levels (C4, C5)
+        all_scales = [len(slices_PAM50)/len(slices_im) for slices_PAM50, slices_im
+                      in zip(level_slices_PAM50, level_slices_im) if len(slices_im) > 0]
+        scale_mean = np.mean(all_scales) if all_scales else np.nan
 
     # Initialize a metrics dict filled by NaN with number of rows equal to number of slices in PAM50 template
     z = im_seg_labeled_PAM50.dim[2]  # z == number of slices

spinalcordtoolbox/scripts/sct_extract_metric.py

@@ -26,6 +26,8 @@
 from spinalcordtoolbox.utils.sys import init_sct, printv, __data_dir__, set_loglevel
 from spinalcordtoolbox.utils.fs import check_file_exist, extract_fname, get_absolute_path, TempFolder
 from spinalcordtoolbox.scripts import sct_maths
+from spinalcordtoolbox.metrics_to_PAM50 import interpolate_metrics
+from spinalcordtoolbox.template import get_vertebral_level_from_slice
 
 
 class Param:
@@ -206,6 +208,19 @@ def get_parser():
         """),
     )
 
+    optional.add_argument(
+        '-normalize-PAM50',
+        metavar=Metavar.int,
+        type=int,
+        choices=[0, 1],
+        default=0,
+        help="Set to 1 to interpolate the extracted metric values into PAM50 anatomical "
+             "dimensions per slice. "
+             "Requires `-vertfile` and `-perslice 1`. "
+             "Inspired by: Valošek J, Bédard S et al. Imaging Neuroscience 2024. "
+             "https://doi.org/10.1162/imag_a_00075"
+    )
+
     advanced = parser.add_argument_group("FOR ADVANCED USERS")
     advanced.add_argument(
         '-param',
@@ -264,6 +279,70 @@ def get_parser():
     return parser
 
 
+def _build_pam50_agg_metric(agg_metric_native, nz_native, label_name, method,
+                            fname_vert_level, fname_vert_level_PAM50):
+    """
+    Interpolate per-slice native-space metrics into PAM50 anatomical dimensions.
+    :param agg_metric_native: per-slice native-space metrics (dict output of extract_metric() with perslice=True)
+    :param nz_native: int: total z-slices in native image
+    :param label_name: str: atlas label name (for 'Label' CSV column), e.g., 'white matter'
+    :param method: str: extraction method ('wa', 'ml', 'map', 'bin', 'median', 'max')
+    :param fname_vert_level: str: native vertebral levels file (centerline-masked)
+    :param fname_vert_level_PAM50: str: PAM50 template PAM50_levels.nii.gz
+    :return: dict keyed by (z,) PAM50 slice tuples, suitable for save_as_csv()
+    """
+    method_key_map = {
+        'wa': 'WA()', 'ml': 'ML()', 'map': 'MAP()',
+        'bin': 'BIN()', 'median': 'MEDIAN()', 'max': 'MAX()'
+    }
+    primary_key = method_key_map[method]
+
+    # Build 1D ndarray with per slice values from the agg_metric_native dict
+    metric_1d = np.full(nz_native, np.nan)
+    for (z,), entry in agg_metric_native.items():
+        val = entry.get(primary_key)
+        if val is not None:
+            metric_1d[z] = val
+
+    # Interpolate to PAM50 space; returns Dict[str, Metric] with 1D data of length z_PAM50
+    metrics_pam50 = interpolate_metrics(
+        {primary_key: Metric(data=metric_1d, label=primary_key)},
+        fname_vert_level_PAM50,
+        fname_vert_level
+    )
+
+    # Make the interpolate_metrics output compatible with the expected input of save_as_csv()
+    # Get 1D ndarray with per-slice metric values in PAM50 space
+    pam50_values = metrics_pam50[primary_key].data
+
+    # Determine which vertebral levels are present in the native data to filter PAM50 output
+    im_native_levels = Image(fname_vert_level).change_orientation('RPI')
+    native_levels = set(
+        int(v) for v in np.unique(im_native_levels.data) if 0 < int(v) < 49
+    )
+
+    # Map each PAM50 z-slice to a vertebral level and build the output agg_metric
+    im_pam50_levels = Image(fname_vert_level_PAM50).change_orientation('RPI')
+
+    agg_metric_pam50 = {}
+    for z_pam50, val in enumerate(pam50_values):
+        # nan means that there was no data in the native space for that slice, so we skip it in the PAM50 space as well
+        if np.isnan(val):
+            continue
+        vert_level = get_vertebral_level_from_slice(im_pam50_levels, z_pam50)
+        if vert_level is None or vert_level not in native_levels:
+            continue
+        entry = {
+            'Label': label_name,
+            'VertLevel': (vert_level,),
+            'DistancePMJ': None,    # required by save_as_csv() but not relevant for PAM50 space
+            primary_key: val,
+        }
+        agg_metric_pam50[(z_pam50,)] = entry
+
+    return agg_metric_pam50
+
+
 def main(argv: Sequence[str]):
     # Ensure that the "-list-labels" argument is always parsed last. That way, if `-f` is passed, then `-list-labels`
     # will see the new location and look there. (https://github.com/spinalcordtoolbox/spinalcordtoolbox/issues/3634)
@@ -289,6 +368,7 @@ def main(argv: Sequence[str]):
     fname_vert_level = arguments.vertfile
     perslice = arguments.perslice
     perlevel = arguments.perlevel
+    normalize_pam50 = arguments.normalize_PAM50
 
     # check if path_label is a file (e.g., single binary mask) instead of a folder (e.g., SCT atlas structure which
     # contains info_label.txt file)
@@ -384,6 +464,8 @@ def main(argv: Sequence[str]):
                f"To use vertebral level information, you may need to run "
                f"`sct_warp_template` to generate the appropriate level file in your working directory.", type=message_type)
         fname_vert_level = None
+        if normalize_pam50:
+            parser.error("Option '-normalize-PAM50' requires a valid '-vertfile'.")
     # Get dimensions of data and labels
     nx, ny, nz = data.data.shape
     nx_atlas, ny_atlas, nz_atlas, nt_atlas = labels.shape
@@ -409,11 +491,32 @@ def main(argv: Sequence[str]):
                                              map_cluster=None)
                 labels_id_user = [99]
 
+    if normalize_pam50 and not perslice:
+        parser.error("Option '-normalize-PAM50' requires option '-perslice 1'.")
+
+    fname_vert_level_PAM50 = os.path.join(__data_dir__, 'PAM50', 'template', 'PAM50_levels.nii.gz')
+
     for id_label in labels_id_user:
         printv('Estimation for label: ' + label_struc[id_label].name, verbose)
-        agg_metric = extract_metric(data, labels=labels, slices=slices, levels=levels, perslice=perslice,
-                                    perlevel=perlevel, fname_vert_level=fname_vert_level, method=method,
-                                    label_struc=label_struc, id_label=id_label, indiv_labels_ids=indiv_labels_ids)
+
+        if normalize_pam50:
+            agg_metric_native = extract_metric(
+                data, labels=labels, slices=slices, levels=levels,
+                perslice=1, perlevel=0,
+                fname_vert_level=fname_vert_level, method=method,
+                label_struc=label_struc, id_label=id_label, indiv_labels_ids=indiv_labels_ids)
+            agg_metric = _build_pam50_agg_metric(
+                agg_metric_native=agg_metric_native,
+                nz_native=nz,
+                label_name=label_struc[id_label].name,
+                method=method,
+                fname_vert_level=fname_vert_level,
+                fname_vert_level_PAM50=fname_vert_level_PAM50)
+        else:
+            agg_metric = extract_metric(
+                data, labels=labels, slices=slices, levels=levels, perslice=perslice,
+                perlevel=perlevel, fname_vert_level=fname_vert_level, method=method,
+                label_struc=label_struc, id_label=id_label, indiv_labels_ids=indiv_labels_ids)
 
         save_as_csv(agg_metric, fname_output, fname_in=fname_data, append=append_csv)
         append_csv = True  # when looping across labels, need to append results in the same file

testing/cli/test_cli_sct_extract_metric.py

@@ -1,5 +1,6 @@
 # pytest unit tests for sct_extract_metric
 
+import csv
 import numpy as np
 
 import pytest
@@ -48,3 +49,50 @@ def test_sct_extract_metric_vertfile_doesnt_exists():
                                       '-vert', '1:12', '-perlevel', '1',
                                       '-method', 'wa', '-l', '51', '-o', fname_out])
         assert e.value.code == 2
+
+
+@pytest.mark.sct_testing
+def test_sct_extract_metric_normalize_pam50(tmp_path):
+    """Verify that -normalize-PAM50 outputs per-slice metrics in PAM50 space with VertLevel populated."""
+    fname_out = str(tmp_path / 'quantif_mtr_pam50.csv')
+    sct_extract_metric.main(argv=[
+        '-i', sct_test_path('mt', 'mtr.nii.gz'),
+        '-f', sct_test_path('mt', 'label/atlas'),
+        '-method', 'wa', '-l', '51',
+        '-vertfile', sct_test_path('mt', 'label', 'template', 'PAM50_levels.nii.gz'),
+        '-perslice', '1', '-normalize-PAM50', '1',
+        '-o', fname_out
+    ])
+    with open(fname_out, 'r') as f:
+        rows = list(csv.DictReader(f))
+    assert len(rows) > 0
+    assert all(r['VertLevel'] != '' for r in rows)
+    assert all(r['WA()'] not in ('', 'nan') for r in rows)
+
+
+def test_sct_extract_metric_normalize_pam50_missing_perslice(tmp_path):
+    """Verify that -normalize-PAM50 1 without -perslice 1 raises an error."""
+    with pytest.raises(SystemExit) as e:
+        sct_extract_metric.main(argv=[
+            '-i', sct_test_path('mt', 'mtr.nii.gz'),
+            '-f', sct_test_path('mt', 'label/atlas'),
+            '-method', 'wa', '-l', '51',
+            '-vertfile', sct_test_path('mt', 'label', 'template', 'PAM50_levels.nii.gz'),
+            '-normalize-PAM50', '1',
+            '-o', str(tmp_path / 'out.csv')
+        ])
+    assert e.value.code == 2
+
+
+def test_sct_extract_metric_normalize_pam50_missing_vertfile(tmp_path):
+    """Verify that -normalize-PAM50 1 with a missing vertfile raises an error."""
+    with pytest.raises(SystemExit) as e:
+        sct_extract_metric.main(argv=[
+            '-i', sct_test_path('mt', 'mtr.nii.gz'),
+            '-f', sct_test_path('mt', 'label/atlas'),
+            '-method', 'wa', '-l', '51',
+            '-vertfile', sct_test_path('mt', 'label', 'template', 'does_not_exist.nii.gz'),
+            '-perslice', '1', '-normalize-PAM50', '1',
+            '-o', str(tmp_path / 'out.csv')
+        ])
+    assert e.value.code == 2