Applying Sigma Filter to Dicom Image Stack -- Extremely Slow

I am having some trouble trying to institute a sigma filter on a dicom image stack (512x512x60). Basically: 'midpixel' is found, if any surrounding pixels in a kernel (5x5x3) are in the same grayscale range (within the sigma value), then average those pixels and log that new pixel value in an output matrix. The code is doing what I want it to, the problem is that it takes about 30 minutes for Matlab to complete this code in a stack this large. I'm sure there is a more elegant way to code this nested set of loops, or some way to better use the image processing toolbox that could seriously improve speed. Any help would be appreciated. See code below:

kernelsize=[5 5 3];
sigmavalue=20;
outimage=zeros(512,512,60); %preallocate the output image with zeroes

[xsize, ysize, zsize]=size(dicomstack);
dicomstack=double(dicomstack);
for zz=1:zsize
    for yy=1:ysize
        for xx=1:xsize
            pixsums=0; numpix=0; %set sum and number of pixels to 0
            midpixel=dicomstack(xx,yy,zz); %midpixel HU value is found
            for z2=-1:1:1 % z kernel size = 3
                zzz=zz+z2;
                if zzz < 1 || zzz > zsize %these check if the new pixel will be outside the image bounds
                    continue
                end
                for y2=-2:1:2 % y kernel size = 5
                    yyy=yy+y2;
                    if yyy < 1 || yyy > ysize
                        continue
                    end
                    for x2=-2:1:2 % x kernel size = 5
                        xxx=xx+x2;
                        if xxx < 1 || xxx > xsize
                            continue
                        end
                        newpixel=dicomstack(xxx,yyy,zzz); %newpixel is one to compare to middle pixel
                        if abs(midpixel - newpixel) < sigmavalue %check to see if difference is less than sigma value
                            pixsums = pixsums + newpixel; %if it is, bin that pixel value into a matrix
                            numpix = numpix + 1;
                        end
                    end
                end
            end
            if numpix > 0
                outimage(xx,yy,zz)=pixsums/numpix;
            else
                outimage(xx,yy,zz)=midpixel;
            end
        end
    end
end

Thanks