function [grid] = prepare_leadfield(cfg, data)

% PREPARE_LEADFIELD computes the forward model for many dipole locations
% on a regular 2D or 3D grid and stores it for efficient inverse modelling
%
% Use as
%   [grid] = prepare_leadfield(cfg, data);
%
% It is neccessary to input the data on which you want to perform the
% inverse computations, since that data generally contain the gradiometer
% information and information about the channels that should be included in
% the forward model computation. The data structure can be either obtained
% from PREPROCESSING, FREQANALYSIS or TIMELOCKANALYSIS. If the data is empty,
% all channels will be included in the forward model.
%
% The positions of the sources are either specified as regular orthogonal grid
%   cfg.xgrid      = something like -20:1:20 or 'auto'
%   cfg.ygrid      = something like -20:1:20 or 'auto'
%   cfg.zgrid      = something like   0:1:20 or 'auto'
%   cfg.resolution = number (e.g. 1 cm) required for automatic grid generation
% or as a selected number of dipoles at locations of interest that can be
% obtained from an anatomical or functional MRI
%   cfg.grid.pos    = N*3 matrix with x,y,z-position
%   cfg.grid.ori    = N*3 matrix with orientation (optional)
%
% The sensor and volume conductor information is specified as
%   cfg.channel    = cell-array (optional, see CHANNELSELECTION)
%   cfg.hdmfile    = string (required)
%   cfg.gradfile   = string (optional, can already be present in data)
%   cfg.elecfile   = string (optional, can already be present in data)
% or alternatively
%   cfg.vol        = structure with volume conduction model
%   cfg.grad       = structure with gradiometer definition
%   cfg.elec       = structure with electrode definition
%
% Optionally, you can modify the leadfields by reducing the rank, or by
% normalizing each column.
%   cfg.reducerank  = 'no', or number (e.g. 2)
%   cfg.normalize   = 'no' (default), or 'yes'
%
% See also SOURCEANALYSIS

% Undocumented options for subspace projection
%   cfg.sel50p    = 'yes' or 'no'
%   cfg.lbex      = 'yes' or 'no'
%   cfg.lbexroi   = number

% Undocumented options for cortex segmentation
%   cfg.mri           = can be filename, MRI structure or segmented volume structure
%   cfg.sourceunits   = 'mm' or 'cm' (default is 'cm')
%   cfg.mriunits      = 'mm' or 'cm' (default is 'mm')
%   cfg.threshold     = 0.1, relative to the maximum value in the segmentation
%   cfg.smooth        = 5, smoothing in voxels

% Copyright (C) 2004, Robert Oostenveld
%
% $Log: prepare_leadfield.m,v $
% Revision 1.10  2005/11/08 11:21:58  roboos
% moved the handling of the normalize and reducerank options towards compute_leadfield
%
% Revision 1.9  2005/10/14 16:23:48  roboos
% fixed typo in fprintf (leedfield->leadfield)
%
% Revision 1.8  2005/09/29 12:52:25  jansch
% changed column-wise normalisation into matrix-wise normalisation
%
% Revision 1.7  2005/09/29 10:37:08  jansch
% fixed bug in improper normalisation
%
% Revision 1.6  2005/06/29 12:46:29  roboos
% the following changes were done on a large number of files at the same time, not all of them may apply to this specific file
% - added try-catch around the inclusion of input data configuration
% - moved cfg.version, cfg.previous and the assignment of the output cfg to the end
% - changed help comments around the configuration handling
% - some changes in whitespace
%
% Revision 1.5  2005/06/08 13:40:33  roboos
% replaced the specific call to either meg_leadfield or eeg_leadfield to the generic compute_leadfield
%
% Revision 1.4  2005/05/17 17:50:37  roboos
% changed all "if" occurences of & and | into && and ||
% this makes the code more compatible with Octave and also seems to be in closer correspondence with Matlab documentation on shortcircuited evaluation of sequential boolean constructs
%
% Revision 1.3  2005/04/25 16:05:11  roboos
% includedupdated options for cortex segmentation based grid in the documentation (not in main documentation)
%
% Revision 1.2  2005/01/10 16:28:10  roboos
% some changes in order of code, no functional change
%
% Revision 1.1  2004/12/24 14:31:09  roboos
% Renamed precompute_leadfield into prepare_leadfield. For backward
% compatibility in peoples scripts, precompute_leadfield still will remain
% available for some time (it prints a warning and calls prepare_leadfield).
% The newly committed prepare_leadfield is a copy of precompute_leadfield
% revision 1.6.
%
% Revision 1.6  2004/12/08 18:00:13  roboos
% implemented consistent method of selecting a subset of channels for
% forward and inverse computations using cfg.channel and updated the
% ducumentation
%
% Revision 1.5  2004/10/27 16:13:53  roboos
% changed progress indication into private progress function
% added sel50p subspace projection
%
% Revision 1.4  2004/09/22 08:49:12  roboos
% added computation of LBEX
%
% Revision 1.3  2004/08/19 06:31:31  roboos
% removed consistency check for cfg.xgrid etc., this is now done in prepare_dipole_grid
%
% Revision 1.2  2004/08/17 08:53:55  roboos
% added cfg to the output, updated help
%
% Revision 1.1  2004/08/16 09:10:47  roboos
% initial version, this replaces the leadfield precomputation in sourceanalysis with some extra options
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set the defaults
if ~isfield(cfg, 'normalize'),        cfg.normalize  = 'no';          end
if ~isfield(cfg, 'reducerank'),       cfg.reducerank = 'no';          end
if ~isfield(cfg, 'lbex'),             cfg.lbex       = 'no';          end
if ~isfield(cfg, 'sel50p'),           cfg.sel50p     = 'no';          end
if ~isfield(cfg, 'feedback'),         cfg.feedback   = 'text';        end

if strcmp(cfg.sel50p, 'yes') && strcmp(cfg.lbex, 'yes')
  error('subspace projection with either lbex or sel50p is mutually exclusive');
end

if nargin<2
  data = [];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% collect and preprocess the electrodes/gradiometer and head model
[vol, sens, cfg] = prepare_vol_sens(cfg, data);

% construct the grid on which the scanning will be done
[grid, cfg] = prepare_dipole_grid(cfg, vol, sens);

progress('init', cfg.feedback, 'computing leadfield');
for i=1:length(grid.inside)
  % compute the leadfield on all grid positions inside the brain
  progress(i/length(grid.inside), 'computing leadfield %d/%d\n', i, length(grid.inside));
  dipindx = grid.inside(i);
  grid.leadfield{dipindx} = compute_leadfield(grid.pos(dipindx,:), sens, vol, 'reducerank', cfg.reducerank, 'normalize', cfg.normalize);

  if isfield(cfg, 'grid') && isfield(cfg.grid, 'ori')
    % multiply with the normalized dipole moment to get the leadfield in the desired orientation
    grid.leadfield{dipindx} = grid.leadfield{dipindx} * grid.ori(dipindx,:)';
  end

end % for all grid locations inside the brain
progress('close');

% fill the positions outside the brain with NaNs
grid.leadfield(grid.outside) = {nan};

% compute the 50 percent channel selection subspace projection
if strcmp(cfg.sel50p, 'yes')
  grid = sel50p(cfg, grid, sens);
end

% compute the local basis function expansion (LBEX) subspace projection
if strcmp(cfg.lbex, 'yes')
  grid = lbex(cfg, grid);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% add version information to the configuration
try
  % get the full name of the function
  cfg.version.name = mfilename('fullpath');
catch
  % required for compatibility with Matlab versions prior to release 13 (6.5)
  [st, i] = dbstack;
  cfg.version.name = st(i);
end
cfg.version.id = '$Id: prepare_leadfield.m,v 1.10 2005/11/08 11:21:58 roboos Exp $';
% remember the configuration details of the input data
try, cfg.previous = data.cfg; end
% remember the exact configuration details in the output 
grid.cfg = cfg;

