cleaned up code and updated documentation, added an example script for testing

Author: robertoostenveld

ni2_activation.m

@@ -1,24 +1,19 @@
 function [data, time] = ni2_activation(varargin)
 
-% NI2_ACTIVATION creates a timecourse of simplified activity.
+% NI2_ACTIVATION creates a timecourse of the activity of a dipole in the brain.
 %
 % Use as
 %   [data, time] = ni2_activation
 %   [data, time] = ni2_activation('key1', value1, 'key2', value2)
 %
-% Input arguments:
-%   key-value pairs determine the shape of the activations.
-%
-%   frequency, scalar value (in Hz): frequency of oscillatory burst, default = 10 Hz   
-%   phase,     scalar value (in radians): phase of the oscillation at the
-%                peak latency, default = 0
-%   latency,   scalar value (in s): latency at which the amplitude of the
-%                oscillation peaks, default = 0.5
-%   length,    scalar value (in s): length of the activation time course,
-%                the sampling interval is fixed to 1 ms/sampl, default = 1
-%   ncycle,    scalar value: number of oscillatory cycles in the burst, default = 5
-%   powerup,    binary value: 1 = power burst (default); 0 = power loss
-%   fsample,   scalar value: sampling rate (Hz) (default 1000 Hz)
+% Options should be specified as key-value pairs and can be
+%   frequency = scalar value (in Hz): frequency of oscillatory burst, default = 10 Hz   
+%   phase     = scalar value (in radians): phase of the oscillation at the peak latency, default = 0
+%   latency   = scalar value (in s): latency at which the amplitude of the oscillation peaks, default = 0.5
+%   length    = scalar value (in s): length of the activation time course, the sampling interval is fixed to 1 ms/sampl, default = 1
+%   ncycle    = scalar value: number of oscillatory cycles in the burst, default = 5
+%   powerup   = binary value: 1 = power burst (default); 0 = power loss
+%   fsample   = scalar value: sampling rate (Hz) (default 1000 Hz)
 %
 % Output arguments:
 %   data = 1xN vector with the activation time course

ni2_example.m

@@ -0,0 +1,68 @@
+[activity, time] = ni2_activation;
+figure; plot(time, activity);
+
+%%
+
+[activity, time] = ni2_activation('frequency', 100);
+figure; plot(time, activity);
+
+%%
+
+[activity, time] = ni2_activation('frequency', 5);
+figure; plot(time, activity);
+
+%%
+
+headmodel = ni2_headmodel('type', 'singleshell');
+figure; ft_plot_headmodel(headmodel, 'axes', true);
+
+%%
+
+headmodel = ni2_headmodel('type', 'spherical', 'nshell', 3);
+figure; ft_plot_headmodel(headmodel, 'axes', true);
+
+%%
+
+headmodel = ni2_headmodel('type', 'spherical', 'nshell', 1);
+figure; ft_plot_headmodel(headmodel, 'axes', true);
+
+%%
+
+sens = ni2_sensors('type', 'eeg');
+figure; ft_plot_sens(sens, 'label', 'number')
+
+%%
+
+sens = ni2_sensors('type', 'ctf151');
+figure; ft_plot_sens(sens, 'label', 'label', 'coil', false)
+figure; ft_plot_sens(sens, 'label', 'label', 'coil', true)
+
+%%
+
+sens = ni2_sensors('type', 'meg');
+figure; ft_plot_sens(sens, 'label', 'number', 'coil', true)
+
+%%
+
+leadfield = ni2_leadfield(sens, headmodel, [0 0 6 1 0 0]);
+figure; ni2_topoplot(sens, leadfield)
+
+%%
+
+data = leadfield * activity;
+
+figure; ni2_topomovie(sens, data, time)
+
+%%
+
+sourcemodel = ni2_sourcemodel('type', 'mesh');
+
+figure; ft_plot_mesh(sourcemodel)
+figure; ft_plot_mesh(sourcemodel.pos(sourcemodel.inside,:))
+
+%%
+
+sourcemodel = ni2_sourcemodel('type', 'grid', 'resolution', 1);
+
+figure; ft_plot_mesh(sourcemodel)
+figure; ft_plot_mesh(sourcemodel.pos(sourcemodel.inside,:))

ni2_headmodel.m

@@ -1,15 +1,17 @@
 function headmodel = ni2_headmodel(varargin)
 
-% NI2_HEADMODEL generates a simple headmodel for forward modeling. 
+% NI2_HEADMODEL generates a simple headmodel for forward modeling.
 %
-% Use as:
-%  headmodel = ni2_headmodel('type', 'spherical', 'nshell', 1);
-%  headmodel = ni2_headmodel('type', 'spherical', 'nshell', 3);
+% Use as
+%   headmodel = ni2_headmodel('type', 'spherical', 'nshell', 1);
+%   headmodel = ni2_headmodel('type', 'spherical', 'nshell', 3);
 %
 % The first call creates a singlesphere volume conduction model that can be
 % used in combination with MEG sensor-arrays. The second call creates a
 % three concentric spheres model that can be used in combination with EEG
 % sensor-arrays.
+%
+% The headmodel is expressed in cm and in the CTF coordinate system.
 
 type = ft_getopt(varargin, 'type', 'spherical');
 n    = ft_getopt(varargin, 'nshell', 1);
@@ -25,6 +27,7 @@
         headmodel.type = 'singlesphere';
         headmodel      = ft_datatype_headmodel(headmodel);
         headmodel.coordsys = 'ctf';
+
       case 3
         headmodel.o    = [0 0 0];
         headmodel.r    = [0.88 0.92 1.00]*10;
@@ -33,15 +36,16 @@
         headmodel.type = 'concentricspheres';
         headmodel      = ft_datatype_headmodel(headmodel);
         headmodel.coordsys = 'ctf';
+
       otherwise
         error('number of spheres other than 1 or 3 is not supported');
     end
 
   case 'singleshell'
     % note that this is in CTF coordinates
-    load('Subject01_headmodel');
+    load('Subject01_headmodel.mat');
     headmodel.coordsys = 'ctf';
 
   otherwise
-    error('unknown type requested');
+    error('unknown headmodel type');
 end

ni2_leadfield.m

@@ -3,40 +3,45 @@
 % NI2_LEADFIELD generates a leadfield for a given source (or set of
 % sources) using a specified sensor array and volume conductor model.
 %
-% Use as:
+% Use as
 %   leadfield = ni2_leadfield(sens, headmodel, dippar)
 %
 % Input arguments:
-%   - sens      = a sensor array, obtained with ni2_sensors
-%   - headmodel = a volume conductor model, obtained with ni2_headmodel
-%   - dippar    = Nx6, or Nx3 matrix with dipole parameters.
+%   sens      = a sensor array, obtained with NI2_SENSORS
+%   headmodel = a volume conductor model, obtained with NI2_HEADMODEL
+%   dippar    = Nx3 or Nx6 matrix with dipole parameters, where N is the number of dipoles
+%
+% Each row in the dippar-matrix represents a source, the first columns 1-3 are the
+% dipole position parameters (x,y,z coordinates in Cartesian space), and the optional
+% columns 4-6 are the dipole moment parameters (x,y,z orientation, these can include
+% the amplitude).
 %
-% Each row in the dippar-matrix represents a source, the first 3 columns
-% are the position parameters (x,y,z coordinates in Cartesian space), and
-% the options 4-6 columns are the dipole moment parameters (x,y,z
-% orientation and amplitude).
-% 
 % Output argument:
-%   - leadfield = MxN or Mx(Nx3) matrix with the leadfield, where M is the
-%                 number of sensors, and N the number of sources
+%   leadfield = MxN or Mx(Nx3) matrix with the leadfield, where M is the number of sensors
+%               and N is the number of dipoles. There are three columns for dipoles with
+%               an unspecified orientation.
 
 ndip = size(dippar, 1);
 
 cfg = [];
+
 if strcmp(sens.chantype{1}, 'eeg')
   cfg.elec = sens;
 else
   cfg.grad = sens;
 end
-cfg.headmodel   = headmodel;
+
+cfg.headmodel          = headmodel;
 cfg.sourcemodel.pos    = dippar(:,1:3);
 cfg.sourcemodel.inside = 1:ndip;
 cfg.sourcemodel.unit   = headmodel.unit;
-cfg.reducerank  = 'no';
+cfg.reducerank         = 'no';
+
 if strcmp(headmodel.type, 'singleshell')
   cfg.singleshell.batchsize = 2500;
 end
-leadf           = ft_prepare_leadfield(cfg);
+
+leadf = ft_prepare_leadfield(cfg);
 
 if size(dippar, 2)==6
   leadfield = zeros(size(leadf.leadfield{1},1), ndip);

ni2_sensors.m

@@ -1,16 +1,23 @@
 function sens = ni2_sensors(varargin)
 
-% NI2_SENSORS create a fictive sensor-array.
+% NI2_SENSORS creates an EEG or MEG sensor array.
 %
 % Use as
-%  sens = ni2_sensors('type', senstype);
+%   sens = ni2_sensors('type', 'eeg')
+%   sens = ni2_sensors('type', 'meg')
+%   sens = ni2_sensors('type', 'ctf151')
+%   sens = ni2_sensors('type', 'ctf275')
+%   sens = ni2_sensors('type', 'bti248')
+%   sens = ni2_sensors('type', 'neuromag306')
+%   sens = ni2_sensors('type', 'eeg1020')
+%   sens = ni2_sensors('type', 'eeg1010')
 %
-% Where senstype can be any of the following:
-%  'eeg' generates a 91-channel eeg sensor array
-%  'meg' generates a 301-channel meg magnetometer sensor array
-%  'meg_grad_axial' generates a 301-channel meg axial gradiometer sensor
-%  array.
-
+% The default EEG sensor array consists of 91 electrodes over the upper half of a
+% sphere. The default MEG sensor array consists of 301 axial magnetometers over the
+% upper half of a sphere.
+%
+% Note that the units and the coordinate systems in which the sensors are described
+% can differ.
 
 type   = ft_getopt(varargin, 'type', 'eeg');
 jitter = ft_getopt(varargin, 'jitter', 0);

ni2_sourcemodel.m

@@ -1,37 +1,38 @@
 function sourcemodel = ni2_sourcemodel(varargin)
 
-% NI2_SOURCEMODEL generates a sourcemodel for forward and inverse modeling. 
+% NI2_SOURCEMODEL generates a sourcemodel for forward and inverse modeling.
 %
-% Use as:
-%  sourcemodel = ni2_sourcemodel('type', 'grid', 'resolution', res)
+% Use as
+%   sourcemodel = ni2_sourcemodel('type', 'mesh');
+%   sourcemodel = ni2_sourcemodel('type', 'grid', 'resolution', res)
 %
-% Where res is a scalar that specifies the spacing between the dipoles (in
-% cm). A regular 3-dimensional grid of dipoles is created.
+% The first call returns a source model that consists of many dipoles that form a
+% triangulated mesh that describes the cortical sheet.
+%
+% The second call creates a regular 3-dimensional grid of dipoles. The parameter res
+% is a scalar (in cm) that specifies the resolution, i.e., the spacing between the
+% dipoles.
 
-type       = ft_getopt(varargin, 'type',      'grid');
+type       = ft_getopt(varargin, 'type', 'grid');
 headmodel  = ft_getopt(varargin, 'headmodel', []);
 resolution = ft_getopt(varargin, 'resolution', 0.5); % in cm
 
 switch type
   case 'grid'
     if isempty(headmodel)
-      resolution = ft_getopt(varargin, 'resolution', 0.5);
       ax         = 0:resolution:9;
       ax         = [-fliplr(ax(2:end)) ax];
       [x,y,z]    = ndgrid(ax, ax, ax(ax>-1));
-      in         = false(size(x));
-      in(sqrt(x.^2+y.^2+z.^2)<=9) = true;
+      inside     = false(size(x));
+      inside(sqrt(x.^2+y.^2+z.^2)<=9) = true;
 
       pos        = [x(:) y(:) z(:)];
-      % inside     = find(in);
-      % outside    = find(in==0);
       clear x y z
 
       sourcemodel.pos     = pos;
-      % sourcemodel.inside  = inside;
-      % sourcemodel.outside = outside;
-      sourcemodel.inside  = in(:);
+      sourcemodel.inside  = inside(:);
       sourcemodel.dim     = [numel(ax) numel(ax) numel(ax(ax>-1))];
+      sourcemodel.unit    = 'cm';
 
     else
       headmodel = ft_convert_units(headmodel);
@@ -40,7 +41,7 @@
 
       cfg = [];
       cfg.headmodel  = headmodel;
-      cfg.resolution = resolution;
+      cfg.resolution = resolution; % in cm
       cfg.inwardshift = -1.*resolution;
       sourcemodel = ft_prepare_sourcemodel(cfg);
 

ni2_topomovie.m

@@ -1,18 +1,30 @@
-function ni2_topomovie(sens, dat, time, cfg)
+function ni2_topomovie(sens, data, time, cfg)
+
+% NI2_TOPOMOVIE makes a movie of the 2D-projection of the spatial EEG or MEG
+% topography as it changes over time.
+%
+% Use as
+%   ni2_topoplot(sens, data, time)
+%
+% Input arguments:
+%   sens = a sensor-array as obtained with NI2_SENSORS
+%   data = a NxT matrix representing a spatial topography over time. The number of rows (N) 
+%          should be equal to the number of sensors in the sens structure, the number of 
+%          columns should be equal to the number of time points.
+
+if nargin<4
+  cfg = [];
+end
 
 % ensure the sensor description to be according to what FieldTrip expects
 sens = ft_datatype_sens(sens);
 
-% construct a timelocked ERP data structure
+% construct a timelocked ERP data structure that FieldTrip can deal with
 timelock = [];
 timelock.time = time;
-timelock.avg = dat;
+timelock.avg = data;
 timelock.label = sens.label;
 
-if nargin<4
-  cfg = [];
-end
-
 if ~isfield(cfg, 'layout')
   % create layout
   if strncmp(sens.chantype{1}, 'meg', 3)