function out = medianFilter(img, rad)
% Apply median filter with radius `rad` to image array `img`.
% Works with both greyscale (2D) and color/multichannel (3D) image arrays
% (filtering is performed per-channel).

% Query size of input image
[nr, nc, np] = size(img);
% Pre-allocate output image of same size
out = uint8(zeros(nr, nc, np));

% Loop over all pixels
for r = 1:nr
    for c = 1:nc
        % Clamp neighborhood to image boundary
        rmin = max(1, r - rad);
        rmax = min(nr, r + rad);
        cmin = max(1, c - rad);
        cmax = min(nc, c + rad);
        neighborhood = img(rmin:rmax,cmin:cmax,:);
        
        % Apply median filter to each channel
        % Vectorized alternative: out(r,c,:)=median(neighborhood, [1,2])
        for p = 1:np
            nbp = neighborhood(:,:,p);
            out(r, c, p) = median(nbp(:));
        end
    end
end
