bem.py

# Authors: Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>
#          Matti Hamalainen <msh@nmr.mgh.harvard.edu>
#          Eric Larson <larson.eric.d@gmail.com>
#          Lorenzo De Santis <lorenzo.de-santis@u-psud.fr>
#
# License: BSD (3-clause)

from functools import partial
import glob
import os
import os.path as op
import shutil
from copy import deepcopy

import numpy as np
from scipy import linalg

from .io.constants import FIFF, FWD
from .digitization.base import _dig_kind_dict, _dig_kind_rev, _dig_kind_ints
from .io.write import (start_file, start_block, write_float, write_int,
                       write_float_matrix, write_int_matrix, end_block,
                       end_file)
from .io.tag import find_tag
from .io.tree import dir_tree_find
from .io.open import fiff_open
from .surface import (read_surface, write_surface, complete_surface_info,
                      _compute_nearest, _get_ico_surface, read_tri,
                      _fast_cross_nd_sum, _get_solids)
from .transforms import _ensure_trans, apply_trans
from .utils import (verbose, logger, run_subprocess, get_subjects_dir, warn,
                    _pl, _validate_type, _TempDir)
from .fixes import einsum


# ############################################################################
# Compute BEM solution

# The following approach is based on:
#
# de Munck JC: "A linear discretization of the volume conductor boundary
# integral equation using analytically integrated elements",
# IEEE Trans Biomed Eng. 1992 39(9) : 986 - 990
#


class ConductorModel(dict):
    """BEM or sphere model."""

    def __repr__(self):  # noqa: D105
        if self['is_sphere']:
            center = ', '.join('%0.1f' % (x * 1000.) for x in self['r0'])
            rad = self.radius
            if rad is None:  # no radius / MEG only
                extra = 'Sphere (no layers): r0=[%s] mm' % center
            else:
                extra = ('Sphere (%s layer%s): r0=[%s] R=%1.f mm'
                         % (len(self['layers']) - 1, _pl(self['layers']),
                            center, rad * 1000.))
        else:
            extra = ('BEM (%s layer%s)' % (len(self['surfs']),
                                           _pl(self['surfs'])))
        return '<ConductorModel  |  %s>' % extra

    def copy(self):
        """Return copy of ConductorModel instance."""
        return deepcopy(self)

    @property
    def radius(self):
        """Sphere radius if an EEG sphere model."""
        if not self['is_sphere']:
            raise RuntimeError('radius undefined for BEM')
        return None if len(self['layers']) == 0 else self['layers'][-1]['rad']


def _calc_beta(rk, rk_norm, rk1, rk1_norm):
    """Compute coefficients for calculating the magic vector omega."""
    rkk1 = rk1[0] - rk[0]
    size = np.linalg.norm(rkk1)
    rkk1 /= size
    num = rk_norm + np.dot(rk, rkk1)
    den = rk1_norm + np.dot(rk1, rkk1)
    res = np.log(num / den) / size
    return res


def _lin_pot_coeff(fros, tri_rr, tri_nn, tri_area):
    """Compute the linear potential matrix element computations."""
    omega = np.zeros((len(fros), 3))

    # we replicate a little bit of the _get_solids code here for speed
    # (we need some of the intermediate values later)
    v1 = tri_rr[np.newaxis, 0, :] - fros
    v2 = tri_rr[np.newaxis, 1, :] - fros
    v3 = tri_rr[np.newaxis, 2, :] - fros
    triples = _fast_cross_nd_sum(v1, v2, v3)
    l1 = np.linalg.norm(v1, axis=1)
    l2 = np.linalg.norm(v2, axis=1)
    l3 = np.linalg.norm(v3, axis=1)
    ss = l1 * l2 * l3
    ss += einsum('ij,ij,i->i', v1, v2, l3)
    ss += einsum('ij,ij,i->i', v1, v3, l2)
    ss += einsum('ij,ij,i->i', v2, v3, l1)
    solids = np.arctan2(triples, ss)

    # We *could* subselect the good points from v1, v2, v3, triples, solids,
    # l1, l2, and l3, but there are *very* few bad points. So instead we do
    # some unnecessary calculations, and then omit them from the final
    # solution. These three lines ensure we don't get invalid values in
    # _calc_beta.
    bad_mask = np.abs(solids) < np.pi / 1e6
    l1[bad_mask] = 1.
    l2[bad_mask] = 1.
    l3[bad_mask] = 1.

    # Calculate the magic vector vec_omega
    beta = [_calc_beta(v1, l1, v2, l2)[:, np.newaxis],
            _calc_beta(v2, l2, v3, l3)[:, np.newaxis],
            _calc_beta(v3, l3, v1, l1)[:, np.newaxis]]
    vec_omega = (beta[2] - beta[0]) * v1
    vec_omega += (beta[0] - beta[1]) * v2
    vec_omega += (beta[1] - beta[2]) * v3

    area2 = 2.0 * tri_area
    n2 = 1.0 / (area2 * area2)
    # leave omega = 0 otherwise
    # Put it all together...
    yys = [v1, v2, v3]
    idx = [0, 1, 2, 0, 2]
    for k in range(3):
        diff = yys[idx[k - 1]] - yys[idx[k + 1]]
        zdots = _fast_cross_nd_sum(yys[idx[k + 1]], yys[idx[k - 1]], tri_nn)
        omega[:, k] = -n2 * (area2 * zdots * 2. * solids -
                             triples * (diff * vec_omega).sum(axis=-1))
    # omit the bad points from the solution
    omega[bad_mask] = 0.
    return omega


def _correct_auto_elements(surf, mat):
    """Improve auto-element approximation."""
    pi2 = 2.0 * np.pi
    tris_flat = surf['tris'].ravel()
    misses = pi2 - mat.sum(axis=1)
    for j, miss in enumerate(misses):
        # How much is missing?
        n_memb = len(surf['neighbor_tri'][j])
        assert n_memb > 0  # should be guaranteed by our surface checks
        # The node itself receives one half
        mat[j, j] = miss / 2.0
        # The rest is divided evenly among the member nodes...
        miss /= (4.0 * n_memb)
        members = np.where(j == tris_flat)[0]
        mods = members % 3
        offsets = np.array([[1, 2], [-1, 1], [-1, -2]])
        tri_1 = members + offsets[mods, 0]
        tri_2 = members + offsets[mods, 1]
        for t1, t2 in zip(tri_1, tri_2):
            mat[j, tris_flat[t1]] += miss
            mat[j, tris_flat[t2]] += miss
    return


def _fwd_bem_lin_pot_coeff(surfs):
    """Calculate the coefficients for linear collocation approach."""
    # taken from fwd_bem_linear_collocation.c
    nps = [surf['np'] for surf in surfs]
    np_tot = sum(nps)
    coeff = np.zeros((np_tot, np_tot))
    offsets = np.cumsum(np.concatenate(([0], nps)))
    for si_1, surf1 in enumerate(surfs):
        rr_ord = np.arange(nps[si_1])
        for si_2, surf2 in enumerate(surfs):
            logger.info("        %s (%d) -> %s (%d) ..." %
                        (_surf_name[surf1['id']], nps[si_1],
                         _surf_name[surf2['id']], nps[si_2]))
            tri_rr = surf2['rr'][surf2['tris']]
            tri_nn = surf2['tri_nn']
            tri_area = surf2['tri_area']
            submat = coeff[offsets[si_1]:offsets[si_1 + 1],
                           offsets[si_2]:offsets[si_2 + 1]]  # view
            for k in range(surf2['ntri']):
                tri = surf2['tris'][k]
                if si_1 == si_2:
                    skip_idx = ((rr_ord == tri[0]) |
                                (rr_ord == tri[1]) |
                                (rr_ord == tri[2]))
                else:
                    skip_idx = list()
                # No contribution from a triangle that
                # this vertex belongs to
                # if sidx1 == sidx2 and (tri == j).any():
                #     continue
                # Otherwise do the hard job
                coeffs = _lin_pot_coeff(surf1['rr'], tri_rr[k], tri_nn[k],
                                        tri_area[k])
                coeffs[skip_idx] = 0.
                submat[:, tri] -= coeffs
            if si_1 == si_2:
                _correct_auto_elements(surf1, submat)
    return coeff


def _fwd_bem_multi_solution(solids, gamma, nps):
    """Do multi surface solution.

    * Invert I - solids/(2*M_PI)
    * Take deflation into account
    * The matrix is destroyed after inversion
    * This is the general multilayer case
    """
    pi2 = 1.0 / (2 * np.pi)
    n_tot = np.sum(nps)
    assert solids.shape == (n_tot, n_tot)
    nsurf = len(nps)
    defl = 1.0 / n_tot
    # Modify the matrix
    offsets = np.cumsum(np.concatenate(([0], nps)))
    for si_1 in range(nsurf):
        for si_2 in range(nsurf):
            mult = pi2 if gamma is None else pi2 * gamma[si_1, si_2]
            slice_j = slice(offsets[si_1], offsets[si_1 + 1])
            slice_k = slice(offsets[si_2], offsets[si_2 + 1])
            solids[slice_j, slice_k] = defl - solids[slice_j, slice_k] * mult
    solids += np.eye(n_tot)
    return linalg.inv(solids, overwrite_a=True)


def _fwd_bem_homog_solution(solids, nps):
    """Make a homogeneous solution."""
    return _fwd_bem_multi_solution(solids, None, nps)


def _fwd_bem_ip_modify_solution(solution, ip_solution, ip_mult, n_tri):
    """Modify the solution according to the IP approach."""
    n_last = n_tri[-1]
    mult = (1.0 + ip_mult) / ip_mult

    logger.info('        Combining...')
    offsets = np.cumsum(np.concatenate(([0], n_tri)))
    for si in range(len(n_tri)):
        # Pick the correct submatrix (right column) and multiply
        sub = solution[offsets[si]:offsets[si + 1], np.sum(n_tri[:-1]):]
        # Multiply
        sub -= 2 * np.dot(sub, ip_solution)

    # The lower right corner is a special case
    sub[-n_last:, -n_last:] += mult * ip_solution

    # Final scaling
    logger.info('        Scaling...')
    solution *= ip_mult
    return


def _check_complete_surface(surf, copy=False, incomplete='raise'):
    surf = complete_surface_info(surf, copy=copy, verbose=False)
    fewer = np.where([len(t) < 3 for t in surf['neighbor_tri']])[0]
    if len(fewer) > 0:
        msg = ('Surface %s has topological defects: %d / %d '
               'vertices have fewer than three neighboring '
               'triangles [%s]' % (_surf_name[surf['id']],
                                   len(fewer), surf['ntri'],
                                   ', '.join(str(f) for f in fewer)))
        if incomplete == 'raise':
            raise RuntimeError(msg)
        else:
            warn(msg)
    return surf


def _fwd_bem_linear_collocation_solution(m):
    """Compute the linear collocation potential solution."""
    # first, add surface geometries
    for surf in m['surfs']:
        _check_complete_surface(surf)

    logger.info('Computing the linear collocation solution...')
    logger.info('    Matrix coefficients...')
    coeff = _fwd_bem_lin_pot_coeff(m['surfs'])
    m['nsol'] = len(coeff)
    logger.info("    Inverting the coefficient matrix...")
    nps = [surf['np'] for surf in m['surfs']]
    m['solution'] = _fwd_bem_multi_solution(coeff, m['gamma'], nps)
    if len(m['surfs']) == 3:
        ip_mult = m['sigma'][1] / m['sigma'][2]
        if ip_mult <= FWD.BEM_IP_APPROACH_LIMIT:
            logger.info('IP approach required...')
            logger.info('    Matrix coefficients (homog)...')
            coeff = _fwd_bem_lin_pot_coeff([m['surfs'][-1]])
            logger.info('    Inverting the coefficient matrix (homog)...')
            ip_solution = _fwd_bem_homog_solution(coeff,
                                                  [m['surfs'][-1]['np']])
            logger.info('    Modify the original solution to incorporate '
                        'IP approach...')
            _fwd_bem_ip_modify_solution(m['solution'], ip_solution, ip_mult,
                                        nps)
    m['bem_method'] = FWD.BEM_LINEAR_COLL
    logger.info("Solution ready.")


@verbose
def make_bem_solution(surfs, verbose=None):
    """Create a BEM solution using the linear collocation approach.

    Parameters
    ----------
    surfs : list of dict
        The BEM surfaces to use (`from make_bem_model`)
    %(verbose)s

    Returns
    -------
    bem : instance of ConductorModel
        The BEM solution.

    Notes
    -----
    .. versionadded:: 0.10.0

    See Also
    --------
    make_bem_model
    read_bem_surfaces
    write_bem_surfaces
    read_bem_solution
    write_bem_solution
    """
    logger.info('Approximation method : Linear collocation\n')
    if isinstance(surfs, str):
        # Load the surfaces
        logger.info('Loading surfaces...')
        surfs = read_bem_surfaces(surfs)
    bem = ConductorModel(is_sphere=False, surfs=surfs)
    _add_gamma_multipliers(bem)
    if len(bem['surfs']) == 3:
        logger.info('Three-layer model surfaces loaded.')
    elif len(bem['surfs']) == 1:
        logger.info('Homogeneous model surface loaded.')
    else:
        raise RuntimeError('Only 1- or 3-layer BEM computations supported')
    _check_bem_size(bem['surfs'])
    _fwd_bem_linear_collocation_solution(bem)
    logger.info('BEM geometry computations complete.')
    return bem


# ############################################################################
# Make BEM model

def _ico_downsample(surf, dest_grade):
    """Downsample the surface if isomorphic to a subdivided icosahedron."""
    n_tri = len(surf['tris'])
    bad_msg = ("Cannot decimate to requested ico grade %d. The provided "
               "BEM surface has %d triangles, which cannot be isomorphic with "
               "a subdivided icosahedron. Consider manually decimating the "
               "surface to a suitable density and then use ico=None in "
               "make_bem_model." % (dest_grade, n_tri))
    if n_tri % 20 != 0:
        raise RuntimeError(bad_msg)
    n_tri = n_tri // 20
    found = int(round(np.log(n_tri) / np.log(4)))
    if n_tri != 4 ** found:
        raise RuntimeError(bad_msg)
    del n_tri

    if dest_grade > found:
        raise RuntimeError('For this surface, decimation grade should be %d '
                           'or less, not %s.' % (found, dest_grade))

    source = _get_ico_surface(found)
    dest = _get_ico_surface(dest_grade, patch_stats=True)
    del dest['tri_cent']
    del dest['tri_nn']
    del dest['neighbor_tri']
    del dest['tri_area']
    if not np.array_equal(source['tris'], surf['tris']):
        raise RuntimeError('The source surface has a matching number of '
                           'triangles but ordering is wrong')
    logger.info('Going from %dth to %dth subdivision of an icosahedron '
                '(n_tri: %d -> %d)' % (found, dest_grade, len(surf['tris']),
                                       len(dest['tris'])))
    # Find the mapping
    dest['rr'] = surf['rr'][_get_ico_map(source, dest)]
    return dest


def _get_ico_map(fro, to):
    """Get a mapping between ico surfaces."""
    nearest, dists = _compute_nearest(fro['rr'], to['rr'], return_dists=True)
    n_bads = (dists > 5e-3).sum()
    if n_bads > 0:
        raise RuntimeError('No matching vertex for %d destination vertices'
                           % (n_bads))
    return nearest


def _order_surfaces(surfs):
    """Reorder the surfaces."""
    if len(surfs) != 3:
        return surfs
    # we have three surfaces
    surf_order = [FIFF.FIFFV_BEM_SURF_ID_HEAD,
                  FIFF.FIFFV_BEM_SURF_ID_SKULL,
                  FIFF.FIFFV_BEM_SURF_ID_BRAIN]
    ids = np.array([surf['id'] for surf in surfs])
    if set(ids) != set(surf_order):
        raise RuntimeError('bad surface ids: %s' % ids)
    order = [np.where(ids == id_)[0][0] for id_ in surf_order]
    surfs = [surfs[idx] for idx in order]
    return surfs


def _assert_complete_surface(surf, incomplete='raise'):
    """Check the sum of solid angles as seen from inside."""
    # from surface_checks.c
    # Center of mass....
    cm = surf['rr'].mean(axis=0)
    logger.info('%s CM is %6.2f %6.2f %6.2f mm' %
                (_surf_name[surf['id']],
                 1000 * cm[0], 1000 * cm[1], 1000 * cm[2]))
    tot_angle = _get_solids(surf['rr'][surf['tris']], cm[np.newaxis, :])[0]
    prop = tot_angle / (2 * np.pi)
    if np.abs(prop - 1.0) > 1e-5:
        msg = ('Surface %s is not complete (sum of solid angles '
               'yielded %g, should be 1.)'
               % (_surf_name[surf['id']], prop))
        if incomplete == 'raise':
            raise RuntimeError(msg)
        else:
            warn(msg)


_surf_name = {
    FIFF.FIFFV_BEM_SURF_ID_HEAD: 'outer skin ',
    FIFF.FIFFV_BEM_SURF_ID_SKULL: 'outer skull',
    FIFF.FIFFV_BEM_SURF_ID_BRAIN: 'inner skull',
    FIFF.FIFFV_BEM_SURF_ID_UNKNOWN: 'unknown    ',
}


def _assert_inside(fro, to):
    """Check one set of points is inside a surface."""
    # this is "is_inside" in surface_checks.c
    tot_angle = _get_solids(to['rr'][to['tris']], fro['rr'])
    if (np.abs(tot_angle / (2 * np.pi) - 1.0) > 1e-5).any():
        raise RuntimeError('Surface %s is not completely inside surface %s'
                           % (_surf_name[fro['id']], _surf_name[to['id']]))


def _check_surfaces(surfs, incomplete='raise'):
    """Check that the surfaces are complete and non-intersecting."""
    for surf in surfs:
        _assert_complete_surface(surf, incomplete=incomplete)
    # Then check the topology
    for surf_1, surf_2 in zip(surfs[:-1], surfs[1:]):
        logger.info('Checking that %s surface is inside %s surface...' %
                    (_surf_name[surf_2['id']], _surf_name[surf_1['id']]))
        _assert_inside(surf_2, surf_1)


def _check_surface_size(surf):
    """Check that the coordinate limits are reasonable."""
    sizes = surf['rr'].max(axis=0) - surf['rr'].min(axis=0)
    if (sizes < 0.05).any():
        raise RuntimeError('Dimensions of the surface %s seem too small '
                           '(%9.5f mm). Maybe the the unit of measure is '
                           'meters instead of mm' %
                           (_surf_name[surf['id']], 1000 * sizes.min()))


def _check_thicknesses(surfs):
    """Compute how close we are."""
    for surf_1, surf_2 in zip(surfs[:-1], surfs[1:]):
        min_dist = _compute_nearest(surf_1['rr'], surf_2['rr'],
                                    return_dists=True)[0]
        min_dist = min_dist.min()
        logger.info('Checking distance between %s and %s surfaces...' %
                    (_surf_name[surf_1['id']], _surf_name[surf_2['id']]))
        logger.info('Minimum distance between the %s and %s surfaces is '
                    'approximately %6.1f mm' %
                    (_surf_name[surf_1['id']], _surf_name[surf_2['id']],
                     1000 * min_dist))


def _surfaces_to_bem(surfs, ids, sigmas, ico=None, rescale=True,
                     incomplete='raise'):
    """Convert surfaces to a BEM."""
    # equivalent of mne_surf2bem
    # surfs can be strings (filenames) or surface dicts
    if len(surfs) not in (1, 3) or not (len(surfs) == len(ids) ==
                                        len(sigmas)):
        raise ValueError('surfs, ids, and sigmas must all have the same '
                         'number of elements (1 or 3)')
    for si, surf in enumerate(surfs):
        if isinstance(surf, str):
            surfs[si] = read_surface(surf, return_dict=True)[-1]
    # Downsampling if the surface is isomorphic with a subdivided icosahedron
    if ico is not None:
        for si, surf in enumerate(surfs):
            surfs[si] = _ico_downsample(surf, ico)
    for surf, id_ in zip(surfs, ids):
        # Do topology checks (but don't save data) to fail early
        surf['id'] = id_
        _check_complete_surface(surf, copy=True, incomplete=incomplete)
        surf['coord_frame'] = surf.get('coord_frame', FIFF.FIFFV_COORD_MRI)
        surf.update(np=len(surf['rr']), ntri=len(surf['tris']))
        if rescale:
            surf['rr'] /= 1000.  # convert to meters

    # Shifting surfaces is not implemented here...

    # Order the surfaces for the benefit of the topology checks
    for surf, sigma in zip(surfs, sigmas):
        surf['sigma'] = sigma
    surfs = _order_surfaces(surfs)

    # Check topology as best we can
    _check_surfaces(surfs, incomplete=incomplete)
    for surf in surfs:
        _check_surface_size(surf)
    _check_thicknesses(surfs)
    logger.info('Surfaces passed the basic topology checks.')
    return surfs


@verbose
def make_bem_model(subject, ico=4, conductivity=(0.3, 0.006, 0.3),
                   subjects_dir=None, verbose=None):
    """Create a BEM model for a subject.

    .. note:: To get a single layer bem corresponding to the --homog flag in
              the command line tool set the ``conductivity`` parameter
              to a list/tuple with a single value (e.g. [0.3]).

    Parameters
    ----------
    subject : str
        The subject.
    ico : int | None
        The surface ico downsampling to use, e.g. 5=20484, 4=5120, 3=1280.
        If None, no subsampling is applied.
    conductivity : array of int, shape (3,) or (1,)
        The conductivities to use for each shell. Should be a single element
        for a one-layer model, or three elements for a three-layer model.
        Defaults to ``[0.3, 0.006, 0.3]``. The MNE-C default for a
        single-layer model would be ``[0.3]``.
    subjects_dir : string, or None
        Path to SUBJECTS_DIR if it is not set in the environment.
    %(verbose)s

    Returns
    -------
    surfaces : list of dict
        The BEM surfaces. Use `make_bem_solution` to turn these into a
        `ConductorModel` suitable for forward calculation.

    Notes
    -----
    .. versionadded:: 0.10.0

    See Also
    --------
    make_bem_solution
    make_sphere_model
    read_bem_surfaces
    write_bem_surfaces
    """
    conductivity = np.array(conductivity, float)
    if conductivity.ndim != 1 or conductivity.size not in (1, 3):
        raise ValueError('conductivity must be 1D array-like with 1 or 3 '
                         'elements')
    subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
    subject_dir = op.join(subjects_dir, subject)
    bem_dir = op.join(subject_dir, 'bem')
    inner_skull = op.join(bem_dir, 'inner_skull.surf')
    outer_skull = op.join(bem_dir, 'outer_skull.surf')
    outer_skin = op.join(bem_dir, 'outer_skin.surf')
    surfaces = [inner_skull, outer_skull, outer_skin]
    ids = [FIFF.FIFFV_BEM_SURF_ID_BRAIN,
           FIFF.FIFFV_BEM_SURF_ID_SKULL,
           FIFF.FIFFV_BEM_SURF_ID_HEAD]
    logger.info('Creating the BEM geometry...')
    if len(conductivity) == 1:
        surfaces = surfaces[:1]
        ids = ids[:1]
    surfaces = _surfaces_to_bem(surfaces, ids, conductivity, ico)
    _check_bem_size(surfaces)
    logger.info('Complete.\n')
    return surfaces


# ############################################################################
# Compute EEG sphere model

def _fwd_eeg_get_multi_sphere_model_coeffs(m, n_terms):
    """Get the model depended weighting factor for n."""
    nlayer = len(m['layers'])
    if nlayer in (0, 1):
        return 1.

    # Initialize the arrays
    c1 = np.zeros(nlayer - 1)
    c2 = np.zeros(nlayer - 1)
    cr = np.zeros(nlayer - 1)
    cr_mult = np.zeros(nlayer - 1)
    for k in range(nlayer - 1):
        c1[k] = m['layers'][k]['sigma'] / m['layers'][k + 1]['sigma']
        c2[k] = c1[k] - 1.0
        cr_mult[k] = m['layers'][k]['rel_rad']
        cr[k] = cr_mult[k]
        cr_mult[k] *= cr_mult[k]

    coeffs = np.zeros(n_terms - 1)
    for n in range(1, n_terms):
        # Increment the radius coefficients
        for k in range(nlayer - 1):
            cr[k] *= cr_mult[k]

        # Multiply the matrices
        M = np.eye(2)
        n1 = n + 1.0
        for k in range(nlayer - 2, -1, -1):
            M = np.dot([[n + n1 * c1[k], n1 * c2[k] / cr[k]],
                        [n * c2[k] * cr[k], n1 + n * c1[k]]], M)
        num = n * (2.0 * n + 1.0) ** (nlayer - 1)
        coeffs[n - 1] = num / (n * M[1, 1] + n1 * M[1, 0])
    return coeffs


def _compose_linear_fitting_data(mu, u):
    """Get the linear fitting data."""
    k1 = np.arange(1, u['nterms'])
    mu1ns = mu[0] ** k1
    # data to be fitted
    y = u['w'][:-1] * (u['fn'][1:] - mu1ns * u['fn'][0])
    # model matrix
    M = u['w'][:-1, np.newaxis] * (mu[1:] ** k1[:, np.newaxis] -
                                   mu1ns[:, np.newaxis])
    uu, sing, vv = linalg.svd(M, full_matrices=False)
    ncomp = u['nfit'] - 1
    uu, sing, vv = uu[:, :ncomp], sing[:ncomp], vv[:ncomp]
    return y, uu, sing, vv


def _compute_linear_parameters(mu, u):
    """Compute the best-fitting linear parameters."""
    y, uu, sing, vv = _compose_linear_fitting_data(mu, u)

    # Compute the residuals
    vec = np.dot(y, uu)
    resi = y - np.dot(uu, vec)
    vec /= sing

    lambda_ = np.zeros(u['nfit'])
    lambda_[1:] = np.dot(vec, vv)
    lambda_[0] = u['fn'][0] - np.sum(lambda_[1:])
    rv = np.dot(resi, resi) / np.dot(y, y)
    return rv, lambda_


def _one_step(mu, u):
    """Evaluate the residual sum of squares fit for one set of mu values."""
    if np.abs(mu).max() > 1.0:
        return 1.0

    # Compose the data for the linear fitting, compute SVD, then residuals
    y, uu, sing, vv = _compose_linear_fitting_data(mu, u)
    resi = y - np.dot(uu, np.dot(y, uu))
    return np.dot(resi, resi)


def _fwd_eeg_fit_berg_scherg(m, nterms, nfit):
    """Fit the Berg-Scherg equivalent spherical model dipole parameters."""
    from scipy.optimize import fmin_cobyla
    assert nfit >= 2
    u = dict(nfit=nfit, nterms=nterms)

    # (1) Calculate the coefficients of the true expansion
    u['fn'] = _fwd_eeg_get_multi_sphere_model_coeffs(m, nterms + 1)

    # (2) Calculate the weighting
    f = (min([layer['rad'] for layer in m['layers']]) /
         max([layer['rad'] for layer in m['layers']]))

    # correct weighting
    k = np.arange(1, nterms + 1)
    u['w'] = np.sqrt((2.0 * k + 1) * (3.0 * k + 1.0) /
                     k) * np.power(f, (k - 1.0))
    u['w'][-1] = 0

    # Do the nonlinear minimization, constraining mu to the interval [-1, +1]
    mu_0 = np.zeros(3)
    fun = partial(_one_step, u=u)
    max_ = 1. - 2e-4  # adjust for fmin_cobyla "catol" that not all scipy have
    cons = list()
    for ii in range(nfit):
        def mycon(x, ii=ii):
            return max_ - np.abs(x[ii])
        cons.append(mycon)
    mu = fmin_cobyla(fun, mu_0, cons, rhobeg=0.5, rhoend=1e-5, disp=0)

    # (6) Do the final step: calculation of the linear parameters
    rv, lambda_ = _compute_linear_parameters(mu, u)
    order = np.argsort(mu)[::-1]
    mu, lambda_ = mu[order], lambda_[order]  # sort: largest mu first

    m['mu'] = mu
    # This division takes into account the actual conductivities
    m['lambda'] = lambda_ / m['layers'][-1]['sigma']
    m['nfit'] = nfit
    return rv


@verbose
def make_sphere_model(r0=(0., 0., 0.04), head_radius=0.09, info=None,
                      relative_radii=(0.90, 0.92, 0.97, 1.0),
                      sigmas=(0.33, 1.0, 0.004, 0.33), verbose=None):
    """Create a spherical model for forward solution calculation.

    Parameters
    ----------
    r0 : array-like | str
        Head center to use (in head coordinates). If 'auto', the head
        center will be calculated from the digitization points in info.
    head_radius : float | str | None
        If float, compute spherical shells for EEG using the given radius.
        If 'auto', estimate an appropriate radius from the dig points in Info,
        If None, exclude shells (single layer sphere model).
    info : instance of Info | None
        Measurement info. Only needed if ``r0`` or ``head_radius`` are
        ``'auto'``.
    relative_radii : array-like
        Relative radii for the spherical shells.
    sigmas : array-like
        Sigma values for the spherical shells.
    %(verbose)s

    Returns
    -------
    sphere : instance of ConductorModel
        The resulting spherical conductor model.

    Notes
    -----
    .. versionadded:: 0.9.0

    See Also
    --------
    make_bem_model
    make_bem_solution
    """
    for name in ('r0', 'head_radius'):
        param = locals()[name]
        if isinstance(param, str):
            if param != 'auto':
                raise ValueError('%s, if str, must be "auto" not "%s"'
                                 % (name, param))
    relative_radii = np.array(relative_radii, float).ravel()
    sigmas = np.array(sigmas, float).ravel()
    if len(relative_radii) != len(sigmas):
        raise ValueError('relative_radii length (%s) must match that of '
                         'sigmas (%s)' % (len(relative_radii),
                                          len(sigmas)))
    if len(sigmas) <= 1 and head_radius is not None:
        raise ValueError('at least 2 sigmas must be supplied if '
                         'head_radius is not None, got %s' % (len(sigmas),))
    if (isinstance(r0, str) and r0 == 'auto') or \
       (isinstance(head_radius, str) and head_radius == 'auto'):
        if info is None:
            raise ValueError('Info must not be None for auto mode')
        head_radius_fit, r0_fit = fit_sphere_to_headshape(info, units='m')[:2]
        if isinstance(r0, str):
            r0 = r0_fit
        if isinstance(head_radius, str):
            head_radius = head_radius_fit
    sphere = ConductorModel(is_sphere=True, r0=np.array(r0),
                            coord_frame=FIFF.FIFFV_COORD_HEAD)
    sphere['layers'] = list()
    if head_radius is not None:
        # Eventually these could be configurable...
        relative_radii = np.array(relative_radii, float)
        sigmas = np.array(sigmas, float)
        order = np.argsort(relative_radii)
        relative_radii = relative_radii[order]
        sigmas = sigmas[order]
        for rel_rad, sig in zip(relative_radii, sigmas):
            # sort layers by (relative) radius, and scale radii
            layer = dict(rad=rel_rad, sigma=sig)
            layer['rel_rad'] = layer['rad'] = rel_rad
            sphere['layers'].append(layer)

        # scale the radii
        R = sphere['layers'][-1]['rad']
        rR = sphere['layers'][-1]['rel_rad']
        for layer in sphere['layers']:
            layer['rad'] /= R
            layer['rel_rad'] /= rR

        #
        # Setup the EEG sphere model calculations
        #

        # Scale the relative radii
        for k in range(len(relative_radii)):
            sphere['layers'][k]['rad'] = (head_radius *
                                          sphere['layers'][k]['rel_rad'])
        rv = _fwd_eeg_fit_berg_scherg(sphere, 200, 3)
        logger.info('\nEquiv. model fitting -> RV = %g %%' % (100 * rv))
        for k in range(3):
            logger.info('mu%d = %g    lambda%d = %g'
                        % (k + 1, sphere['mu'][k], k + 1,
                           sphere['layers'][-1]['sigma'] *
                           sphere['lambda'][k]))
        logger.info('Set up EEG sphere model with scalp radius %7.1f mm\n'
                    % (1000 * head_radius,))
    return sphere


# #############################################################################
# Sphere fitting

@verbose
def fit_sphere_to_headshape(info, dig_kinds='auto', units='m', verbose=None):
    """Fit a sphere to the headshape points to determine head center.

    Parameters
    ----------
    info : instance of Info
        Measurement info.
    dig_kinds : list of str | str
        Kind of digitization points to use in the fitting. These can be any
        combination of ('cardinal', 'hpi', 'eeg', 'extra'). Can also
        be 'auto' (default), which will use only the 'extra' points if
        enough (more than 10) are available, and if not, uses 'extra' and
        'eeg' points.
    units : str
        Can be "m" (default) or "mm".

        .. versionadded:: 0.12

    %(verbose)s

    Returns
    -------
    radius : float
        Sphere radius.
    origin_head: ndarray, shape (3,)
        Head center in head coordinates.
    origin_device: ndarray, shape (3,)
        Head center in device coordinates.

    Notes
    -----
    This function excludes any points that are low and frontal
    (``z < 0 and y > 0``) to improve the fit.
    """
    if not isinstance(units, str) or units not in ('m', 'mm'):
        raise ValueError('units must be a "m" or "mm"')
    radius, origin_head, origin_device = _fit_sphere_to_headshape(
        info, dig_kinds)
    if units == 'mm':
        radius *= 1e3
        origin_head *= 1e3
        origin_device *= 1e3
    return radius, origin_head, origin_device


@verbose
def get_fitting_dig(info, dig_kinds='auto', verbose=None):
    """Get digitization points suitable for sphere fitting.

    Parameters
    ----------
    info : instance of Info
        The measurement info.
    dig_kinds : list of str | str
        Kind of digitization points to use in the fitting. These can be any
        combination of ('cardinal', 'hpi', 'eeg', 'extra'). Can also
        be 'auto' (default), which will use only the 'extra' points if
        enough (more than 10) are available, and if not, uses 'extra' and
        'eeg' points.
    %(verbose)s

    Returns
    -------
    dig : array, shape (n_pts, 3)
        The digitization points (in head coordinates) to use for fitting.

    Notes
    -----
    This will exclude digitization locations that have ``z < 0 and y > 0``,
    i.e. points on the nose and below the nose on the face.

    .. versionadded:: 0.14
    """
    _validate_type(info, "info")
    if info['dig'] is None:
        raise RuntimeError('Cannot fit headshape without digitization '
                           ', info["dig"] is None')
    if isinstance(dig_kinds, str):
        if dig_kinds == 'auto':
            # try "extra" first
            try:
                return get_fitting_dig(info, 'extra')
            except ValueError:
                pass
            return get_fitting_dig(info, ('extra', 'eeg'))
        else:
            dig_kinds = (dig_kinds,)
    # convert string args to ints (first make dig_kinds mutable in case tuple)
    dig_kinds = list(dig_kinds)
    for di, d in enumerate(dig_kinds):
        dig_kinds[di] = _dig_kind_dict.get(d, d)
        if dig_kinds[di] not in _dig_kind_ints:
            raise ValueError('dig_kinds[#%d] (%s) must be one of %s'
                             % (di, d, sorted(list(_dig_kind_dict.keys()))))

    # get head digization points of the specified kind(s)
    hsp = [p['r'] for p in info['dig'] if p['kind'] in dig_kinds]
    if any(p['coord_frame'] != FIFF.FIFFV_COORD_HEAD for p in info['dig']):
        raise RuntimeError('Digitization points not in head coordinates, '
                           'contact mne-python developers')

    # exclude some frontal points (nose etc.)
    hsp = np.array([p for p in hsp if not (p[2] < -1e-6 and p[1] > 1e-6)])

    if len(hsp) <= 10:
        kinds_str = ', '.join(['"%s"' % _dig_kind_rev[d]
                               for d in sorted(dig_kinds)])
        msg = ('Only %s head digitization points of the specified kind%s (%s,)'
               % (len(hsp), _pl(dig_kinds), kinds_str))
        if len(hsp) < 4:
            raise ValueError(msg + ', at least 4 required')
        else:
            warn(msg + ', fitting may be inaccurate')
    return hsp


@verbose
def _fit_sphere_to_headshape(info, dig_kinds, verbose=None):
    """Fit a sphere to the given head shape."""
    hsp = get_fitting_dig(info, dig_kinds)
    radius, origin_head = _fit_sphere(np.array(hsp), disp=False)
    # compute origin in device coordinates
    head_to_dev = _ensure_trans(info['dev_head_t'], 'head', 'meg')
    origin_device = apply_trans(head_to_dev, origin_head)
    logger.info('Fitted sphere radius:'.ljust(30) + '%0.1f mm'
                % (radius * 1e3,))
    # 99th percentile on Wikipedia for Giabella to back of head is 21.7cm,
    # i.e. 108mm "radius", so let's go with 110mm
    # en.wikipedia.org/wiki/Human_head#/media/File:HeadAnthropometry.JPG
    if radius > 0.110:
        warn('Estimated head size (%0.1f mm) exceeded 99th '
             'percentile for adult head size' % (1e3 * radius,))
    # > 2 cm away from head center in X or Y is strange
    if np.linalg.norm(origin_head[:2]) > 0.02:
        warn('(X, Y) fit (%0.1f, %0.1f) more than 20 mm from '
             'head frame origin' % tuple(1e3 * origin_head[:2]))
    logger.info('Origin head coordinates:'.ljust(30) +
                '%0.1f %0.1f %0.1f mm' % tuple(1e3 * origin_head))
    logger.info('Origin device coordinates:'.ljust(30) +
                '%0.1f %0.1f %0.1f mm' % tuple(1e3 * origin_device))
    return radius, origin_head, origin_device


def _fit_sphere(points, disp='auto'):
    """Fit a sphere to an arbitrary set of points."""
    from scipy.optimize import fmin_cobyla
    if isinstance(disp, str) and disp == 'auto':
        disp = True if logger.level <= 20 else False
    # initial guess for center and radius
    radii = (np.max(points, axis=1) - np.min(points, axis=1)) / 2.
    radius_init = radii.mean()
    center_init = np.median(points, axis=0)

    # optimization
    x0 = np.concatenate([center_init, [radius_init]])

    def cost_fun(center_rad):
        d = np.linalg.norm(points - center_rad[:3], axis=1) - center_rad[3]
        d *= d
        return d.sum()

    def constraint(center_rad):
        return center_rad[3]  # radius must be >= 0

    x_opt = fmin_cobyla(cost_fun, x0, constraint, rhobeg=radius_init,
                        rhoend=radius_init * 1e-6, disp=disp)

    origin = x_opt[:3]
    radius = x_opt[3]
    return radius, origin


def _check_origin(origin, info, coord_frame='head', disp=False):
    """Check or auto-determine the origin."""
    if isinstance(origin, str):
        if origin != 'auto':
            raise ValueError('origin must be a numerical array, or "auto", '
                             'not %s' % (origin,))
        if coord_frame == 'head':
            R, origin = fit_sphere_to_headshape(info, verbose=False,
                                                units='m')[:2]
            logger.info('    Automatic origin fit: head of radius %0.1f mm'
                        % (R * 1000.,))
            del R
        else:
            origin = (0., 0., 0.)
    origin = np.array(origin, float)
    if origin.shape != (3,):
        raise ValueError('origin must be a 3-element array')
    if disp:
        origin_str = ', '.join(['%0.1f' % (o * 1000) for o in origin])
        msg = ('    Using origin %s mm in the %s frame'
               % (origin_str, coord_frame))
        if coord_frame == 'meg' and info['dev_head_t'] is not None:
            o_dev = apply_trans(info['dev_head_t'], origin)
            origin_str = ', '.join('%0.1f' % (o * 1000,) for o in o_dev)
            msg += ' (%s mm in the head frame)' % (origin_str,)
        logger.info(msg)
    return origin


# ############################################################################
# Create BEM surfaces

@verbose
def make_watershed_bem(subject, subjects_dir=None, overwrite=False,
                       volume='T1', atlas=False, gcaatlas=False, preflood=None,
                       show=False, copy=False, verbose=None):
    """Create BEM surfaces using the FreeSurfer watershed algorithm.

    Parameters
    ----------
    subject : str
        Subject name (required)
    subjects_dir : str
        Directory containing subjects data. If None use
        the Freesurfer SUBJECTS_DIR environment variable.
    overwrite : bool
        Write over existing files
    volume : str
        Defaults to T1
    atlas : bool
        Specify the --atlas option for mri_watershed
    gcaatlas : bool
        Use the subcortical atlas
    preflood : int
        Change the preflood height
    show : bool
        Show surfaces to visually inspect all three BEM surfaces (recommended).

        .. versionadded:: 0.12

    copy : bool
        If True (default False), use copies instead of symlinks for surfaces
        (if they do not already exist).

        .. versionadded:: 0.18
    %(verbose)s

    Notes
    -----
    .. versionadded:: 0.10
    """
    from .viz.misc import plot_bem
    env, mri_dir = _prepare_env(subject, subjects_dir,
                                requires_freesurfer=True)[:2]
    tempdir = _TempDir()  # fsl and Freesurfer create some random junk in CWD
    run_subprocess_env = partial(run_subprocess, env=env,
                                 cwd=tempdir)

    subjects_dir = env['SUBJECTS_DIR']
    subject_dir = op.join(subjects_dir, subject)
    mri_dir = op.join(subject_dir, 'mri')
    T1_dir = op.join(mri_dir, volume)
    T1_mgz = op.join(mri_dir, volume + '.mgz')
    bem_dir = op.join(subject_dir, 'bem')
    ws_dir = op.join(subject_dir, 'bem', 'watershed')
    if not op.isdir(bem_dir):
        os.makedirs(bem_dir)
    if not op.isdir(T1_dir) and not op.isfile(T1_mgz):
        raise RuntimeError('Could not find the MRI data')
    if op.isdir(ws_dir):
        if not overwrite:
            raise RuntimeError('%s already exists. Use the --overwrite option'
                               ' to recreate it.' % ws_dir)
        else:
            shutil.rmtree(ws_dir)
    # put together the command
    cmd = ['mri_watershed']
    if preflood:
        cmd += ["-h", "%s" % int(preflood)]

    if gcaatlas:
        cmd += ['-atlas', '-T1', '-brain_atlas', env['FREESURFER_HOME'] +
                '/average/RB_all_withskull_2007-08-08.gca',
                subject_dir + '/mri/transforms/talairach_with_skull.lta']
    elif atlas:
        cmd += ['-atlas']
    if op.exists(T1_mgz):
        cmd += ['-useSRAS', '-surf', op.join(ws_dir, subject), T1_mgz,
                op.join(ws_dir, 'ws')]
    else:
        cmd += ['-useSRAS', '-surf', op.join(ws_dir, subject), T1_dir,
                op.join(ws_dir, 'ws')]
    # report and run
    logger.info('\nRunning mri_watershed for BEM segmentation with the '
                'following parameters:\n\n'
                'SUBJECTS_DIR = %s\n'
                'SUBJECT = %s\n'
                'Results dir = %s\n' % (subjects_dir, subject, ws_dir))
    os.makedirs(op.join(ws_dir, 'ws'))
    run_subprocess_env(cmd)

    if op.isfile(T1_mgz):
        new_info = _extract_volume_info(T1_mgz)
        if new_info is None:
            warn('nibabel is required to replace the volume info. Volume info'
                 'not updated in the written surface.')
            new_info = dict()
        surfs = ['brain', 'inner_skull', 'outer_skull', 'outer_skin']
        for s in surfs:
            surf_ws_out = op.join(ws_dir, '%s_%s_surface' % (subject, s))

            rr, tris, volume_info = read_surface(surf_ws_out,
                                                 read_metadata=True)
            volume_info.update(new_info)  # replace volume info, 'head' stays

            write_surface(s, rr, tris, volume_info=volume_info)
            # Create symbolic links
            surf_out = op.join(bem_dir, '%s.surf' % s)
            if not overwrite and op.exists(surf_out):
                skip_symlink = True
            else:
                if op.exists(surf_out):
                    os.remove(surf_out)
                _symlink(surf_ws_out, surf_out, copy)
                skip_symlink = False

        if skip_symlink:
            logger.info("Unable to create all symbolic links to .surf files "
                        "in bem folder. Use --overwrite option to recreate "
                        "them.")
            dest = op.join(bem_dir, 'watershed')
        else:
            logger.info("Symbolic links to .surf files created in bem folder")
            dest = bem_dir

    logger.info("\nThank you for waiting.\nThe BEM triangulations for this "
                "subject are now available at:\n%s." % dest)

    # Write a head file for coregistration
    fname_head = op.join(bem_dir, subject + '-head.fif')
    if op.isfile(fname_head):
        os.remove(fname_head)

    surf = _surfaces_to_bem([op.join(ws_dir, subject + '_outer_skin_surface')],
                            [FIFF.FIFFV_BEM_SURF_ID_HEAD], sigmas=[1])
    write_bem_surfaces(fname_head, surf)

    # Show computed BEM surfaces
    if show:
        plot_bem(subject=subject, subjects_dir=subjects_dir,
                 orientation='coronal', slices=None, show=True)

    logger.info('Created %s\n\nComplete.' % (fname_head,))


def _extract_volume_info(mgz, raise_error=True):
    """Extract volume info from a mgz file."""
    try:
        import nibabel as nib
    except ImportError:
        return  # warning raised elsewhere
    header = nib.load(mgz).header
    vol_info = dict()
    version = header['version']
    if version == 1:
        version = '%s  # volume info valid' % version
    else:
        raise ValueError('Volume info invalid.')
    vol_info['valid'] = version
    vol_info['filename'] = mgz
    vol_info['volume'] = header['dims'][:3]
    vol_info['voxelsize'] = header['delta']
    vol_info['xras'], vol_info['yras'], vol_info['zras'] = header['Mdc'].T
    vol_info['cras'] = header['Pxyz_c']
    return vol_info


# ############################################################################
# Read

@verbose
def read_bem_surfaces(fname, patch_stats=False, s_id=None, verbose=None):
    """Read the BEM surfaces from a FIF file.

    Parameters
    ----------
    fname : string
        The name of the file containing the surfaces.
    patch_stats : bool, optional (default False)
        Calculate and add cortical patch statistics to the surfaces.
    s_id : int | None
        If int, only read and return the surface with the given s_id.
        An error will be raised if it doesn't exist. If None, all
        surfaces are read and returned.
    %(verbose)s

    Returns
    -------
    surf: list | dict
        A list of dictionaries that each contain a surface. If s_id
        is not None, only the requested surface will be returned.

    See Also
    --------
    write_bem_surfaces, write_bem_solution, make_bem_model
    """
    # Default coordinate frame
    coord_frame = FIFF.FIFFV_COORD_MRI
    # Open the file, create directory
    f, tree, _ = fiff_open(fname)
    with f as fid:
        # Find BEM
        bem = dir_tree_find(tree, FIFF.FIFFB_BEM)
        if bem is None or len(bem) == 0:
            raise ValueError('BEM data not found')

        bem = bem[0]
        # Locate all surfaces
        bemsurf = dir_tree_find(bem, FIFF.FIFFB_BEM_SURF)
        if bemsurf is None:
            raise ValueError('BEM surface data not found')

        logger.info('    %d BEM surfaces found' % len(bemsurf))
        # Coordinate frame possibly at the top level
        tag = find_tag(fid, bem, FIFF.FIFF_BEM_COORD_FRAME)
        if tag is not None:
            coord_frame = tag.data
        # Read all surfaces
        if s_id is not None:
            surf = [_read_bem_surface(fid, bsurf, coord_frame, s_id)
                    for bsurf in bemsurf]
            surf = [s for s in surf if s is not None]
            if not len(surf) == 1:
                raise ValueError('surface with id %d not found' % s_id)
        else:
            surf = list()
            for bsurf in bemsurf:
                logger.info('    Reading a surface...')
                this = _read_bem_surface(fid, bsurf, coord_frame)
                surf.append(this)
                logger.info('[done]')
            logger.info('    %d BEM surfaces read' % len(surf))
        for this in surf:
            if patch_stats or this['nn'] is None:
                _check_complete_surface(this)
    return surf[0] if s_id is not None else surf


def _read_bem_surface(fid, this, def_coord_frame, s_id=None):
    """Read one bem surface."""
    # fid should be open as a context manager here
    res = dict()
    # Read all the interesting stuff
    tag = find_tag(fid, this, FIFF.FIFF_BEM_SURF_ID)

    if tag is None:
        res['id'] = FIFF.FIFFV_BEM_SURF_ID_UNKNOWN
    else:
        res['id'] = int(tag.data)

    if s_id is not None and res['id'] != s_id:
        return None

    tag = find_tag(fid, this, FIFF.FIFF_BEM_SIGMA)
    res['sigma'] = 1.0 if tag is None else float(tag.data)

    tag = find_tag(fid, this, FIFF.FIFF_BEM_SURF_NNODE)
    if tag is None:
        raise ValueError('Number of vertices not found')

    res['np'] = int(tag.data)

    tag = find_tag(fid, this, FIFF.FIFF_BEM_SURF_NTRI)
    if tag is None:
        raise ValueError('Number of triangles not found')
    res['ntri'] = int(tag.data)

    tag = find_tag(fid, this, FIFF.FIFF_MNE_COORD_FRAME)
    if tag is None:
        tag = find_tag(fid, this, FIFF.FIFF_BEM_COORD_FRAME)
        if tag is None:
            res['coord_frame'] = def_coord_frame
        else:
            res['coord_frame'] = tag.data
    else:
        res['coord_frame'] = tag.data

    # Vertices, normals, and triangles
    tag = find_tag(fid, this, FIFF.FIFF_BEM_SURF_NODES)
    if tag is None:
        raise ValueError('Vertex data not found')

    res['rr'] = tag.data.astype(np.float)  # XXX : double because of mayavi bug
    if res['rr'].shape[0] != res['np']:
        raise ValueError('Vertex information is incorrect')

    tag = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_NORMALS)
    if tag is None:
        tag = find_tag(fid, this, FIFF.FIFF_BEM_SURF_NORMALS)
    if tag is None:
        res['nn'] = None
    else:
        res['nn'] = tag.data.copy()
        if res['nn'].shape[0] != res['np']:
            raise ValueError('Vertex normal information is incorrect')

    tag = find_tag(fid, this, FIFF.FIFF_BEM_SURF_TRIANGLES)
    if tag is None:
        raise ValueError('Triangulation not found')

    res['tris'] = tag.data - 1  # index start at 0 in Python
    if res['tris'].shape[0] != res['ntri']:
        raise ValueError('Triangulation information is incorrect')

    return res


@verbose
def read_bem_solution(fname, verbose=None):
    """Read the BEM solution from a file.

    Parameters
    ----------
    fname : string
        The file containing the BEM solution.
    %(verbose)s

    Returns
    -------
    bem : instance of ConductorModel
        The BEM solution.

    See Also
    --------
    read_bem_surfaces
    write_bem_surfaces
    make_bem_solution
    write_bem_solution
    """
    # mirrors fwd_bem_load_surfaces from fwd_bem_model.c
    logger.info('Loading surfaces...')
    bem_surfs = read_bem_surfaces(fname, patch_stats=True, verbose=False)
    if len(bem_surfs) == 3:
        logger.info('Three-layer model surfaces loaded.')
        needed = np.array([FIFF.FIFFV_BEM_SURF_ID_HEAD,
                           FIFF.FIFFV_BEM_SURF_ID_SKULL,
                           FIFF.FIFFV_BEM_SURF_ID_BRAIN])
        if not all(x['id'] in needed for x in bem_surfs):
            raise RuntimeError('Could not find necessary BEM surfaces')
        # reorder surfaces as necessary (shouldn't need to?)
        reorder = [None] * 3
        for x in bem_surfs:
            reorder[np.where(x['id'] == needed)[0][0]] = x
        bem_surfs = reorder
    elif len(bem_surfs) == 1:
        if not bem_surfs[0]['id'] == FIFF.FIFFV_BEM_SURF_ID_BRAIN:
            raise RuntimeError('BEM Surfaces not found')
        logger.info('Homogeneous model surface loaded.')

    # convert from surfaces to solution
    bem = ConductorModel(is_sphere=False, surfs=bem_surfs)
    logger.info('\nLoading the solution matrix...\n')
    f, tree, _ = fiff_open(fname)
    with f as fid:
        # Find the BEM data
        nodes = dir_tree_find(tree, FIFF.FIFFB_BEM)
        if len(nodes) == 0:
            raise RuntimeError('No BEM data in %s' % fname)
        bem_node = nodes[0]

        # Approximation method
        tag = find_tag(f, bem_node, FIFF.FIFF_BEM_APPROX)
        if tag is None:
            raise RuntimeError('No BEM solution found in %s' % fname)
        method = tag.data[0]
        if method not in (FIFF.FIFFV_BEM_APPROX_CONST,
                          FIFF.FIFFV_BEM_APPROX_LINEAR):
            raise RuntimeError('Cannot handle BEM approximation method : %d'
                               % method)

        tag = find_tag(fid, bem_node, FIFF.FIFF_BEM_POT_SOLUTION)
        dims = tag.data.shape
        if len(dims) != 2:
            raise RuntimeError('Expected a two-dimensional solution matrix '
                               'instead of a %d dimensional one' % dims[0])

        dim = 0
        for surf in bem['surfs']:
            if method == FIFF.FIFFV_BEM_APPROX_LINEAR:
                dim += surf['np']
            else:  # method == FIFF.FIFFV_BEM_APPROX_CONST
                dim += surf['ntri']

        if dims[0] != dim or dims[1] != dim:
            raise RuntimeError('Expected a %d x %d solution matrix instead of '
                               'a %d x %d one' % (dim, dim, dims[1], dims[0]))
        sol = tag.data
        nsol = dims[0]

    bem['solution'] = sol
    bem['nsol'] = nsol
    bem['bem_method'] = method

    # Gamma factors and multipliers
    _add_gamma_multipliers(bem)
    kind = {
        FIFF.FIFFV_BEM_APPROX_CONST: 'constant collocation',
        FIFF.FIFFV_BEM_APPROX_LINEAR: 'linear_collocation',
    }[bem['bem_method']]
    logger.info('Loaded %s BEM solution from %s', kind, fname)
    return bem


def _add_gamma_multipliers(bem):
    """Add gamma and multipliers in-place."""
    bem['sigma'] = np.array([surf['sigma'] for surf in bem['surfs']])
    # Dirty trick for the zero conductivity outside
    sigma = np.r_[0.0, bem['sigma']]
    bem['source_mult'] = 2.0 / (sigma[1:] + sigma[:-1])
    bem['field_mult'] = sigma[1:] - sigma[:-1]
    # make sure subsequent "zip"s work correctly
    assert len(bem['surfs']) == len(bem['field_mult'])
    bem['gamma'] = ((sigma[1:] - sigma[:-1])[np.newaxis, :] /
                    (sigma[1:] + sigma[:-1])[:, np.newaxis])


_surf_dict = {'inner_skull': FIFF.FIFFV_BEM_SURF_ID_BRAIN,
              'outer_skull': FIFF.FIFFV_BEM_SURF_ID_SKULL,
              'head': FIFF.FIFFV_BEM_SURF_ID_HEAD}


def _bem_find_surface(bem, id_):
    """Find surface from already-loaded BEM."""
    if isinstance(id_, str):
        name = id_
        id_ = _surf_dict[id_]
    else:
        name = _surf_name[id_]
    idx = np.where(np.array([s['id'] for s in bem['surfs']]) == id_)[0]
    if len(idx) != 1:
        raise RuntimeError('BEM model does not have the %s triangulation'
                           % name.replace('_', ' '))
    return bem['surfs'][idx[0]]


# ############################################################################
# Write

def write_bem_surfaces(fname, surfs):
    """Write BEM surfaces to a fiff file.

    Parameters
    ----------
    fname : str
        Filename to write.
    surfs : dict | list of dict
        The surfaces, or a single surface.
    """
    if isinstance(surfs, dict):
        surfs = [surfs]
    with start_file(fname) as fid:
        start_block(fid, FIFF.FIFFB_BEM)
        write_int(fid, FIFF.FIFF_BEM_COORD_FRAME, surfs[0]['coord_frame'])
        _write_bem_surfaces_block(fid, surfs)
        end_block(fid, FIFF.FIFFB_BEM)
        end_file(fid)


def _write_bem_surfaces_block(fid, surfs):
    """Write bem surfaces to open file handle."""
    for surf in surfs:
        start_block(fid, FIFF.FIFFB_BEM_SURF)
        write_float(fid, FIFF.FIFF_BEM_SIGMA, surf['sigma'])
        write_int(fid, FIFF.FIFF_BEM_SURF_ID, surf['id'])
        write_int(fid, FIFF.FIFF_MNE_COORD_FRAME, surf['coord_frame'])
        write_int(fid, FIFF.FIFF_BEM_SURF_NNODE, surf['np'])
        write_int(fid, FIFF.FIFF_BEM_SURF_NTRI, surf['ntri'])
        write_float_matrix(fid, FIFF.FIFF_BEM_SURF_NODES, surf['rr'])
        # index start at 0 in Python
        write_int_matrix(fid, FIFF.FIFF_BEM_SURF_TRIANGLES,
                         surf['tris'] + 1)
        if 'nn' in surf and surf['nn'] is not None and len(surf['nn']) > 0:
            write_float_matrix(fid, FIFF.FIFF_BEM_SURF_NORMALS, surf['nn'])
        end_block(fid, FIFF.FIFFB_BEM_SURF)


def write_bem_solution(fname, bem):
    """Write a BEM model with solution.

    Parameters
    ----------
    fname : str
        The filename to use.
    bem : instance of ConductorModel
        The BEM model with solution to save.

    See Also
    --------
    read_bem_solution
    """
    _check_bem_size(bem['surfs'])
    with start_file(fname) as fid:
        start_block(fid, FIFF.FIFFB_BEM)
        # Coordinate frame (mainly for backward compatibility)
        write_int(fid, FIFF.FIFF_BEM_COORD_FRAME,
                  bem['surfs'][0]['coord_frame'])
        # Surfaces
        _write_bem_surfaces_block(fid, bem['surfs'])
        # The potential solution
        if 'solution' in bem:
            if bem['bem_method'] != FWD.BEM_LINEAR_COLL:
                raise RuntimeError('Only linear collocation supported')
            write_int(fid, FIFF.FIFF_BEM_APPROX, FIFF.FIFFV_BEM_APPROX_LINEAR)
            write_float_matrix(fid, FIFF.FIFF_BEM_POT_SOLUTION,
                               bem['solution'])
        end_block(fid, FIFF.FIFFB_BEM)
        end_file(fid)


# #############################################################################
# Create 3-Layers BEM model from Flash MRI images

def _prepare_env(subject, subjects_dir, requires_freesurfer):
    """Prepare an env object for subprocess calls."""
    env = os.environ.copy()
    if requires_freesurfer and not os.environ.get('FREESURFER_HOME'):
        raise RuntimeError('I cannot find freesurfer. The FREESURFER_HOME '
                           'environment variable is not set.')

    _validate_type(subject, "str")

    subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
    if not op.isdir(subjects_dir):
        raise RuntimeError('Could not find the MRI data directory "%s"'
                           % subjects_dir)
    subject_dir = op.join(subjects_dir, subject)
    if not op.isdir(subject_dir):
        raise RuntimeError('Could not find the subject data directory "%s"'
                           % (subject_dir,))
    env['SUBJECT'] = subject
    env['SUBJECTS_DIR'] = subjects_dir
    mri_dir = op.join(subject_dir, 'mri')
    bem_dir = op.join(subject_dir, 'bem')
    return env, mri_dir, bem_dir


@verbose
def convert_flash_mris(subject, flash30=True, convert=True, unwarp=False,
                       subjects_dir=None, verbose=None):
    """Convert DICOM files for use with make_flash_bem.

    Parameters
    ----------
    subject : str
        Subject name.
    flash30 : bool
        Use 30-degree flip angle data.
    convert : bool
        Assume that the Flash MRI images have already been converted
        to mgz files.
    unwarp : bool
        Run grad_unwarp with -unwarp option on each of the converted
        data sets. It requires FreeSurfer's MATLAB toolbox to be properly
        installed.
    subjects_dir : string, or None
        Path to SUBJECTS_DIR if it is not set in the environment.
    %(verbose)s

    Notes
    -----
    Before running this script do the following:
    (unless convert=False is specified)

        1. Copy all of your FLASH images in a single directory <source> and
           create a directory <dest> to hold the output of mne_organize_dicom
        2. cd to <dest> and run
           $ mne_organize_dicom <source>
           to create an appropriate directory structure
        3. Create symbolic links to make flash05 and flash30 point to the
           appropriate series:
           $ ln -s <FLASH 5 series dir> flash05
           $ ln -s <FLASH 30 series dir> flash30
           Some partition formats (e.g. FAT32) do not support symbolic links.
           In this case, copy the file to the appropriate series:
           $ cp <FLASH 5 series dir> flash05
           $ cp <FLASH 30 series dir> flash30
        4. cd to the directory where flash05 and flash30 links are
        5. Set SUBJECTS_DIR and SUBJECT environment variables appropriately
        6. Run this script

    This function assumes that the Freesurfer segmentation of the subject
    has been completed. In particular, the T1.mgz and brain.mgz MRI volumes
    should be, as usual, in the subject's mri directory.
    """
    env, mri_dir = _prepare_env(subject, subjects_dir,
                                requires_freesurfer=True)[:2]
    tempdir = _TempDir()  # fsl and Freesurfer create some random junk in CWD
    run_subprocess_env = partial(run_subprocess, env=env,
                                 cwd=tempdir)
    # Step 1a : Data conversion to mgz format
    if not op.exists(op.join(mri_dir, 'flash', 'parameter_maps')):
        os.makedirs(op.join(mri_dir, 'flash', 'parameter_maps'))
    echos_done = 0
    if convert:
        logger.info("\n---- Converting Flash images ----")
        echos = ['001', '002', '003', '004', '005', '006', '007', '008']
        if flash30:
            flashes = ['05', '30']
        else:
            flashes = ['05']
        #
        missing = False
        for flash in flashes:
            for echo in echos:
                if not op.isdir(op.join('flash' + flash, echo)):
                    missing = True
        if missing:
            echos = ['002', '003', '004', '005', '006', '007', '008', '009']
            for flash in flashes:
                for echo in echos:
                    if not op.isdir(op.join('flash' + flash, echo)):
                        raise RuntimeError("Directory %s is missing."
                                           % op.join('flash' + flash, echo))
        #
        for flash in flashes:
            for echo in echos:
                if not op.isdir(op.join('flash' + flash, echo)):
                    raise RuntimeError("Directory %s is missing."
                                       % op.join('flash' + flash, echo))
                sample_file = glob.glob(op.join('flash' + flash, echo, '*'))[0]
                dest_file = op.join(mri_dir, 'flash',
                                    'mef' + flash + '_' + echo + '.mgz')
                # do not redo if already present
                if op.isfile(dest_file):
                    logger.info("The file %s is already there")
                else:
                    cmd = ['mri_convert', sample_file, dest_file]
                    run_subprocess_env(cmd)
                    echos_done += 1
    # Step 1b : Run grad_unwarp on converted files
    flash_dir = op.join(mri_dir, "flash")
    template = op.join(flash_dir, "mef*.mgz")
    files = glob.glob(template)
    if len(files) == 0:
        raise ValueError('No suitable source files found (%s)' % template)
    if unwarp:
        logger.info("\n---- Unwarp mgz data sets ----")
        for infile in files:
            outfile = infile.replace(".mgz", "u.mgz")
            cmd = ['grad_unwarp', '-i', infile, '-o', outfile, '-unwarp',
                   'true']
            run_subprocess_env(cmd)
    # Clear parameter maps if some of the data were reconverted
    pm_dir = op.join(flash_dir, 'parameter_maps')
    if echos_done > 0 and op.exists(pm_dir):
        shutil.rmtree(pm_dir)
        logger.info("\nParameter maps directory cleared")
    if not op.exists(pm_dir):
        os.makedirs(pm_dir)
    # Step 2 : Create the parameter maps
    if flash30:
        logger.info("\n---- Creating the parameter maps ----")
        if unwarp:
            files = glob.glob(op.join(flash_dir, "mef05*u.mgz"))
        if len(os.listdir(pm_dir)) == 0:
            cmd = (['mri_ms_fitparms'] +
                   files +
                   [op.join(flash_dir, 'parameter_maps')])
            run_subprocess_env(cmd)
        else:
            logger.info("Parameter maps were already computed")
        # Step 3 : Synthesize the flash 5 images
        logger.info("\n---- Synthesizing flash 5 images ----")
        if not op.exists(op.join(pm_dir, 'flash5.mgz')):
            cmd = ['mri_synthesize', '20', '5', '5',
                   op.join(pm_dir, 'T1.mgz'),
                   op.join(pm_dir, 'PD.mgz'),
                   op.join(pm_dir, 'flash5.mgz')
                   ]
            run_subprocess_env(cmd)
            os.remove(op.join(pm_dir, 'flash5_reg.mgz'))
        else:
            logger.info("Synthesized flash 5 volume is already there")
    else:
        logger.info("\n---- Averaging flash5 echoes ----")
        template = op.join(flash_dir,
                           "mef05*u.mgz" if unwarp else "mef05*.mgz")
        files = glob.glob(template)
        if len(files) == 0:
            raise ValueError('No suitable source files found (%s)' % template)
        cmd = (['mri_average', '-noconform'] +
               files +
               [op.join(pm_dir, 'flash5.mgz')])
        run_subprocess_env(cmd)
        if op.exists(op.join(pm_dir, 'flash5_reg.mgz')):
            os.remove(op.join(pm_dir, 'flash5_reg.mgz'))
    del tempdir  # finally done running subprocesses
    assert op.isfile(op.join(pm_dir, 'flash5.mgz'))


@verbose
def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None,
                   flash_path=None, copy=False, verbose=None):
    """Create 3-Layer BEM model from prepared flash MRI images.

    Parameters
    ----------
    subject : str
        Subject name.
    overwrite : bool
        Write over existing .surf files in bem folder.
    show : bool
        Show surfaces to visually inspect all three BEM surfaces (recommended).
    subjects_dir : string, or None
        Path to SUBJECTS_DIR if it is not set in the environment.
    flash_path : str | None
        Path to the flash images. If None (default), mri/flash/parameter_maps
        within the subject reconstruction is used.

        .. versionadded:: 0.13.0
    copy : bool
        If True (default False), use copies instead of symlinks for surfaces
        (if they do not already exist).

        .. versionadded:: 0.18
    %(verbose)s

    Notes
    -----
    This program assumes that FreeSurfer is installed and sourced properly.

    This function extracts the BEM surfaces (outer skull, inner skull, and
    outer skin) from multiecho FLASH MRI data with spin angles of 5 and 30
    degrees, in mgz format.

    See Also
    --------
    convert_flash_mris
    """
    from .viz.misc import plot_bem

    is_test = os.environ.get('MNE_SKIP_FS_FLASH_CALL', False)

    env, mri_dir, bem_dir = _prepare_env(subject, subjects_dir,
                                         requires_freesurfer=True)
    tempdir = _TempDir()  # fsl and Freesurfer create some random junk in CWD
    run_subprocess_env = partial(run_subprocess, env=env,
                                 cwd=tempdir)

    if flash_path is None:
        flash_path = op.join(mri_dir, 'flash', 'parameter_maps')
    else:
        flash_path = op.abspath(flash_path)
    subjects_dir = env['SUBJECTS_DIR']

    logger.info('\nProcessing the flash MRI data to produce BEM meshes with '
                'the following parameters:\n'
                'SUBJECTS_DIR = %s\n'
                'SUBJECT = %s\n'
                'Result dir = %s\n' % (subjects_dir, subject,
                                       op.join(bem_dir, 'flash')))
    # Step 4 : Register with MPRAGE
    logger.info("\n---- Registering flash 5 with MPRAGE ----")
    flash5 = op.join(flash_path, 'flash5.mgz')
    flash5_reg = op.join(flash_path, 'flash5_reg.mgz')
    if not op.exists(flash5_reg):
        if op.exists(op.join(mri_dir, 'T1.mgz')):
            ref_volume = op.join(mri_dir, 'T1.mgz')
        else:
            ref_volume = op.join(mri_dir, 'T1')
        cmd = ['fsl_rigid_register', '-r', ref_volume, '-i', flash5,
               '-o', flash5_reg]
        run_subprocess_env(cmd)
    else:
        logger.info("Registered flash 5 image is already there")
    # Step 5a : Convert flash5 into COR
    logger.info("\n---- Converting flash5 volume into COR format ----")
    flash5_dir = op.join(mri_dir, 'flash5')
    shutil.rmtree(flash5_dir, ignore_errors=True)
    os.makedirs(flash5_dir)
    if not is_test:  # CIs don't have freesurfer, skipped when testing.
        cmd = ['mri_convert', flash5_reg, op.join(mri_dir, 'flash5')]
        run_subprocess_env(cmd)
    # Step 5b and c : Convert the mgz volumes into COR
    convert_T1 = False
    T1_dir = op.join(mri_dir, 'T1')
    if not op.isdir(T1_dir) or len(glob.glob(op.join(T1_dir, 'COR*'))) == 0:
        convert_T1 = True
    convert_brain = False
    brain_dir = op.join(mri_dir, 'brain')
    if not op.isdir(brain_dir) or \
            len(glob.glob(op.join(brain_dir, 'COR*'))) == 0:
        convert_brain = True
    logger.info("\n---- Converting T1 volume into COR format ----")
    if convert_T1:
        T1_fname = op.join(mri_dir, 'T1.mgz')
        if not op.isfile(T1_fname):
            raise RuntimeError("Both T1 mgz and T1 COR volumes missing.")
        os.makedirs(T1_dir)
        cmd = ['mri_convert', T1_fname, T1_dir]
        run_subprocess_env(cmd)
    else:
        logger.info("T1 volume is already in COR format")
    logger.info("\n---- Converting brain volume into COR format ----")
    if convert_brain:
        brain_fname = op.join(mri_dir, 'brain.mgz')
        if not op.isfile(brain_fname):
            raise RuntimeError("Both brain mgz and brain COR volumes missing.")
        os.makedirs(brain_dir)
        cmd = ['mri_convert', brain_fname, brain_dir]
        run_subprocess_env(cmd)
    else:
        logger.info("Brain volume is already in COR format")
    # Finally ready to go
    if not is_test:  # CIs don't have freesurfer, skipped when testing.
        logger.info("\n---- Creating the BEM surfaces ----")
        cmd = ['mri_make_bem_surfaces', subject]
        run_subprocess_env(cmd)
    del tempdir  # ran our last subprocess; clean up directory

    logger.info("\n---- Converting the tri files into surf files ----")
    flash_bem_dir = op.join(bem_dir, 'flash')
    if not op.exists(flash_bem_dir):
        os.makedirs(flash_bem_dir)
    surfs = ['inner_skull', 'outer_skull', 'outer_skin']
    for surf in surfs:
        out_fname = op.join(flash_bem_dir, surf + '.tri')
        shutil.move(op.join(bem_dir, surf + '.tri'), out_fname)
        nodes, tris = read_tri(out_fname, swap=True)
        # Do not write volume info here because the tris are already in
        # standard Freesurfer coords
        write_surface(op.splitext(out_fname)[0] + '.surf', nodes, tris)

    # Cleanup section
    logger.info("\n---- Cleaning up ----")
    os.remove(op.join(bem_dir, 'inner_skull_tmp.tri'))
    # os.chdir(mri_dir)
    if convert_T1:
        shutil.rmtree(T1_dir)
        logger.info("Deleted the T1 COR volume")
    if convert_brain:
        shutil.rmtree(brain_dir)
        logger.info("Deleted the brain COR volume")
    shutil.rmtree(flash5_dir)
    logger.info("Deleted the flash5 COR volume")
    # Create symbolic links to the .surf files in the bem folder
    logger.info("\n---- Creating symbolic links ----")
    # os.chdir(bem_dir)
    for surf in surfs:
        surf = op.join(bem_dir, surf + '.surf')
        if not overwrite and op.exists(surf):
            skip_symlink = True
        else:
            if op.exists(surf):
                os.remove(surf)
            _symlink(op.join(flash_bem_dir, op.basename(surf)), surf, copy)
            skip_symlink = False
    if skip_symlink:
        logger.info("Unable to create all symbolic links to .surf files "
                    "in bem folder. Use --overwrite option to recreate them.")
        dest = op.join(bem_dir, 'flash')
    else:
        logger.info("Symbolic links to .surf files created in bem folder")
        dest = bem_dir
    logger.info("\nThank you for waiting.\nThe BEM triangulations for this "
                "subject are now available at:\n%s.\nWe hope the BEM meshes "
                "created will facilitate your MEG and EEG data analyses."
                % dest)
    # Show computed BEM surfaces
    if show:
        plot_bem(subject=subject, subjects_dir=subjects_dir,
                 orientation='coronal', slices=None, show=True)


def _check_bem_size(surfs):
    """Check bem surface sizes."""
    if len(surfs) > 1 and surfs[0]['np'] > 10000:
        warn('The bem surfaces have %s data points. 5120 (ico grade=4) '
             'should be enough. Dense 3-layer bems may not save properly.' %
             surfs[0]['np'])


def _symlink(src, dest, copy=False):
    """Create a relative symlink (or just copy)."""
    if not copy:
        src_link = op.relpath(src, op.dirname(dest))
        try:
            os.symlink(src_link, dest)
        except OSError:
            warn('Could not create symbolic link %s. Check that your '
                 'partition handles symbolic links. The file will be copied '
                 'instead.' % dest)
            copy = True
    if copy:
        shutil.copy(src, dest)