Efficient way to apply function to every submatrix of a given size in MATLAB

Question

In MATLAB, I want to apply a function to each submatrix of the given size of a matrix.

My current code works but is very slow.

%inputs
%    f: function to apply to submatrices
%    M: main matrix
%    h: height of submatrices
%    w: width of submatrices
%output
%    out: matrix of results of applying f to submatrices of M
function out = metricmap(f,M,h,w)
    [r,c] = size(M);
    h = h-1;w = w-1; %to make these values deltas 
    out = zeros([r,c]-[h,w]);
    
    for i = 1:(r - h)
        for j = 1:(c - w)
            subM = M(i:(i+h),j:(j+w));
            out(i,j) = f(subM);
        end
    end
end

Any suggestions are greatly appreciated!

What function do you need to apply? If it is a (weighted) average, or a max or min operation, then this can be very efficient. Otherwise it’s not going to be better than what you have now. The answer below encapsulates your code in a neat little function, but it’s the same thing and not faster. — Cris Luengo, Mar 16 '21 at 20:59
It is corr2(sliding box,some other matrix). Possible solution I thought of: 1. construct a 4 dimensional array newM where newM(i,j,:,:) = M(i:i+h,j,j+h). 2. use arrayfun to apply f to newM along the first 2 dimensions. I am not sure this is possible in MATLAB. That's how I would do it in mathematica. — user2752148, Mar 16 '21 at 22:17
If you can do without subtracting the mean, then you have a trivial convolution. That is very efficient. I don't think you can make `corr2` into a form that would be easy to apply to each window. — Cris Luengo, Mar 16 '21 at 23:02

score 2 · Answer 1 · answered Mar 16 '21 at 20:45

2

If you have the image processing toolbox, you can use blockproc to apply f on a [h w] sliding window of M:

out = blockproc(M, [h w], f);

answered Mar 16 '21 at 20:45

tdy

36,675
19
86
83

2

`blockproc` does exactly the same as OP’s function, it might be a bit more clever with boundary extension and so on, I don’t know. But it’s not any faster. – Cris Luengo Mar 16 '21 at 21:01
1

Thank you! If nothing else, this could make my code look cleaner. – user2752148 Mar 16 '21 at 22:14
Please correct me if I'm wrong, but it seems that blockproc operates on non-overlapping blocks of the source image. I want to operate on all subimages of a given size, including overlapping ones. – user2752148 Mar 16 '21 at 22:34
Just tested with `blockproc(eye(4),[2,2],@(x) mean2(x.data))` which returns a 2 by 2 matrix. I want it to return a 3 by 3 matrix, since there are 9 distinct 2 by 2 submatrices of eye(4), not just 4. – user2752148 Mar 16 '21 at 22:38
It looks like you have to do a bit of [hacking to make `blockproc` use overlapping windows](https://stackoverflow.com/a/29041292/13138364). – tdy Mar 16 '21 at 22:42
2

Actually I think [`nlfilter`](https://www.mathworks.com/help/images/ref/nlfilter.html) is the function you're looking for, not `blockproc`. – Cris Luengo Mar 16 '21 at 23:02

score 0 · Answer 2 · answered Mar 17 '21 at 18:23

As corr2 is the function that to be applied to submatrices you can use the following vectorized method of computing sliding_corr2(A,B) where size of B is smaller than A:

B_B = B - mean(B(:));
mA = conv2(A, ones(h, w) / (h * w), 'valid');
numerator = conv2(A, flipud(fliplr(B_B)), 'valid') - mA * sum(B_B(:));
denominator = sqrt(                        ...
    (                                      ...
      conv2(A .^ 2, ones(h, w), 'valid') - ...
      (h * w) * (mA .^ 2)                  ...
    )                                      ...
    * (B_B(:).' * B_B(:))                  ...
);
out = numerator ./ denominator;

Assuming that a is i,j sub-matrix of A the corr2 function is formulated as:

It can be rewritten as:

Now it can be written as convolution operations:

Where * is the convolution operator.

Efficient way to apply function to every submatrix of a given size in MATLAB

2 Answers2