Use convn() with a kernel like this:
kernel = ones(3,3,3);
kernel(2,2,2) = 0; % Exclude the central pixel
out = convn(inputImage, kernel, 'same');
Now that gives the mean at each 3x3x3 volume. Now to get 26 of those for each voxel, you will have to call convn() again with kernels that are big enough to contain a single voxel out at the center of the subvolumes. I trust you can figure that out on your own.
