Reassigning RGB values to parcellation

[amsinha9]

I am trying to convert the ROI label colors from hcp.mmp to those that correspond to a figure I have. I generated a colormap.csv file that contains the new RGB values and a reference map called roimap that lists the Glasser parcellation ROIs with the same ROI label and network ID as before. See code below:

import nibabel as nib
import nilearn.plotting as plotting
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import pandas
import csv as csv
import os
import copy
%matplotlib inline
import hcp_utils as hcp

colormap = "colormap.csv"
roimap = "ROI_Mapping_Python_csv.csv"
cmap = pandas.read_csv(colormap, header=None)
rmap = pandas.read_csv(roimap)

myparc = copy.deepcopy(hcp.mmp)

for l in range(1,361):
    # get network id
    networkID = (rmap['Network #'].loc[rmap['Region Label']==hcp.mmp.labels[l]]).values
    # get rgb
    newR, newG, newB = cmap[0][networkID].values[0], cmap[1][networkID].values[0], cmap[2][networkID].values[0]
    rgbval = [float(newR)/255.0, float(newG)/255.0, float(newB)/255.0]
    # saturate
    hsv = matplotlib.colors.rgb_to_hsv(rgbval) 
    hsv[1] = hsv[1] + 0.1
    newrgbval = matplotlib.colors.hsv_to_rgb(hsv)

    # set new rgb map
    myparc.rgba[l][0:3] = newrgbval

Although I have 379 ROIs as in hcp.mmp, I can't access the final 361-379. Moreover, the new RGB values in myparc.rgba are not being converted correctly. Could you provide some help with this?

colormap.xlsx ROI_Mapping_Python.xlsx

[amsinha9]

Actually, I was able to resolve the above issue. The only thing now is when I try to make a subset parcellation, I get the error:

ValueError Traceback (most recent call last)

in () 38 mesh_sub = hcp.load_surfaces() 39 # hcp.view_parcellation(hcp.mesh.inflated, submmp) ---> 40 hcp.view_parcellation(mesh_sub.inflated, submmp) 41 # hcp.parcellation_labels(submmp) 42 # hcp.parcellation_labels(myparc) ~\Anaconda3\lib\site-packages\hcp_utils\hcp_utils.py in view_parcellation(meshLR, parcellation) 274 275 cmap = matplotlib.colors.ListedColormap(rgba) --> 276 return plotting.view_surf(meshLR, normalized_cortex_map, symmetric_cmap=False, cmap=cmap) 277 278 def parcellation_labels(parcellation): ~\AppData\Roaming\Python\Python36\site-packages\nilearn\plotting\html_surface.py in view_surf(surf_mesh, surf_map, bg_map, threshold, cmap, black_bg, vmax, vmin, symmetric_cmap, colorbar, colorbar_height, colorbar_fontsize, title, title_fontsize) 322 surf_map=surf_map, surf_mesh=surf_mesh, threshold=threshold, 323 cmap=cmap, black_bg=black_bg, bg_map=bg_map, --> 324 symmetric_cmap=symmetric_cmap, vmax=vmax, vmin=vmin) 325 info['colorbar'] = colorbar 326 info['cbar_height'] = colorbar_height ~\AppData\Roaming\Python\Python36\site-packages\nilearn\plotting\html_surface.py in one_mesh_info(surf_map, surf_mesh, threshold, cmap, black_bg, bg_map, symmetric_cmap, vmax, vmin) 52 colors = colorscale( 53 cmap, surf_map, threshold, symmetric_cmap=symmetric_cmap, ---> 54 vmax=vmax, vmin=vmin) 55 info['inflated_left'] = mesh_to_plotly(surf_mesh) 56 info['vertexcolor_left'] = _get_vertexcolor( ~\AppData\Roaming\Python\Python36\site-packages\nilearn\plotting\js_plotting_utils.py in colorscale(cmap, values, threshold, symmetric_cmap, vmax, vmin) 106 'Custom cmap', cmaplist, cmap.N) 107 x = np.linspace(0, 1, 100) --> 108 rgb = our_cmap(x, bytes=True)[:, :3] 109 rgb = np.array(rgb, dtype=int) 110 colors = [] ~\AppData\Roaming\Python\Python36\site-packages\matplotlib\colors.py in __call__(self, X, alpha, bytes) 480 # See class docstring for arg/kwarg documentation. 481 if not self._isinit: --> 482 self._init() 483 mask_bad = None 484 if not cbook.iterable(X): ~\AppData\Roaming\Python\Python36\site-packages\matplotlib\colors.py in _init(self) 686 self._lut = np.ones((self.N + 3, 4), float) 687 self._lut[:-3, 0] = makeMappingArray( --> 688 self.N, self._segmentdata['red'], self._gamma) 689 self._lut[:-3, 1] = makeMappingArray( 690 self.N, self._segmentdata['green'], self._gamma) ~\AppData\Roaming\Python\Python36\site-packages\matplotlib\colors.py in makeMappingArray(N, data, gamma) 400 if x[0] != 0. or x[-1] != 1.0: 401 raise ValueError( --> 402 "data mapping points must start with x=0 and end with x=1") 403 if (np.diff(x) < 0).any(): 404 raise ValueError("data mapping points must have x in increasing order") ValueError: data mapping points must start with x=0 and end with x=1 Below is my code to generate a subset of the entire parcellation to generate a brain figure: ``` # subset parcellation from sklearn.utils import Bunch def make_subset_parcellation(parcellation, subset): """ Takes the given parcellation and produces a new one where only parcels in the subset are retained. Remaining part of the brain is set to 0 (unassigned). """ ids = np.arange(len(subset)+1) nontrivial_ids = np.arange(1, len(subset)+1) # labels labels = dict() labels[0] = 'None' for i, c in enumerate(subset): labels[i+1] = parcellation.labels[c] # colors rgba = dict() rgba[0] = np.array([0.9,0.9,0.9,1]) for i, c in enumerate(subset): rgba[i+1] = parcellation.rgba[c] # mapping map_all = np.zeros_like(parcellation.map_all) for i, c in enumerate(subset): map_all[parcellation.map_all==c] = i+1 new_parcellation = Bunch() new_parcellation.map_all = map_all new_parcellation.labels = labels new_parcellation.ids = ids new_parcellation.nontrivial_ids = nontrivial_ids new_parcellation.rgba = rgba return new_parcellation # submmp = make_subset_parcellation(myparc, [100, 105, 107, 110, 200, 250, 260, 300]) submmp = make_subset_parcellation(myparc, [370]) mesh_sub = hcp.load_surfaces() # hcp.view_parcellation(hcp.mesh.inflated, submmp) hcp.view_parcellation(mesh_sub.inflated, submmp) ````

[rmldj]

Hi,

Could you attach the code for constructing myparc? (as you corrected the one from the first post) Could you attach also the csv files (instead of xlsx). Romuald

[amsinha9]

Can we switch this over to email instead? I tried attaching the csv files, but it says "We don't support that file type". If so, let me know what email address to contact you at.

[rmldj]

Send me the csv's by e-mail to romuald.janik@gmail.com. But best put the corrected code for constructing myparc here..

[amsinha9]

Here is the full code for generating myparc to take the existing RGB values, replace them with new RGB values and generate a figure that shows a subset of the parcellation rather than the entire parcellation.

import nibabel as nib
import nilearn.plotting as plotting
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import pandas
import csv as csv
import os
import copy
import hcp_utils as hcp
%matplotlib inline

colormap = "colormap.csv"
roimap = "ROI_Mapping_Python_csv.csv"

cmap = pandas.read_csv(colormap, header=None)
rmap = pandas.read_csv(roimap)

myparc = copy.deepcopy(hcp.mmp)

for l in range(1,380):
    # get network id
    # networkID = (rmap['Network #'].loc[rmap['Region Label']==hcp.mmp.labels[l]]).values
    networkID = (rmap['Network #'].loc[rmap['Region Label']==hcp.mmp.labels[l]])-1
    # get rgb
    newR, newG, newB = cmap[0][networkID], cmap[1][networkID], cmap[2][networkID]
    rgbval = [float(newR)/255.0, float(newG)/255.0, float(newB)/255.0]
    # saturate
    hsv = matplotlib.colors.rgb_to_hsv(rgbval) 
    hsv[1] = hsv[1] + 0.1
    newrgbval = matplotlib.colors.hsv_to_rgb(hsv)

    # set new rgb map
    myparc.rgba[l][0:3] = newrgbval

def view_parcellation(meshLR, parcellation, subset=None, color=(0.9,0.9,0.9,1)):
    """
    View the given parcellation on an a whole brain surface mesh.
    """
    # for some parcellations the numerical ids need not be consecutive
    cortex_map = hcp.cortex_data(parcellation.map_all)
    ids = np.unique(cortex_map)
    normalized_cortex_map = np.zeros_like(cortex_map)
    rgba = np.zeros((len(ids), 4))
    for i in range(len(ids)):
        ind = cortex_map==ids[i]
        normalized_cortex_map[ind] = i
        if subset is None or i in subset:
            rgba[i,:] = parcellation.rgba[ids[i]]
        else:
            rgba[i,:] = color

    cmap = matplotlib.colors.ListedColormap(rgba)
    return plotting.view_surf(meshLR, normalized_cortex_map, symmetric_cmap=False, cmap=cmap)
from sklearn.utils import Bunch

def make_subset_parcellation(parcellation, subset):
    """
    Takes the given parcellation and produces a new one where only parcels in the subset are retained.
    Remaining part of the brain is set to 0 (unassigned).
    """
    ids = np.arange(len(subset)+1)
    nontrivial_ids = np.arange(1, len(subset)+1)
    # labels
    labels = dict()
    labels[0] = 'None'
    for i, c in enumerate(subset):
        labels[i+1] = parcellation.labels[c]
    # colors
    rgba = dict()
    rgba[0] = np.array([0.9,0.9,0.9,1])
    for i, c in enumerate(subset):
        rgba[i+1] = parcellation.rgba[c]
    # mapping
    map_all = np.zeros_like(parcellation.map_all)
    for i, c in enumerate(subset):
        map_all[parcellation.map_all==c] = i+1

    new_parcellation = Bunch()
    new_parcellation.map_all = map_all
    new_parcellation.labels = labels
    new_parcellation.ids = ids
    new_parcellation.nontrivial_ids = nontrivial_ids
    new_parcellation.rgba = rgba

    return new_parcellation

submmp = make_subset_parcellation(myparc, [361])
mesh_sub = hcp.load_surfaces()
hcp.view_parcellation(mesh_sub.inflated, submmp)

For some reason, it won't generate the subset parcellation for the remaining 361-379 parcels. Could you provide some help with this? The line submmp = make_subset_parcellation(myparc, [ ]) works for all 1-360 parcels, but no indices larger than 360. I get the error, "ValueError: data mapping points must start with x=0 and end with x=1"

[rmldj]

The MMP cortical parcellation has 360 cortical regions. The regions numbered 361-379 are the standard subcortical structures available in cifti. If you display myparc.labels you can see at the end:

 359: 'R_a32pr',
 360: 'R_p24',
 361: 'accumbens_left',
 362: 'accumbens_right',
 363: 'amygdala_left',
 364: 'amygdala_right',
 365: 'caudate_left',
 366: 'caudate_right',
 367: 'cerebellum_left',
 368: 'cerebellum_right',
 369: 'diencephalon_left',
 370: 'diencephalon_right',
 371: 'hippocampus_left',
 372: 'hippocampus_right',
 373: 'pallidum_left',
 374: 'pallidum_right',
 375: 'putamen_left',
 376: 'putamen_right',
 377: 'thalamus_left',
 378: 'thalamus_right',
 379: 'brainStem'

The line

submmp = make_subset_parcellation(myparc, [361])

works ok. The error came from

hcp.view_parcellation(mesh_sub.inflated, submmp)

because there was no nontrivial region in the cortex (361 is 'accumbens_left').

[amsinha9]

Right, and if I do

hcp.parcellation_labels(hcp.mmp)

I get the full 379 parcels (360 + 19 subcortical). I just changed the RGB values for all of the parcels and want to plot a subset of them in a brain figure. Is there a way to resolve the error in using hcp.view_parcellation(mesh_sub.inflated, submmp)? Shuld I be using something else instead of mesh_sub.inflated?

I was able to generate figures with subsets of the 360 parcels with the above code, but just need a way to generate the same figure with the subset of 361-379.

[rmldj]

hcp.view_parcellation shows only the cortical part of the parcellation. I.e. this is a surface plot on the surface of the cortex. You cannot see subcortical structures there by definition.

[amsinha9]

Okay so hcp.view_parcellation(mesh_sub.inflated, hcp.mmp) only shows the 360 cortical structures. Do you have any suggestions on how I could visualize the subcortical structures in any way using hcp-utils or a function I could write to do so?

[rmldj]

Unfortunately not - I do not know of a good way to plot them in 3D. But these 19 subcortical structures are fairly standard - do you really need to visualize them?

[amsinha9]

Yes, I was trying to see if I could work with one of the other preloaded parcellations in hcp-utils, but encountered the same issue.

In terms of visualizing the subcortical network, yes our results include identifying subcortical regions, so we would like to include a similar brain figure showing the subcortical regions in a brain figure for a publication. We are citing your toolbox in our work as a tool for generating the figures because it has been immensely helpful, but we just need some sort of similar figure that highlights the 19 subcortical regions.

[rmldj]

Indeed it would be nice to have something like this but I am not aware of a library which would be doing that...

[amsinha9]

I will try to find a FreeSurfer subcortical image to use instead. Thank you for your time and help with this; hcp-utils has been a very valuable resource for us, and the support has been excellent.