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,recursion_limit);
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.
Optional Inputs:
recursion_limit - The number of recursive embeddings you want MATLAB to
allow. MATLAB limits this and, since clusters are
found using recursion, it is possible that you will
exceed MATLAB's limit with large clusters. The
default behavior is to estimate the largest number of
embeddings, set this as the limit, and then re-set
the original limit after the clustering is done.
However, if this doesn't work, use a higher number
(MATLAB's default is 500 in our lab).
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.
used_rec_lim - The recursion limit used by the function. This may
differ from the specified recursion limit since this
function will automatically double the recursion
limit whenever it is too small. This output is
useful when you need to re-run this function on a
similarly sized data set (e.g., when one is
performing a cluster-based permutation test) since
having to automatically increase the recursion limit
wastes time.
Notes:
-A global variable, clust_ids, is created to keep track of which cluster
each channel/time point pair belongs to. This is done to minimize memory
demands during the recusion. It is erased after all the clusters are
established.
Author:
David Groppe
Kutaslab, 5/20110001 % 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,recursion_limit); 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 % Optional Inputs: 0028 % recursion_limit - The number of recursive embeddings you want MATLAB to 0029 % allow. MATLAB limits this and, since clusters are 0030 % found using recursion, it is possible that you will 0031 % exceed MATLAB's limit with large clusters. The 0032 % default behavior is to estimate the largest number of 0033 % embeddings, set this as the limit, and then re-set 0034 % the original limit after the clustering is done. 0035 % However, if this doesn't work, use a higher number 0036 % (MATLAB's default is 500 in our lab). 0037 % 0038 % Outputs: 0039 % clust_membership - A 2D matrix (channel x time point) indicating which 0040 % cluster each channel/time point pair belongs to. For 0041 % example, if clust_membership(2,10)=3, then the t-score 0042 % at the 2nd channel and 10th time point belongs to the 0043 % 3rd cluster. 0 indicates that the t-score was not 0044 % included in any cluster. 0045 % n_clust - The number of clusters found. 0046 % used_rec_lim - The recursion limit used by the function. This may 0047 % differ from the specified recursion limit since this 0048 % function will automatically double the recursion 0049 % limit whenever it is too small. This output is 0050 % useful when you need to re-run this function on a 0051 % similarly sized data set (e.g., when one is 0052 % performing a cluster-based permutation test) since 0053 % having to automatically increase the recursion limit 0054 % wastes time. 0055 % 0056 % Notes: 0057 % -A global variable, clust_ids, is created to keep track of which cluster 0058 % each channel/time point pair belongs to. This is done to minimize memory 0059 % demands during the recusion. It is erased after all the clusters are 0060 % established. 0061 % 0062 % Author: 0063 % David Groppe 0064 % Kutaslab, 5/2011 0065 0066 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%% 0067 % ?/?/??- 0068 0069 function [clust_membership, n_clust, used_rec_lim]=find_clusters(tscores,thresh,chan_hood,thresh_sign,recursion_limit) 0070 0071 [n_chan, n_tpt]=size(tscores); 0072 0073 global clust_ids; % channel x time point matrix indicated which cluster each "voxel" belongs to 0074 0075 clust_ids=zeros(n_chan,n_tpt); 0076 0077 time_matrix=repmat(1:n_tpt,n_chan,1); 0078 chan_matrix=repmat([1:n_chan]',1,n_tpt); 0079 0080 if thresh_sign>0 0081 %looking for positive clusters 0082 above_thresh_ids=find(tscores>=thresh); 0083 else 0084 %looking for negative clusters 0085 above_thresh_ids=find(tscores<=thresh); 0086 end 0087 above_thresh_times=time_matrix(above_thresh_ids); 0088 above_thresh_chans=chan_matrix(above_thresh_ids); 0089 0090 %Clear chan & time matrix to save memory (just in case) 0091 clear chan_matrix time_matrix 0092 0093 n_above=length(above_thresh_ids); 0094 0095 orig_limit=get(0,'RecursionLimit'); 0096 if nargin<5 || isempty(recursion_limit), 0097 if orig_limit<n_above, 0098 set(0,'RecursionLimit',n_above); 0099 end 0100 else 0101 set(0,'RecursionLimit',recursion_limit); 0102 end 0103 0104 n_clust=0; 0105 for a=1:n_above, 0106 voxel_id=above_thresh_ids(a); 0107 if ~clust_ids(voxel_id) 0108 %this "voxel" isn't in a cluster yet 0109 n_clust=n_clust+1; 0110 clust_ids(voxel_id)=n_clust; 0111 0112 %go through all the remaining voxels and find all the above 0113 %threshold voxels this voxel is neighbors with and give them this 0114 %cluster # 0115 did_it=0; 0116 clust_ids_copy=clust_ids; 0117 while ~did_it, 0118 try 0119 follow_clust(a+1,n_above,a,n_clust,above_thresh_ids,above_thresh_times,above_thresh_chans,chan_hood,n_chan); 0120 did_it=1; 0121 catch me 0122 did_it=0; 0123 clust_ids=clust_ids_copy; %reset clust_ids to what it was before the attempt 0124 if strcmpi(me.identifier,'MATLAB:recursionLimit') 0125 current_limit=get(0,'RecursionLimit'); 0126 %fprintf('Doubling RecursionLimit to %d\n',current_limit*2); 0127 set(0,'RecursionLimit',current_limit*2); 0128 end 0129 end 0130 end 0131 end 0132 end 0133 0134 %reset RecursionLimit to its original value 0135 used_rec_lim=get(0,'RecursionLimit'); 0136 set(0,'RecursionLimit',orig_limit); 0137 clust_membership=clust_ids; 0138 clear global clust_ids 0139 0140 %% End of Main Function 0141 0142 function follow_clust(start_voxel,n_above,current_voxel_id,current_clust_num,above_thresh_ids,above_thresh_times,above_thresh_chans,chan_hood,n_chan) 0143 %Function for recursively finding all the members of cluster based on 0144 %single "voxel" seed 0145 % 0146 % start_voxel - All voxels above this might belong to the cluster 0147 % (i.e., all voxels before this have already been 0148 % assigned to clusters) 0149 % n_above - The total number of above threshold voxels 0150 % current_voxel_id - The index of the seed voxel into 0151 % above_thresh_ids, above_thresh_chans, above_thresh_times 0152 % current_clust_num - The # of the cluster the seed voxel belongs to 0153 % above_thresh_ids - A vector of indices of all the above threshold 0154 % voxels into the original 2d data matrix (channel x 0155 % time point) 0156 % above_thresh_times - The time points corresponding to above_thresh_ids 0157 % above_thresh_chans - The time points corresponding to above_thresh_chans 0158 % chan_hood - A symmetric 2d matrix indicating which channels are 0159 % neighbors with other channels. If chan_hood(a,b)=1, 0160 % then Channel A and B are neighbors. 0161 % n_chan - The total number of channels 0162 % 0163 0164 global clust_ids; % channel x time point matrix indicated which cluster each "voxel" belongs to 0165 0166 new_members=zeros(1,n_chan*3); %pre-allocate memory 0167 new_members_ct=0; 0168 for b=start_voxel:n_above, 0169 if ~clust_ids(above_thresh_ids(b)) 0170 temp_dist=abs(above_thresh_times(b)-above_thresh_times(current_voxel_id)); 0171 if above_thresh_chans(current_voxel_id)==above_thresh_chans(b) 0172 %voxels are at same channel 0173 chan_dist=0; 0174 elseif chan_hood(above_thresh_chans(current_voxel_id),above_thresh_chans(b)) 0175 %channels are neighbors 0176 chan_dist=1; 0177 else 0178 %voxels aren't spatially compatible 0179 chan_dist=2; 0180 end 0181 if (temp_dist+chan_dist)<=1, 0182 %if voxels are at same time point and neighboring channels OR 0183 %if voxels are at same channel and neighboring time points, 0184 %merge them into the same cluster 0185 0186 clust_ids(above_thresh_ids(b))=current_clust_num; 0187 %keep track of which other voxels are joined to this 0188 %cluster 0189 new_members_ct=new_members_ct+1; 0190 new_members(new_members_ct)=b; 0191 end 0192 end 0193 end 0194 0195 %recursively search for above threshold neighbors of the other 0196 %voxels that were just joined to this cluster 0197 for c=1:new_members_ct, 0198 follow_clust(start_voxel,n_above,new_members(c),current_clust_num,above_thresh_ids,above_thresh_times,above_thresh_chans,chan_hood,n_chan); 0199 end 0200