find_clusters() - Given a 2D matrix of t-scores, this function bundles
neighboring above threshold t-scores into clusters
and returns a 2D matrix of cluster assignments. For
use with cluster-based permutation tests.
Usage:
>>clust_membership=find_clusters(tscores,thresh,chan_hood,thresh_sign);
Required Inputs:
tscores - A 2D matrix of t-scores (channel x time point). Note,
t-scores SHOULD be signed so that negative and positive
deviations from the null hypothesis are NOT clustered
together.
thresh - The thershold for cluster inclusion (i.e., only t-scores
more extreme than the threshold will be included in clusters).
If thresh is positive, clusters of positive t-scores will
be returned. Otherwise clusters of negative t-scores will
be returned.
chan_hood - A symmetric 2d matrix indicating which channels are
neighbors with other channels. If chan_hood(a,b)=1,
then Channel A and B are neighbors. This is produced by the
function spatial_neighbors.m.
thresh_sign - If greater than zero, t-scores greater than thresh will be
included in clusters. Otherwise, t-scores less than thresh
will be included in clusters.
Outputs:
clust_membership - A 2D matrix (channel x time point) indicating which
cluster each channel/time point pair belongs to. For
example, if clust_membership(2,10)=3, then the t-score
at the 2nd channel and 10th time point belongs to the
3rd cluster. 0 indicates that the t-score was not
included in any cluster.
n_clust - The number of clusters found.
Notes:
-A global variable, clust_ids, is created to keep track of which cluster
each channel/time point pair belongs to. This speeds things up because
of the subfunction this variable has to be passed to.
Authors:
Amy Guthormsen
Applied Modern Physics Group
Los Alamos National Laboratory
David M. Groppe
Kutaslab
University of California, San Diego0001 % find_clusters() - Given a 2D matrix of t-scores, this function bundles 0002 % neighboring above threshold t-scores into clusters 0003 % and returns a 2D matrix of cluster assignments. For 0004 % use with cluster-based permutation tests. 0005 % 0006 % Usage: 0007 % >>clust_membership=find_clusters(tscores,thresh,chan_hood,thresh_sign); 0008 % 0009 % Required Inputs: 0010 % tscores - A 2D matrix of t-scores (channel x time point). Note, 0011 % t-scores SHOULD be signed so that negative and positive 0012 % deviations from the null hypothesis are NOT clustered 0013 % together. 0014 % thresh - The thershold for cluster inclusion (i.e., only t-scores 0015 % more extreme than the threshold will be included in clusters). 0016 % If thresh is positive, clusters of positive t-scores will 0017 % be returned. Otherwise clusters of negative t-scores will 0018 % be returned. 0019 % chan_hood - A symmetric 2d matrix indicating which channels are 0020 % neighbors with other channels. If chan_hood(a,b)=1, 0021 % then Channel A and B are neighbors. This is produced by the 0022 % function spatial_neighbors.m. 0023 % thresh_sign - If greater than zero, t-scores greater than thresh will be 0024 % included in clusters. Otherwise, t-scores less than thresh 0025 % will be included in clusters. 0026 % 0027 % 0028 % Outputs: 0029 % clust_membership - A 2D matrix (channel x time point) indicating which 0030 % cluster each channel/time point pair belongs to. For 0031 % example, if clust_membership(2,10)=3, then the t-score 0032 % at the 2nd channel and 10th time point belongs to the 0033 % 3rd cluster. 0 indicates that the t-score was not 0034 % included in any cluster. 0035 % n_clust - The number of clusters found. 0036 % 0037 % 0038 % Notes: 0039 % -A global variable, clust_ids, is created to keep track of which cluster 0040 % each channel/time point pair belongs to. This speeds things up because 0041 % of the subfunction this variable has to be passed to. 0042 % 0043 % Authors: 0044 % Amy Guthormsen 0045 % Applied Modern Physics Group 0046 % Los Alamos National Laboratory 0047 % 0048 % David M. Groppe 0049 % Kutaslab 0050 % University of California, San Diego 0051 0052 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%% 0053 % 12/11/2011-Made reursiveless by Amy Guthormsen 0054 0055 function [clust_membership, n_clust]=find_clusters(tscores,thresh,chan_hood,thresh_sign) 0056 0057 [n_chan, n_tpt]=size(tscores); 0058 0059 global clust_ids; % channel x time point matrix indicated which cluster each "voxel" belongs to 0060 0061 clust_ids=zeros(n_chan,n_tpt); 0062 0063 time_matrix=repmat(1:n_tpt,n_chan,1); 0064 chan_matrix=repmat([1:n_chan]',1,n_tpt); 0065 0066 if thresh_sign>0 0067 %looking for positive clusters 0068 above_thresh_ids=find(tscores>=thresh); 0069 else 0070 %looking for negative clusters 0071 above_thresh_ids=find(tscores<=thresh); 0072 end 0073 above_thresh_times=time_matrix(above_thresh_ids); 0074 above_thresh_chans=chan_matrix(above_thresh_ids); 0075 0076 %Clear chan & time matrix to save memory (just in case) 0077 clear chan_matrix time_matrix 0078 0079 n_above=length(above_thresh_ids); 0080 0081 n_clust=0; 0082 for a=1:n_above, 0083 0084 voxel_id=above_thresh_ids(a); 0085 0086 if ~clust_ids(voxel_id) 0087 %this "voxel" isn't in a cluster yet 0088 n_clust=n_clust+1; 0089 0090 clust_ids(voxel_id)=n_clust; 0091 0092 %go through all the remaining voxels and find all the above 0093 %threshold voxels this voxel is neighbors with and give them this 0094 %cluster # 0095 voxels_not_checked = ones(size(above_thresh_ids)); 0096 check_me = zeros(size(above_thresh_ids)); 0097 0098 check_me(a) = 1; 0099 0100 while sum(check_me>0) 0101 first = find(check_me,1); 0102 new = follow_clust(n_above,first,n_clust,above_thresh_ids,above_thresh_times,above_thresh_chans,chan_hood,n_chan); 0103 check_me(new)=1; 0104 voxels_not_checked(first)=0; 0105 check_me = check_me&voxels_not_checked; 0106 end 0107 end 0108 end 0109 0110 clust_membership=clust_ids; 0111 clear global clust_ids 0112 0113 %% End of Main Function 0114 function [new_members] = follow_clust(n_above,current_voxel_id,current_clust_num,above_thresh_ids,above_thresh_times,above_thresh_chans,chan_hood,n_chan) 0115 %Function for finding all the members of cluster based on 0116 %single "voxel" seed. If it finds new members of a cluster, it indicates 0117 %which cluster they are a member of in the global variable clust_ids. 0118 % 0119 % Inputs: 0120 % clust_ids - Matrix indicating which cluster each 0121 % voxel belongs to (if any) 0122 % n_above - The total number of above threshold voxels 0123 % current_voxel_id - The index of the seed voxel into 0124 % above_thresh_ids, above_thresh_chans, above_thresh_times 0125 % current_clust_num - The # of the cluster the seed voxel belongs to 0126 % above_thresh_ids - A vector of indices of all the above threshold 0127 % voxels into the original 2d data matrix (channel x 0128 % time point) 0129 % above_thresh_times - The time points corresponding to above_thresh_ids 0130 % above_thresh_chans - The time points corresponding to above_thresh_chans 0131 % chan_hood - A symmetric 2d matrix indicating which channels are 0132 % neighbors with other channels. If chan_hood(a,b)=1, 0133 % then Channel A and B are neighbors. 0134 % n_chan - The total number of channels 0135 % 0136 % 0137 % Output: 0138 % new_members - above_thresh_ids indices to new members of the cluster 0139 0140 % global clust_ids; % channel x time point matrix indicated which cluster each "voxel" belongs to 0141 0142 new_members=zeros(1,n_chan*3); %pre-allocate memory 0143 new_members_ct=0; 0144 0145 for b=current_clust_num:n_above, 0146 if ~clust_ids(above_thresh_ids(b)) 0147 temp_dist=abs(above_thresh_times(b)-above_thresh_times(current_voxel_id)); 0148 if above_thresh_chans(current_voxel_id)==above_thresh_chans(b) 0149 %voxels are at same channel 0150 chan_dist=0; 0151 elseif chan_hood(above_thresh_chans(current_voxel_id),above_thresh_chans(b)) 0152 %channels are neighbors 0153 chan_dist=1; 0154 else 0155 %voxels aren't spatially compatible 0156 chan_dist=2; 0157 end 0158 if (temp_dist+chan_dist)<=1, 0159 %if voxels are at same time point and neighboring channels OR 0160 %if voxels are at same channel and neighboring time points, 0161 %merge them into the same cluster 0162 0163 clust_ids(above_thresh_ids(b))=current_clust_num; 0164 %keep track of which other voxels are joined to this 0165 %cluster 0166 new_members_ct=new_members_ct+1; 0167 new_members(new_members_ct)=b; 0168 end 0169 end 0170 end 0171 0172 new_members=nonzeros(new_members); 0173 end 0174 end