ENH - updates and fixes

Author: schoffelen

ni2_headmodel.m

@@ -22,12 +22,14 @@
         headmodel.r    = 10;
         headmodel.c    = 1;
         headmodel.unit = 'cm';
+        headmodel.type = 'singlesphere';
         headmodel      = ft_datatype_headmodel(headmodel);
       case 3
         headmodel.o    = [0 0 0];
         headmodel.r    = [0.88 0.92 1.00]*10;
         headmodel.c    = [1 1/80 1];
         headmodel.unit = 'cm';
+        headmodel.type = 'concentricspheres';
         headmodel      = ft_datatype_headmodel(headmodel);
       otherwise
         error('number of spheres other than 1 or 3 is not supported');

ni2_leadfield.m

@@ -29,9 +29,9 @@
   cfg.grad = sens;
 end
 cfg.headmodel   = headmodel;
-cfg.grid.pos    = dippar(:,1:3);
-cfg.grid.inside = 1:ndip;
-cfg.grid.unit   = headmodel.unit;
+cfg.sourcemodel.pos    = dippar(:,1:3);
+cfg.sourcemodel.inside = 1:ndip;
+cfg.sourcemodel.unit   = headmodel.unit;
 cfg.reducerank  = 'no';
 if strcmp(headmodel.type, 'singleshell')
   cfg.singleshell.batchsize = 2500;

ni2_sensors.m

@@ -12,15 +12,24 @@
 %  array.
 
 
-type = ft_getopt(varargin, 'type', 'eeg');
+type   = ft_getopt(varargin, 'type', 'eeg');
+jitter = ft_getopt(varargin, 'jitter', 0);
+n      = ft_getopt(varargin, 'n', 162); % number of vertices for the sphere, determines the number of electrodes, this is the old default
 
 switch type
   case 'eeg'
     % create an eeg electrode array
-    [chanpos, tri] = icosahedron162;
+    [chanpos, tri] = mesh_sphere(n);
     chanpos        = chanpos*10;
     chanpos(chanpos(:,3)<0,:) = [];
 
+    if jitter
+      [th,phi,r] = cart2sph(chanpos(:,1), chanpos(:,2), chanpos(:,3));
+      shift1 = 2*jitter*(rand(numel(th),1)  - 1);
+      shift2 = 2*jitter*(rand(numel(phi),1) - 1);
+      [chanpos(:,1),chanpos(:,2),chanpos(:,3)] = sph2cart(th+shift1(:),phi+shift2(:),r);
+    end
+      
     sens.chanpos = chanpos;
     sens.elecpos = chanpos;
     for k = 1:size(chanpos,1)
@@ -29,6 +38,7 @@
     end
     sens = ft_datatype_sens(sens);
 
+    
   case {'meg_mag' 'meg'}
     % create an meg magnetometer array
     [chanpos, tri] = icosahedron642;
@@ -78,7 +88,7 @@
     % not implemented yet
   case 'meg_ctf275'
     load('ctf275');
-    case 'meg_ctf151'
+  case 'meg_ctf151'
     load('ctf151');
   case 'meg_bti248'
     load('bti248');

not_in_zip/ni2_dipolefitting_leadfielderror.m

@@ -0,0 +1,88 @@
+% create some data
+[S,time] = ni2_activation;
+sens        = ni2_sensors('type','eeg','n',64);
+
+% in order to fool around with the reference of the leadfield, we need to
+% add an explicit 'tra', otherwise ft_compute_leadfield will subtract the
+% mean by default. By adding a tra as Identity matrix, the leadfields will
+% be returned 'without' reference. Here we add a reference on the left, 'eeg74'
+ref            = 25;
+sens.tra       = eye(numel(sens.label));
+sens.tra(:,ref) = sens.tra(:,ref)-1;
+
+% create a volume conductor model
+headmodel   = ni2_headmodel('type', 'spherical', 'nshell', 3);
+
+% create a 3D grid for the initial gridsearch
+sourcemodel = ni2_sourcemodel('type','grid','resolution',1);
+
+% simulate some sensor level data
+dippar      = [-8./sqrt(2) 0.5 8./sqrt(2) [1 1 1].*(3^(-1/3))];
+leadfield   = ni2_leadfield(sens, headmodel, dippar);
+noise       = sens.tra*randn(numel(sens.label),1000)*1e-3;
+sensordata  = leadfield*S+noise;
+assert(all(sensordata(ref,:)==0));
+
+% visualize it
+topo_observed = sensordata(:,500);
+ni2_topoplot(sens, topo_observed);
+
+data        = []; 
+data.avg    = sensordata; 
+data.time   = time; 
+data.label  = sens.label; 
+data.elec   = sens; 
+data.dimord ='chan_time';
+
+% this is what I think Sarang means, i.e. data with an implicitref present,
+% and then svd'ed in order to remove the last component
+[u,s,v]  = svd(data.avg, 'econ');
+montage2.tra = u(:,1:end-1)';
+montage2.labelold = data.label;
+for k = 1:size(montage2.tra,1)
+  montage2.labelnew{k,1} = sprintf('comp%02d',k);
+end
+data_svd      = ft_apply_montage(data, montage2);
+data_svd.elec = ft_apply_montage(data.elec, montage2);
+
+% this is the data with the reference channel removed
+cfg         = [];
+cfg.channel = data.label;
+cfg.channel(ref) = [];
+data        = ft_selectdata(cfg, data);
+
+% do some referencing by hand to obtain the average referenced data
+montage1.tra  = eye(numel(data.label))- ones(numel(data.label))./numel(data.label);
+montage1.labelold = data.label;
+montage1.labelnew = data.label;
+data_avg      = ft_apply_montage(data, montage1);
+data_avg.elec = ft_apply_montage(data.elec, montage1);
+
+cfg            = []; 
+cfg.gridsearch = 'yes'; 
+cfg.model      = 'regional'; 
+cfg.headmodel  = headmodel; 
+cfg.sourcemodel = sourcemodel;
+cfg.latency    = [0.49 0.51]; 
+cfg.nonlinear  = 'yes'; 
+cfg.numdipoles = 1;
+dip            = ft_dipolefitting(cfg, data);
+dip_avg        = ft_dipolefitting(cfg, data_avg);
+dip_svd        = ft_dipolefitting(cfg, data_svd);
+
+%ni2_topoplot(sens,dip.Vmodel(:,6));
+%ni2_topoplot(sens,dip.Vdata(:,6));
+
+sens_jittered = ni2_sensors('type', 'eeg', 'jitter', 0.1, 'n', 64);
+
+dataj = data;
+dataj_avg = data_avg;
+dataj_svd = data_svd;
+
+dataj.elec.elecpos = sens_jittered.chanpos;
+dipj = ft_dipolefitting(cfg, dataj);
+dataj_avg.elec.elecpos = sens_jittered.chanpos;
+dipj_avg = ft_dipolefitting(cfg, dataj_avg);
+dataj_svd.elec.elecpos = sens_jittered.chanpos;
+dipj_svd = ft_dipolefitting(cfg, dataj_svd);
+