function result = extentThreshold(image, indices, k)
% function result = extentThreshold(image, indices, k)
% Applies the extent threshold function of SPM5 on 'image' with given
% 'indices' and cluster size 'k' 

[x, y, z] = ind2sub(size(image), indices);

XYZ = [x y z];

% from spm_getSPM line 676:

%-Calculate extent threshold filtering (from spm_getSPM, line 676)
%-------------------------------------------------------------------
A     = spm_clusters(XYZ');
Q     = [];
for i = 1:max(A)
    j = find(A == i);
    if length(j) >= k; Q = [Q j]; end
end

% ...eliminate voxels
%-------------------------------------------------------------------
XYZ   = XYZ(Q,:);

result = zeros(size(image));
inds = sub2ind(size(image), XYZ(:,1), XYZ(:,2), XYZ(:,3));
result(inds) = image(inds);