def test_io_bem_surfaces():
    """Testing reading of bem surfaces
    """
    surf = read_bem_surfaces(fname, add_geom=True)
    surf = read_bem_surfaces(fname, add_geom=False)
    print "Number of surfaces : %d" % len(surf)

    write_bem_surface('bem_surf.fif', surf[0])
    surf_read = read_bem_surfaces('bem_surf.fif', add_geom=False)

    for key in surf[0].keys():
        assert_array_almost_equal(surf[0][key], surf_read[0][key])

def test_io_bem_surfaces():
    """Test reading of bem surfaces
    """
    tempdir = _TempDir()
    surf = read_bem_surfaces(fname, add_geom=True)
    surf = read_bem_surfaces(fname, add_geom=False)
    print("Number of surfaces : %d" % len(surf))

    write_bem_surface(op.join(tempdir, 'bem_surf.fif'), surf[0])
    surf_read = read_bem_surfaces(op.join(tempdir, 'bem_surf.fif'),
                                  add_geom=False)

    for key in surf[0].keys():
        assert_array_almost_equal(surf[0][key], surf_read[0][key])

def test_bem_model():
    """Test BEM model creation from Python with I/O"""
    tempdir = _TempDir()
    fname_temp = op.join(tempdir, "temp-bem.fif")
    for kwargs, fname in zip((dict(), dict(conductivity=[0.3])), [fname_bem_3, fname_bem_1]):
        model = make_bem_model("sample", ico=2, subjects_dir=subjects_dir, **kwargs)
        model_c = read_bem_surfaces(fname)
        _compare_bem_surfaces(model, model_c)
        write_bem_surfaces(fname_temp, model)
        model_read = read_bem_surfaces(fname_temp)
        _compare_bem_surfaces(model, model_c)
        _compare_bem_surfaces(model_read, model_c)
    assert_raises(
        ValueError, make_bem_model, "sample", conductivity=[0.3, 0.006], subjects_dir=subjects_dir  # bad conductivity
    )

def test_volume_source_space():
    """Test setting up volume source spaces."""
    tempdir = _TempDir()
    src = read_source_spaces(fname_vol)
    temp_name = op.join(tempdir, 'temp-src.fif')
    surf = read_bem_surfaces(fname_bem, s_id=FIFF.FIFFV_BEM_SURF_ID_BRAIN)
    surf['rr'] *= 1e3  # convert to mm
    # The one in the testing dataset (uses bem as bounds)
    for bem, surf in zip((fname_bem, None), (None, surf)):
        src_new = setup_volume_source_space(
            'sample', pos=7.0, bem=bem, surface=surf, mri='T1.mgz',
            subjects_dir=subjects_dir)
        write_source_spaces(temp_name, src_new, overwrite=True)
        src[0]['subject_his_id'] = 'sample'  # XXX: to make comparison pass
        _compare_source_spaces(src, src_new, mode='approx')
        del src_new
        src_new = read_source_spaces(temp_name)
        _compare_source_spaces(src, src_new, mode='approx')
    pytest.raises(IOError, setup_volume_source_space, 'sample',
                  pos=7.0, bem=None, surface='foo',  # bad surf
                  mri=fname_mri, subjects_dir=subjects_dir)
    assert repr(src) == repr(src_new)
    assert src.kind == 'volume'
    # Spheres
    sphere = make_sphere_model(r0=(0., 0., 0.), head_radius=0.1,
                               relative_radii=(0.9, 1.0), sigmas=(0.33, 1.0))
    src = setup_volume_source_space(pos=10)
    src_new = setup_volume_source_space(pos=10, sphere=sphere)
    _compare_source_spaces(src, src_new, mode='exact')
    pytest.raises(ValueError, setup_volume_source_space, sphere='foo')
    # Need a radius
    sphere = make_sphere_model(head_radius=None)
    pytest.raises(ValueError, setup_volume_source_space, sphere=sphere)

def test_make_bem_model(tmpdir):
    """Test BEM model creation from Python with I/O."""
    tempdir = str(tmpdir)
    fname_temp = op.join(tempdir, 'temp-bem.fif')
    for kwargs, fname in zip((dict(), dict(conductivity=[0.3])),
                             [fname_bem_3, fname_bem_1]):
        model = make_bem_model('sample', ico=2, subjects_dir=subjects_dir,
                               **kwargs)
        model_c = read_bem_surfaces(fname)
        _compare_bem_surfaces(model, model_c)
        write_bem_surfaces(fname_temp, model)
        model_read = read_bem_surfaces(fname_temp)
        _compare_bem_surfaces(model, model_c)
        _compare_bem_surfaces(model_read, model_c)
    # bad conductivity
    with pytest.raises(ValueError, match='conductivity must be'):
        make_bem_model('sample', 4, [0.3, 0.006], subjects_dir=subjects_dir)

def test_bem_solution():
    """Test making a BEM solution from Python with I/O."""
    # test degenerate conditions
    surf = read_bem_surfaces(fname_bem_1)[0]
    pytest.raises(RuntimeError, _ico_downsample, surf, 10)  # bad dec grade
    s_bad = dict(tris=surf['tris'][1:], ntri=surf['ntri'] - 1, rr=surf['rr'])
    pytest.raises(RuntimeError, _ico_downsample, s_bad, 1)  # not isomorphic
    s_bad = dict(tris=surf['tris'].copy(), ntri=surf['ntri'],
                 rr=surf['rr'])  # bad triangulation
    s_bad['tris'][0] = [0, 0, 0]
    pytest.raises(RuntimeError, _ico_downsample, s_bad, 1)
    s_bad['id'] = 1
    pytest.raises(RuntimeError, _assert_complete_surface, s_bad)
    s_bad = dict(tris=surf['tris'], ntri=surf['ntri'], rr=surf['rr'].copy())
    s_bad['rr'][0] = 0.
    pytest.raises(RuntimeError, _get_ico_map, surf, s_bad)

    surfs = read_bem_surfaces(fname_bem_3)
    pytest.raises(RuntimeError, _assert_inside, surfs[0], surfs[1])  # outside
    surfs[0]['id'] = 100  # bad surfs
    pytest.raises(RuntimeError, _order_surfaces, surfs)
    surfs[1]['rr'] /= 1000.
    pytest.raises(RuntimeError, _check_surface_size, surfs[1])

    # actually test functionality
    tempdir = _TempDir()
    fname_temp = op.join(tempdir, 'temp-bem-sol.fif')
    # use a model and solution made in Python
    conductivities = [(0.3,), (0.3, 0.006, 0.3)]
    fnames = [fname_bem_sol_1, fname_bem_sol_3]
    for cond, fname in zip(conductivities, fnames):
        for model_type in ('python', 'c'):
            if model_type == 'python':
                model = make_bem_model('sample', conductivity=cond, ico=2,
                                       subjects_dir=subjects_dir)
            else:
                model = fname_bem_1 if len(cond) == 1 else fname_bem_3
        solution = make_bem_solution(model)
        solution_c = read_bem_solution(fname)
        _compare_bem_solutions(solution, solution_c)
        write_bem_solution(fname_temp, solution)
        solution_read = read_bem_solution(fname_temp)
        _compare_bem_solutions(solution, solution_c)
        _compare_bem_solutions(solution_read, solution_c)

def test_make_scalp_surfaces():
    """Test mne make_scalp_surfaces."""
    check_usage(mne_make_scalp_surfaces)
    # Copy necessary files to avoid FreeSurfer call
    tempdir = _TempDir()
    surf_path = op.join(subjects_dir, 'sample', 'surf')
    surf_path_new = op.join(tempdir, 'sample', 'surf')
    os.mkdir(op.join(tempdir, 'sample'))
    os.mkdir(surf_path_new)
    subj_dir = op.join(tempdir, 'sample', 'bem')
    os.mkdir(subj_dir)
    shutil.copy(op.join(surf_path, 'lh.seghead'), surf_path_new)

    orig_fs = os.getenv('FREESURFER_HOME', None)
    if orig_fs is not None:
        del os.environ['FREESURFER_HOME']
    cmd = ('-s', 'sample', '--subjects-dir', tempdir)
    os.environ['_MNE_TESTING_SCALP'] = 'true'
    dense_fname = op.join(subj_dir, 'sample-head-dense.fif')
    medium_fname = op.join(subj_dir, 'sample-head-medium.fif')
    try:
        with ArgvSetter(cmd, disable_stdout=False, disable_stderr=False):
            pytest.raises(RuntimeError, mne_make_scalp_surfaces.run)
            os.environ['FREESURFER_HOME'] = tempdir  # don't actually use it
            mne_make_scalp_surfaces.run()
            assert op.isfile(dense_fname)
            assert op.isfile(medium_fname)
            pytest.raises(IOError, mne_make_scalp_surfaces.run)  # no overwrite
    finally:
        if orig_fs is not None:
            os.environ['FREESURFER_HOME'] = orig_fs
        else:
            del os.environ['FREESURFER_HOME']
        del os.environ['_MNE_TESTING_SCALP']
    # actually check the outputs
    head_py = read_bem_surfaces(dense_fname)
    assert_equal(len(head_py), 1)
    head_py = head_py[0]
    head_c = read_bem_surfaces(op.join(subjects_dir, 'sample', 'bem',
                                       'sample-head-dense.fif'))[0]
    assert_allclose(head_py['rr'], head_c['rr'])

def test_io_bem():
    """Test reading and writing of bem surfaces and solutions."""
    tempdir = _TempDir()
    temp_bem = op.join(tempdir, 'temp-bem.fif')
    pytest.raises(ValueError, read_bem_surfaces, fname_raw)
    pytest.raises(ValueError, read_bem_surfaces, fname_bem_3, s_id=10)
    surf = read_bem_surfaces(fname_bem_3, patch_stats=True)
    surf = read_bem_surfaces(fname_bem_3, patch_stats=False)
    write_bem_surfaces(temp_bem, surf[0])
    surf_read = read_bem_surfaces(temp_bem, patch_stats=False)
    _compare_bem_surfaces(surf, surf_read)

    pytest.raises(RuntimeError, read_bem_solution, fname_bem_3)
    temp_sol = op.join(tempdir, 'temp-sol.fif')
    sol = read_bem_solution(fname_bem_sol_3)
    assert 'BEM' in repr(sol)
    write_bem_solution(temp_sol, sol)
    sol_read = read_bem_solution(temp_sol)
    _compare_bem_solutions(sol, sol_read)
    sol = read_bem_solution(fname_bem_sol_1)
    pytest.raises(RuntimeError, _bem_find_surface, sol, 3)

def test_make_scalp_surfaces():
    """Test mne make_scalp_surfaces."""
    check_usage(mne_make_scalp_surfaces)
    # Copy necessary files to avoid FreeSurfer call
    tempdir = _TempDir()
    surf_path = op.join(subjects_dir, "sample", "surf")
    surf_path_new = op.join(tempdir, "sample", "surf")
    os.mkdir(op.join(tempdir, "sample"))
    os.mkdir(surf_path_new)
    subj_dir = op.join(tempdir, "sample", "bem")
    os.mkdir(subj_dir)
    shutil.copy(op.join(surf_path, "lh.seghead"), surf_path_new)

    orig_fs = os.getenv("FREESURFER_HOME", None)
    if orig_fs is not None:
        del os.environ["FREESURFER_HOME"]
    cmd = ("-s", "sample", "--subjects-dir", tempdir)
    os.environ["_MNE_TESTING_SCALP"] = "true"
    dense_fname = op.join(subj_dir, "sample-head-dense.fif")
    medium_fname = op.join(subj_dir, "sample-head-medium.fif")
    try:
        with ArgvSetter(cmd, disable_stdout=False, disable_stderr=False):
            assert_raises(RuntimeError, mne_make_scalp_surfaces.run)
            os.environ["FREESURFER_HOME"] = tempdir  # don't actually use it
            mne_make_scalp_surfaces.run()
            assert_true(op.isfile(dense_fname))
            assert_true(op.isfile(medium_fname))
            assert_raises(IOError, mne_make_scalp_surfaces.run)  # no overwrite
    finally:
        if orig_fs is not None:
            os.environ["FREESURFER_HOME"] = orig_fs
        else:
            del os.environ["FREESURFER_HOME"]
        del os.environ["_MNE_TESTING_SCALP"]
    # actually check the outputs
    head_py = read_bem_surfaces(dense_fname)
    assert_equal(len(head_py), 1)
    head_py = head_py[0]
    head_c = read_bem_surfaces(op.join(subjects_dir, "sample", "bem", "sample-head-dense.fif"))[0]
    assert_allclose(head_py["rr"], head_c["rr"])

def test_make_scalp_surfaces(tmpdir):
    """Test mne make_scalp_surfaces."""
    check_usage(mne_make_scalp_surfaces)
    has = 'SUBJECTS_DIR' in os.environ
    # Copy necessary files to avoid FreeSurfer call
    tempdir = str(tmpdir)
    surf_path = op.join(subjects_dir, 'sample', 'surf')
    surf_path_new = op.join(tempdir, 'sample', 'surf')
    os.mkdir(op.join(tempdir, 'sample'))
    os.mkdir(surf_path_new)
    subj_dir = op.join(tempdir, 'sample', 'bem')
    os.mkdir(subj_dir)
    shutil.copy(op.join(surf_path, 'lh.seghead'), surf_path_new)

    cmd = ('-s', 'sample', '--subjects-dir', tempdir)
    with modified_env(**{'_MNE_TESTING_SCALP': 'true'}):
        dense_fname = op.join(subj_dir, 'sample-head-dense.fif')
        medium_fname = op.join(subj_dir, 'sample-head-medium.fif')
        with ArgvSetter(cmd, disable_stdout=False, disable_stderr=False):
            with modified_env(FREESURFER_HOME=None):
                pytest.raises(RuntimeError, mne_make_scalp_surfaces.run)
            with modified_env(FREESURFER_HOME=tempdir):
                mne_make_scalp_surfaces.run()
                assert op.isfile(dense_fname)
                assert op.isfile(medium_fname)
                with pytest.raises(IOError, match='overwrite'):
                    mne_make_scalp_surfaces.run()
    # actually check the outputs
    head_py = read_bem_surfaces(dense_fname)
    assert_equal(len(head_py), 1)
    head_py = head_py[0]
    head_c = read_bem_surfaces(op.join(subjects_dir, 'sample', 'bem',
                                       'sample-head-dense.fif'))[0]
    assert_allclose(head_py['rr'], head_c['rr'])
    if not has:
        assert 'SUBJECTS_DIR' not in os.environ

    def __init__(self, bem, unit='m'):
        if isinstance(bem, basestring):
            bem = mne.read_bem_surfaces(bem)[0]

        pts = bem['rr']
        tri = bem['tris']

        if unit == 'mm':
            pts *= 1000
        elif unit == 'm':
            pass
        else:
            raise ValueError('Unit: %r' % unit)

        super(geom_bem, self).__init__(pts, tri)

def check_bem(sbj_id, sbj_dir, raw_info, trans_fname, report):
    import os.path as op
    import mne
    from mne.viz import plot_bem, plot_trans
    
    from mayavi import mlab
    import numpy as np
    
    ### plot bem surfs to MRI in the 3 different views
    fig1 = plot_bem(subject=sbj_id, subjects_dir=sbj_dir, orientation='axial', show=False)
    fig2 = plot_bem(subject=sbj_id, subjects_dir=sbj_dir, orientation='sagittal', show=False)
    fig3 = plot_bem(subject=sbj_id, subjects_dir=sbj_dir, orientation='coronal', show=False)
    report.add_figs_to_section([fig1, fig2, fig3], captions=['axial view', 'sagittal view', 'coronal view'], section='BEM surfaces')
    
    ### plot bem surf and source space
    bem_fname = op.join(sbj_dir, sbj_id + '/bem/%s-5120-bem-sol.fif' % sbj_id)
    surf = mne.read_bem_surfaces(bem_fname, patch_stats=True)
    print 'Number of surfaces : %d' % len(surf)
    brain_col = (0.95, 0.83, 0.83)
    cortex_col = (0.91, 0.89, 0.67)
    points = surf[0]['rr']
    faces = surf[0]['tris']
    fig4 = mlab.figure(size=(600, 600), bgcolor=(0, 0, 0))
    mlab.triangular_mesh(points[:, 0], points[:, 1], points[:, 2], faces, color=brain_col, opacity=0.3)
    
    src_fname = op.join(sbj_dir, sbj_id + '/bem/%s-ico-5-src.fif' % sbj_id)
    src = mne.read_source_spaces(src_fname)
    lh_points = src[0]['rr']
    lh_faces = src[0]['tris']
    rh_points = src[1]['rr']
    rh_faces = src[1]['tris']
    mlab.triangular_mesh(lh_points[:, 0], lh_points[:, 1], lh_points[:, 2], lh_faces, color=cortex_col, opacity=0.8)
    mlab.triangular_mesh(rh_points[:, 0], rh_points[:, 1], rh_points[:, 2], rh_faces, color=cortex_col, opacity=0.8)
    picks = mne.pick_types(raw_info, meg=True, ref_meg=False, eeg=False, stim=False, eog=False, exclude='bads')

    ### plot sensors    
    sens_loc = [ raw_info['chs'][picks[i]]['loc'][:3] for i in range(len(picks)) ]
    sens_loc = np.array(sens_loc)
    mlab.points3d(sens_loc[:, 0], sens_loc[:, 1], sens_loc[:, 2], color=(1, 1, 1), opacity=1, scale_factor=0.005)
    
    report.add_figs_to_section(fig4, captions=['source space'], section='BEM cortex sensors')
    fig5 = plot_trans(raw_info, trans_fname, subject=sbj_id, subjects_dir=sbj_dir)
    report.add_figs_to_section(fig5, captions=['MEG <-> MRI coregistration quality'], section='MEG <-> MRI')

def test_volume_source_space():
    """Test setting up volume source spaces
    """
    tempdir = _TempDir()
    src = read_source_spaces(fname_vol)
    temp_name = op.join(tempdir, 'temp-src.fif')
    surf = read_bem_surfaces(fname_bem, s_id=FIFF.FIFFV_BEM_SURF_ID_BRAIN)
    surf['rr'] *= 1e3  # convert to mm
    # The one in the testing dataset (uses bem as bounds)
    for bem, surf in zip((fname_bem, None), (None, surf)):
        src_new = setup_volume_source_space('sample', temp_name, pos=7.0,
                                            bem=bem, surface=surf,
                                            mri=fname_mri,
                                            subjects_dir=subjects_dir)
        _compare_source_spaces(src, src_new, mode='approx')
        del src_new
        src_new = read_source_spaces(temp_name)
        _compare_source_spaces(src, src_new, mode='approx')
    assert_raises(IOError, setup_volume_source_space, 'sample', temp_name,
                  pos=7.0, bem=None, surface='foo',  # bad surf
                  mri=fname_mri, subjects_dir=subjects_dir)

def plot_BEM_surface(subject, fnout_img=None):
    """"Reads in and plots the BEM surface."""

    # get name of BEM-file
    subjects_dir = os.environ.get('SUBJECTS_DIR')
    fn_bem = os.path.join(subjects_dir,
                          subject,
                          'bem',
                          subject + "-5120-5120-5120-bem-sol.fif")

    surfaces = mne.read_bem_surfaces(fn_bem, add_geom=True)

    print "Number of surfaces : %d" % len(surfaces)

    # Show result
    head_col = (0.95, 0.83, 0.83)  # light pink
    skull_col = (0.91, 0.89, 0.67)
    brain_col = (0.67, 0.89, 0.91)  # light blue
    colors = [head_col, skull_col, brain_col]

    # create figure and plot results
    fig_BEM = mlab.figure(size=(1200, 1200),
                          bgcolor=(0, 0, 0))
    for c, surf in zip(colors, surfaces):
        points = surf['rr']
        faces = surf['tris']
        mlab.triangular_mesh(points[:, 0],
                             points[:, 1],
                             points[:, 2],
                             faces,
                             color=c,
                             opacity=0.3)

    # save result
    if fnout_img is not None:
        mlab.savefig(fnout_img,
                     figure=fig_BEM,
                     size=(1200, 1200))

    mlab.close(all=True)

def test_volume_source_space():
    """Test setting up volume source spaces."""
    tempdir = _TempDir()
    src = read_source_spaces(fname_vol)
    temp_name = op.join(tempdir, 'temp-src.fif')
    surf = read_bem_surfaces(fname_bem, s_id=FIFF.FIFFV_BEM_SURF_ID_BRAIN)
    surf['rr'] *= 1e3  # convert to mm
    # The one in the testing dataset (uses bem as bounds)
    for bem, surf in zip((fname_bem, None), (None, surf)):
        src_new = setup_volume_source_space(
            'sample', pos=7.0, bem=bem, surface=surf, mri='T1.mgz',
            subjects_dir=subjects_dir)
        write_source_spaces(temp_name, src_new, overwrite=True)
        src[0]['subject_his_id'] = 'sample'  # XXX: to make comparison pass
        _compare_source_spaces(src, src_new, mode='approx')
        del src_new
        src_new = read_source_spaces(temp_name)
        _compare_source_spaces(src, src_new, mode='approx')
    assert_raises(IOError, setup_volume_source_space, 'sample',
                  pos=7.0, bem=None, surface='foo',  # bad surf
                  mri=fname_mri, subjects_dir=subjects_dir)
    assert_equal(repr(src), repr(src_new))
    assert_equal(src.kind, 'volume')

def _run(subjects_dir, subject, force, overwrite, verbose=None):
    this_env = copy.copy(os.environ)
    this_env['SUBJECTS_DIR'] = subjects_dir
    this_env['SUBJECT'] = subject

    if 'SUBJECTS_DIR' not in this_env:
        raise RuntimeError('The environment variable SUBJECTS_DIR should '
                           'be set')

    if not op.isdir(subjects_dir):
        raise RuntimeError('subjects directory %s not found, specify using '
                           'the environment variable SUBJECTS_DIR or '
                           'the command line option --subjects-dir')

    if 'MNE_ROOT' not in this_env:
        raise RuntimeError('MNE_ROOT environment variable is not set')

    if 'FREESURFER_HOME' not in this_env:
        raise RuntimeError('The FreeSurfer environment needs to be set up '
                           'for this script')
    force = '--force' if force else '--check'
    subj_path = op.join(subjects_dir, subject)
    if not op.exists(subj_path):
        raise RuntimeError('%s does not exits. Please check your subject '
                           'directory path.' % subj_path)

    if op.exists(op.join(subj_path, 'mri', 'T1.mgz')):
        mri = 'T1.mgz'
    else:
        mri = 'T1'

    logger.info('1. Creating a dense scalp tessellation with mkheadsurf...')

    def check_seghead(surf_path=op.join(subj_path, 'surf')):
        for k in ['/lh.seghead', '/lh.smseghead']:
            surf = surf_path + k if op.exists(surf_path + k) else None
            if surf is not None:
                break
        return surf

    my_seghead = check_seghead()
    if my_seghead is None:
        run_subprocess(['mkheadsurf', '-subjid', subject, '-srcvol', mri],
                       env=this_env)

    surf = check_seghead()
    if surf is None:
        raise RuntimeError('mkheadsurf did not produce the standard output '
                           'file.')

    dense_fname = '{0}/{1}/bem/{1}-head-dense.fif'.format(subjects_dir,
                                                          subject)
    logger.info('2. Creating %s ...' % dense_fname)
    _check_file(dense_fname, overwrite)
    run_subprocess(['mne_surf2bem', '--surf', surf, '--id', '4', force,
                    '--fif', dense_fname], env=this_env)
    levels = 'medium', 'sparse'
    my_surf = mne.read_bem_surfaces(dense_fname)[0]
    tris = [30000, 2500]
    if os.getenv('_MNE_TESTING_SCALP', 'false') == 'true':
        tris = [len(my_surf['tris'])]  # don't actually decimate
    for ii, (n_tri, level) in enumerate(zip(tris, levels), 3):
        logger.info('%i. Creating %s tessellation...' % (ii, level))
        logger.info('%i.1 Decimating the dense tessellation...' % ii)
        points, tris = mne.decimate_surface(points=my_surf['rr'],
                                            triangles=my_surf['tris'],
                                            n_triangles=n_tri)
        other_fname = dense_fname.replace('dense', level)
        logger.info('%i.2 Creating %s' % (ii, other_fname))
        _check_file(other_fname, overwrite)
        tempdir = _TempDir()
        surf_fname = tempdir + '/tmp-surf.surf'
        # convert points to meters, make mne_analyze happy
        mne.write_surface(surf_fname, points * 1e3, tris)
        # XXX for some reason --check does not work here.
        try:
            run_subprocess(['mne_surf2bem', '--surf', surf_fname, '--id', '4',
                            '--force', '--fif', other_fname], env=this_env)
        finally:
            del tempdir

# Author: Alexandre Gramfort <[email protected]>
#
# License: BSD (3-clause)

print(__doc__)

import numpy as np
import mne
from mne.datasets import sample

data_path = sample.data_path()
fwd_fname = data_path + '/MEG/sample/sample_audvis-meg-oct-6-fwd.fif'
dip_fname = data_path + '/MEG/sample/sample_audvis_set1.dip'
bem_fname = data_path + '/subjects/sample/bem/sample-5120-bem-sol.fif'

brain_surface = mne.read_bem_surfaces(bem_fname, add_geom=True)[0]
points = brain_surface['rr']
faces = brain_surface['tris']

fwd = mne.read_forward_solution(fwd_fname)
src = fwd['src']

# read dipoles
time, pos, amplitude, ori, gof = mne.read_dip(dip_fname)

print("Time (ms): %s" % time)
print("Amplitude (nAm): %s" % amplitude)
print("GOF (%%): %s" % gof)

# only plot those for which GOF is above 50%
pos = pos[gof > 50.]

fsaverage_path = op.join(op.dirname(mne.__file__), 'data', 'fsaverage')
data_path = mne.datasets.sample.data_path()
subjects_dir = op.join(data_path, 'subjects')

###############################################################################
# Load the source digitization, destination surfaces, and source surfaces
# (with everything in head coordinates):

# Digitization
info = mne.io.read_info(op.join(data_path, 'MEG', 'sample',
                                'sample_audvis_raw.fif'))
hsp = mne.bem.get_fitting_dig(info, ('cardinal', 'extra'))

# Destination head surface
fsaverage_surfs = [mne.read_bem_surfaces(op.join(fsaverage_path,
                                                 'fsaverage-%s.fif' % kind))[0]
                   for kind in ('head', 'inner_skull-bem')]
fsaverage_trans = mne.read_trans(op.join(fsaverage_path,
                                         'fsaverage-trans.fif'))
for surf in fsaverage_surfs:
    mne.surface.transform_surface_to(surf, 'head', fsaverage_trans, copy=False)

# Some source surfaces to transform as examples
sample_surfs = mne.read_bem_surfaces(
    op.join(subjects_dir, 'sample', 'bem', 'sample-5120-bem.fif'))
sample_trans = mne.read_trans(op.join(data_path, 'MEG', 'sample',
                                      'sample_audvis_raw-trans.fif'))
for surf in sample_surfs:
    mne.surface.transform_surface_to(surf, 'head', sample_trans, copy=False)

###############################################################################

def test_io_bem_surfaces():
    """Testing reading of bem surfaces
    """
    surf = mne.read_bem_surfaces(fname, add_geom=False)
    surf = mne.read_bem_surfaces(fname, add_geom=True)
    print "Number of surfaces : %d" % len(surf)

using a cloud of digitization points for coordinate alignment
instead of e.g. EEG-cap positions.

"""
print(__doc__)

# Authors: Denis Engemann <[email protected]>
#          Alexandre Gramfort <[email protected]>
#
# License: BSD (3-clause)

import mne
from mne.surface import decimate_surface

path = mne.datasets.sample.data_path()
surf = mne.read_bem_surfaces(path + '/subjects/sample/bem/sample-head.fif')[0]

points, triangles = surf['rr'], surf['tris']

# reduce to 30000 meshes equaling ${SUBJECT}-head-medium.fif output from
# mne_make_scalp_surfaces.py and mne_make_scalp_surfaces
points_dec, triangles_dec = decimate_surface(points, triangles,
                                             n_triangles=30000)

try:
    from enthought.mayavi import mlab
except:
    from mayavi import mlab

head_col = (0.95, 0.83, 0.83)  # light pink