def test_io_surface():
    """Test reading and writing of Freesurfer surface mesh files
    """
    fname = op.join(data_path, 'subjects', 'fsaverage', 'surf', 'lh.inflated')
    pts, tri = read_surface(fname)
    write_surface(op.join(tempdir, 'tmp'), pts, tri)
    c_pts, c_tri = read_surface(op.join(tempdir, 'tmp'))
    assert_array_equal(pts, c_pts)
    assert_array_equal(tri, c_tri)

def test_io_surface():
    """Test reading and writing of Freesurfer surface mesh files
    """
    tempdir = _TempDir()
    fname_quad = op.join(data_path, 'subjects', 'bert', 'surf',
                         'lh.inflated.nofix')
    fname_tri = op.join(data_path, 'subjects', 'fsaverage', 'surf',
                        'lh.inflated')
    for fname in (fname_quad, fname_tri):
        pts, tri = read_surface(fname)
        write_surface(op.join(tempdir, 'tmp'), pts, tri)
        c_pts, c_tri = read_surface(op.join(tempdir, 'tmp'))
        assert_array_equal(pts, c_pts)
        assert_array_equal(tri, c_tri)

def test_flash_bem():
    """Test mne flash_bem."""
    check_usage(mne_flash_bem, force_help=True)
    # Using the sample dataset
    subjects_dir = op.join(sample.data_path(download=False), 'subjects')
    # Copy necessary files to tempdir
    tempdir = _TempDir()
    mridata_path = op.join(subjects_dir, 'sample', 'mri')
    subject_path_new = op.join(tempdir, 'sample')
    mridata_path_new = op.join(subject_path_new, 'mri')
    os.makedirs(op.join(mridata_path_new, 'flash'))
    os.makedirs(op.join(subject_path_new, 'bem'))
    shutil.copyfile(op.join(mridata_path, 'T1.mgz'),
                    op.join(mridata_path_new, 'T1.mgz'))
    shutil.copyfile(op.join(mridata_path, 'brain.mgz'),
                    op.join(mridata_path_new, 'brain.mgz'))
    # Copy the available mri/flash/mef*.mgz files from the dataset
    flash_path = op.join(mridata_path_new, 'flash')
    for kind in (5, 30):
        in_fname = op.join(mridata_path, 'flash', 'mef%02d.mgz' % kind)
        shutil.copyfile(in_fname, op.join(flash_path, op.basename(in_fname)))
    # Test mne flash_bem with --noconvert option
    # (since there are no DICOM Flash images in dataset)
    out_fnames = list()
    for kind in ('outer_skin', 'outer_skull', 'inner_skull'):
        out_fnames.append(op.join(subject_path_new, 'bem', 'outer_skin.surf'))
    assert not any(op.isfile(out_fname) for out_fname in out_fnames)
    with ArgvSetter(('-d', tempdir, '-s', 'sample', '-n'),
                    disable_stdout=False, disable_stderr=False):
        mne_flash_bem.run()
    # do they exist and are expected size
    for out_fname in out_fnames:
        _, tris = read_surface(out_fname)
        assert len(tris) == 5120

def plot_coregistration(subject, subjects_dir, hcp_path, recordings_path,
                        info_from=(('data_type', 'rest'), ('run_index', 0)),
                        view_init=(('azim', 0), ('elev', 0))):
    """A diagnostic plot to show the HCP coregistration

    Parameters
    ----------
    subject : str
        The subject
    subjects_dir : str
        The path corresponding to MNE/freesurfer SUBJECTS_DIR (to be created)
    hcp_path : str
        The path where the HCP files can be found.
    recordings_path : str
        The path to converted data (including the head<->device transform).
    info_from : tuple of tuples | dict
        The reader info concerning the data from which sensor positions
        should be read.
        Must not be empty room as sensor positions are in head
        coordinates for 4D systems, hence not available in that case.
        Note that differences between the sensor positions across runs
        are smaller than 12 digits, hence negligible.
    view_init : tuple of tuples | dict
        The initival view, defaults to azimuth and elevation of 0,
        a simple lateral view

    Returns
    -------
    fig : matplotlib.figure.Figure
        The figure object.
    """
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D  #  noqa

    if isinstance(info_from, tuple):
        info_from = dict(info_from)
    if isinstance(view_init, tuple):
        view_init = dict(view_init)

    head_mri_t = read_trans(
        op.join(recordings_path, subject,
                '{}-head_mri-trans.fif'.format(subject)))

    info = read_info(subject=subject, hcp_path=hcp_path, **info_from)

    info = pick_info(info, _pick_data_channels(info, with_ref_meg=False))
    sens_pnts = np.array([c['loc'][:3] for c in info['chs']])
    sens_pnts = apply_trans(head_mri_t, sens_pnts)
    sens_pnts *= 1e3  # put in mm scale

    pnts, tris = read_surface(
        op.join(subjects_dir, subject, 'bem', 'inner_skull.surf'))

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(*sens_pnts.T, color='purple', marker='o')
    ax.scatter(*pnts.T, color='green', alpha=0.3)
    ax.view_init(**view_init)
    fig.tight_layout()
    return fig

def dissolve_label(labels, source, targets, subjects_dir=None,
                   hemi='both'):
    """
    Assign every point from source to the target that is closest to it.

    Parameters
    ----------
    labels : list
        List of labels (as returned by mne.read_annot).
    source : str
        Name of the source label (without hemi affix).
    targets : list of str
        List of target label names (without hemi affix).
    subjects_dir : str
        subjects_dir.
    hemi : 'both', 'lh', 'rh'
        Hemisphere(s) for which to dissolve the label.

    Notes
    -----
    Modifies ``labels`` in-place, returns None.
    """
    subjects_dir = get_subjects_dir(subjects_dir)
    subject = labels[0].subject
    if hemi == 'both':
        hemis = ('lh', 'rh')
    elif hemi == 'lh' or hemi == 'rh':
        hemis = (hemi,)
    else:
        raise ValueError("hemi=%r" % hemi)

    idx = {l.name: i for i, l in enumerate(labels)}

    rm = set()
    for hemi in hemis:
        fpath = os.path.join(subjects_dir, subject, 'surf', hemi + '.inflated')
        points, _ = mne.read_surface(fpath)

        src_name = '-'.join((source, hemi))
        src_idx = idx[src_name]
        rm.add(src_idx)
        src_label = labels[src_idx]
        tgt_names = ['-'.join((name, hemi)) for name in targets]
        tgt_idxs = [idx[name] for name in tgt_names]
        tgt_labels = [labels[i] for i in tgt_idxs]
        tgt_points = [points[label.vertices] for label in tgt_labels]

        vert_by_tgt = {i: [] for i in xrange(len(targets))}
        for src_vert in src_label.vertices:
            point = points[src_vert:src_vert + 1]
            dist = [cdist(point, pts).min() for pts in tgt_points]
            tgt = np.argmin(dist)
            vert_by_tgt[tgt].append(src_vert)

        for i, label in enumerate(tgt_labels):
            new_vertices = vert_by_tgt[i]
            label.vertices = np.union1d(label.vertices, new_vertices)

    for i in sorted(rm, reverse=True):
        del labels[i]

def test_io_surface():
    """Test reading and writing of Freesurfer surface mesh files
    """
    tempdir = _TempDir()
    fname_quad = op.join(data_path, 'subjects', 'bert', 'surf',
                         'lh.inflated.nofix')
    fname_tri = op.join(data_path, 'subjects', 'fsaverage', 'surf',
                        'lh.inflated')
    for fname in (fname_quad, fname_tri):
        pts, tri, vol_info = read_surface(fname, read_metadata=True)
        write_surface(op.join(tempdir, 'tmp'), pts, tri, volume_info=vol_info)
        c_pts, c_tri, c_vol_info = read_surface(op.join(tempdir, 'tmp'),
                                                read_metadata=True)
        assert_array_equal(pts, c_pts)
        assert_array_equal(tri, c_tri)
        assert_true(_is_equal_dict([vol_info, c_vol_info]))

def test_watershed_bem():
    """Test mne watershed bem."""
    check_usage(mne_watershed_bem)
    # Copy necessary files to tempdir
    tempdir = _TempDir()
    mridata_path = op.join(subjects_dir, 'sample', 'mri')
    subject_path_new = op.join(tempdir, 'sample')
    mridata_path_new = op.join(subject_path_new, 'mri')
    os.mkdir(op.join(tempdir, 'sample'))
    os.mkdir(mridata_path_new)
    if op.exists(op.join(mridata_path, 'T1')):
        shutil.copytree(op.join(mridata_path, 'T1'), op.join(mridata_path_new,
                                                             'T1'))
    if op.exists(op.join(mridata_path, 'T1.mgz')):
        shutil.copyfile(op.join(mridata_path, 'T1.mgz'),
                        op.join(mridata_path_new, 'T1.mgz'))
    out_fnames = list()
    for kind in ('outer_skin', 'outer_skull', 'inner_skull'):
        out_fnames.append(op.join(subject_path_new, 'bem', 'inner_skull.surf'))
    assert not any(op.isfile(out_fname) for out_fname in out_fnames)
    with ArgvSetter(('-d', tempdir, '-s', 'sample', '-o'),
                    disable_stdout=False, disable_stderr=False):
        mne_watershed_bem.run()
    for out_fname in out_fnames:
        _, tris = read_surface(out_fname)
        assert len(tris) == 20480

def test_io_surface():
    """Test reading and writing of Freesurfer surface mesh files."""
    tempdir = _TempDir()
    fname_quad = op.join(data_path, 'subjects', 'bert', 'surf',
                         'lh.inflated.nofix')
    fname_tri = op.join(data_path, 'subjects', 'fsaverage', 'surf',
                        'lh.inflated')
    for fname in (fname_quad, fname_tri):
        with pytest.warns(None):  # no volume info
            pts, tri, vol_info = read_surface(fname, read_metadata=True)
        write_surface(op.join(tempdir, 'tmp'), pts, tri, volume_info=vol_info)
        with pytest.warns(None):  # no volume info
            c_pts, c_tri, c_vol_info = read_surface(op.join(tempdir, 'tmp'),
                                                    read_metadata=True)
        assert_array_equal(pts, c_pts)
        assert_array_equal(tri, c_tri)
        assert_equal(object_diff(vol_info, c_vol_info), '')

def test_read_curv():
    """Test reading curvature data."""
    fname_curv = op.join(data_path, 'subjects', 'fsaverage', 'surf', 'lh.curv')
    fname_surf = op.join(data_path, 'subjects', 'fsaverage', 'surf',
                         'lh.inflated')
    bin_curv = read_curvature(fname_curv)
    rr = read_surface(fname_surf)[0]
    assert_true(len(bin_curv) == len(rr))
    assert_true(np.logical_or(bin_curv == 0, bin_curv == 1).all())

def test_io_surface():
    """Test reading and writing of Freesurfer surface mesh files."""
    tempdir = _TempDir()
    fname_quad = op.join(data_path, 'subjects', 'bert', 'surf',
                         'lh.inflated.nofix')
    fname_tri = op.join(data_path, 'subjects', 'fsaverage', 'surf',
                        'lh.inflated')
    for fname in (fname_quad, fname_tri):
        with warnings.catch_warnings(record=True) as w:
            pts, tri, vol_info = read_surface(fname, read_metadata=True)
        assert_true(all('No volume info' in str(ww.message) for ww in w))
        write_surface(op.join(tempdir, 'tmp'), pts, tri, volume_info=vol_info)
        with warnings.catch_warnings(record=True) as w:  # No vol info
            c_pts, c_tri, c_vol_info = read_surface(op.join(tempdir, 'tmp'),
                                                    read_metadata=True)
        assert_array_equal(pts, c_pts)
        assert_array_equal(tri, c_tri)
        assert_equal(object_diff(vol_info, c_vol_info), '')

def plot_stc_time_point(stc, subject, limits=[5, 10, 15], time_index=0,
                        surf='inflated', measure='dSPM', subjects_dir=None):
    """Plot a time instant from a SourceEstimate using matplotlib

    The same could be done with mayavi using proper 3D.

    Parameters
    ----------
    stc : instance of SourceEstimate
        The SourceEstimate to plot.
    subject : string
        The subject name (only needed if surf is a string).
    time_index : int
        Time index to plot.
    surf : str, or instance of surfaces
        Surface to use (e.g., 'inflated' or 'white'), or pre-loaded surfaces.
    measure : str
        The label for the colorbar. None turns the colorbar off.
    subjects_dir : str, or None
        Path to the SUBJECTS_DIR. If None, the path is obtained by using
        the environment variable SUBJECTS_DIR.
    """
    subjects_dir = get_subjects_dir(subjects_dir)
    pl.figure(facecolor='k', figsize=(8, 5))
    hemis = ['lh', 'rh']
    if isinstance(surf, str):
        surf = [read_surface(op.join(subjects_dir, subject, 'surf',
                                     '%s.%s' % (h, surf))) for h in hemis]
    my_cmap = mne_analyze_colormap(limits)
    for hi, h in enumerate(hemis):
        coords = surf[hi][0][stc.vertno[hi]]
        if hi == 0:
            vals = stc_all_cluster_vis.lh_data[:, time_index]
        else:
            vals = stc_all_cluster_vis.rh_data[:, time_index]
        ax = pl.subplot(1, 2, 1 - hi, axis_bgcolor='none')
        pl.tick_params(labelbottom='off', labelleft='off')
        flipper = -1 if hi == 1 else 1
        sc = ax.scatter(flipper * coords[:, 1], coords[:, 2], c=vals,
                        vmin=-limits[2], vmax=limits[2], cmap=my_cmap,
                        edgecolors='none', s=5)
        ax.set_aspect('equal')
        pl.axis('off')
    try:
        pl.tight_layout(0)
    except:
        pass
    if measure is not None:
        cax = pl.axes([0.85, 0.15, 0.025, 0.15], axisbg='k')
        cb = pl.colorbar(sc, cax, ticks=[-limits[2], 0, limits[2]])
        cb.set_label(measure, color='w')
    pl.setp(pl.getp(cb.ax, 'yticklabels'), color='w')
    pl.draw()
    pl.show()

def test_sparse_morph():
    """Test sparse morphing."""
    rng = np.random.RandomState(0)
    vertices_fs = [np.sort(rng.permutation(np.arange(10242))[:4]),
                   np.sort(rng.permutation(np.arange(10242))[:6])]
    data = rng.randn(10, 1)
    stc_fs = SourceEstimate(data, vertices_fs, 1, 1, 'fsaverage')
    spheres_fs = [mne.read_surface(op.join(
        subjects_dir, 'fsaverage', 'surf', '%s.sphere.reg' % hemi))[0]
        for hemi in ('lh', 'rh')]
    spheres_sample = [mne.read_surface(op.join(
        subjects_dir, 'sample', 'surf', '%s.sphere.reg' % hemi))[0]
        for hemi in ('lh', 'rh')]
    morph_fs_sample = compute_source_morph(
        stc_fs, 'fsaverage', 'sample', sparse=True, spacing=None,
        subjects_dir=subjects_dir)
    stc_sample = morph_fs_sample.apply(stc_fs)
    offset = 0
    orders = list()
    for v1, s1, v2, s2 in zip(stc_fs.vertices, spheres_fs,
                              stc_sample.vertices, spheres_sample):
        dists = cdist(s1[v1], s2[v2])
        order = np.argmin(dists, axis=-1)
        assert_array_less(dists[np.arange(len(order)), order], 1.5)  # mm
        orders.append(order + offset)
        offset += len(order)
    assert_allclose(stc_fs.data, stc_sample.data[np.concatenate(orders)])
    # Return
    morph_sample_fs = compute_source_morph(
        stc_sample, 'sample', 'fsaverage', sparse=True, spacing=None,
        subjects_dir=subjects_dir)
    stc_fs_return = morph_sample_fs.apply(stc_sample)
    offset = 0
    orders = list()
    for v1, s, v2 in zip(stc_fs.vertices, spheres_fs, stc_fs_return.vertices):
        dists = cdist(s[v1], s[v2])
        order = np.argmin(dists, axis=-1)
        assert_array_less(dists[np.arange(len(order)), order], 1.5)  # mm
        orders.append(order + offset)
        offset += len(order)
    assert_allclose(stc_fs.data, stc_fs_return.data[np.concatenate(orders)])

def test_annot_io():
    """Test I/O from and to *.annot files."""
    # copy necessary files from fsaverage to tempdir
    tempdir = _TempDir()
    subject = 'fsaverage'
    label_src = os.path.join(subjects_dir, 'fsaverage', 'label')
    surf_src = os.path.join(subjects_dir, 'fsaverage', 'surf')
    label_dir = os.path.join(tempdir, subject, 'label')
    surf_dir = os.path.join(tempdir, subject, 'surf')
    os.makedirs(label_dir)
    os.mkdir(surf_dir)
    shutil.copy(os.path.join(label_src, 'lh.PALS_B12_Lobes.annot'), label_dir)
    shutil.copy(os.path.join(label_src, 'rh.PALS_B12_Lobes.annot'), label_dir)
    shutil.copy(os.path.join(surf_src, 'lh.white'), surf_dir)
    shutil.copy(os.path.join(surf_src, 'rh.white'), surf_dir)

    # read original labels
    with pytest.raises(IOError, match='\nPALS_B12_Lobes$'):
        read_labels_from_annot(subject, 'PALS_B12_Lobesey',
                               subjects_dir=tempdir)
    labels = read_labels_from_annot(subject, 'PALS_B12_Lobes',
                                    subjects_dir=tempdir)

    # test saving parcellation only covering one hemisphere
    parc = [l for l in labels if l.name == 'LOBE.TEMPORAL-lh']
    write_labels_to_annot(parc, subject, 'myparc', subjects_dir=tempdir)
    parc1 = read_labels_from_annot(subject, 'myparc', subjects_dir=tempdir)
    parc1 = [l for l in parc1 if not l.name.startswith('unknown')]
    assert_equal(len(parc1), len(parc))
    for l1, l in zip(parc1, parc):
        assert_labels_equal(l1, l)

    # test saving only one hemisphere
    parc = [l for l in labels if l.name.startswith('LOBE')]
    write_labels_to_annot(parc, subject, 'myparc2', hemi='lh',
                          subjects_dir=tempdir)
    annot_fname = os.path.join(tempdir, subject, 'label', '%sh.myparc2.annot')
    assert os.path.isfile(annot_fname % 'l')
    assert not os.path.isfile(annot_fname % 'r')
    parc1 = read_labels_from_annot(subject, 'myparc2',
                                   annot_fname=annot_fname % 'l',
                                   subjects_dir=tempdir)
    parc_lh = [l for l in parc if l.name.endswith('lh')]
    for l1, l in zip(parc1, parc_lh):
        assert_labels_equal(l1, l)

    # test that the annotation is complete (test Label() support)
    rr = read_surface(op.join(surf_dir, 'lh.white'))[0]
    label = sum(labels, Label(hemi='lh', subject='fsaverage')).lh
    assert_array_equal(label.vertices, np.arange(len(rr)))

def run():
    from mne.commands.utils import get_optparser

    parser = get_optparser(__file__)

    parser.add_option("-s", "--surf", dest="surf",
                      help="Surface in Freesurfer format", metavar="FILE")
    parser.add_option("-f", "--fif", dest="fif",
                      help="FIF file produced", metavar="FILE")
    parser.add_option("-i", "--id", dest="id", default=4,
                      help=("Surface Id (e.g. 4 sur head surface)"))

    options, args = parser.parse_args()

    if options.surf is None:
        parser.print_help()
        sys.exit(1)

    print("Converting %s to BEM FIF file." % options.surf)
    points, tris = mne.read_surface(options.surf)
    points *= 1e-3
    surf = dict(coord_frame=5, id=int(options.id), nn=None, np=len(points),
                ntri=len(tris), rr=points, sigma=1, tris=tris)
    mne.write_bem_surface(options.fif, surf)

def test_random_parcellation():
    """Test generation of random cortical parcellation."""
    hemi = 'both'
    n_parcel = 50
    surface = 'sphere.reg'
    subject = 'sample_ds'
    rng = np.random.RandomState(0)

    # Parcellation
    labels = random_parcellation(subject, n_parcel, hemi, subjects_dir,
                                 surface=surface, random_state=rng)

    # test number of labels
    assert_equal(len(labels), n_parcel)
    if hemi == 'both':
        hemi = ['lh', 'rh']
    hemis = np.atleast_1d(hemi)
    for hemi in set(hemis):
        vertices_total = []
        for label in labels:
            if label.hemi == hemi:
                # test that labels are not empty
                assert (len(label.vertices) > 0)

                # vertices of hemi covered by labels
                vertices_total = np.append(vertices_total, label.vertices)

        # test that labels don't intersect
        assert_equal(len(np.unique(vertices_total)), len(vertices_total))

        surf_fname = op.join(subjects_dir, subject, 'surf', hemi + '.' +
                             surface)
        vert, _ = read_surface(surf_fname)

        # Test that labels cover whole surface
        assert_array_equal(np.sort(vertices_total), np.arange(len(vert)))

centers_file='nki_region_centers.txt'
surface='pial'

new_names_file='nki_region_names_fix.txt'
annot_name='nki'


names=[]

with open(names_file,'r') as fd:
	for ln in fd:
		names.append(ln.strip())

centers=np.loadtxt(centers_file)

lhsurf=mne.read_surface('lh.%s'%surface)
rhsurf=mne.read_surface('rh.%s'%surface)

labs_lh=[]
labs_rh=[]
region_nums={}

with open(new_names_file,'w') as fd:
	#find correct number for this label
	for i,(c,n) in enumerate(zip(centers,names)):
		if n in region_nums:
			nr=region_nums[n]+1
		else:
			nr=1
		region_nums[n]=nr

def _stc_to_label(stc, src, smooth, subjects_dir=None):
    """Compute a label from the non-zero sources in an stc object.

    Parameters
    ----------
    stc : SourceEstimate
        The source estimates.
    src : SourceSpaces | str | None
        The source space over which the source estimates are defined.
        If it's a string it should the subject name (e.g. fsaverage).
        Can be None if stc.subject is not None.
    smooth : int
        Number of smoothing iterations.
    subjects_dir : str | None
        Path to SUBJECTS_DIR if it is not set in the environment.

    Returns
    -------
    labels : list of Labels | list of list of Labels
        The generated labels. If connected is False, it returns
        a list of Labels (one per hemisphere). If no Label is available
        in a hemisphere, None is returned. If connected is True,
        it returns for each hemisphere a list of connected labels
        ordered in decreasing order depending of the maximum value in the stc.
        If no Label is available in an hemisphere, an empty list is returned.
    """
    src = stc.subject if src is None else src

    if isinstance(src, string_types):
        subject = src
    else:
        subject = stc.subject

    if isinstance(src, string_types):
        subjects_dir = get_subjects_dir(subjects_dir)
        surf_path_from = op.join(subjects_dir, src, 'surf')
        rr_lh, tris_lh = read_surface(op.join(surf_path_from,
                                      'lh.white'))
        rr_rh, tris_rh = read_surface(op.join(surf_path_from,
                                      'rh.white'))
        rr = [rr_lh, rr_rh]
        tris = [tris_lh, tris_rh]
    else:
        if not isinstance(src, SourceSpaces):
            raise TypeError('src must be a string or a set of source spaces')
        if len(src) != 2:
            raise ValueError('source space should contain the 2 hemispheres')
        rr = [1e3 * src[0]['rr'], 1e3 * src[1]['rr']]
        tris = [src[0]['tris'], src[1]['tris']]

    labels = []
    cnt = 0
    for hemi_idx, (hemi, this_vertno, this_tris, this_rr) in enumerate(
            zip(['lh', 'rh'], stc.vertices, tris, rr)):
        this_data = stc.data[cnt:cnt + len(this_vertno)]
        e = mesh_edges(this_tris)
        e.data[e.data == 2] = 1
        n_vertices = e.shape[0]
        e = e + sparse.eye(n_vertices, n_vertices)

        clusters = [this_vertno[np.any(this_data, axis=1)]]

        cnt += len(this_vertno)

        clusters = [c for c in clusters if len(c) > 0]

        if len(clusters) == 0:
            this_labels = None
        else:
            this_labels = []
            colors = _n_colors(len(clusters))
            for c, color in zip(clusters, colors):
                idx_use = c
                for k in range(smooth):
                    e_use = e[:, idx_use]
                    data1 = e_use * np.ones(len(idx_use))
                    idx_use = np.where(data1)[0]

                label = Label(idx_use, this_rr[idx_use], None, hemi,
                              'Label from stc', subject=subject,
                              color=color)

                this_labels.append(label)

            this_labels = this_labels[0]

        labels.append(this_labels)

    return labels

def labels_from_mni_coords(seeds, extent=30., subject='fsaverage',
                           surface='white', mask=None, subjects_dir=None,
                           parc=None):
    """Create a parcellation from seed coordinates in MNI space

    Parameters
    ----------
    seeds : dict
        Seed coordinates. Keys are label names, including -hemi tags. values
        are seeds (array_like of shape (3,) or (3, n_seeds)).
    extent : scalar
        Extent of the label in millimeters (maximum distance from the seed).
    subject : str
        MRI-subject to use (default 'fsaverage').
    surface : str
        Surface to use (default 'white').
    mask : None | str
        A parcellation used to mask the parcellation under construction.
    subjects_dir : str
        SUBJECTS_DIR.
    parc : None | str
        Name of the parcellation under construction (only used for error
        messages).
    """
    name_re = re.compile("\w+-(lh|rh)$")
    if not all(name.endswith(('lh', 'rh')) for name in seeds):
        err = ("Names need to end in 'lh' or 'rh' so that the proper "
               "hemisphere can be selected")
        raise ValueError(err)

    # load surfaces
    subjects_dir = get_subjects_dir(subjects_dir)
    fpath = os.path.join(subjects_dir, subject, 'surf', '.'.join(('%s', surface)))
    surfs = {hemi: mne.read_surface(fpath % hemi) for hemi in ('lh', 'rh')}

    # prepare seed properties for mne.grow_labels
    vertices = []
    names = []
    hemis = []
    for name, coords_ in seeds.iteritems():
        m = name_re.match(name)
        if not m:
            raise ValueError("Invalid seed name in %r parc: %r. Names must "
                             "conform to the 'xxx-lh' or 'xxx-rh' scheme."
                             % (parc, name))
        coords = np.atleast_2d(coords_)
        if coords.ndim != 2 or coords.shape[1] != 3:
            raise ValueError("Invalid coordinate specification for seed %r in "
                             "parc %r: %r. Seeds need to be specified as "
                             "arrays with shape (3,) or (n_seeds, 3)."
                             % (name, parc, coords_))
        hemi = m.group(1)
        seed_verts = []
        for coord in coords:
            dist = np.sqrt(np.sum((surfs[hemi][0] - coord) ** 2, axis=1))
            seed_verts.append(np.argmin(dist))
        vertices.append(seed_verts)
        names.append(name)
        hemis.append(hemi == 'rh')

    # grow labels
    labels = mne.grow_labels(subject, vertices, extent, hemis, subjects_dir,
                             1, False, names, surface)

    # apply mask
    if mask is not None:
        mlabels = mne.read_labels_from_annot(subject, mask,
                                             subjects_dir=subjects_dir)
        unknown = {l.hemi: l for l in mlabels if l.name.startswith('unknown-')}

        for label in labels:
            rm = unknown[label.hemi]
            if np.any(np.in1d(label.vertices, rm.vertices)):
                label.vertices = np.setdiff1d(label.vertices, rm.vertices, True)

    return labels

Find distances between vertices and plot vertices in a small region.

mainly using functions from within mne.label.grow_labels
'''

import mne
from mne.datasets import sample

data_path = sample.data_path()
subjects_dir = data_path + '/subjects'

tris, vert, dist = {}, {}, {}
hemi = 0  # lh

# read the surface
vert[hemi], tris[hemi] = mne.read_surface(subjects_dir + '/fsaverage/surf/lh.inflated')

# obtain distance matrix
dist[hemi] = mne.label.mesh_dist(tris[hemi], vert[hemi])

# choose seed vertex as 20 and plot vertices within 5mm radius around it
# obtain neighbouring vertices within 5mm distance
my_verts, my_dist = mne.label._verts_within_dist(dist[hemi], [20], 5)

# number of vertices in a given radius
print len(my_verts)

from surfer import Brain
brain = Brain('fsaverage', hemi='lh', surf='inflated',
              subjects_dir='/Users/psripad/sciebo/resting_state_analysis/')

import sys

import mne

if __name__ == '__main__':

    from mne.commands.utils import get_optparser

    parser = get_optparser(__file__)

    parser.add_option("-s", "--surf", dest="surf",
                    help="Surface in Freesurfer format", metavar="FILE")
    parser.add_option("-f", "--fif", dest="fif",
                    help="FIF file produced", metavar="FILE")
    parser.add_option("-i", "--id", dest="id", default=4,
                    help=("Surface Id (e.g. 4 sur head surface)"))

    options, args = parser.parse_args()

    if options.surf is None:
        parser.print_help()
        sys.exit(1)

    print("Converting %s to BEM FIF file." % options.surf)

    points, tris = mne.read_surface(options.surf)
    points *= 1e-3
    surf = dict(coord_frame=5, id=int(options.id), nn=None, np=len(points),
                ntri=len(tris), rr=points, sigma=1, tris=tris)
    mne.write_bem_surface(options.fif, surf)