Home > matlabmk > find_clustersRECURSIVE.m

find_clustersRECURSIVE

PURPOSE ^

find_clusters() - Given a 2D matrix of t-scores, this function bundles

SYNOPSIS ^

function [clust_membership, n_clust, used_rec_lim]=find_clusters(tscores,thresh,chan_hood,thresh_sign,recursion_limit)

DESCRIPTION ^

 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/2011

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % 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

Generated on Tue 10-May-2016 16:37:45 by m2html © 2005