[BUG] Potential wrong edge_mask calculation and out-of-range indexing issues with QuickShear defacing

[uturkbey]

Description of the potential edge_mask calculation issue

@sezginer and I might have found an issue in brainles_preprocessing/defacing/quickshear/nipy_quickshear.py where the edge_mask computation can become incorrect when the brain area in the brain masks touch the array boundary.

In particular, the current implementation of edge_mask uses np.roll to compare shifted versions of the mask:

preprocessing/brainles_preprocessing/defacing/quickshear/nipy_quickshear.py

Lines 86 to 93 in 6878ce1

	edgemask = (
	4 * brain
	- np.roll(brain, 1, 0)
	- np.roll(brain, -1, 0)
	- np.roll(brain, 1, 1)
	- np.roll(brain, -1, 1)
	!= 0
	)

When the brain mask touches an array edge, np.roll introduces wrap-around: voxels rolled off from one side, appear on the opposite side. This creates an artificial edge on the opposite boundary. As a result, edgemask contains spurious edges and can propagate incorrect geometry into subsequent steps.

An example case is visualized in the top row of the first figure below. In RPS orientation, posterior part of the brain mask touches the posterior array boundary. During edge_mask calculation, 'fake' edge voxels appear on the anterior side even though the real boundary contact is posterior. And these 'fake' edges propagate through convex_hull, and slope calculation steps. Second figure zooms in to 'fake' edge area.

Description of the potential out-of-range indexing during defaced_mask calculation issue

Additionally (related, but not necessarily dependent), a separate error can occur during the calculation of the defacing mask, if the computed defacing plane intersects the array edge posterior to the brain (e.g., as in the top row right plot of the first figure above).

Based on the ys calculation with mask_RPS.shape[2] and the for loop iterating over the positive ys indices, x values might go beyond defaced_mask_RPS.shape[1] range and cause index out-of-range errors for such special cases.

preprocessing/brainles_preprocessing/defacing/quickshear/nipy_quickshear.py

Lines 162 to 166 in 6878ce1

	ys = np.arange(0, mask_RPS.shape[2]) * slope + yint
	defaced_mask_RPS = np.ones(mask_RPS.shape, dtype="bool")
	for x, y in zip(np.nonzero(ys > 0)[0], ys.astype(int)):
	defaced_mask_RPS[:, x, :y] = 0

In our case, the incorrect edgemask triggered this later crash, but we believe the out-of-range indexing can happen independently depending on geometry of brain masks and plane placement (although we acknowledge this should be a very rare case with regular brain segmentations).

To Reproduce

Steps to reproduce the behavior:

1. Download this brain segmentation mask: native__t1-der-cor-nc_brain_mask.nii.gz
2. Install 'BrainLes-Preprocessing' as suggested
3. Run defacing with the example code below:

   from brainles_preprocessing.defacing import QuickshearDefacer
   
   defacer = QuickshearDefacer(
       buffer=10.0,
       force_atlas_registration=False,  # We don't need atlas registration
   )
   
   defacer.deface(
       input_image_path=Path(your_path_to_native__t1-der-cor-nc_brain_mask),
       mask_image_path=Path(your_path_to_output_defacing_mask)
   )

4. You should see an error message similar to this:

   Traceback (most recent call last):
     File "/home/user/deface_test.py", line 112, in <module>
       defacer.deface(
     File "/home/user/preprocessing/brainles_preprocessing/defacing/quickshear/quickshear.py", line 68, in deface
       mask = run_quickshear(bet_img=bet_img, buffer=self.buffer)
     File "/home/user/preprocessing/brainles_preprocessing/defacing/quickshear/nipy_quickshear.py", line 166, in run_quickshear
       defaced_mask_RPS[:, x, :y] = 0
   IndexError: index 180 is out of bounds for axis 1 with size 180

5. As intermediate values such as edgemask are implicit, you would need to extract those arrays manually.

Expected behavior

Second row of the first figure above shows our expected behavior. To get these results:

1. We replace the np.roll implementation in edge_mask calculation with convolution as below:

   kernel = np.array([
       [0, -1,  0],
       [-1, 4, -1],
       [0, -1,  0]
   ])
   edgemask = ndimage.convolve(brain, kernel, mode='constant', cval=0)

operating system and version

Debian GNU/Linux 11 (bullseye)

Python environment and version?

miniforge3 environment with Python 3.10.19

version of brainles_preprocessing

brainles_preprocessing==0.6.8
numpy==2.2.6