function [grid, cfg] = ft_prepare_sourcemodel(cfg, headmodel, sens)

% FT_PREPARE_SOURCEMODEL constructs a source model, for example a 3D grid or a
% cortical sheet. The source model that can be used for source reconstruction,
% beamformer scanning, linear estimation and MEG interpolation.
%
% Use as
%   grid = ft_prepare_sourcemodel(cfg)
%
% where the configuration structure contains the details on how the source
% model should be constructed.
%
% A source model can be constructed based on
%   - regular 3D grid with explicit specification
%   - regular 3D grid with specification of the resolution
%   - regular 3D grid, based on segmented MRI, restricted to gray matter
%   - regular 3D grid, based on a warped template grid, based on the MNI brain
%   - surface grid based on the brain surface from the volume conduction model
%   - surface grid based on the head surface from an external file
%   - cortical sheet that was created in MNE or Freesurfer
%   - using user-supplied grid positions, which can be regular or irregular
% The approach that will be used depends on the configuration options that
% you specify.
%
% Configuration options for generating a regular 3D grid
%   cfg.grid.xgrid      = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
%   cfg.grid.ygrid      = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
%   cfg.grid.zgrid      = vector (e.g.   0:1:20) or 'auto' (default = 'auto')
%   cfg.grid.resolution = number (e.g. 1 cm) for automatic grid generation
%
% Configuration options for a predefined grid
%   cfg.grid.pos        = N*3 matrix with position of each source
%   cfg.grid.inside     = N*1 vector with boolean value whether grid point is inside brain (optional)
%   cfg.grid.dim        = [Nx Ny Nz] vector with dimensions in case of 3D grid (optional)
%
% The following fields are not used in this function, but will be copied along to the output
%   cfg.grid.leadfield
%   cfg.grid.filter or alternatively cfg.grid.avg.filter
%   cfg.grid.subspace
%   cfg.grid.lbex
%
% Configuration options for a warped MNI grid
%   cfg.mri             = can be filename or MRI structure, containing the individual anatomy
%   cfg.grid.warpmni    = 'yes'
%   cfg.grid.resolution = number (e.g. 6) of the resolution of the
%                         template MNI grid, defined in mm
%   cfg.grid.template   = specification of a template grid (grid structure), or a
%                         filename of a template grid (defined in MNI space),
%                         either cfg.grid.resolution or cfg.grid.template needs
%                         to be defined. If both are defined cfg.grid.template
%                         prevails
%   cfg.grid.nonlinear  = 'no' (or 'yes'), use non-linear normalization
%
% Configuration options for cortex segmentation, i.e. for placing dipoles in grey matter
%   cfg.mri           = can be filename, MRI structure or segmented MRI structure
%   cfg.threshold     = 0.1, relative to the maximum value in the segmentation
%   cfg.smooth        = 5, smoothing in voxels
%
% Configuration options for reading a cortical sheet from file
%   cfg.headshape     = string, should be a *.fif file
%
% The EEG or MEG sensor positions can be present in the data or can be specified as
%   cfg.elec          = structure with electrode positions, see FT_DATATYPE_SENS
%   cfg.grad          = structure with gradiometer definition, see FT_DATATYPE_SENS
% or alternatively
%   cfg.elecfile      = name of file containing the electrode positions, see FT_READ_SENS
%   cfg.gradfile      = name of file containing the gradiometer definition, see FT_READ_SENS
%
% The headmodel or volume conduction model can be specified as
%   cfg.headmodel       = structure with volume conduction model, see FT_PREPARE_HEADMODEL
%
% Other configuration options
%   cfg.grid.unit       = string, can be 'mm', 'cm', 'm' (default is automatic)
%   cfg.grid.tight   = 'yes' or 'no' (default is automatic)
%   cfg.inwardshift  = number, how much should the innermost surface be moved inward to constrain
%                      sources to be considered inside the source compartment (default = 0)
%   cfg.moveinward   = number, move dipoles inward to ensure a certain distance to the innermost
%                      surface of the source compartment (default = 0)
%   cfg.spherify     = 'yes' or 'no', scale the source model so that it fits inside a sperical
%                      volume conduction model (default = 'no')
%   cfg.symmetry     = 'x', 'y' or 'z' symmetry for two dipoles, can be empty (default = [])
%   cfg.headshape    = a filename for the headshape, a structure containing a single surface,
%                      or a Nx3 matrix with headshape surface points (default = [])
%   cfg.spmversion   = string, 'spm2', 'spm8', 'spm12' (default = 'spm8')
%
% See also FT_PREPARE_LEADFIELD, FT_PREPARE_HEADMODEL, FT_SOURCEANALYSIS,
% FT_DIPOLEFITTING, FT_MEGREALIGN

% Copyright (C) 2004-2013, Robert Oostenveld
%
% 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 debug
ft_preamble provenance headmodel sens
ft_preamble trackconfig

% 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', 'mriunits');
cfg = ft_checkconfig(cfg, 'renamed', {'hdmfile', 'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed', {'vol',     'headmodel'});

% put the low-level options pertaining to the dipole grid in their own field
cfg = ft_checkconfig(cfg, 'renamed', {'tightgrid', 'tight'}); % this is moved to cfg.grid.tight by the subsequent createsubcfg
cfg = ft_checkconfig(cfg, 'renamed', {'sourceunits', 'unit'}); % this is moved to cfg.grid.unit by the subsequent createsubcfg
cfg = ft_checkconfig(cfg, 'createsubcfg',  {'grid'});

% set the defaults
cfg.moveinward = ft_getopt(cfg, 'moveinward',  []); % the default is automatic and depends on a triangulation being present
cfg.spherify   = ft_getopt(cfg, 'spherify',  'no');
cfg.headshape  = ft_getopt(cfg, 'headshape',  []);
cfg.symmetry   = ft_getopt(cfg, 'symmetry',   []);
cfg.spmversion = ft_getopt(cfg, 'spmversion', 'spm8');
cfg.grid       = ft_getopt(cfg, 'grid',       []);
cfg.grid.unit  = ft_getopt(cfg.grid, 'unit',  'auto');

% this code expects the inside to be represented as a logical array
if isfield(cfg, 'grid')
  cfg.grid = ft_checkconfig(cfg.grid, 'renamed',  {'pnt' 'pos'});
  if isfield(cfg.grid, 'template')
    cfg.grid.template = ft_checkconfig(cfg.grid.template, 'renamed',  {'pnt' 'pos'});
  end
end
cfg = ft_checkconfig(cfg, 'inside2logical', 'yes');

if ~isfield(cfg, 'headmodel') && nargin>1
  % put it in the configuration structure
  % this is for backward compatibility, 13 Januari 2011
  cfg.headmodel = headmodel;
end

if ~isfield(cfg, 'grad') && ~isfield(cfg, 'elec') && nargin>2
  % put it in the configuration structure
  % this is for backward compatibility, 13 Januari 2011
  cfg.grad = sens;
end

if isfield(cfg.grid, 'resolution') && isfield(cfg.grid, 'xgrid') && ~ischar(cfg.grid.xgrid)
  ft_error('You cannot specify cfg.grid.resolution and an explicit cfg.grid.xgrid simultaneously');
end
if isfield(cfg.grid, 'resolution') && isfield(cfg.grid, 'ygrid') && ~ischar(cfg.grid.ygrid)
  ft_error('You cannot specify cfg.grid.resolution and an explicit cfg.grid.ygrid simultaneously');
end
if isfield(cfg.grid, 'resolution') && isfield(cfg.grid, 'zgrid') && ~ischar(cfg.grid.zgrid)
  ft_error('You cannot specify cfg.grid.resolution and an explicit cfg.grid.zgrid simultaneously');
end

% the source model can be constructed in a number of ways
basedongrid       = isfield(cfg.grid, 'xgrid') && ~ischar(cfg.grid.xgrid);                              % regular 3D grid with explicit specification
basedonpos        = isfield(cfg.grid, 'pos');                                                           % using user-supplied grid positions, which can be regular or irregular
basedonshape      = ~isempty(cfg.headshape);                                                            % surface grid based on inward shifted head surface from external file
basedonmri        = isfield(cfg, 'mri') && ~(isfield(cfg.grid, 'warpmni') && istrue(cfg.grid.warpmni)); % regular 3D grid, based on segmented MRI, restricted to gray matter
basedonmni        = isfield(cfg, 'mri') &&  (isfield(cfg.grid, 'warpmni') && istrue(cfg.grid.warpmni)); % regular 3D grid, based on warped MNI template
basedonvol        = false;                                                                              % surface grid based on inward shifted brain surface from volume conductor
basedoncortex     = isfield(cfg, 'headshape') && (iscell(cfg.headshape) || any(ft_filetype(cfg.headshape, {'neuromag_fif', 'freesurfer_triangle_binary', 'caret_surf', 'gifti'}))); % cortical sheet from external software such as Caret or FreeSurfer, can also be two separate hemispheres
basedonresolution = isfield(cfg.grid, 'resolution') && ~basedonmri && ~basedonmni;                      % regular 3D grid with specification of the resolution

if basedonshape && basedoncortex
  % treating it as cortical sheet has preference
  basedonshape = false;
end

if basedongrid && basedonpos
  % fall back to default behaviour, in which the pos overrides the grid
  basedongrid = false;
end

if ~any([basedonresolution basedongrid basedonpos basedonshape basedonmri basedoncortex basedonmni]) && ~isempty(cfg.headmodel)
  % fall back to default behaviour, which is to create a surface grid (e.g. used in MEGREALIGN)
  basedonvol = 1;
end

% these are mutually exclusive, but printing all requested methods here
% facilitates debugging of weird configs. Also specify the defaults here to
% keep the overview
if basedonresolution
  fprintf('creating dipole grid based on automatic 3D grid with specified resolution\n');
  cfg.grid.xgrid  = ft_getopt(cfg.grid, 'xgrid',  'auto');
  cfg.grid.ygrid  = ft_getopt(cfg.grid, 'ygrid',  'auto');
  cfg.grid.zgrid  = ft_getopt(cfg.grid, 'zgrid',  'auto');
  cfg.inwardshift = ft_getopt(cfg,      'inwardshift', 0); %in this case for inside detection, FIXME move to cfg.grid
  cfg.grid.tight  = ft_getopt(cfg.grid, 'tight',   'yes');
end

if basedongrid
  fprintf('creating dipole grid based on user specified 3D grid\n');
  cfg.inwardshift = ft_getopt(cfg,      'inwardshift', 0); %in this case for inside detection, FIXME move to cfg.grid
  cfg.grid.tight  = ft_getopt(cfg.grid, 'tight',   'yes');
end

if basedonpos
  fprintf('creating dipole grid based on user specified dipole positions\n');
  cfg.inwardshift = ft_getopt(cfg,      'inwardshift', 0); %in this case for inside detection, FIXME move to cfg.grid
  cfg.grid.tight  = ft_getopt(cfg.grid, 'tight',    'no');
end

if basedonshape
  fprintf('creating dipole grid based on inward-shifted head shape\n');
  cfg.inwardshift = ft_getopt(cfg,      'inwardshift',  0); %in this case for inside detection, FIXME move to cfg.grid
  cfg.spheremesh  = ft_getopt(cfg,      'spheremesh', 642); % FIXME move spheremesh to cfg.grid
  cfg.grid.tight  = ft_getopt(cfg.grid, 'tight',    'yes');
end

if basedoncortex
  cfg.grid.tight  = ft_getopt(cfg.grid, 'tight', 'yes');
end

if basedonmri
  fprintf('creating dipole grid based on an anatomical volume\n');
  cfg.threshold   = ft_getopt(cfg,      'threshold', 0.1); % relative
  cfg.smooth      = ft_getopt(cfg,      'smooth',      5);   % in voxels
  cfg.grid.tight  = ft_getopt(cfg.grid, 'tight',   'yes');
end

if basedonvol
  fprintf('creating dipole grid based on inward-shifted brain surface from volume conductor model\n');
  cfg.inwardshift = ft_getopt(cfg,      'inwardshift',   0); %in this case for inside detection, FIXME move to cfg.grid
  cfg.spheremesh  = ft_getopt(cfg,      'spheremesh',  642); % FIXME move spheremesh to cfg.grid
  cfg.grid.tight  = ft_getopt(cfg.grid, 'tight',      'no');
end

if basedonmni
  cfg.grid.tight     = ft_getopt(cfg.grid, 'tight',       'no');
  cfg.grid.nonlinear = ft_getopt(cfg.grid, 'nonlinear',   'no');
end

% these are mutually exclusive
if sum([basedonresolution basedongrid basedonpos basedonshape basedonmri basedonvol basedoncortex basedonmni])~=1
  ft_error('incorrect cfg specification for constructing a dipole grid');
end

if (isfield(cfg, 'smooth') && ~strcmp(cfg.smooth, 'no')) || basedonmni
  % check that the preferred SPM version is on the path
  ft_hastoolbox(cfg.spmversion, 1);
end

% start with an empty grid
grid = [];

% get the volume conduction model
try
  headmodel = ft_fetch_vol(cfg);
catch
  headmodel = [];
end

% get the gradiometer or electrode definition
try
  sens = ft_fetch_sens(cfg);
catch
  sens = [];
end

if strcmp(cfg.grid.unit, 'auto')
  if isfield(cfg.grid, 'pos') && size(cfg.grid.pos,1)>10
    % estimate the units based on the existing source positions
    cfg.grid = rmfield(cfg.grid, 'unit'); % remove 'auto' and have ft_determine_units determine it properly
    cfg.grid = ft_determine_units(cfg.grid);
  elseif ~isempty(sens)
    % copy the units from the sensor array
    cfg.grid.unit = sens.unit;
  elseif ~isempty(headmodel)
    % copy the units from the volume conduction model
    cfg.grid.unit = headmodel.unit;
  else
    ft_warning('assuming "cm" as default source units');
    cfg.grid.unit = 'cm';
  end
end

% convert the sensor array to the desired units for the source model
if ~isempty(sens)
  sens = ft_convert_units(sens, cfg.grid.unit);
end

% convert the head model to the desired units for the source model
if ~isempty(headmodel)
  headmodel = ft_convert_units(headmodel, cfg.grid.unit);
end

if basedonresolution
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % construct a regular 3D grid that spans a box encompassing all electrode
  % or gradiometer coils, this will typically also cover the complete brain
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if ~isempty(sens) && isfield(headmodel, 'chanpos')
    % determine the bounding box of the sensor array
    minpos = min(sens.chanpos,[],1);
    maxpos = max(sens.chanpos,[],1);
  elseif ~isempty(headmodel)
    % determine the bounding box of the volume conduction model
    if isfield(headmodel, 'bnd') && isfield(headmodel.bnd, 'pos')
      pos = cat(1, headmodel.bnd(:).pos);
    elseif isfield(headmodel, 'bnd') && isfield(headmodel.bnd, 'pnt')
      pos = cat(1, headmodel.bnd(:).pnt);
    elseif isfield(headmodel, 'pos')
      pos = headmodel.pos;
    elseif ft_voltype(headmodel, 'localspheres')
      pos = headsurface(headmodel, sens);
    elseif ft_voltype(headmodel, 'singlesphere')
      pos = [
        headmodel.o - headmodel.r
        headmodel.o + headmodel.r
        ];
    elseif ft_voltype(headmodel, 'concentricspheres')
      pos = [
        headmodel.o - max(headmodel.r)
        headmodel.o + max(headmodel.r)
        ];
    end
    minpos = min(pos,[],1);
    maxpos = max(pos,[],1);
    % add a few percent on either side
    minpos(minpos<0) = minpos(minpos<0).*1.08;
    maxpos(maxpos>0) = maxpos(maxpos>0).*1.08;
    minpos(minpos>0) = minpos(minpos>0).*0.92;
    maxpos(maxpos<0) = maxpos(maxpos<0).*0.92;
  else
    ft_error('creating a 3D grid based on resolution requires either sensor positions or a headmodel to estimate the extent');
  end
  
  fprintf('creating 3D grid with %g %s resolution\n', cfg.grid.resolution, cfg.grid.unit);
  
  % round the bounding box limits to the nearest cm
  switch cfg.grid.unit
    case 'm'
      minpos = floor(minpos*100)/100;
      maxpos = ceil(maxpos*100)/100;
    case 'cm'
      minpos = floor(minpos);
      maxpos = ceil(maxpos);
    case 'mm'
      minpos = floor(minpos/10)*10;
      maxpos = ceil(maxpos/10)*10;
  end
  
  if ischar(cfg.grid.xgrid) && strcmp(cfg.grid.xgrid, 'auto')
    grid.xgrid = minpos(1):cfg.grid.resolution:maxpos(1);
  end
  if ischar(cfg.grid.ygrid) && strcmp(cfg.grid.ygrid, 'auto')
    grid.ygrid = minpos(2):cfg.grid.resolution:maxpos(2);
  end
  if ischar(cfg.grid.zgrid) && strcmp(cfg.grid.zgrid, 'auto')
    grid.zgrid = minpos(3):cfg.grid.resolution:maxpos(3);
  end
  grid.dim   = [length(grid.xgrid) length(grid.ygrid) length(grid.zgrid)];
  [X, Y, Z]  = ndgrid(grid.xgrid, grid.ygrid, grid.zgrid);
  grid.pos   = [X(:) Y(:) Z(:)];
  grid.unit  = cfg.grid.unit;
  fprintf('initial 3D grid dimensions are [%d %d %d]\n', grid.dim(1), grid.dim(2), grid.dim(3));
end

if basedongrid
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % a detailed xgrid/ygrid/zgrid has been specified, the other details
  % still need to be determined
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  grid.xgrid = cfg.grid.xgrid;
  grid.ygrid = cfg.grid.ygrid;
  grid.zgrid = cfg.grid.zgrid;
  grid.dim   = [length(grid.xgrid) length(grid.ygrid) length(grid.zgrid)];
  [X, Y, Z]  = ndgrid(grid.xgrid, grid.ygrid, grid.zgrid);
  grid.pos   = [X(:) Y(:) Z(:)];
  grid.unit  = cfg.grid.unit;
end

if basedonpos
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % a grid is already specified in the configuration, reuse as much of the
  % prespecified grid as possible (but only known objects)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  grid = keepfields(cfg.grid, {'pos', 'unit', 'xgrid', 'ygrid', 'zgrid', 'mom', 'tri', 'dim', 'transform', 'inside', 'lbex', 'subspace', 'leadfield', 'filter', 'label', 'leadfielddimord'});
end

if basedonmri
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % construct a grid based on the segmented MRI that is provided in the
  % configuration, only voxels in gray matter will be used
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if ischar(cfg.mri)
    mri = ft_read_mri(cfg.mri);
  else
    mri = cfg.mri;
  end

  % ensure the mri to have units
  if ~isfield(mri, 'unit')
    mri = ft_determine_units(mri);
  end

  if ~isfield(cfg.grid, 'resolution')
    switch cfg.grid.unit
      case 'mm'
        cfg.grid.resolution = 10;
      case 'cm'
        cfg.grid.resolution = 1;
      case 'dm'
        cfg.grid.resolution = 0.1;
      case 'm'
        cfg.grid.resolution = 0.01;
    end
  end

  issegmentation = false;
  if isfield(mri, 'gray')
    % this is not a boolean segmentation, but based on tissue probability
    % maps, being the original implementation here.
    dat = double(mri.gray);

    % apply a smoothing of a certain amount of voxels
    if ~strcmp(cfg.smooth, 'no');
      dat = volumesmooth(dat, cfg.smooth, 'MRI gray matter');
    end

  elseif isfield(mri, 'anatomy')
    % this could be a tpm stored on disk, i.e. the result of
    % ft_volumesegment. Reading it in always leads to the field 'anatomy'.
    % Note this could be any anatomical mask
    dat = double(mri.anatomy);

    % apply a smoothing of a certain amount of voxels
    if ~strcmp(cfg.smooth, 'no');
      dat = volumesmooth(dat, cfg.smooth, 'anatomy');
    end

  elseif ft_datatype(mri, 'segmentation')
    % this is a proper segmentation, where a set of boolean masks is in the
    % input, or and indexed volume, along with labels. FIXME for now still
    % only works for boolean volumes.
    issegmentation = true;
    fn = booleanfields(mri);
    if isempty(fn)
      % convert indexed segmentation into probabilistic
      mri = ft_datatype_segmentation(mri, 'segmentationstyle', 'probabilistic');
      fn  = booleanfields(mri);
    end

    dat = false(mri.dim);
    for i=1:numel(fn)
      if ~strcmp(cfg.smooth, 'no')
        mri.(fn{i}) = volumesmooth(double(mri.(fn{i})), cfg.smooth, fn{i}) > cfg.threshold;
      end
      dat = dat | mri.(fn{i});
    end
    dat = double(dat);
  else
    ft_error('cannot determine the format of the segmentation in cfg.mri');
  end


  % determine for each voxel whether it belongs to the grey matter
  fprintf('thresholding MRI data at a relative value of %f\n', cfg.threshold);
  head = dat./max(dat(:)) > cfg.threshold;

  % convert the source/functional data into the same units as the anatomical MRI
  scale = ft_scalingfactor(cfg.grid.unit, mri.unit);

  ind                 = find(head(:));
  fprintf('%d from %d voxels in the segmentation are marked as ''inside'' (%.0f%%)\n', length(ind), numel(head), 100*length(ind)/numel(head));
  [X,Y,Z]             = ndgrid(1:mri.dim(1), 1:mri.dim(2), 1:mri.dim(3));  % create the grid in MRI-coordinates
  posmri              = [X(ind) Y(ind) Z(ind)];                            % take only the inside voxels
  poshead             = ft_warp_apply(mri.transform, posmri);                 % transform to head coordinates
  resolution          = cfg.grid.resolution*scale;                                        % source and mri can be expressed in different units (e.g. cm and mm)
  xgrid               = floor(min(poshead(:,1))):resolution:ceil(max(poshead(:,1)));      % create the grid in head-coordinates
  ygrid               = floor(min(poshead(:,2))):resolution:ceil(max(poshead(:,2)));      % with 'consistent' x,y,z definitions
  zgrid               = floor(min(poshead(:,3))):resolution:ceil(max(poshead(:,3)));
  [X,Y,Z]             = ndgrid(xgrid,ygrid,zgrid);
  pos2head            = [X(:) Y(:) Z(:)];
  pos2mri             = ft_warp_apply(inv(mri.transform), pos2head);        % transform to MRI voxel coordinates
  pos2mri             = round(pos2mri);
  inside              = getinside(pos2mri, head);                           % use helper subfunction

  grid.pos            = pos2head/scale;                                     % convert to source units
  grid.xgrid          = xgrid/scale;                                        % convert to source units
  grid.ygrid          = ygrid/scale;                                        % convert to source units
  grid.zgrid          = zgrid/scale;                                        % convert to source units
  grid.dim            = [length(grid.xgrid) length(grid.ygrid) length(grid.zgrid)];
  grid.inside         = inside(:);
  grid.unit           = cfg.grid.unit;

  if issegmentation
    % pass on the segmentation information on the grid points, the
    % individual masks have been smoothed above
    fn = booleanfields(mri);
    for i=1:numel(fn)
      grid.(fn{i}) = getinside(pos2mri, mri.(fn{i}));
    end
    % convert back is not in general possible because the masks can be
    % overlapping due to smoothing
    % grid = ft_datatype_segmentation(grid, 'segmentationstyle', segstyle);
  end
  fprintf('the full grid contains %d grid points\n', numel(grid.inside));
  fprintf('%d grid points are mared as inside the brain\n',  sum(grid.inside));
end

if basedoncortex
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % read it from a *.fif file that was created using Freesurfer and MNE
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if iscell(cfg.headshape)
    % FIXME loop over all files, this should be two hemispheres
    keyboard
  else
    shape = ft_read_headshape(cfg.headshape);
  end
  % ensure that the headshape is in the same units as the source
  shape     = ft_convert_units(shape, cfg.grid.unit);
  % return both the vertices and triangles from the cortical sheet
  grid.pos  = shape.pos;
  grid.tri  = shape.tri;
  grid.unit = shape.unit;
end

if basedonshape
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % use the headshape  to make a superficial dipole layer (e.g.
  % for megrealign). Assume that all points are inside the volume.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % get the surface describing the head shape
  if isstruct(cfg.headshape) && isfield(cfg.headshape, 'pos')
    % use the headshape surface specified in the configuration
    headshape = cfg.headshape;
  elseif isnumeric(cfg.headshape) && size(cfg.headshape,2)==3
    % use the headshape points specified in the configuration
    headshape.pos = cfg.headshape;
  elseif ischar(cfg.headshape)
    % read the headshape from file
    headshape = ft_read_headshape(cfg.headshape);
  else
    ft_error('cfg.headshape is not specified correctly')
  end
  % ensure that the headshape is in the same units as the source
  headshape = ft_convert_units(headshape, cfg.grid.unit);
  if ~isfield(headshape, 'tri')
    % generate a closed triangulation from the surface points
    headshape.pos = unique(headshape.pos, 'rows');
    headshape.tri = projecttri(headshape.pos);
  end
  % please note that cfg.inwardshift should be expressed in the units consistent with cfg.grid.unit
  grid.pos     = headsurface([], [], 'headshape', headshape, 'inwardshift', cfg.inwardshift, 'npnt', cfg.spheremesh);
  grid.tri     = headshape.tri;
  grid.unit    = headshape.unit;
  grid.inside  = true(size(grid.pos,1),1);
end

if basedonvol
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % use the volume conduction model to make a superficial dipole layer (e.g.
  % for megrealign). Assume that all points are inside the volume.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % please note that cfg.inwardshift should be expressed in the units consistent with cfg.grid.unit
  grid.pos     = headsurface(headmodel, sens, 'inwardshift', cfg.inwardshift, 'npnt', cfg.spheremesh);
  grid.unit    = cfg.grid.unit;
  grid.inside  = true(size(grid.pos,1),1);
end

if basedonmni
  if ~isfield(cfg.grid, 'template') && ~isfield(cfg.grid, 'resolution')
    ft_error('you either need to specify the filename of a template grid in cfg.grid.template, or a resolution in cfg.grid.resolution');
  elseif isfield(cfg.grid, 'template')
    % let the template filename prevail
    fname = cfg.grid.template;
  elseif isfield(cfg.grid, 'resolution') && cfg.grid.resolution==round(cfg.grid.resolution)
    % use one of the templates that are in Fieldtrip, this requires a
    % resolution
    fname = ['standard_sourcemodel3d',num2str(cfg.grid.resolution),'mm.mat'];
  elseif isfield(cfg.grid, 'resolution') && cfg.grid.resolution~=round(cfg.grid.resolution)
    fname = ['standard_sourcemodel3d',num2str(floor(cfg.grid.resolution)),'point',num2str(10*(cfg.grid.resolution-floor(cfg.grid.resolution))),'mm.mat'];
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % check whether the mni template grid exists for the specified resolution
  % if not create it: FIXME (this needs to be done still)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % get the mri
  if ischar(cfg.mri)
    if ~exist(fname, 'file')
      ft_error('the MNI template grid based on the specified resolution does not exist');
    end
    mri = ft_read_mri(cfg.mri);
  else
    mri = cfg.mri;
  end

  % get the template grid
  if ischar(fname)
    mnigrid = load(fname, 'sourcemodel');
    mnigrid = mnigrid.sourcemodel;
  else
    mnigrid = cfg.grid.template;
  end

  % ensure these to have units in mm, the conversion of the source model is done further down
  mri     = ft_convert_units(mri,     'mm');
  mnigrid = ft_convert_units(mnigrid, 'mm');

  % ensure that it is specified with logical inside
  mnigrid = fixinside(mnigrid);

  % spatial normalisation of mri and construction of subject specific dipole grid positions
  tmpcfg           = keepfields(cfg, {'spmversion', 'spmmethod'});
  tmpcfg.nonlinear = cfg.grid.nonlinear;
  if isfield(cfg.grid, 'templatemri')
    tmpcfg.template = cfg.grid.templatemri;
  end
  normalise = ft_volumenormalise(tmpcfg, mri);

  if ~isfield(normalise, 'params') && ~isfield(normalise, 'initial')
    fprintf('applying an inverse warp based on a linear transformation only\n');
    grid.pos = ft_warp_apply(inv(normalise.cfg.final), mnigrid.pos);
  else
    grid.pos = ft_warp_apply(inv(normalise.initial), ft_warp_apply(normalise.params, mnigrid.pos, 'sn2individual'));
  end
  if isfield(mnigrid, 'dim')
    grid.dim   = mnigrid.dim;
  end
  if isfield(mnigrid, 'tri')
    grid.tri   = mnigrid.tri;
  end
  grid.unit    = mnigrid.unit;
  grid.inside  = mnigrid.inside;
  grid.params  = normalise.params;
  grid.initial = normalise.initial;
  if ft_datatype(mnigrid, 'parcellation')
    % copy the boolean fields over
    grid = copyfields(mnigrid, grid, booleanfields(mnigrid));
  end

end

% in most cases the source model will already be in the desired units, but e.g. for "basedonmni" it will be in 'mm'
% convert to the requested units
grid = ft_convert_units(grid, cfg.grid.unit);

if strcmp(cfg.spherify, 'yes')
  if ~ft_voltype(headmodel, 'singlesphere') && ~ft_voltype(headmodel, 'concentricspheres')
    ft_error('this only works for spherical volume conduction models');
  end
  % deform the cortex so that it fits in a unit sphere
  pos = mesh_spherify(grid.pos, [], 'shift', 'range');
  % scale it to the radius of the innermost sphere, make it a tiny bit smaller to
  % ensure that the support point with the exact radius 1 is still inside the sphere
  pos = pos*min(headmodel.r)*0.999;
  pos(:,1) = pos(:,1) + headmodel.o(1);
  pos(:,2) = pos(:,2) + headmodel.o(2);
  pos(:,3) = pos(:,3) + headmodel.o(3);
  grid.pos = pos;
end

if ~isempty(cfg.moveinward)
  % construct a triangulated boundary of the source compartment
  [pos1, tri1] = headsurface(headmodel, [], 'inwardshift', cfg.moveinward, 'surface', 'brain');
  inside = bounding_mesh(grid.pos, pos1, tri1);
  if ~all(inside)
    pos2 = grid.pos(~inside,:);
    [dum, pos3] = project_elec(pos2, pos1, tri1);
    grid.pos(~inside,:) = pos3;
  end
  if cfg.moveinward>cfg.inwardshift
    grid.inside  = true(size(grid.pos,1),1);
  end
end

% determine the dipole locations that are inside the source compartment of the
% volume conduction model, i.e. inside the brain
if ~isfield(grid, 'inside')
  grid.inside = ft_inside_vol(grid.pos, headmodel, 'grad', sens, 'headshape', cfg.headshape, 'inwardshift', cfg.inwardshift); % this returns a boolean vector
end

if strcmp(cfg.grid.tight, 'yes')
  fprintf('%d dipoles inside, %d dipoles outside brain\n', sum(grid.inside), sum(~grid.inside));
  fprintf('making tight grid\n');
  xmin = min(grid.pos(grid.inside,1));
  ymin = min(grid.pos(grid.inside,2));
  zmin = min(grid.pos(grid.inside,3));
  xmax = max(grid.pos(grid.inside,1));
  ymax = max(grid.pos(grid.inside,2));
  zmax = max(grid.pos(grid.inside,3));
  xmin_indx = find(grid.xgrid==xmin);
  ymin_indx = find(grid.ygrid==ymin);
  zmin_indx = find(grid.zgrid==zmin);
  xmax_indx = find(grid.xgrid==xmax);
  ymax_indx = find(grid.ygrid==ymax);
  zmax_indx = find(grid.zgrid==zmax);
  sel =       (grid.pos(:,1)>=xmin & grid.pos(:,1)<=xmax); % select all grid positions inside the tight box
  sel = sel & (grid.pos(:,2)>=ymin & grid.pos(:,2)<=ymax); % select all grid positions inside the tight box
  sel = sel & (grid.pos(:,3)>=zmin & grid.pos(:,3)<=zmax); % select all grid positions inside the tight box
  % update the grid locations that are marked as inside the brain
  grid.pos   = grid.pos(sel,:);
  % update the boolean fields, this requires the original dim
  fn = booleanfields(grid);
  for i=1:numel(fn)
    grid.(fn{i}) = grid.(fn{i})(sel);
  end
  grid.xgrid   = grid.xgrid(xmin_indx:xmax_indx);
  grid.ygrid   = grid.ygrid(ymin_indx:ymax_indx);
  grid.zgrid   = grid.zgrid(zmin_indx:zmax_indx);
  grid.dim     = [length(grid.xgrid) length(grid.ygrid) length(grid.zgrid)];
end
fprintf('%d dipoles inside, %d dipoles outside brain\n', sum(grid.inside), sum(~grid.inside));

% apply the symmetry constraint, i.e. add a symmetric dipole for each location that was defined sofar
if ~isempty(cfg.symmetry)
  if size(grid.pos,2)>3
    % sanity check, see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3119
    ft_warning('the construction of a symmetric dipole model requires to start with a Nx3 description of the dipole positions, discarding subsequent columns');
    grid.pos = grid.pos(:,1:3);
  end
  if strcmp(cfg.symmetry, 'x')
    reduce = [1 2 3];         % select the parameters [x1 y1 z1]
    expand = [1 2 3 1 2 3];   % repeat them as [x1 y1 z1 x1 y1 z1]
    mirror = [1 1 1 -1 1 1];  % multiply each of them with 1 or -1, resulting in [x1 y1 z1 -x1 y
  elseif strcmp(cfg.symmetry, 'y')
    reduce = [1 2 3];         % select the parameters [x1 y1 z1]
    expand = [1 2 3 1 2 3];   % repeat them as [x1 y1 z1 x1 y1 z1]
    mirror = [1 1 1 1 -1 1];  % multiply each of them with 1 or -1, resulting in [x1 y1 z1 x1 -y
  elseif strcmp(cfg.symmetry, 'z')
    reduce = [1 2 3];         % select the parameters [x1 y1 z1]
    expand = [1 2 3 1 2 3];   % repeat them as [x1 y1 z1 x1 y1 z1]
    mirror = [1 1 1 1 1 -1];  % multiply each of them with 1 or -1, resulting in [x1 y1 z1 x1 y1
  else
    ft_error('unrecognized symmetry constraint');
  end
  fprintf('each source describes two dipoles with symmetry along %s axis\n', cfg.symmetry);
  % expand the number of parameters from one (3) to two dipoles (6)
  grid.pos = grid.pos(:,expand) .* repmat(mirror, size(grid.pos,1), 1);
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble provenance grid
ft_postamble history    grid

%--------------------------------------------------------------
% helper function for basedonmri method to determine the inside
% returns a boolean vector
function inside = getinside(pos, mask)

% it might be that the box with the points does not completely fit into the
% mask
dim = size(mask);
sel = find(pos(:,1)<1 |  pos(:,1)>dim(1) | ...
  pos(:,2)<1 |  pos(:,2)>dim(2) | ...
  pos(:,3)<1 |  pos(:,3)>dim(3));
if isempty(sel)
  % use the efficient implementation
  inside = mask(sub2ind(dim, pos(:,1), pos(:,2), pos(:,3)));
else
  % only loop over the points that can be dealt with
  inside = false(size(pos,1), 1);
  for i=setdiff(1:size(pos,1), sel(:)')
    inside(i) = mask(pos(i,1), pos(i,2), pos(i,3));
  end
end

%--------------------------------------------------------------------------
% helper function to return the fieldnames of the boolean fields in a
% segmentation, should work both for volumetric and for source
function fn = booleanfields(mri)

fn = fieldnames(mri);
isboolean = false(1,numel(fn));
for i=1:numel(fn)
  if islogical(mri.(fn{i})) && isequal(numel(mri.(fn{i})),prod(mri.dim))
    isboolean(i) = true;
  end
end
fn  = fn(isboolean);
