Print

---

Following up on this thread, in case people find it useful: here is some code to find MEG sensors nearest to an MNI region of interest 


%% jg_find_sensors_nearest_mni_coords.m
%
% (thanks to vladimir litvak for showing 
%  where to get the coordinate info from; 
%  ...c.f. our discussion on SPM e-mail list)
%
% This works for Elekta neuromag data
% 
% Note that nearest sensors are only relevant/useful
% in this context if they are gradiometers. 
% 
% Data needs to be coregistered to MNI space
%
% Abbreviations: 'ms' = mni space; 'ss' = sensor space
%
%                                 JG 02/052012
%


%% Specify params etc:


clear all close all


addpath('/home/jg03/code/spm8_fil_devel/')
spm eeg


roi_coords =             [  -42, -22,  7  ;... % LA1  
   	                            -61, -32,  8  ;... % LSTG
                                     -8,  -93, -6  ;... % V2 BA18 L  
                                     -9,  -96,  2  ;... % V1 BA17 L 
			             46, -14,  8  ;... % RA1                         
			             59, -25,  8  ;... % RSTG
                                     15, -84, -7  ;... % V2 BA18 R  
                                     15, -92,  4  ];   % V1 BA17 R    
                                                         
file = <filename> 
n_nearest_sensors = 5;


%% Do it. 

D = spm_eeg_load(file);
             
ms_verts = D.inv{1}.mesh.tess_mni.vert;

ss_verts = D.inv{1}.forward.mesh.vert;
ss_coilpos = D.inv{1}.datareg.sensors.coilpos; 
ss_chanpos = D.inv{1}.datareg.sensors.chanpos; 
 
for r = 1:size(roi_coords,1)
 
  % List vertices in ascending order of Euclidean distance from ROI
  ms_eds = sqrt(sum((repmat(roi_coords(r,:)',1,length(ms_verts))-ms_verts').^2));
  [ms_eds_sorted{r}, ms_eds_sorted_idx{r}] = sort(ms_eds, 'ascend');
 
  % Find the n nearest sensors to the vertex nearest to roi r
  ss_eds = sqrt(sum((repmat(ss_verts(ms_eds_sorted_idx{r}(1),:)',1,length(ss_chanpos))-ss_chanpos').^2));
  [ss_eds_sorted{r}, ss_eds_sorted_idx{r}] = sort(ss_eds, 'ascend');

  % add to list of nearest channels
   si = ss_eds_sorted_idx{r}(1:n_nearest_sensors);
   nearest_sensors_inds(1:n_nearest_sensors,r) = si; % ss_eds_sorted_idx{r}(1:n_nearest_sensors);
   nearest_sensors_coords(1:n_nearest_sensors,:,r) = ss_chanpos(si,:); %ss_eds_sorted_idx{r (1:n_nearest_sensors),:);
  

end

% Get the channel types of the identified sensors. 
nearest_sensors_chantypes = reshape(D.chantype(nearest_sensors_inds(:)), size(nearest_sensors_inds));


%% Visualize results to check:

% Vertex closest tot he mnui roi is a red cross+circle
% Sensor(s) nearest to this vertex are black cross+circles

% Find the n nearest sensors to the vertex nearest to roi r
dm = 2; % decimation number for brain (so can see the roi inside) 
for r = 1:size(roi_coords,1)

   fig = figure; 
  
   % brain (decimated)
   scatter3(ss_verts(1:dm:end,1), ss_verts(1:dm:end,2), ss_verts(1:dm:end,3), 'cyan.');
   hold on; % (if you put hold on before scatter3 it doesn't add the grid lines...)
  
   % sensors
   plot3(ss_chanpos(:,1), ss_chanpos(:,2), ss_chanpos(:,3), 'go')
  
   % roi vertex
   plot3(ss_verts(ms_eds_sorted_idx{r}(1),1),ss_verts(ms_eds_sorted_idx{r}(1),2), ss_verts(ms_eds_sorted_idx{r}(1),3), 'rX', 'Markersize', 30); 
   plot3(ss_verts(ms_eds_sorted_idx{r}(1),1),ss_verts(ms_eds_sorted_idx{r}(1),2), ss_verts(ms_eds_sorted_idx{r}(1),3), 'ro', 'Markersize', 30);   % (add the 'o' for clarity)       

   % n nearest sensors
   plot3(nearest_sensors_coords(:,1,r), nearest_sensors_coords(:,2,r), nearest_sensors_coords(:,3,r), 'kX', 'MarkerSize', 30);
   plot3(nearest_sensors_coords(:,1,r), nearest_sensors_coords(:,2,r), nearest_sensors_coords(:,3,r), 'ko', 'MarkerSize', 30);
  
   title(sprintf('Sensors nearest to ROI %d, MNI coords: %d %d %d',...
        r, roi_coords(r,1), roi_coords(r,2), roi_coords(r,3)),...
        'FontWeight', 'Bold');
       
end