function Ib = gauss_blur3d(I, sigma)


%convolve along the x-axis
Gx = gauss1d(sigma(1));
Ix = convn(I, Gx, 'same');

%convolve along the y-axis
Gy = gauss1d(sigma(2))';
Iy = convn(Ix, Gy, 'same');

%convolve along the z-axis
Gz_prime = gauss1d(sigma(3));
Gz = reshape(Gz_prime, [1, 1, length(Gz_prime)]);
Iz = convn(Iy, Gz, 'same');

Ib = Iz;