function [ img_out ] = resamplegrid ( img_inDat, resample_factor )
%RESAMPLEGRID resample a matrix by interpolation
%   Usage:
%   [ M_resample ] = resamplegrid ( M, resample_factor )
%       M_resample: resampled matrix
%       M: input matrix
%       resample_factor: floating point vector, size of resampled grid is
%           round(resample_factor*dimX), round(resample_factor*dimY), dimZ
%
%   Development status: under development
%   Copyleft (c) Xiaohu Mo
%   Version: 0.1
%
%   TODO
%   3d resample
%   different resample factor in each dimension
%   interpolation method

if isinteger(img_inDat)
    warning('Converted integer input array to float');
    img_inDat = double(img_inDat);
end

is_input_flat = false;
if ~isstruct(img_inDat)
    is_input_flat = true;
    img_in.data = img_inDat;
    img_in.start = [0 0 0];
    img_in.spacing = [1 1 1];
end

if numel(resample_factor) == 1
    % resample_factor = repmat(resample_factor, size(size(img_in.data)));
    resample_factor = [resample_factor resample_factor 1];
end

xg1 = img_in.start(1) + img_in.spacing(1) * (0:size(img_in.data,1)-1);
xg2 = img_in.start(2) + img_in.spacing(2) * (0:size(img_in.data,2)-1);
xg3 = img_in.start(3) + img_in.spacing(3) * (0:size(img_in.data,3)-1);
ext_out = round(size(img_in.data) .* resample_factor);
img_out.data = zeros(ext_out);
img_out.spacing = img_in.spacing .* size(img_in.data) ./ size(img_out.data);
img_out.start = img_in.start - img_in.spacing/2 + img_out.spacing/2;
xq1 = img_out.start(1) + img_out.spacing(1) * (0:size(img_out.data,1)-1);
xq2 = img_out.start(2) + img_out.spacing(2) * (0:size(img_out.data,2)-1);
xq3 = img_out.start(3) + img_out.spacing(3) * (0:size(img_out.data,3)-1);

[Xg1, Xg2, Xg3] = ndgrid(xg1, xg2, xg3);
[Xq1, Xq2, Xq3] = ndgrid(xq1, xq2, xq3);
% note: 'linear' results in all NaN in case of extrapolation
img_out.data = interpn(Xg1,Xg2,Xg3,img_in.data, Xq1, Xq2, Xq3, 'spline');
% F = griddedInterpolant(Xg1, Xg2, Xg3, img_in.data, 'spline');
% img_out.data = F(Xq1, Xq2, Xq3);

if is_input_flat
    img_out = img_out.data;
end