hcpmmp_lut.py

# In fsaverage/mri ran:
#
# $ mri_aparc2aseg --s fsaverage --volmask --annot HCPMMP1
# $ mri_aparc2aseg --s fsaverage --volmask --annot HCPMMP1_combined
#
# Then this script can be used to create the lookup table for atlas_ids.

import os.path as op
import numpy as np
import mne

mri_dir = op.join(mne.get_config('SUBJECTS_DIR'), 'fsaverage', 'mri')

fs_lut = op.join(op.dirname(mne.__file__), 'data', 'FreeSurferColorLUT.txt')
for parc in ('HCPMMP1', 'HCPMMP1_combined'):
    labels = mne.read_labels_from_annot('fsaverage', parc, 'both', sort=False)
    header = """
# Below is the color table for the cortical labels of the seg volume
# created by mri_aparc2aseg in which the aseg cortex label is replaced
# by the labels in the %s aparc.

#No.    Label Name:                         R   G   B   A
""" % (parc,)
    ids = dict()
    with open(op.join(mri_dir, '%sColorLUT.txt' % (parc,)), 'w') as fid:
        with open(fs_lut, 'r') as fid_:
            for line in fid_:
                fid.write(line)
                if line.strip() and not line.strip().startswith('#'):
                    id_, name = line.strip().split()[:2]
                    id_ = int(id_)
                    assert name not in ids
                    ids[name] = id_
                if 'SUSPICIOUS' in line:
                    break
        fid.write(header)
        n_left = sum(label.hemi == 'lh' for label in labels)
        for li, label in enumerate(labels):
            offset = 1000 if label.hemi == 'lh' else 2000 - n_left
            id_ = li + offset
            assert label.name not in ids
            ids[label.name] = id_
            name = label.name.replace(' ', '_').ljust(35)
            color = np.round(np.array(label.color) * 255).astype(int)
            color[-1] = 255 - color[-1]
            color = ' '.join('%3d' % x for x in color)
            line = '%d    %s %s\n' % (id_, name, color)
            fid.write(line)
    # Eventually we could remove the hippocampus labels, but people
    # can easily do this just by popping from the LUT
    # (also it's nice to keep the labels complete)
    """
    if parc == 'HCPMMP1':
        import nibabel as nib
        hip_names = (('L_H_ROI-lh', 'Left-Hippocampus'),
                     ('R_H_ROI-rh', 'Right-Hippocampus'))
        assert all(h[0] in ids for h in hip_names)
        assert all(h[1] in ids for h in hip_names)
        atlas = nib.load(op.join(mri_dir, '%s+aseg.mgz' % (parc,)))
        data = atlas.get_fdata()
        raise RuntimeError
    """