I am trying to use parfor for a large matrix operation. I am getting Out of Memory during deserialization error. Is there a way to minimize the memory? Below is the minimal example code:

clc; clear;
warpedImages = num2cell(uint8(randi([0,255], 1654, 6288, 3, 35)),1:3);
warpedImages = reshape(warpedImages,1,[]);
% Initialze
n = length(warpedImages);
sigmaN = 10;
sigmag = 0.1;
panoramasize = size(warpedImages{1});
Amat = cell(n);
Bvec = zeros(n,1);   
IuppeIdx = nonzeros(triu(reshape(1:numel(Amat), size(Amat))));
Amat_temp = cell(1,length(IuppeIdx));    
matSize = size(Amat);
    
% 4D warped images (Slicing)
wim_4d = cell2mat(reshape(warpedImages,1,1,1,[]));
% Get the Ibarijs and Nijs
parfor i = 1:length(IuppeIdx)
    % Index to subscripts
    [ii,jj] = ind2sub(matSize, IuppeIdx(i));                 
    if ii == jj
        diag_val_1 = 0;
        diag_val_2 = 0;
       
        Z = 1:n;
        Z(Z==ii) = [];
        for d = Z
            [Ibarij, Ibarji, Nij] = getIbarNij(panoramasize, wim_4d(:,:,:,ii), wim_4d(:,:,:,d));
            diag_val_1 = diag_val_1 + ( (Nij + Nij) .*  Ibarij.^2 );
            diag_val_2 = diag_val_2 + Nij;
        end
        diag_val = diag_val_1 + (sigmaN^2/sigmag^2) * diag_val_2;
        B_val = (sigmaN^2/sigmag^2) * diag_val_2;
        Amat_temp{i} = diag_val;
        Bvec(i) = B_val
    end       
    if ii ~= jj
        [Ibarij,Ibarji,Nij] = getIbarNij(panoramasize, wim_4d(:,:,:,ii), wim_4d(:,:,:,jj));
        Amat_temp{i} = -(Nij+Nij) .* (Ibarij .* Ibarji);
    end
end
function [Ibarij,Ibarji,Nij] = getIbarNij(panoramasize, Imij, Imji)
    Ibarij = zeros(panoramasize,'uint8');
    Ibarji = zeros(panoramasize,'uint8');
    
    % Overlay the warpedImage onto the panorama.
    maski = imbinarize(rgb2gray(255 * Imij));  
    maskj = imbinarize(rgb2gray(255 * Imji));
    
    % Find the overlap mask
    Nij_im  = maski & maskj;
    Nij_im  = imfill(Nij_im, 'holes');    
    Nijidx = repmat(Nij_im, 1, 1, size(Imij,3));
    
    % Get the overlapping region RGB values for two images
    Ibarij(Nijidx) = Imij(Nijidx);
    Ibarji(Nijidx) = Imji(Nijidx);
    
    % Convert to double
    Ibarij_double = double(Ibarij);
    Ibarji_double = double(Ibarji);
    
    % Nij
    Nij = sum(sum(Nij_im));
    
    % Ibar ijs
    Ibarij = reshape(sum(sum(Ibarij_double)) ./ Nij, 1, 3);
    Ibarji = reshape(sum(sum(Ibarji_double)) ./ Nij, 1, 3); 
    
    % Replace NaNs by zeros
    Ibarij(isnan(Ibarij)) = 0;
    Ibarji(isnan(Ibarji)) = 0;
end

Line [Ibarij, Ibarji, Nij] = getIbarNij(panoramasize, wim_4d(:,:,:,ii), wim_4d(:,:,:,d)); throws a warning message: The entire array or structure ‘wim_4d’ is a broadcast variable. This might result in unnecessary communication overhead. I have used ind2sub for getting the subscripts as it is easy to work. However, wim_4d(:,:,:,ii) and others cannot be sliced. Any other suggestion and help is appreciated!