Rendering problems with add_contours when binary images have rotations set in the affine matrix

[vnckppl]

The problem
When I plot a binary image that has some rotation to it in the affine header using Nilearn's add_contours function, the resulting image includes many additional lines that are not part of the surface of the ROI in the binary image. I don't encounter these issues when I use the same code with the same image that has not rotation set in the affine matrix.

Reproduction of the issue
I have attached example data and code to reproduce this problem:
contour_issue.zip

My environment
Python 3.8.3 (default, May 27 2020, 20:53:40)
matplotlib==3.2.2
nilearn==0.6.2

[bthirion]

Thx for sharing the problem !
Clearly the interpolation involved in rotation application has a disastrous effect on the edge display.
Can you identify more rpecisely where in the code the issue takes place ?

[vnckppl]

In the example that I attached, that would be here:
display.add_contours(rr, colors="#00ff00")

add_contours() calls _map_show. I think that's where the real issue may be, because _map_show applies an affine matrix transformation.

[NicolasGensollen]

I had a look at this today and I was able to reproduce the issue with the code and data you provided @vnckppl

import numpy as np
import matplotlib.pyplot as plt
from nilearn import plotting
from nilearn.image import load_img, get_data

ddir = '/home/nicolas/Downloads/contour_issue/contour_issue/data'

# File names
template_default = ddir + '/SUIT.nii.gz'
template_rotated = ddir + '/SUIT_rot.nii.gz'
binROI_default = ddir + '/suitROI.nii.gz'
binROI_rotated = ddir + '/suitROI_rot.nii.gz'

# Load files
td = load_img(template_default)
tr = load_img(template_rotated)
rd = load_img(binROI_default)
rr = load_img(binROI_rotated)

# Check the data
assert np.all(td.affine == rd.affine)
assert np.all(tr.affine == rr.affine)
assert np.all(get_data(td) == get_data(tr))
assert np.all(get_data(rd) == get_data(rr))
assert np.all(np.unique(get_data(rd)) == [0,1])
assert np.all(np.unique(get_data(rr)) == [0,1])

# Plot td with rd contours
display = plotting.plot_anat(td)
display.add_contours(rd, colors="#00ff00")

# Plot tr with rr contours with same cut for comparison
display = plotting.plot_anat(tr, cut_coords = (-1,-36,-24))
display.add_contours(rr, colors="#00ff00")

After investigation, I believe the problem comes indeed from the function _map_show when called from add_contours with rr as argument. More precisely, at this line:

nilearn/nilearn/plotting/displays.py

Line 753 in 0c6388a

img = reorder_img(img, resample=resampling_interpolation)

rr (variable img in the code) is resampled with a continuous interpolation, such that the mask computed here:

nilearn/nilearn/plotting/displays.py

Line 780 in 0c6388a

not_mask = np.logical_or(data > thr, data < -thr)

isn't correct anymore which results in the weird looking output.

As a fix, I propose to call reorder_img only when type isn't contour or contourf:

if type not in ('contour', 'contourf'):
        img = reorder_img(img, resample=resampling_interpolation)

If you add this fix and re-run the code above, you should get the following output:

[NicolasGensollen]

Hmmm... After looking more into this I realize that my previous message and the proposed fix are wrong.

I believe the problem comes from the fact that, since the affine is not considered diagonal in the code below, the binary image is resampled with continuous interpolation:

nilearn/nilearn/image/resampling.py

Lines 694 to 706 in bb3cd50

	if not np.all((np.abs(A) > 0.001).sum(axis=0) == 1):
	# The affine is not nearly diagonal
	if resample is None:
	raise ValueError('Cannot reorder the axes: '
	'the image affine contains rotations')
	else:
	# Identify the voxel size using a QR decomposition of the
	# affine
	Q, R = np.linalg.qr(affine[:3, :3])
	target_affine = np.diag(np.abs(np.diag(R))[
	np.abs(Q).argmax(axis=1)])
	return resample_img(img, target_affine=target_affine,
	interpolation=resample)

[bthirion]

So the solution would be to turn to nearest interpolation for binary images ?

Actually, what makes more sense would be to interpolate the image from which the edges are estimated, and not interpolate an image (possibly raise a warning if this is done)...

[NicolasGensollen]

The problem is happening here:

nilearn/nilearn/image/resampling.py

Lines 270 to 275 in bb3cd50

	ndimage.affine_transform(data, A,
	offset=b,
	output_shape=target_shape,
	output=out,
	cval=fill_value,
	order=interpolation_order)

Before this line data is binary while the output of affine_transform with interpolation_order=3 ('continuous`) isn't.

So the solution would be to turn to nearest interpolation for binary images ?

This is indeed a possible solution. For example, doing something like this:

 def _map_show(self, img, type='imshow',
              resampling_interpolation='continuous',
              threshold=None, **kwargs):

    # This function would move somewhere more appropriate...
    def _is_binary_image(img):
        img = load_img(img)
        data = _utils.niimg._safe_get_data(img, ensure_finite=True)
        unique_values = np.unique(data)
        if len(unique_values) != 2:
            return False
        return sorted(list(unique_values)) == [0,1]

    if _is_binary_image(img):
        img = reorder_img(img, resample='nearest')
    else:
        img = reorder_img(img, resample=resampling_interpolation)
    threshold = float(threshold) if threshold is not None else None
....

gives this output with the example code of my first message:

Actually, what makes more sense would be to interpolate the image from which the edges are estimated, and not interpolate an image (possibly raise a warning if this is done)...

I'm not sure I understand what you mean here. In the code above, the resampling operation is applied to both tr when calling plot_anat, and rr when calling add_countours, since they both exhibit the same non diagonal affine. My previous wrong fix in #2656 was simply avoiding the transformation for rr which was not visible here but is clearly visible for this example:

• example in the online doc
• same example with wrong fix

[bthirion]

My point is that to avoid surprises, you should avoid reinterpolating a binary image. Since the binary image is always obtained from another initial image, I'd resample this original image rather than the binary (edge) image.
Anyway, your fix is right, so you can forget this comment.

[bthirion]

So please go ahead with this solution. thx !