MATLAB: Inevitable broadcast variable in parfor

broadcast variableimage processingParallel Computing Toolboxparfor

Hi,
I want to perform some calculations in a parfor loop using the neighbouring data in a 3D matrix. Currently I am using the following code, which has the problem that data is a 'broadcast variable'. This slows down the calculations and causes memory problems. (This is a small scale example where 18.5 MB is send to the workers.The actual data matrix will be around 100x100x65336 (i.e. 2.6 GB vs 3.6 MB in this example).)
rng('default');
height = 60;
width  = 60;
ntimes = 2^8;
radius = 10;
block_size = 2*radius+1;
data   = randn(height, width, ntimes, 'single'); % create some random data

output = zeros(height, width, ntimes, 'single');
result = zeros(height-2*radius, width-2*radius, block_size, block_size, ntimes, 'single');
ticBytes(gcp); tic;
parfor y = radius+1:height-radius   
    % create temporary array for parfor loop
    temp_r = zeros(width-2*radius, block_size, block_size, ntimes, 'single'); 
                
    for x = radius+1:width-radius                       
        % here I calculate the fft of a block surrounding the pixel of interest.
        % Since it violates rules of indexing, data is a broadcast variable
        a = fft(data(y-radius:y+radius, x-radius:x+radius, :), [], 3);   
        
        % For simplicity sake I store the fft (In reality I do more calculations).
        temp_r(x-radius, :, :, :) = real(a); 
    end
    
    result(y-radius, :, :, :, :) = temp_r;
end
toc; tocBytes(gcp);
To overcome this problem I considered using linear indexing. However, I then cannot figure out how to overcome the broadcasting. Due to the large size of the actual data, I think that separating the data in blocks (2*radius+1) for every single pixel is not feasible.
rng('default');
height = 60;
width  = 60;
ntimes = 2^8;
radius = 10;
data = randn(height, width, ntimes, 'single'); % create some random data
data = reshape(data, [height*width, ntimes]);
% determination of linear indices
[cols, rows] = meshgrid(radius+1:height-radius, radius+1:width-radius);
linear_ids = sub2ind([height, width], rows, cols);
% create an array with linear ids of neighbouring pixels
[x,y] = meshgrid(radius*(-height-1):height:radius*height-radius, 0:radius*2);
mask = x + y;
% the neighbour matrix has 2*radius+1 columns (the neighbours) 
% and the rows are the different pixels under study
neighbour_matrix = bsxfun(@plus, linear_ids(:), mask(:)');
output = zeros(length(cols) * length(rows),(2*radius+1)^2, ntimes);
ticBytes(gcp); tic;
parfor i = 1:length(cols) * length(rows) 
    neighbours = neighbour_matrix(i,:);
    
    % calculate fft using the neighbours here, but data is still a broadcast variable
    % but how do I make data a sliced variable
    a = fft(data(neighbours, :), [], 2); 
    output(i, :, :) = a;
end
toc; tocBytes(gcp);
I would be very grateful if someone could point me in the right direction. Thank you in advance!
tl;dr
Is it possible to avoid having a broadcast variable when a vector of indices is being used?

Best Answer

I doubt you need any for-loops at all. Because your fft is along the 3rd dimension in all the sub-blocks, you could just all the FFT effort with,
dataFFT= real( fft(data,[],3) );
The rest of the operations involved in deriving "result" consist of extracting sliding 2D windows from dataFFT and copying them into a larger matrix. You should probably not be doing that as it entails a lot of duplicate storage. What are you planning to do with the sub-blocks one you have them?
Related Question