Savitzky Golay filter on 4D images

I have stacks of MRI images for a number of patients, and I want to use the savitzky golay filter on each slice for all the patients. My problem is that I cant figure out how to apply the filter on every slice in the stack.

Here is my code so far:

 [Files, Path, OK] = uigetfile('*.*',...
        'Choose the images (stacks)to analyze', 'MultiSelect', 'on');
    if ~OK % Sjekkar at alt er OK så langt.
        return 
    end

    NumFiles = numel(Files); % Numer of files
    XPas = cell(NumFiles,1); % Cell array for X-matrix for each patient
    YPas = cell(NumFiles,1); % Cell array for Y-vector for each patient
    Tsteps = nan(NumFiles,1); % Number of time steps for each stack. This can be different for each patient 
    indsPas = nan(NumFiles,2); % start and end-index for each patient
    ImSize = nan(NumFiles,2); % Imagesize for each patient
       k = 1;
      for i = 1:NumFiles

        % Import stack no i 
        im = importdata([Path,Files{i}]);

        [m,n,t,NumFrames] = size(im);

        Tsteps(i) = t; % number of timesteps

       ImSavitzky = zeros(size(im));

                for l=1:NumFrames
%                 [row columns numberOfColorBands] = size(imageArray);
%                 imageArray = imread('cameraman.tif');
                imageArray = double(im);
                    % Apply the Savitzky-Golay filter.
                    % First apply it in the vertical (row) direction.
                    p = 2; % Order of the polynomial
                    verticallySmoothedImage = sgolayfilt(imageArray, p, 3);

                    % Apply the Savitzky-Golay filter.
                    % First apply it in the vertical (row) direction.
%                     p = 2; % Order of the polynomial
%                     horizontallySmoothedImage = sgolayfilt(imageArray, p, 3);

                    imSavitzky(:,:,t,l) = sgolayfilt(verticallySmoothedImage, p, 3);         
                end

        indsPas(i,1) = k;
        indsPas(i,2) = k + size(YPas{i},1) - 1;

        ImSize(i,1) = m;
        ImSize(i,2) = n;

        k = k + size(YPas{i},1);

        implay(ImSavitzky)

      end