Efficient way to extract neighbours pixel of thousands locations simultaneously

Question

I tried to extract the 8 neighbors pixel of an image. However I would like to perform this operation simultaneously for hundreds of locations (row and column in the function). That mean the output Block is a 3D matrix, each slice of the 3D matrix is corresponding to single location.

function Block = getNeighbours(image, row, column)
% Create a 3x3 matrix that contains the neighbors of the point (row,column)
    row=round(row(:));
    column=round(column(:));
    neighbors_x = [row(:,1)-1 row(:,1) row(:,1)+1];
    neighbors_y = [column(:,1)-1 column(:,1) column(:,1)+1];
    Block = zeros(3,3,size(row,1));
    for i=1:size(row,1)
    Block(:,:,i) = image(neighbors_x(i,:), neighbors_y(i,:)); %can I avoid this loop?
    end
end

From the code above, I need to loop to finish the job. It seems to be the bottle neck of this function, and for thousands of locations, this is definitely not efficient. Any way to avoid this?

You can try this function by:

image=randi([1 255], [300 300]);
row=randi([1 200], [1 1000]);
column=randi([1 200], [1 1000]);

block=getNeighbours(image, row, column);

Answer 1

If you have the Image Processing Toolbox, one option is to use im2col to preprocess your image into blocks, generating all of the neighborhoods at once. This would be particularly effective if you need to access blocks numerous times.

Given an image:

>> img = reshape(1:25,5,5)

img =

    1    6   11   16   21
    2    7   12   17   22
    3    8   13   18   23
    4    9   14   19   24
    5   10   15   20   25

im2col converts each valid block (this is important...) into a column of a matrix:

C = im2col(img,[3,3])

C =

    1    2    3    6    7    8   11   12   13
    2    3    4    7    8    9   12   13   14
    3    4    5    8    9   10   13   14   15
    6    7    8   11   12   13   16   17   18
    7    8    9   12   13   14   17   18   19
    8    9   10   13   14   15   18   19   20
   11   12   13   16   17   18   21   22   23
   12   13   14   17   18   19   22   23   24
   13   14   15   18   19   20   23   24   25

As you can see, the first column of C has the values of the neighborhood of index 7, or img(2,2). The next one is index 8, or img(3,2), and so on. The offset is because im2col only gives you valid blocks. If this is what you want, you just have to remember to subtract 1 from the row/column subscripts before finding the column number in C using sub2ind:

C_idx = sub2ind(size(img), 2-1, 2-1)

C_idx =  1

(If you pad the image, the indices and subscripts are the same as in the original image.)

So the column index of the neighborhood of img(2,2) is 1.

You could actually stop right here. If you're applying the same kernel to all of the queried neighborhoods, you could just make the kernel into a row vector and multiply, choosing which columns of C you want to consider:

kernel = ones(3)/9
kernel =

   0.11111   0.11111   0.11111
   0.11111   0.11111   0.11111
   0.11111   0.11111   0.11111

kernel(:).' * C
ans =

    7.0000    8.0000    9.0000   12.0000   13.0000   14.0000   17.0000   18.0000   19.0000

C_idx = [1:9];   % or randi(9, 1, 4), or any valid list of column numbers
k_vec = kernel(:).';   % turn kernel into a row vector
result = k_vec * C(:, C_idx);

result =

    7.0000    8.0000    9.0000   12.0000   13.0000   14.0000   17.0000   18.0000   19.0000

---

If you're not happy with the neighborhoods being column vectors and really want them to be 3x3 blocks, you can reshape the C matrix into a 3D matrix:

D = reshape(C, 3, 3, []);

D =

ans(:,:,1) =

    1    6   11
    2    7   12
    3    8   13

ans(:,:,2) =

    2    7   12
    3    8   13
    4    9   14

ans(:,:,3) =

    3    8   13
    4    9   14
    5   10   15

ans(:,:,4) =

    6   11   16
    7   12   17
    8   13   18

ans(:,:,5) =

    7   12   17
    8   13   18
    9   14   19

ans(:,:,6) =

    8   13   18
    9   14   19
   10   15   20

ans(:,:,7) =

   11   16   21
   12   17   22
   13   18   23

ans(:,:,8) =

   12   17   22
   13   18   23
   14   19   24

ans(:,:,9) =

   13   18   23
   14   19   24
   15   20   25

From here, you use the same indices as you used before for the columns, but now they're indices into the third dimension, e.g.:

>> D(:, :, 1)
ans =

    1    6   11
    2    7   12
    3    8   13

gives you the same neighborhood as C(:, 1).

Answer 2

you can use sub2ind on the row-column indexes and then add a general 9-neighborhood indexes to form block indexes. For large number of row-column pairs it's about 10 times faster than using a loop:

image=randi([1 255], [300 300]);
row=randi([2 200], [1 10000]);
column=randi([2 200], [1 10000]);
tic;
sz = size(image);
% general 9-neighberhood indexes
hoodIdxs = [-sz(2)-1,-1,sz(2)-1;-sz(2),0,sz(2);-sz(2)+1,1,sz(2)+1];
% linear indexes of requested pixels
idxs = permute(sub2ind(sz,row,column),[1 3 2]);
% block indexes
blockIdxs = bsxfun(@plus,idxs,hoodIdxs);
% generate blocks
blocks = image(blockIdxs);
toc;
% OP's function
tic;
blocks_loop = getNeighbours(image, row, column); % ~10 times slower
toc
% check
isequal(blocks, blocks_loop) % yes

note that when generating random row and columns I set the lower end of the randi interval to 2, otherwise a row-column on the image edge will throw an indexing error when trying to call blocks = image(blockIdxs);. You can easily overcome this by padding the image (here padded with nan's):

% pad image edges with nans
padded = padarray(image,[1 1],nan,'both');
sz = size(padded);
% general 9-neighberhood indexes
hoodIdxs = [-sz(2)-1,-1,sz(2)-1;-sz(2),0,sz(2);-sz(2)+1,1,sz(2)+1];
% linear indexes of requested pixels (+1 because of the padding)
idxs = permute(sub2ind(sz,row + 1,column + 1),[1 3 2]);
% block indexes
blockIdxs = bsxfun(@plus,idxs,hoodIdxs);
% generate blocks
blocks = padded(blockIdxs);