Working with ROIs (Freesurfer)
The .label files created during autorecon3, or by saving an overlay are plaintext lists of vertices. It should be straightforward to find the intersection or union of the vertices in two or more .label files using common Unix shell commands.
.label Syntax
Below is the first few lines of the lh.banksts.label file for a participant. This file was produced during the Lausanne parcellation procedure.
#!ascii label , from subject FS_0043 vox2ras=TkReg 1457 32992 -47.867 -62.891 33.007 0.0000000000 32993 -48.881 -62.907 32.963 0.0000000000 34014 -47.298 -61.696 33.538 0.0000000000 34015 -47.972 -61.259 33.648 0.0000000000 34024 -47.634 -62.264 33.251 0.0000000000 34025 -48.286 -61.775 33.303 0.0000000000 ...etc
The first line is a header line. The second line indicates the number of vertices in the label file. The remaining lines indicate the vertex number, x,y,z coordinates of the vertex, and I have no idea what that last column indicates but it's not terribly important for our purposes.
Preparation of .label Files
The group analysis used the fsaverage surface, and since mris_divide_parcellation (which you may also be using, or use later) assumes we're subdividing a particular subject's annotation, it makes sense to copy the source .annot files into the fsaverage label/ folder (giving the appropriate ?h prefix)
cp cache.th30.abs.sig.ocn.annot \ $SUBJECTS_DIR\fsaverage\label\rh.cache.th30.abs.sig.ocn.annot
The steps described below do not work on .annot files because they are not plaintext and easily manipulated. However, mri_annotation2label will convert a .annot file to a series of .label files.
e.g. mri_annotation2label --annotation cache.th30.abs.sig.ocn --hemi rh --subject fsaverage --surf pial --labelbase rh.T1
Set Operations on .label Files
Intersecting .label Files (Intersection)
Given what we know about the contents of a label file, the unix comm command should suffice to find the common entries between two label files, which is the set intersection of vertices appearing in each.
https://mail.nmr.mgh.harvard.edu/pipermail//freesurfer/2009-August/011743.html
Merging .label Files (Union)
A union operation will find all vertices appearing in either .label file. You will want to make sure that only unique values are retained (i.e., you don't want to list the same vertex multiple times in the output .label file).
Or check out:
mri_mergelabels -i label1 -i label2 ... -o outputlabel
This is a simple script that will merge two or more label files. It does this by catting the label files together (after removing the first two lines). It inserts a new header (the first two lines). The number of entries in the new file (ie, the number on the second line), is computed by summing those from the input files.
ROI 'VennDiagram' Matlab Function
ROIVennDiagram.m uses intersect and setdiff functions to extract intersecting and unique vertices. Variable names suggest that it's used for T1 and T2 data, but you can use it for any paired data. The label files generated above (using mri_annotation2label) need to be copied into a lh and rh directory within a parent dir. A simple modification would allow these files to remain in your fsaverage or other subject's dir, if you wanted to clutter those up.
function ROIvenndiagram = ROIvenndiagram(hemi)
%% function ROIvenndiagram = ROIvenndiagram('hemi');
%Pull out intersecting and different vertices for FS label files,
%and store in one of 3 text files (T1, T1T2, T2)
%
%Currently takes one string input indicating lh or rh.
%
%You'll need to adjust the file names on lines 42 and 43 to match your
%file names.
%
%Note - The headers in the output text files have quotes around them, as
%the headers begin with matlab unfriendly characters (#). So, you'll need to
%go into those 3 output text files and delete the quotes.
%
%Intended for use in parent dir containing lh and rh dirs, which
%contain label files.
% _____________________________________
% G.J. Smith
% gjsmith4@buffalo.edu
% _____________________________________
parentdir=pwd;
hemidir = sprintf('%s/%s', pwd, hemi);
cd(hemidir);
%% Pre-initialize array for cating labels
T1Agg = cell(1,5);
T2Agg = cell(1,5);
% Begin cycle through label files
for labeln = 1:3000 %arbirary number higher than any label number
%% Apply appropriate number of 0's before number - eg. 001, 010, 100
numcorrection=num2str(labeln);
if labeln < 10
labelnum=strcat('00',numcorrection);
elseif labeln < 100
labelnum=strcat('0',numcorrection);
else
labelnum=strcat(numcorrection);
end
%% T1 and T2 label file names -- Edit these to match your file naming convention
T1_filename = sprintf('%s.T1-%s.label', hemi, labelnum);
T2_filename = sprintf('%s.T2-%s.label', hemi, labelnum);
%% Build one large aggregate file with all vertices from all label files
header = 2; % lines 1 and 2 of label file stored seperately
delimiter = ' ';
if exist([pwd filesep T1_filename], 'file')
%load data
T1file = importdata(T1_filename,delimiter,header);
%Store in one array
T1Agg = vertcat(T1Agg, num2cell(T1file.data));
end
if exist([pwd filesep T2_filename], 'file')
%load data
T2file = importdata(T2_filename,delimiter,header);
%Store in one array
T2Agg = vertcat(T2Agg, num2cell(T2file.data));
end
end
%% Check that label files were found
exist T1file;
if ans == 1
display('Labels found, pulling out intersects and setdiffs');
else
display('No label files found, check your dir or file names');
return;
end
%% Find intersection and differences
T1Agg = cell2mat(T1Agg);
T2Agg = cell2mat(T2Agg);
%compare
%union = union(T1file.data(:,1), T2file.data(:,1));
intersection = intersect(T1Agg, T2Agg, 'rows'); %'rows' works with same # of columns
T1_diffs = setdiff(T1Agg, T2Agg, 'rows');
T2_diffs = setdiff(T2Agg, T1Agg, 'rows');
%% Build new label files
%T1
T1_differences = cell(50000, 5); %hard coded arbitrarily high number of rows for pre-init
T1_differences(1,1) = T1file.textdata(1,1); % header
T1_differences{2,1} = length(T1_diffs); % header line 2, new number of rows
T1_differences(3:length(T1_diffs) + 2,:) = num2cell(T1_diffs); % data
%T2
T2_differences = cell(50000, 5); %hard coded arbitrarily high number of rows for pre-init
T2_differences(1,1) = T2file.textdata(1,1); % header
T2_differences{2,1} = length(T2_diffs); % header line 2, new number of rows
T2_differences(3:length(T2_diffs) + 2,:) = num2cell(T2_diffs); % data
%Match
Match_labelfile = cell(50000, 5); %hard coded arbitrarily high number of rows for pre-init
Match_labelfile(1,1) = T1file.textdata(1,1); % header T1 and T2 for same subject is the same
Match_labelfile{2,1} = length(intersection); % header line 2, new number of rows
Match_labelfile(3:length(intersection) + 2, :) = num2cell(intersection); % data
%% Remove blank rows from oversized pre-alo arrays
T1_differences(all(cellfun(@isempty,T1_differences),2),:) = [];
T2_differences(all(cellfun(@isempty,T2_differences),2),:) = [];
Match_labelfile(all(cellfun(@isempty,Match_labelfile),2),:) = [];
%% Write files
cd(parentdir);
T1_fname = sprintf('%s.T1.label.txt', hemi);
writetable(cell2table(T1_differences), T1_fname,'Delimiter',' ','WriteVariableNames',false);
%T2
T2_fname = sprintf('%s.T2.label.txt', hemi);
writetable(cell2table(T2_differences), T2_fname,'Delimiter',' ','WriteVariableNames',false);
%Match
Match_fname = sprintf('%s.T1T2.label.txt', hemi);
writetable(cell2table(Match_labelfile), Match_fname,'Delimiter',' ','WriteVariableNames',false);
end
Visualize and Cluster (Intersected Label File)
tksurfer fsaverage ?h pial