function [downSampledData, start, stepSize, idx] = downSampleData(data, origStepSize, desiredStepSize, downSamplingType)

if nargin < 4
    downSamplingType = 'noInterp';
end

downSamplingType = lower(downSamplingType);

d0 = size(data);

switch downSamplingType
    case 'nointerp'
        % (From SPM): Discrete number of voxels to skip when downsampling
        sk = max([1, 1, 1], round(desiredStepSize * [1, 1, 1] ./ origStepSize));
        
        % downsample
        downSampledData = data(1:sk(1):d0(1), 1:sk(2):d0(2), 1:sk(3):d0(3));
        
        % voxels that are kept in original image
        idx = zeros(size(data));
        idx(1:sk(1):d0(1), 1:sk(2):d0(2), 1:sk(3):d0(3)) = 1;
        idx = idx > 0;
        
        % First voxel starts in 0,0,0 (center)
        start = [0,0,0];
        
        % Even if we specify a desiredStepSize, we're not guaranteed to get it if we do not interpolate
        stepSize = sk .* origStepSize;
        
    case 'nearestneighbor'
        d = ceil((d0-1)/desiredStepSize)+1;
        
        % The resampling assumes that we start and end at the outermost voxels, effectively moving the first and last
        % center to center -/+ 0.5*step respectively
        % round gives us the index that's closest to the voxel we're sampling for
        idxX = round(linspace(1,d0(1),d(1)));
        idxY = round(linspace(1,d0(2),d(2)));
        idxZ = round(linspace(1,d0(3),d(3)));              
        downSampledData = data(idxX, idxY, idxZ);
        
        % First voxel starts in 0,0,0 (center)
        start = [0,0,0];
        
        % This is actually not that precise...
        stepSize = desiredStepSize * [1,1,1];
    otherwise
        error('downsampling type not defined.');
end

end