Conjunction Map: overlap not working for brain with both hemispheres

[mattvan83]

Dear PySurfer,

I tried to save different brain views of data overlap inspired from the example below:

https://pysurfer.github.io/auto_examples/plot_fmri_conjunction.html#sphx-glr-auto-examples-plot-fmri-conjunction-py

Buth when comparing the different views' results I realized that conjunction map was not working correctly in some of the cases on the views with both hemispheres. Indeed, you can see that on the right hemisphere views conjunction map is well noticeable in purple colour, whereas on the views with both hemipsheres the conjunction map is not visible in the right hemisphere.

Please find below the code I used and attached the raw data and images resulting:

"""
Load fdg overlay arrays.
"""
## Load TFCE p-corrected maps
overlay_file_fdg_lh = os.path.join("fdg", "lh.palm.fdr0.CorrPos_tfce_tstat_fwep.mgz")
overlay_file_fdg_rh = os.path.join("fdg", "rh.palm.fdr0.CorrPos_tfce_tstat_fwep.mgz")

overlay_array_fdg_lh = nib.load(overlay_file_fdg_lh)
overlay_array_fdg_lh = overlay_array_fdg_lh.get_fdata()
overlay_array_fdg_lh = np.squeeze(overlay_array_fdg_lh)

overlay_array_fdg_rh = nib.load(overlay_file_fdg_rh)
overlay_array_fdg_rh = overlay_array_fdg_rh.get_fdata()
overlay_array_fdg_rh = np.squeeze(overlay_array_fdg_rh)

"""
Load eAV45 overlay arrays.
"""
## Load TFCE p-corrected maps
overlay_file_eAV45_lh = os.path.join("eAV45", "lh.palm.fdr0.CorrPos_tfce_tstat_fwep.mgz")
overlay_file_eAV45_rh = os.path.join("eAV45", "rh.palm.fdr0.CorrPos_tfce_tstat_fwep.mgz")

overlay_array_eAV45_lh = nib.load(overlay_file_eAV45_lh)
overlay_array_eAV45_lh = overlay_array_eAV45_lh.get_fdata()
overlay_array_eAV45_lh = np.squeeze(overlay_array_eAV45_lh)

overlay_array_eAV45_rh = nib.load(overlay_file_eAV45_rh)
overlay_array_eAV45_rh = overlay_array_eAV45_rh.get_fdata()
overlay_array_eAV45_rh = np.squeeze(overlay_array_eAV45_rh)

"""
Define conjunction overlay arrays.
"""
conjunct_lh = np.min(np.vstack((overlay_array_fdg_lh, overlay_array_eAV45_lh)), axis=0)
conjunct_rh = np.min(np.vstack((overlay_array_fdg_rh, overlay_array_eAV45_rh)), axis=0)

for hemi in ["lh", "rh", "both"]:
         """
         Initialize the visualization.
         """
         brain = Brain("fsaverage", hemi, "white, background="black")

         if hemi == "rh":
               brain.add_data(overlay_array_fdg_rh, min=1.3, mid=1.3+0.01, max=20, thresh=1.3, colormap="Blues", colorbar=False, hemi="rh", alpha=1)
               brain.add_data(overlay_array_eAV45_rh, min=1.3, mid=1.3+0.01, max=20, thresh=1.3, colormap="Reds", colorbar=False, hemi="rh", alpha=1)
               brain.add_data(conjunct_rh, min=1.3, mid=1.3+0.01, max=20, thresh=1.3, colormap="Purples", colorbar=False, hemi="rh", alpha=1)

               brain.save_montage(os.path.join(outdir,"rh." + filename + ".tiff"), order=['medial','lateral'], orientation='h', border_size=15, colorbar=None, row=-1, col=-1)
           elif hemi == "both":
                """
                Load rh overlays.
                """
                brain.add_data(overlay_array_fdg_rh, min=1.3, mid=1.3+0.01, max=20, thresh=1.3, colormap="Blues", colorbar=False, hemi="rh", alpha=1)
                brain.add_data(overlay_array_eAV45_rh, min=1.3, mid=1.3+0.01, max=20, thresh=1.3, colormap="Reds", colorbar=False, hemi="rh", alpha=1)
                brain.add_data(conjunct_rh, min=1.3, mid=1.3+0.01, max=20, thresh=1.3, colormap="Purples", colorbar=False, hemi="rh", alpha=1)
                """
                Load lh overlays.
                """
                brain.add_data(overlay_array_fdg_lh, min=1.3, mid=1.3+0.01, max=20, thresh=1.3, colormap="Blues", colorbar=False, hemi="lh", alpha=1)
                brain.add_data(overlay_array_eAV45_lh, min=1.3, mid=1.3+0.01, max=20, thresh=1.3, colormap="Reds", colorbar=False, hemi="lh", alpha=1)
                brain.add_data(conjunct_lh, min=1.3, mid=1.3+0.01, max=20, thresh=1.3, colormap="Purples", colorbar=False, hemi="lh", alpha=1)
               """
               Save lh/rh montage.
               """
               brain.save_montage(os.path.join(outdir,"both." + filename + ".tiff"), order=['ventral', 'dorsal'], orientation='v', border_size=15, colorbar=None, row=-1, col=-1)

Could anyone help me with this problem?

Best,
Matthieu
ConjunctionOverlays.zip

[mattvan83]

Any suggestions please ?

[mwaskom]

I think I can confirm this one, but I don't understand it yet.

Here is a self-contained example:

b = Brain("fsaverage", "both", "white", cortex=None, background="black")
b.show_view("dor")
n = 163842
for hemi in ["lh", "rh"]:
    for colormap in ["Reds", "Blues"]:
        data = np.random.randn(n)
        b.add_data(data, min=0, max=1, thresh=0, colormap=colormap, colorbar=False, hemi=hemi)

With PySurfer 0.9.0, this renders:

Watching the plot appear, it's clear that the right hemisphere has two separate overlays but they both appear with the Blues colormap.

[mwaskom]

Hm, this issue is very flaky and didn't persist across a notebook kernel restart.

[mwaskom]

Doh, forgot that I had activated the master branch to start a bisect before restarting the kernel.

So whatever this is, it seems fixed on master.

[mwaskom]

@mattvan83 it also looks like your original example runs fine on master.

[mattvan83]

@mwaskom ok I will test it. So what do you call master, the last release of PySurfer?
How is the proper way to install it on OSX?

[mwaskom]

“Master” means the current status of the master branch on this github repository. It implies a newer version of the library than what pip or conda give you by default, although you can use pip to install it pretty easily (I’m not sure if that’s true for conda as well).

By the way @larsoner we haven’t released in a while and may want to do that after we sort out the other issues that @mattvan83 has run into.

[mattvan83]

I tried with the master (0.10.dev0) and it still didn't work.

Please find the exact code and raw data I used attached:

#### Overlap both hemispheres FDG/eAV45 Figure BEM DEBUG ####
WD = "/Users/matthieu/Google Drive/Code/Python/PySurfer_DEBUG_overlap/"
min_thr_color = 1.3

"""
Load fdg overlay arrays.
"""
## Load TFCE p-corrected maps
overlay_file_fdg_lh = os.path.join(WD, "fdg", "lh.palm.fdr0.CorrPos_tfce_tstat_fwep.mgz")
overlay_file_fdg_rh = os.path.join(WD, "fdg", "rh.palm.fdr0.CorrPos_tfce_tstat_fwep.mgz")

overlay_array_fdg_lh = nib.load(overlay_file_fdg_lh)
overlay_array_fdg_lh = overlay_array_fdg_lh.get_fdata()
overlay_array_fdg_lh = np.squeeze(overlay_array_fdg_lh)

overlay_array_fdg_rh = nib.load(overlay_file_fdg_rh)
overlay_array_fdg_rh = overlay_array_fdg_rh.get_fdata()
overlay_array_fdg_rh = np.squeeze(overlay_array_fdg_rh)

"""
Check if one of both hemispheres contains values above the minimum threshold (p=0.05 --> -log10(p)=1.3).
"""
lh_size_fdg_over_thr = overlay_array_fdg_lh[np.where(overlay_array_fdg_lh >= min_thr_color)].size
rh_size_fdg_over_thr = overlay_array_fdg_rh[np.where(overlay_array_fdg_rh >= min_thr_color)].size

"""
Load eAV45 overlay arrays.
"""
## Load TFCE p-corrected maps
overlay_file_eAV45_lh = os.path.join(WD, "eAV45", "lh.palm.fdr0.CorrPos_tfce_tstat_fwep.mgz")
overlay_file_eAV45_rh = os.path.join(WD, "eAV45", "rh.palm.fdr0.CorrPos_tfce_tstat_fwep.mgz")

overlay_array_eAV45_lh = nib.load(overlay_file_eAV45_lh)
overlay_array_eAV45_lh = overlay_array_eAV45_lh.get_fdata()
overlay_array_eAV45_lh = np.squeeze(overlay_array_eAV45_lh)

overlay_array_eAV45_rh = nib.load(overlay_file_eAV45_rh)
overlay_array_eAV45_rh = overlay_array_eAV45_rh.get_fdata()
overlay_array_eAV45_rh = np.squeeze(overlay_array_eAV45_rh)

"""
Check if one of both hemispheres contains values above the minimum threshold (p=0.05 --> -log10(p)=1.3).
"""
lh_size_eAV45_over_thr = overlay_array_eAV45_lh[np.where(overlay_array_eAV45_lh >= min_thr_color)].size
rh_size_eAV45_over_thr = overlay_array_eAV45_rh[np.where(overlay_array_eAV45_rh >= min_thr_color)].size

"""
Define filename of the overlap between fdg and eAV45 patterns.
"""
filename = "test"

"""
Define conjunction overlay arrays.
"""
conjunct_lh = np.min(np.vstack((overlay_array_fdg_lh, overlay_array_eAV45_lh)), axis=0)
conjunct_rh = np.min(np.vstack((overlay_array_fdg_rh, overlay_array_eAV45_rh)), axis=0)

brain = Brain("fsaverage", "both", "white", background="black", subjects_dir="/Applications/freesurfer/subjects/")

"""
Load rh overlays.
"""
if rh_size_fdg_over_thr > 0:
    brain.add_data(overlay_array_fdg_rh, min=1.3, mid=1.3+0.01, max=20, thresh=1.3, colormap="Blues", colorbar=False, hemi="rh", alpha=1)
if rh_size_eAV45_over_thr > 0:
    brain.add_data(overlay_array_eAV45_rh, min=1.3, mid=1.3+0.01, max=20, thresh=1.3, colormap="Reds", colorbar=False, hemi="rh", alpha=1)
if rh_size_fdg_over_thr > 0 and rh_size_eAV45_over_thr > 0:
    brain.add_data(conjunct_rh, min=1.3, mid=1.3+0.01, max=20, thresh=1.3, colormap="Purples", colorbar=False, hemi="rh", alpha=1)
"""
Load lh overlays.
"""
if lh_size_fdg_over_thr > 0:
    brain.add_data(overlay_array_fdg_lh, min=1.3, mid=1.3+0.01, max=20, thresh=1.3, colormap="Blues", colorbar=False, hemi="lh", alpha=1)
if lh_size_eAV45_over_thr > 0:
    brain.add_data(overlay_array_eAV45_lh, min=1.3, mid=1.3+0.01, max=20, thresh=1.3, colormap="Reds", colorbar=False, hemi="lh", alpha=1)
if lh_size_fdg_over_thr > 0 and lh_size_eAV45_over_thr > 0:
    brain.add_data(conjunct_lh, min=1.3, mid=1.3+0.01, max=20, thresh=1.3, colormap="Purples", colorbar=False, hemi="lh", alpha=1)
"""
Save lh/rh montage.
"""
brain.save_montage(os.path.join(WD,"both." + filename + ".tiff"), order=['ventral','dorsal'], orientation='v', border_size=15, colorbar=None, row=-1, col=-1)
[ConjunctionMaps.zip](https://github.com/nipy/PySurfer/files/4083289/ConjunctionMaps.zip)

Do you have any idea?

[mattvan83]

ConjunctionMaps.zip

Sorry with the attached file.

[mwaskom]

Can you please specify what "it didn't work" means?

[mattvan83]

Sorry, it gave this image.

[mwaskom]

Can you please tell me what happens when you run the self-contained example I posted in this thread? It's really a lot easier to work with simple examples that don't require downloading external data, if possible.

[mattvan83]

Please find below the image it gave with your simple self-contained example.

It seems to work, doesn't it ? The problem may be in the save_montage method ?

[mwaskom]

Unfortunately, I can confirm that there is still a major bug in add_data. Here is another self-contained example that shows what's going on:

import os
import numpy as np
import nibabel as nib
import matplotlib.pyplot as plt
from surfer import Brain

b = Brain("fsaverage", "both", "white", background="black")
b.show_view("dorsal")

screenshots = []

colormaps = iter(["Reds", "Blues", "Purples", "Greens"])

for hemi in ["lh", "rh"]:
    coords, _ = nib.freesurfer.read_geometry(
        os.path.join(os.environ["SUBJECTS_DIR"], f"fsaverage/surf/{hemi}.white")
    )
    for sign in [-1, +1]:
        colormap = next(colormaps)
        b.add_data(np.sign(coords[:, 1]) == sign,
                   colormap=colormap, min=0, max=2, thresh=.5,
                   colorbar=False, verbose=False, hemi=hemi)
        screenshots.append(b.screenshot())
b.close()

f, axes = plt.subplots(1, len(screenshots), figsize=(len(screenshots) * 3, 3))
for ax, img in zip(axes, screenshots):
    ax.imshow(img)
    ax.set_axis_off()
f.tight_layout()
f.savefig("add_data_bug.png")

You can see that adding data on the right hemisphere is changing the properties of the last data array added to the left hemisphere.

I suspect that it has to do with scale_data_colormap looping over both hemispheres, even when add_data was called with a specific hemisphere. The fix may be simple. But I didn't write the code and am a little confused by it. @larsoner do you know what's going on here?

[larsoner]

Unfortunately I don't remember -- I always have to re-learn what those loops are doing when I look at them, too :(

[mattvan83]

Any trick to subvert the problem? This is pretty common when mapping on the cortical surface overlap between two patterns.

[mwaskom]

Plenty of tricks. The issue arises when you're not plotting the same set of arrays on the left and right hemispheres. If you can work around having to conditionally add overlays depending on whether they pass the threshold (a constraint imposed by mayavi that i find annoying), you won't see the issue.

One idea would be to use the fact that adding transparent=True will make the colormap go from 0% to 100% opacity between mid and max.

So instead of min=thresh, mid=thresh + .01, thresh=thresh, you should be able to do min=thresh - .01, mid=thresh.

In a quick test, this produces a figure that looks right.

[larsoner]

@mwaskom do you think we should just release now? We can take up add_data tweaks separately if need be, though it sounds like there might already be some suitable workaround.

[mwaskom]

No, I think this issue is a blocker because it produces incorrect figures.

[larsoner]

@mwaskom pushed a release to PyPi and updated docs, would you be up for sending a release email (if you think it's necessary/useful) whenever there's time?