fmri_clust_filt

 [Z,XYZ,C] = clust_filt(Y,k);

 Finds clusters of k contiguous voxels of non-zero data in the data volume specified by Y.
 The function is based on SPM5 code (spm_getSPM.m) and calls spm_clusters().

0001 function [Y,XYZ,C] = fmri_clust_filt(Y,k)
0002 % [Z,XYZ,C] = clust_filt(Y,k);
0003 %
0004 % Finds clusters of k contiguous voxels of non-zero data in the data volume specified by Y.
0005 % The function is based on SPM5 code (spm_getSPM.m) and calls spm_clusters().
0006 %
0007 
0008 % 09/14/07 Petr Janata
0009 % 2010.08.20 FB - copied from the autobio/fmri project directory into the
0010 %   private database fmri folder, added "fmri_" prefix to the name to match
0011 %   the convention for fmri analysis helper functions
0012 
0013 XYZ = {};
0014 
0015 % Find all non-NaN or non-zero voxels in the volume
0016 if any(isnan(Y))
0017   good_vox_idxs = find(~isnan(Y));
0018   multiplier = NaN;
0019 else
0020   good_vox_idxs = find(Y);
0021   multiplier = 0;
0022 end
0023 
0024 fprintf('Finding clusters of size >=%d in %d good voxels\n', k, length(good_vox_idxs));
0025 
0026 % Create an XYZ list
0027 [XYZ{1:3}] = ind2sub(size(Y),good_vox_idxs);
0028 XYZ = cat(2,XYZ{:})';
0029 
0030 % Get a list of cluster memberships
0031 A = spm_clusters(XYZ);
0032 
0033 % Create a list of indices that exceed cluster threshold
0034 Q = [];
0035 C = [];
0036 nc = 0;
0037 for i = 1:max(A)
0038   j = find(A == i);
0039   if length(j) >= k
0040     Q = [Q j]; 
0041     nc = nc+1;
0042     C(nc).size = length(j);
0043     C(nc).XYZ = XYZ(:,j);
0044     curr_good_idxs = good_vox_idxs(j);
0045     [C(nc).maxvoxval maxidxs] = max(Y(curr_good_idxs));
0046     C(nc).maxvoxidx_in_vol = curr_good_idxs(maxidxs);
0047     C(nc).maxvoxidx_in_clust = find(Y(curr_good_idxs)==C(nc).maxvoxval);
0048     C(nc).voxidx_in_vol = good_vox_idxs(j);
0049   end
0050 end
0051 
0052 % Weed out voxels that didn't pass the threshold
0053 fprintf('Found %d voxels belonging to %d clusters\n', length(Q), nc);
0054 XYZ = XYZ(:,Q);
0055 good_vox_idxs = good_vox_idxs(Q);
0056 
0057 % Filter the data volume
0058 Y1 = zeros(size(Y))*multiplier;
0059 Y1(good_vox_idxs) = 1;
0060 
0061 Y = Y.*Y1;
0062 
0063 return