How do I efficiently extract image patches at specified locations?

Question

I need to extract image patches of size s x s x 3 around specified 2D locations from an image (3 channels).

How can I do this efficiently without a for loop? I know I can extract one patch around (x,y) location as:

apatch = I(y-s/2:y+s/2, x-s/2:x+s/2, :)

How can I do this for many patches? I know I can use MATLAB's function blockproc but I can't specify the locations.

Answer 1

You can use im2col from the image processing toolbox to transform each pixel neighbourhood into a single column. The pixel neighbourhoods are selected such that each block is chose on a column-basis, which means that the blocks are constructed by traversing down the rows first, then proceeding to the next column and getting the neighbourhoods there.

You call im2col this way:

B = im2col(A, [M N]);

I'm assuming you'll want sliding / overlapping neighbourhoods and not distinct neighbourhoods, which are what is normally used when performing any kind of image filtering. A is your image and you want to find M x N pixel neighbourhoods transformed as columns. B would be the output where each neighbourhood is a single column and horizontally-tiled together. However, you'll probably want to handle the case where you want to grab pixel neighbourhoods along the borders of the image. In this case, you'll want to pad the image first. We're going to assume that M and N are odd to allow the padding to be easier. Specifically, you want to be sure that there are floor(M/2) rows padded on top of the image as well as the bottom as well as floor(N/2) columns padded to the left of the image as well as the right. As such, we should pad A first by using padarray. Let's assume that the border pixels will be replicated, which means that the padded rows and columns will simply be those grabbed from the top or bottom row, or the left and right column, depending on where we need to pad. Therefore:

Apad = padarray(A, floor([M N]/2), 'replicate');

For the next part, if you want to choose specify neighbourhoods, you can use sub2ind to convert your 2D co-ordinates into linear indices so you can select the right columns to get the right pixel blocks. However, because you have a colour image, you'll want to perform im2col on each colour channel. Unfortunately, im2col only works on grayscale images, and so you'd have to repeat this for each channel in your image.

As such, to get ready for patch sampling, do something like this:

B = arrayfun(@(x) im2col(Apad(:,:,x), [M N]), 1:size(A,3), 'uni', 0);
B = cat(3, B{:});

The above code will create a 3D version of im2col, where each 3D slice would be what im2col produces for each colour channel. Now, we can use sub2ind to convert your (x,y) co-ordinates into linear indices so that we can choose which pixel neighbourhoods we want. Therefore, assuming your positions are stored in vectors x and y, you would do something like this:

%// Generate linear indices
ind = sub2ind([size(A,1) size(A,2)], y, x);

%// Select neighbourhoods
%// Should be shaped as a MN x len(ind) x 3 matrix
neigh = B(:,ind,:);

%// Create cell arrays for each patch
patches = arrayfun(@(x) reshape(B(:,x,:), [M N 3]), 1:numel(ind), 'uni', 0);

patches will be a cell array where each element contains your desired patch at each location of (x,y) that you specify. Therefore, patches{1} would be the patch located at (x(1), y(1)), patches{2} would be the patch located at (x(2), y(2)), etc. For your copying and pasting pleasure, this is what we have:

%// Define image, M and N here
%//...
%//...

Apad = padarray(A, floor([M N]/2), 'replicate');
B = arrayfun(@(x) im2col(Apad(:,:,x), [M N]), 1:size(A,3), 'uni', 0);
B = cat(3, B{:});

ind = sub2ind([size(A,1) size(A,2)], y, x);
neigh = B(:,ind,:);
patches = arrayfun(@(x) reshape(neigh(:,x,:), [M N 3]), 1:numel(ind), 'uni', 0);

Answer 2

As unexpected as this may seem, but for me the naive for-loop is actually the fastest. This might depend on your version of MATLAB though, as with newer versions they keep on improving the JIT compiler.

Common data:

A = rand(30, 30, 3); % Image
I = [5,2,3,21,24]; % I = y 
J = [3,7,5,20,22]; % J = x
s = 3; % Block size

Naive approach: (faster than im2col and arrayfun!)

Patches = cell(size(I));
steps = -(s-1)/2:(s-1)/2;
for k = 1:numel(Patches);
    Patches{k} = A(I(k)+steps, ...
                   J(k)+steps, ...
                   :);
end

Approach using arrayfun: (slower than the loop)

steps = -(s-1)/2:(s-1)/2;
Patches = arrayfun(@(ii,jj) A(ii+steps,jj+steps,:), I, J, 'UniformOutput', false);