function [headmodel, cfg] = ft_prepare_headmodel(cfg, data)

% FT_PREPARE_HEADMODEL constructs a volume conduction model from the geometry
% of the head. The volume conduction model specifies how currents that are
% generated by sources in the brain, e.g. dipoles, are propagated through the
% tissue and how these result in externally measureable EEG potentials or MEG
% fields.
%
% FieldTrip implements a variety of forward solutions, partially with internal
% code and some of them using external toolboxes or executables. Each of the
% forward solutions requires a set of configuration options which are listed
% below. This function takes care of all the preparatory steps in the
% construction of the volume conduction model and sets it up so that
% subsequent computations are efficient and fast.
%
% Use as
%   headmodel = ft_prepare_headmodel(cfg)               or
%   headmodel = ft_prepare_headmodel(cfg, mesh)         with the output of FT_PREPARE_MESH or FT_READ_HEADSHAPE
%   headmodel = ft_prepare_headmodel(cfg, seg)          with the output of FT_VOLUMESEGMENT
%   headmodel = ft_prepare_headmodel(cfg, elec)         with the output of FT_READ_SENS
%   headmodel = ft_prepare_headmodel(cfg, sourcemodel)  with the output of FT_PREPARE_LEADFIELD
%
% In general the input to this function is a geometrical description of the
% shape of the head and a description of the electrical conductivity. The
% geometrical description can be a set of surface points obtained from
% fT_READ_HEADSHAPE, a surface mesh that was obtained from FT_PREPARE_MESH or
% a segmented anatomical MRI that was obtained from FT_VOLUMESEGMENT.
%
% The cfg argument is a structure that can contain:
%   cfg.method         = string that specifies the forward solution, see below
%   cfg.conductivity   = a number or a vector containing the conductivities of the compartments
%   cfg.tissue         = a string or integer, to be used in combination with a 'seg' for the
%                          second intput. If 'brain', 'skull', and 'scalp' are fields
%                          present in 'seg', then cfg.tissue need not be specified, as
%                          these are defaults, depending on cfg.method. Otherwise,
%                          cfg.tissue should refer to which field(s) of seg should be used.
%
% For EEG the following methods are available:
%   singlesphere       analytical single sphere model
%   concentricspheres  analytical concentric sphere model with up to 4 spheres
%   openmeeg           boundary element method, based on the OpenMEEG software
%   bemcp              boundary element method, based on the implementation from Christophe Phillips
%   dipoli             boundary element method, based on the implementation from Thom Oostendorp
%   hbf                boundary element method, based on the implementation from Matti Stenroos
%   asa                boundary element method, based on the (commercial) ASA software
%   simbio             finite element method, based on the SimBio software
%   duneuro            finite element method, based on the DUNEuro software
%   fns                finite difference method, based on the FNS software
%   infinite           electric dipole in an infinite homogenous medium
%   halfspace          infinite homogenous medium on one side, vacuum on the other
%   besa               finite element leadfield matrix from BESA
%   interpolate        interpolate the precomputed leadfield
%
% For MEG the following methods are available:
%   openmeeg           boundary element method, based on the OpenMEEG software
%   hbf                boundary element method, based on the implementation from Matti Stenroos
%   singlesphere       analytical single sphere model
%   localspheres       local spheres model for MEG, one sphere per channel
%   singleshell        realisically shaped single shell approximation, based on the implementation from Guido Nolte
%   infinite           magnetic dipole in an infinite vacuum
%
% Each specific method has its own specific configuration options which are listed below.
%
% BEMCP, DIPOLI, OPENMEEG
%   cfg.tissue            see above; in combination with 'seg' input
%   cfg.isolatedsource    (optional)
%   cfg.tempdir           (optional)
%   cfg.tempname          (optional)
%
% CONCENTRICSPHERES
%   cfg.tissue            see above; in combination with 'seg' input
%   cfg.order             (optional)
%   cfg.fitind            (optional)
%
% LOCALSPHERES
%   cfg.grad
%   cfg.tissue            see above; in combination with 'seg' input; default options are 'brain' or 'scalp'
%   cfg.feedback          (optional)
%   cfg.radius            (optional)
%   cfg.maxradius         (optional)
%   cfg.baseline          (optional)
%
% SIMBIO
%   cfg.conductivity
%
% DUNEURO
%   cfg.conductivity      An array with the conductivities must be provided. (see above)
%   cfg.grid_filename     Alternatively,  a filename for the grid and a filename for the conductivities can be passed.
%   cfg.tensors_filename  "
%   cfg.duneuro           (optional) Additional settings can be provided for duneuro (see http://www.duneuro.org, and duneuro_defaults.m).
%
% SINGLESHELL
%   cfg.tissue            see above; in combination with 'seg' input; default options are 'brain' or 'scalp'
%   cfg.order             (optional)
%
% SINGLESPHERE
%   cfg.tissue            see above; in combination with 'seg' input; default options are 'brain' or 'scalp'; must be only 1 value
%
% INTERPOLATE
%    cfg.outputfile       (required) string, filename prefix for the output files
%
% BESA
%   cfg.headmodel         (required) string, filename of precomputed FEM leadfield
%   cfg.elec              (required) structure with electrode positions or filename, see FT_READ_SENS
%   cfg.outputfile        (required) string, filename prefix for the output files
%
% FNS
%   cfg.tissue
%   cfg.tissueval
%   cfg.conductivity
%   cfg.elec
%   cfg.grad
%   cfg.transform
%
% HALFSPACE
%   cfg.point
%   cfg.submethod         (optional)
%
% HBF
%   cfg.conductivity      (required) [2 x n_boundaries] array with conductivity values for inside and outside of a boundary
%   cfg.isolatedsource    (optional) set to 1 to apply isolated source approach on innermost boundary
%   cfg.checkmesh         (optional) ['yes | 'no'] check the integrity and ordering of boundaries
%
%
% More details for each of the specific methods can be found in the corresponding
% low-level function which is called FT_HEADMODEL_XXX where XXX is the method
% of choise.
%
% See also FT_PREPARE_MESH, FT_PREPARE_SOURCEMODEL, FT_PREPARE_LEADFIELD,
% FT_HEADMODEL_BEMCP, FT_HEADMODEL_ASA, FT_HEADMODEL_DIPOLI,
% FT_HEADMODEL_SIMBIO, FT_HEADMODEL_FNS, FT_HEADMODEL_HALFSPACE,
% FT_HEADMODEL_INFINITE, FT_HEADMODEL_OPENMEEG, FT_HEADMODEL_SINGLESPHERE,
% FT_HEADMODEL_CONCENTRICSPHERES, FT_HEADMODEL_LOCALSPHERES,
% FT_HEADMODEL_SINGLESHELL, FT_HEADMODEL_INTERPOLATE, FT_HEADMODEL_DUNEURO,
% FT_HEADMODEL_HBF

% Copyright (C) 2011, Cristiano Micheli
% Copyright (C) 2011-2012, Jan-Mathijs Schoffelen, Robert Oostenveld
% Copyright (C) 2013-2018, Robert Oostenveld, Johanna Zumer
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init

ft_preamble provenance data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'deprecated', 'geom');
cfg = ft_checkconfig(cfg, 'renamed', {'geom', 'headshape'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'bem_openmeeg', 'openmeeg'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'bem_dipoli', 'dipoli'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'bem_cp', 'bemcp'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'nolte', 'singleshell'});
cfg = ft_checkconfig(cfg, 'renamed', {'duneuro_settings', 'duneuro'});
cfg = ft_checkconfig(cfg, 'renamed', {'hdmfile', 'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed', {'vol',     'headmodel'});

if isfield(cfg, 'headmodel') && ischar(cfg.headmodel) && ~endsWith(cfg.headmodel, '.lft')
  cfg.method = 'file'; % FIXME this is not documented, note that it does not apply to BESA headmodels
elseif isfield(cfg, 'headmodel') && isstruct(cfg.headmodel)
  cfg.method = 'existing'; % FIXME this is not documented
end

cfg = ft_checkconfig(cfg, 'required', 'method');

% set the general defaults
cfg.headshape       = ft_getopt(cfg, 'headshape');
cfg.conductivity    = ft_getopt(cfg, 'conductivity');

% volume related options
cfg.tissue          = ft_getopt(cfg, 'tissue');
cfg.smooth          = ft_getopt(cfg, 'smooth');
cfg.threshold       = ft_getopt(cfg, 'threshold');

% other options
cfg.isolatedsource  = ft_getopt(cfg, 'isolatedsource');   % used for dipoli and openmeeg
cfg.tempdir         = ft_getopt(cfg, 'tempdir');          % used for dipoli
cfg.tempname        = ft_getopt(cfg, 'tempname');         % used for dipoli
cfg.point           = ft_getopt(cfg, 'point');            % used for halfspace
cfg.submethod       = ft_getopt(cfg, 'submethod');        % used for halfspace
cfg.feedback        = ft_getopt(cfg, 'feedback');
cfg.radius          = ft_getopt(cfg, 'radius');
cfg.maxradius       = ft_getopt(cfg, 'maxradius');
cfg.baseline        = ft_getopt(cfg, 'baseline');
cfg.singlesphere    = ft_getopt(cfg, 'singlesphere');
cfg.grid_filename   = ft_getopt(cfg, 'grid_filename');    % used for duneuro
cfg.tensors_filename= ft_getopt(cfg, 'tensors_filename'); % used for duneuro
cfg.duneuro         = ft_getopt(cfg, 'duneuro');
cfg.tissueval       = ft_getopt(cfg, 'tissueval');        % used for simbio
cfg.transform       = ft_getopt(cfg, 'transform');
cfg.siunits         = ft_getopt(cfg, 'siunits', 'no');    % yes/no, convert the input and continue with SI units
cfg.unit            = ft_getopt(cfg, 'unit');
cfg.smooth          = ft_getopt(cfg, 'smooth');           % used for interpolate
cfg.headmodel       = ft_getopt(cfg, 'headmodel');        % can contain CTF localspheres model
cfg.checkmesh       = ft_getopt(cfg, 'checkmesh');        % used for hbf

% the data can be passed as input arguments or can be read from disk
hasdata = exist('data', 'var');

if hasdata
  % the data should describe the geometrical mesh
  if isfield(data, 'bnd')
    data = data.bnd;
  end
  % check if the input data is valid for this function and ensure that it has the units specified
  data = ft_checkdata(data, 'hasunit', 'yes');
  % replace pnt by pos
  data = fixpos(data);
else
  data = [];
end

% convert to SI units
if istrue(cfg.siunits)
  cfg.unit = 'm';
end

% convert the geometrical data to the desired units for the source model
if ~isempty(cfg.unit)
  if ~isempty(data)
    data = ft_convert_units(data, cfg.unit);
  end
  if isfield(cfg, 'grad') && ~isempty(cfg.grad)
    cfg.grad = ft_convert_units(cfg.grad, cfg.unit);
  end
  if isfield(cfg, 'elec') && ~isempty(cfg.elec)
    cfg.elec = ft_convert_units(cfg.elec, cfg.unit);
  end
end

% the conductivity in the data overrules cfg.conductivity
if hasdata && isfield(data, 'cond')
  cfg.conductivity = data.cond;
end

% boolean variables to manages the different geometrical input data objects
input_mesh  = ft_datatype(data, 'mesh');
input_seg   = ft_datatype(data, 'segmentation');
input_elec  = ft_datatype(data, 'sens');
input_pos   = ~input_mesh && isfield(data, 'pos'); % surface points without triangulation

% the construction of the volume conductor model is performed below
switch cfg.method
  case 'file'
    % read it from file
    headmodel = ft_read_headmodel(cfg.headmodel);
    cfg = rmfield(cfg, 'method'); % FIXME this is not documented

  case 'existing'
    % return an existing one
    headmodel = cfg.headmodel;
    cfg = rmfield(cfg, 'method'); % FIXME this is not documented

  case 'interpolate'
    % the "data" here represents the output of FT_PREPARE_LEADFIELD, i.e. a regular dipole
    % grid with pre-computed leadfields
    sens = ft_fetch_sens(cfg, data);
    headmodel = ft_headmodel_interpolate(cfg.outputfile, sens, data, 'smooth', cfg.smooth);

  case 'besa'
    % cfg.headmodel points to the filename of the FEM solution that was computed in BESA
    % cfg.elec points to the filename of the corresponding electrode specification
    sens = ft_fetch_sens(cfg, data);
    headmodel = ft_headmodel_interpolate(cfg.outputfile, sens, cfg.headmodel, 'smooth', cfg.smooth);

  case 'asa'
    if ~ft_filetype(cfg.headmodel, 'asa_vol')
      ft_error('You must supply a valid cfg.headmodel for use with ASA headmodel')
    end
    headmodel = ft_headmodel_asa(cfg.headmodel);

    case {'bemcp' 'dipoli' 'hbf' 'openmeeg'}
    % the low-level functions all need a mesh
    if isfield(data, 'pos') && isfield(data, 'tri')
      if ~isfield(cfg, 'numvertices') || isempty(cfg.numvertices) || isequal(cfg.numvertices, arrayfun(@(x) size(x.pos, 1), data))
        % copy the input data
        geometry = data;
      else
        % retriangulate the input to the user-specified number of vertices
        tmpcfg.method = 'headshape';
        tmpcfg.headshape = data;
        tmpcfg.numvertices = cfg.numvertices;
        geometry = ft_prepare_mesh(tmpcfg);
      end
    elseif isfield(data, 'transform') && isfield(data, 'dim')
      tmpcfg   = keepfields(cfg, {'numvertices', 'tissue', 'spmversion'});
      geometry = ft_prepare_mesh(tmpcfg, data);
    else
      ft_error('Either a segmented MRI or data with closed triangulated mesh is required as data input for the bemcp, dipoli or openmeeg method');
    end

    if strcmp(cfg.method, 'bemcp')
      headmodel = ft_headmodel_bemcp(geometry, 'conductivity', cfg.conductivity);
      if any(isnan(headmodel.mat(:)))
        % HACK add a little bit of noise, with the NatMEG tutorial data, I discovered that this prevents the warning
        % Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.
        geometry(1).pos = geometry(1).pos + randn(size(geometry(1).pos))*ft_scalingfactor('um', geometry(1).unit);
        geometry(2).pos = geometry(2).pos + randn(size(geometry(2).pos))*ft_scalingfactor('um', geometry(2).unit);
        geometry(3).pos = geometry(3).pos + randn(size(geometry(3).pos))*ft_scalingfactor('um', geometry(3).unit);
        ft_warning('NaN detected, trying once more with slightly different vertex positions');
        headmodel = ft_headmodel_bemcp(geometry, 'conductivity', cfg.conductivity);
      end
    elseif strcmp(cfg.method, 'dipoli')
      headmodel = ft_headmodel_dipoli(geometry, 'conductivity', cfg.conductivity, 'isolatedsource', cfg.isolatedsource, 'tempdir', cfg.tempdir, 'tempname', cfg.tempname);
    elseif strcmp(cfg.method, 'hbf')
      % coonvert meshes to m
      geometry =  ft_convert_units(geometry,'m');
      headmodel = ft_headmodel_hbf(geometry, 'conductivity', cfg.conductivity, 'isolatedsource', cfg.isolatedsource, 'checkmesh', cfg.checkmesh);
    else
      headmodel = ft_headmodel_openmeeg(geometry, 'conductivity', cfg.conductivity, 'isolatedsource', cfg.isolatedsource, 'tissue', cfg.tissue);
    end

  case 'concentricspheres'
    cfg.fitind = ft_getopt(cfg, 'fitind');
    cfg.order  = ft_getopt(cfg, 'order');

    % the low-level functions needs surface points, triangles are not needed
    if input_mesh || input_pos
      geometry = data;
    elseif input_seg
      tmpcfg   = keepfields(cfg, {'numvertices', 'tissue', 'spmversion'});
      geometry = ft_prepare_mesh(tmpcfg, data);
    elseif input_elec
      geometry.pos = data.chanpos;
      geometry.unit = data.unit;
    elseif ~isempty(cfg.headshape) && ischar(cfg.headshape)
      geometry = ft_read_headshape(cfg.headshape);
    elseif ~isempty(cfg.headshape)
      geometry = fixpos(cfg.headshape);
    else
      ft_error('You must give a mesh, segmented MRI, sensor data type, or cfg.headshape');
    end

    headmodel = ft_headmodel_concentricspheres(geometry, 'conductivity', cfg.conductivity, 'fitind', cfg.fitind, 'order', cfg.order);

  case 'halfspace'
    if input_mesh || input_pos
      geometry = data;
    else
      ft_error('a surface mesh is required as input for the halfspace method');
    end
    if isempty(cfg.point)
      ft_error('cfg.point is required for halfspace method');
    end

    headmodel = ft_headmodel_halfspace(geometry, cfg.point, 'conductivity', cfg.conductivity, 'sourcemodel', cfg.submethod);

  case 'infinite'
    % this takes no input arguments
    headmodel = ft_headmodel_infinite();

  case {'localspheres' 'singlesphere' 'singleshell'}
    cfg.grad = ft_getopt(cfg, 'grad');           % used for localspheres

    % these three methods all require a single set of surface points
    if input_mesh || input_pos
      geometry = data;
    elseif input_seg
      tmpcfg = keepfields(cfg, {'numvertices' 'spmversion'});
      if ~isempty(cfg.tissue)
        % extract the specified surface
        tmpcfg.tissue = cfg.tissue;
        geometry = ft_prepare_mesh(tmpcfg, data);
      else
        % try to extract either the brain or scalp surface
        geometry = [];
        if isempty(geometry)
          try
            tmpcfg.tissue = 'brain';
            geometry = ft_prepare_mesh(tmpcfg, data);
          catch
            me = lasterror;
            if isequal(me.identifier, 'MATLAB:mex:ErrInvalidMEXFile')
              % SPM8 mex file issues are common on macOS, these should not remain invisible
              rethrow(me);
            end
          end
        end
        if isempty(geometry)
          try
            tmpcfg.tissue = 'scalp';
            geometry = ft_prepare_mesh(tmpcfg, data);
          catch
            me = lasterror;
            if isequal(me.identifier, 'MATLAB:mex:ErrInvalidMEXFile')
              % SPM8 mex file issues are common on macOS, these should not remain invisible
              rethrow(me);
            end
          end
        end
        if isempty(geometry)
          ft_error('please specify cfg.tissue and pass an appropriate segmented MRI as input data')
        end
      end
    elseif input_elec
      geometry.pos = data.chanpos;
      geometry.unit = data.unit;
    elseif ~isempty(cfg.headshape)
      % get the surface describing the head shape
      [geometry.pos, geometry.tri] = headsurface([], [], 'headshape', cfg.headshape);
    elseif ~isempty(cfg.headmodel)
      % the CTF *.hdm file will be read further down
    else
      ft_error('this requires a mesh, set of surface points or a segmented mri');
    end

    switch cfg.method
      case 'singlesphere'
        if ~isempty(cfg.headmodel)
          % read the volume conduction model from a CTF *.hdm file
          tmp = ft_read_headmodel(cfg.headmodel);
          try
            % the single sphere is contained in the "orig" field
            headmodel = [];
            headmodel.r =  tmp.orig.MEG_Sphere.RADIUS;
            headmodel.o = [tmp.orig.MEG_Sphere.ORIGIN_X tmp.orig.MEG_Sphere.ORIGIN_Y tmp.orig.MEG_Sphere.ORIGIN_Z];
            headmodel.unit = 'cm';
          catch
            ft_error('the volume conduction model in "%s" is invalid', cfg.headmodel);
          end
        else
          % construct the volume conduction model
          headmodel = ft_headmodel_singlesphere(geometry, 'conductivity', cfg.conductivity);
        end % headmodel

      case 'localspheres'
        if ~isempty(cfg.headmodel)
          % read the volume conduction model from a CTF *.hdm file
          tmp = ft_read_headmodel(cfg.headmodel);
          try
            headmodel = [];
            headmodel.label = tmp.label;
            headmodel.r = tmp.r;
            headmodel.o = tmp.o;
            headmodel.unit = 'cm';
          catch
            ft_error('the volume conduction model in "%s" is invalid', cfg.headmodel);
          end
        else
          % construct the volume conduction model
          cfg.grad = ft_getopt(cfg, 'grad');
          if isempty(cfg.grad)
            ft_error('for cfg.method = %s, you must also supply cfg.grad', cfg.method);
          end
          headmodel = ft_headmodel_localspheres(geometry, cfg.grad, 'feedback', cfg.feedback, 'radius', cfg.radius, 'maxradius', cfg.maxradius, 'baseline', cfg.baseline, 'singlesphere', cfg.singlesphere);
        end % headmodel

      case 'singleshell'
        cfg.order  = ft_getopt(cfg, 'order');
        if ~isfield(geometry, 'tri')
          tmpcfg = [];
          tmpcfg.headshape = geometry;
          geometry = ft_prepare_mesh(tmpcfg);
        end
        headmodel = ft_headmodel_singleshell(geometry, 'order', cfg.order);

      otherwise
        ft_error('unsupported method %s', cfg.method);
    end % switch method

  case {'simbio'}
    if input_elec || isfield(data, 'pos') || input_mesh
      geometry = data; % more serious checks of validity of the mesh occur inside ft_headmodel_simbio
    else
      ft_error('you must provide a mesh with tetrahedral or hexahedral elements, where each element has a scalar or tensor conductivity');
    end
    headmodel = ft_headmodel_simbio(geometry, 'conductivity', cfg.conductivity);

  case {'duneuro'}
    if input_mesh
      geometry = data; % more serious checks of validity of the mesh occur inside ft_headmodel_duneuro
    else
      error('You must provide a mesh with tetrahedral or hexahedral elements, where each element has a scalar or tensor conductivity');
    end
    headmodel = ft_headmodel_duneuro(geometry, 'grid_filename', cfg.grid_filename, 'tensors_filename', cfg.tensors_filename,...
      'conductivity', cfg.conductivity, 'duneuro_settings', cfg.duneuro);
  case {'fns'}
    if input_seg
      data = ft_datatype_segmentation(data, 'segmentationstyle', 'indexed');
    else
      ft_error('segmented MRI must be given as data input')
    end
    sens = ft_fetch_sens(cfg, data);
    headmodel = ft_headmodel_fns(data.seg, 'tissue', cfg.tissue, 'tissueval', cfg.tissueval, 'tissuecond', cfg.conductivity, 'sens', sens, 'transform', cfg.transform);

  otherwise
    ft_error('unsupported method "%s"', cfg.method);
end % switch method

% ensure that the geometrical units are specified
if ~ft_headmodeltype(headmodel, 'infinite')
  headmodel = ft_determine_units(headmodel);
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble provenance
ft_postamble previous data
ft_postamble history headmodel
