Neighbor Indexing Example
From MATLAB Techniques for Image Processing by Steve Eddins.
Many image processing operations can be classified as recursive neighborhood operations. You do something to each neighbor of a pixel, then you do something to each neighbor of those neighbors, and so on. Flood filling is a good example of this type of operation.
Linear indexing can be used to implement such operations efficiently. I call the technique neighbor indexing.
Let's make a sample 5-by-5 matrix and consider three seed locations, or starting locations.
A = magic(5)
A = 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9
Starting locations at (2,2), (1,4), (3,4):
rows = [2 1 3]; cols = [2 4 4];
As you'll see, the neighbor indexing technique is based on linear indexing, so let's convert the starting locations to linear indices:
idx = sub2ind(size(A),rows,cols)
idx = 7 16 18
How do you find the linear indices of the eastern neighbor of each of the starting locations?
Answer: Add the number of rows
M = size(A,1); east_neighbors = idx + M
east_neighbors = 12 21 23
M is called the linear offset for the eastern direction.
A(east_neighbors) gives you the pixel values of the eastern neighbors:
A(east_neighbors)
ans = 7 15 22
What is the linear offset for the northern direction?
northern_offset = -1; northern_neighbors = idx + northern_offset
northern_neighbors = 6 15 17
OOPS! The (1,4) pixel doesn't have a northern neighbor!
The neighbor indexing technique DOES NOT WORK FOR PIXELS ON THE BORDERS!
Therefore some form of padding is usually necessary.
Here are the offsets corresponding to all eight neighbor directions:
East M Southeast M + 1 South 1 Southwest -M + 1 West -M Northwest -M - 1 North -1 Northeast M - 1
The function bsxfun gives you a way to compute all the neighbors of a set of pixels in a single, vectorized expression.
For example, suppose the linear indices of our seed pixels are:
idx = [275 2348 8100]
idx = 275 2348 8100
And we are interested in finding the north, east, south, and west neighbors of all these pixels:
% Assume a 100-by-100 image for this example.
M = 100;
neighbor_offsets = [-1, M, 1, -M]';
Add every neighbor_offset to every linear index:
neighbors = bsxfun(@plus,idx,neighbor_offsets)
neighbors = 274 2347 8099 375 2448 8200 276 2349 8101 175 2248 8000
Now we have enough indexing tools to implement a flood-fill algorithm. Let's try it with this image and starting location:
bw = bwperim(imread('circles.png')); imshow(bw) r = 100; c = 100; hold on plot(c,r,'g*') hold off

Initialize a linear index array containing the "active pixels." This is the set of pixels to be changed from black to white. Keep processing until the list of active pixels becomes empty.
active_pixels = sub2ind(size(bw),r,c); bw_out = bw; M = size(bw,1); neighbor_offsets = [-1, M, 1, -M]; h = imshow(bw_out); while ~isempty(active_pixels) % Set the active pixels to 1. bw_out(active_pixels) = 1; % The new active pixels list is the set of neighbors of the % current list. active_pixels = bsxfun(@plus,active_pixels,neighbor_offsets); active_pixels = active_pixels(:); % Remove from the active_pixels list pixels that are already % 1. This is how we keep the flood fill operation "inside" % the white boundaries. active_pixels(bw_out(active_pixels)) = []; % Remove duplicates from the list. active_pixels = unique(active_pixels); set(h,'CData',bw_out); drawnow; end imshow(bw_out)

How fast is the computation? I have made a function called fillWithNeighborIndexing that contains the code above without the graphic display.
A simple way to measure the execution time of something in MATLAB is to use tic and toc.
tic, bw_out = fillWithNeighborIndexing(bw,100,100); toc
Elapsed time is 0.029368 seconds.
A more reliable way to measure the execution time of things that take less than a second is to use the MATLAB function timeit. To use timeit, make a function handle that does exactly the thing you want to time, and then pass that function handle to timeit.
f = @() fillWithNeighborIndexing(bw,100,100); timeit(f)
ans = 0.0247
How does that time compare with the Image Processing Toolbox function imfill?
idx = sub2ind(size(bw),100,100); g = @() imfill(bw,idx,4); timeit(g)
ans = 0.0012
Note that this result is identical to what is obtained from the morphological operator IMFILL (when the starting point is specified and the connectivity is set to 4).
bw = bwperim(imread('circles.png'));
imshow(bw)

r = 100; c = 100; active_pixels = sub2ind(size(bw),r,c); imshow(imfill(bw,active_pixels,4))
