function R = rtsGaussianBlurND(input, sigma)
%This function convolves an ND array with a gaussian kernel specified by
%sigma.  This takes advantage of the separability of the kernel and is
%therefore very fast.

%disp('begin***********************************************');

%get the number of dimensions for the convolution
dim = ndims(input);

%initialize the result to the source
R = input;


for d=1:dim
    %if the dimension is greater than 1 and the sigma not a dirac
    if size(input, d) > 1 && sigma(d) > 0
        %create a 1D kernel
        kernelshape = ones(1, dim);
        kernelshape(d) = sigma(d)*6;
        kernel = zeros(kernelshape);
        kernel(:) = ndgauss(sigma(d)*6, sigma(d));

        %perform the convolution
        R = convn(R, kernel, 'same');
    end
end