Trying to vectorize a column-wise image processing step

I have a single-precision dataset specdata of size nPix x nObs. Each column 1..nObs is a (diffcols x diffrows) image stored unfolded. (nPix=diffrows*diffcols)

I can batch process them easily like this:

blurredData = zeros( size( specdata ), 'single');
for a=1:nObs
    blurredData(:,a) = reshape(imgaussfilt( reshape(specdata(:,a),...
        diffcols,diffrows), 0.75 ), nPix,1);
end

But, I would like to cleverly vectorize this loop if possible. This:

f = @( D, diffrows, diffcols, nPix ) reshape(imgaussfilt( reshape(D,diffcols,diffrows), 0.75 ), nPix,1);
ii = 1:nObs;
S(:, ii ) = f(specdata(:,ii),diffrows,diffcols,nPix);

blows up ("To RESHAPE the number of elements must not change."), apparently because I'm passing all of specdata to the reshape.

Any thoughts? Is there some obvious vectorization I am missing?

Thanks.

Answer 1

[Guillaume]

Your code is very easy to vectorise because imgaussfilt is very happy to operate at once on a stack of 2D images (i.e. a 3D matrix). All you have to do is reshape your whole matrix into that stack:

blurreddata = imgaussfilt(reshape(specdata, diffrows, diffcols, []), 0.75);

This gives back a stack of filtered images that you could indeed reshape back into a 2D array, but in my opinion it makes more sense to keep the array as a 3D matrix.

Note that this is only possible because imgaussfilt can work on the stack all at once. For other functions that can only operate on one image, you don't have a choice but to loop. If you keep the data as a 3D matrix, you don't have to reshape inside the loop, just iterate over the 3rd dimension. On the other hand reshape is a very fast operation (just rewrite the matrix header) so it's not very critical.

Answer 2

[Chad]

Excellent, thank you! As a bonus, the version you suggested runs about 2x as fast on my GPU (Quadro P5000) than on my CPU (Xeon 8-core).