function matOut = voxelDownsample(matIn,scaleFac,dataType)
% voxelDownsample Make a voxel image smaller
% B = voxelDownsample(A,s) makes the voxel image defined by matrix matIn 
% smaller by a factor of scaleFac and returns the resulting matrix in
% matOut. 
%
% Example: 20x20x20 input image with scaling factor of 2 will result
% in a 10x10x10 output image.
%
% Notes:
% 1. Slight shifting might occur for images scaled with non-integer
%    factors. 
% 2. Upscaling (scaleFac < 1) is currently not supported, but
%    sometimes might produce reasonable results.
% 3. Currently scaling by the same factor in all dimensions is the only
%    supported mode. While it might not make much sense to do, one may
%    eventually wish to scale in different dimensions by a different 
%    factor. This might be added.
%
% +----------------------------------------------------------------------+ 
% | *** INPUTS ***                                                       |
% | matIn - n*m*p matrix with 0 for the matrix phase and 1 for the other |
% |         phase                                                        |
% | scaleFac - value by which to scale the image (tested for 2 and 4)    |
% |                                                                      |
% | *** OUTPUTS ***                                                      |
% | matOut - n/sizeRatio * m/sizeRatio * p/sizeRatio (all rounded down)  |
% |          matrix with 0 for matrix phase and 1 for other phase.       |
% +----------------------------------------------------------------------+
% 
% Orion L. Kafka, Northwestern University
%
% File history:
% 1  Dec 2016 - Written
% 29 Apr 2017 - Added datatype control, changed SumA to sumA

% Check for some easy mistakes
assert(ndims(matIn) == 3,...
    'Invalid entry: matIn. Number of dimensions must = 3')
assert(size(matIn,1) > scaleFac || ...
       size(matIn,2) > scaleFac || ...
       size(matIn,3) > scaleFac, ...
    ['Invalde entry: matIn OR scaleFac.',...
     ' Matrix size in each dim must be > scaleFac'])
assert(scaleFac>0,...
    'Invalid entry: scaleFac = %f Scaling factor must be > zero.',...
    scaleFac)

if scaleFac < 1
    warning('scaleFac = %f will make the image BIGGER.',scaleFac)
    warning('Upscaling is not well supported; proceeding anyway.')
end

% Preallocate output array
matOut = zeros(floor(size(matIn,1)/scaleFac),...
               floor(size(matIn,2)/scaleFac),...
               floor(size(matIn,3)/scaleFac),dataType);

% In a region around the current point of interest, add up all the voxel
% values. If the region is more than half filled with 1's make the voxel a
% 1. If less than half, make it a 0. Region size is defined by the image
% size and scaling factor.
for I = 1:size(matIn,1)/scaleFac
    for J = 1:size(matIn,2)/scaleFac
        for K = 1:size(matIn,3)/scaleFac
            
            sumA = sum(sum(sum(matIn(...
                floor(I*scaleFac-(scaleFac-1)):ceil(I*scaleFac),...
                floor(J*scaleFac-(scaleFac-1)):ceil(J*scaleFac),...
                floor(K*scaleFac-(scaleFac-1)):ceil(K*scaleFac)))));            
            
            if sumA >= (scaleFac^3)/2 + 1
                matOut(I,J,K) = 1;
            else
                matOut(I,J,K) = 0;
            end
            
        end
    end
end