nipy / PySurfer

Cortical neuroimaging visualization in Python
https://pysurfer.github.io/
BSD 3-Clause "New" or "Revised" License
240 stars 97 forks source link

Medial view non-"region" is being mapped onto #200

Closed YSanchezAraujo closed 7 years ago

YSanchezAraujo commented 7 years ago

In the picture below, I think the large area in the middle of the medial views is usually greyed out in images and in this use case (I think it should just be white)I'm not sure why values are being mapped onto it.

xb_brain

That result is gotten using:

import os
import sys
import numpy as np
from custom import utils
import nibabel as nb
from surfer import Brain

result_name = sys.argv[1]
header_name = sys.argv[2]
save_name = sys.argv[3]

brain = Brain("fsaverage", "split", "inflated",
               views=['lat', 'med'], background="white")

hp = "/storage/gablab001/data/genus/fs_cog/exp1/column_headers"
rp = "/storage/gablab001/data/genus/fs_cog/exp1/lg/results"
with open(os.path.join(hp, header_name), "r") as tmp:
    headers = np.array(tmp.read().strip("\n").split("\n"))

res = utils.read_pickle(os.path.join(rp,"{}".format(result_name)))
coefs = [val for key, val in res.items() if "coef" in key]
nz = lambda x: np.nonzero(x)[0]
features_selected_lh = []
features_selected_rh = []
features_all = []

from collections import Counter
import pandas as pd

for coef in coefs:
    idx = nz(coef)
    features_all.extend(headers[idx].tolist())
    for col in headers[idx].tolist():
        if 'lh' in col:
            features_selected_lh.append(col)
        elif 'rh' in col:
            features_selected_rh.append(col)

left_count = dict(Counter(features_selected_lh))
right_count = dict(Counter(features_selected_rh))

left_df = pd.DataFrame({
    'col': [key.replace('lh_','').replace('.','-').replace('_thickness_D','')
            for key, val in left_count.items()],
    'val': [val for key, val in left_count.items()]
})

right_df = pd.DataFrame({
    'col': [key.replace('rh_','').replace('.','-').replace('_thickness_D','')
            for key, val in right_count.items()],
    'val': [val for key, val in right_count.items()]
})

label_dir = '/cm/shared/openmind/freesurfer/5.3.0/subjects/fsaverage/label'
left_label_file = 'lh.aparc.a2009s.annot'
right_label_file = 'rh.aparc.a2009s.annot'
lh_aparc_file = os.path.join(label_dir, left_label_file)
rh_aparc_file = os.path.join(label_dir, right_label_file)
lh_labels, lh_ctab, lh_names = nb.freesurfer.read_annot(lh_aparc_file)
rh_labels, rh_ctab, rh_names = nb.freesurfer.read_annot(rh_aparc_file)
left_df = left_df.set_index('col').loc[lh_names].reset_index().fillna(0)
right_df = right_df.set_index('col').loc[rh_names].reset_index().fillna(0)
vtx_lh = left_df.val.values[lh_labels]
vtx_rh = right_df.val.values[rh_labels]

brain.add_data(vtx_lh,
               0,
               400,
               colormap="Reds",
               alpha=.8,
               hemi='lh')

brain.add_annotation(lh_aparc_file, hemi='lh')

brain.add_data(vtx_rh,
               0,
               400,
               colormap="Reds",
               alpha=.8,
               hemi='rh')

brain.add_annotation(rh_aparc_file, hemi='rh', remove_existing=False)

If I look at left_df for values equal to or less than 70 I get:

                       col   val
0                  Unknown   0.0
7   G_and_S_cingul-Mid-Ant  51.0
13     G_front_inf-Orbital  57.0
24               G_orbital  69.0
26   G_pariet_inf-Supramar  46.0
30             G_precuneus  60.0
34      G_temp_sup-Lateral  36.0
42             Medial_wall   0.0
43          Pole_occipital  57.0
46               S_central  51.0
47     S_cingul-Marginalis  52.0
49   S_circular_insula_inf  63.0
51     S_collat_transv_ant  61.0
53             S_front_inf  54.0

I thought maybe something was happening with "Unknown" and it doesn't seem to be in the map from the names to the labels:

# these returns False
"Unknown" in np.array(lh_names)[lh_labels]
"Unknown" in left_df.col.values[lh_labels]

# this gives a 1. 
(np.array(lh_names)[lh_labels] == left_df.col.values[lh_labels]).mean()

Any ideas on what's happening?

mwaskom commented 7 years ago

It's a little hard to say without being able to run your script, but I think the problem is here:

vtx_lh = left_df.val.values[lh_labels]

The medial wall vertices are -1 in the lh_labels array, so the vtx_lh array is going to put the final value (i.e. the value in the -1 positional index) in the medial wall. Does that make sense?

YSanchezAraujo commented 7 years ago

@mwaskom thanks for the info. So I don't think I understand - I guess I have the wrong idea here but if I do:

lh_labels, lh_ctab, lh_names = nb.freesurfer.read_annot(lh_aparc_file)
rh_labels, rh_ctab, rh_names = nb.freesurfer.read_annot(rh_aparc_file)
zip(lh_labels, lh_names)

then -1 in lh_lables shows up like this:

(-1, 'G_and_S_cingul-Mid-Post')
(-1, 'G_temp_sup-Plan_tempo')
(-1, 'G_temporal_middle')

and Medial_wall in lh_names is (46, 'Medial_wall')

or should I not expect the order of the values in those two lists to correspond to one another?

mwaskom commented 7 years ago

zip(lh_labels, lh_names)

These don't positionally correspond. Look at the lengths. The labels array has as many entries as there are vertices, the names list has one fewer entries than there are unique labels in the labels array, because the undefined (-1) vertices don't have a name.

YSanchezAraujo commented 7 years ago

I think I see what you're saying. Then my question becomes how do I order lh_names so the index representations in lh_labels map to the correct index in lh_names ? I ordered my data based on the order nibabel returns the names list since I was going off https://pysurfer.github.io/examples/plot_parc_values.html

mwaskom commented 7 years ago

I think all you need to do is account separately for the vertices that are undefined in the parcellation. I.e.

vtx_lh[lh_labels == -1] = 0

Then the medial wall will be white in your colormap. Or you can set it to some negative number and do thresh=0 and it will be gray.

In the parcellation example, the Yeo atlas uses a label for the medial vertices (rather than leaving them undefined) so this problem doesn't arise.

YSanchezAraujo commented 7 years ago

Oh ok that works. Thank you!