Python mne.read_surface函数代码示例

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)))