Resample lesion mask into in_file's voxel grid

niworkflows/interfaces/norm.py

@@ -518,14 +518,16 @@ def create_cfm(in_file, lesion_mask=None, global_mask=True, out_path=None):
 
     Notes
     -----
-    in_file and lesion_mask must be in the same
-    image space and have the same dimensions
+    The lesion mask is resampled into in_file's voxel grid before being
+    combined, so the two only need to overlap in world (scanner) space; they
+    do not need to share the same orientation or dimensions.
 
     """
     import os
 
     import nibabel as nb
     import numpy as np
+    from nibabel.processing import resample_from_to
     from nipype.utils.filemanip import fname_presuffix
 
     if out_path is None:
@@ -549,11 +551,13 @@ def create_cfm(in_file, lesion_mask=None, global_mask=True, out_path=None):
 
     # If a lesion mask was provided, combine it with the secondary mask.
     if lesion_mask is not None:
-        # Reorient the lesion mask and get the data.
-        lm_img = nb.as_closest_canonical(nb.load(lesion_mask))
+        # Resample the lesion into in_file's voxel grid so the subtraction is
+        # spatially correct regardless of the lesion's stored orientation or
+        # grid. Nearest-neighbor (order=0) keeps the mask binary.
+        lm_img = resample_from_to(nb.load(lesion_mask), (data.shape, in_img.affine), order=0)
 
         # Subtract lesion mask from secondary mask, set negatives to 0
-        data = np.fmax(data - lm_img.dataobj, 0)
+        data = np.fmax(data - np.asanyarray(lm_img.dataobj), 0)
         # Cost function mask will be created from subtraction
     # Otherwise, CFM will be created from global mask
 

niworkflows/interfaces/tests/test_norm.py

@@ -20,7 +20,12 @@
 #
 #     https://www.nipreps.org/community/licensing/
 #
-from ..norm import SpatialNormalization
+import nibabel as nb
+import numpy as np
+from nibabel.affines import apply_affine
+from nibabel.orientations import axcodes2ornt, io_orientation, ornt_transform
+
+from ..norm import SpatialNormalization, create_cfm
 
 
 def test_get_settings():
@@ -32,3 +37,37 @@ def test_get_settings():
     norm = SpatialNormalization(moving='T1w', flavor='testing')
     settings = norm._get_settings()
     assert len(settings) == 3
+
+
+def test_create_cfm_lesion_orientation(tmp_path):
+    """A lesion stored in RAS must be excluded at the correct world location
+    even when in_file is stored in LPS (regression test for the lesion mask
+    orientation bug)."""
+    shape = (4, 5, 6)
+    ras_affine = np.eye(4)
+
+    # in_file: all-ones brain mask, stored in LPS orientation.
+    ras_img = nb.Nifti1Image(np.ones(shape, dtype=np.uint8), ras_affine)
+    xfm = ornt_transform(io_orientation(ras_img.affine), axcodes2ornt(('L', 'P', 'S')))
+    in_img = ras_img.as_reoriented(xfm)
+    in_file = str(tmp_path / 'in_lps.nii.gz')
+    in_img.to_filename(in_file)
+
+    # lesion: single voxel in RAS at world coordinate (1, 2, 3).
+    lesion_data = np.zeros(shape, dtype=np.uint8)
+    lesion_data[1, 2, 3] = 1
+    lesion_file = str(tmp_path / 'lesion_ras.nii.gz')
+    nb.Nifti1Image(lesion_data, ras_affine).to_filename(lesion_file)
+
+    out = create_cfm(
+        in_file,
+        lesion_mask=lesion_file,
+        global_mask=True,
+        out_path=str(tmp_path / 'cfm.nii.gz'),
+    )
+
+    cfm = nb.load(out)
+    zeros = np.argwhere(np.asanyarray(cfm.dataobj) == 0)
+    # Exactly one voxel excluded, at world coordinate (1, 2, 3).
+    assert zeros.shape[0] == 1
+    assert np.allclose(apply_affine(cfm.affine, zeros[0]), [1, 2, 3])