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
|