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).

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 by changing the base file names and write file names. The label files generated above (using mri_annotation2label) currently need to be copied into a lh and rh directory within a parent dir. Though, a simple modification would allow these files to remain in your fsaverage or other subject's dir.

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