Gallery example code produces speckled/wrong overlays

[mmbannert]

Hi,

I would like to display conjunction results in a similar manner as plot_fmri_conjunction.py. However, the example code plot_fmri_conjunction.py already fails to produce the image shown in the gallery. It produces speckled overlays and the overlap region, which is supposed to be shown in purple, does not show up at all. You can see the result in the attached screenshot on the right.

Is anyone familiar with this problem?

Also, I noticed that it produces the warnings shown on the left. I'm unsure if it matters. I tried it with an older version of vtk as well but it did not seem to make any difference.

I'm trying this on Ubuntu 22.04.1 LTS.

My environment.yml for the conda environment looks as follows:

channels:
  - conda-forge
dependencies:
  - python=3.9
  - ipython
  - pip
  - pip:
    - pysurfer

It would be great if you could help.

Thanks,
Michael

PS: The same problem arises with plot_fmri_activation.py

[mmbannert]

I have done some further testing and I found out that the code also fails on macOS 12.6 and shows the same artifacts.

[mmbannert]

For what it's worth, this is how I circumnavigated PySurfer's problem. I simply plotted both overlays independently, saved them to a PIL image and then used PIL's blend function to overlay them. It is cheeky. But it gets the job done. Here is the code:

import os
import numpy as np
from surfer import io
import mne
from tempfile import NamedTemporaryFile
from PIL import Image
import matplotlib.pyplot as plt

op = os.path

# You'll need to change this based on where your PySurfer examples are stored
# and where your Freesurfer subejcts folder is located
base_dir = op.expanduser('~/experiments/PySurfer/examples')
os.environ['SUBJECTS_DIR'] = op.expanduser(
    '~/group/software/freesurfer/freesurfer60/subjects')

# Use brain class from MNE
Brain = mne.viz.get_brain_class()

"""
Read both of the activation maps in using
surfer's io functions.
"""
sig1 = io.read_scalar_data(op.join(base_dir, 'example_data', "lh.sig.nii.gz"))
sig2 = io.read_scalar_data(op.join(base_dir, 'example_data', \
    "lh.alt_sig.nii.gz"))

"""
Zero out the vertices that do not meet a threshold.
"""
thresh = 4
sig1[sig1 < thresh] = 0
sig2[sig2 < thresh] = 0

"""
A conjunction is defined as the minimum significance
value between the two maps at each vertex.
"""
conjunct = np.min(np.vstack((sig1, sig2)), axis=0)

# This plots some surface overlay on a hemisphere, saves it to disk in a
# temporary file, then reads it as a PIL image. The temporary file is deleted
# afterwards.
def plot_overlay_on_brain(overlay: np.ndarray, hemisphere: str, \
    subject_id='fsaverage', surface='inflated', **kwargs):

    with NamedTemporaryFile(suffix='.png') as ntf:

        brain = Brain(subject_id, hemisphere, surface, offscreen=True,
            cortex='low_contrast', background='white')
        brain.add_data(overlay, **kwargs)
        brain.save_image(ntf.name)
        brain.close()
        brain_img = Image.open(ntf.name)

    return brain_img

# Create image of first overlay on brain
img1 = plot_overlay_on_brain(sig1, 'lh', fmin=4, fmid=5, fmax=30,
    transparent=True, colormap="Reds", colorbar=False, alpha=1)

plt.figure()
plt.imshow(np.asarray(img1))

# Create image of second overlay on brain
img2 = plot_overlay_on_brain(sig2, 'lh', fmin=4, fmid=5, fmax=30,
    transparent=True, colormap='Blues', colorbar=False, alpha=1)

plt.figure()
plt.imshow(np.asarray(img2))


# Combine the two
conj_img = Image.blend(img1, img2, .5)

# Et voilà!
plt.figure()
plt.imshow(np.asarray(conj_img))

plt.show()

[larsoner]

I'd recommend using mne.viz.Brain for this stuff nowadays, overlays should be plotted correctly and the code is actively maintained

[mmbannert]

Dear Eric,

Thanks for your comment. I fully agree that overlays should be plotted correctly and it does appear that MNE sees more activity than PySurfer these days.

After reading your comment on a similar issue a while ago, I in fact posted a question on how to achieve the same result with MNE on their message board. Unfortunately, there were no answers.

Cheers,
Michael

[larsoner]

Can you open an enhancement request on MNE's issue tracker? To me this is a pretty general request and something we could probably add to brain.add_data without too much work (hopefully!):

https://github.com/mne-tools/mne-python/issues/new?assignees=&labels=ENH&template=feature_request.yml

[mmbannert]

Done. The feature request can be seen here. Thanks! Oops, there it already is :) GitHub has added it for me.