Every so often someone asks in comp.soft-sys.matlab for a "function like Image Processing Toolbox function foobar that works in 3-D." If I happen to see the post, I will follow up, explaining that foobar itself works just fine in 3-D.

So in case some blog readers have been looking for a "3-D version of foobar," I thought I'd mention it here.

In the late 1990s we started hearing about requests for three-dimensional image processing support. Then on some customer visits in 2000 we heard some demand for four-dimensional support! As a result, the major Image Processing Toolbox upgrade in 2001 (version 3.0) added arbitrary dimension support to the toolbox.

Many Image Processing Toolbox functions implement mathematical operations (such as convolution) that are perfectly well defined in multidimensional spaces. So our strategy was to go ahead and implement multidimensional support for all such functions.

The meaning of the higher dimensions varies with the application and the data. With volumetric data (pixels become "voxels"), the 3rd dimension is spatial, just like the first two. Sometimes a higher dimension is time, or sequence number. And sometimes a higher dimension is spectral, as in color images or multispectral imagery.

At some point after 2001, we began to realize that users weren't discovering the multidimensional capability of the toolbox. Our first attempt to rectify that problem came when we added documentation sections explaining how to use toolbox functions to process sequences. (See "Working with Sequences" in the Users Guide.) Through a kind of happy mathematical coincidence, it turns out that image processing operators working in higher dimensions often just do "the right thing" when you have an image sequence represented in MATLAB as a three-dimensional array.

Let me explain that statement. Suppose you have a 640-by-480-by-100 array, representing 100 images in a sequence. And suppose you want to filter each image in the sequence with the same filter. Well, it turns out that multidimensional convolution of a 640-by-480-by-100 input with a "flat" filter (that is, an M-by-N-by-1 filter) is exactly equivalent to filtering each image in the sequence independently.

For many users, though, this wasn't at all obvious, even if they did discover the multidimensional capability of many of the functions. So that's what led to the "Working with Sequences" section in the doc. But the multidimensional functions can also handle volumetric processing, or other multidimensional applications.

Still, multidimensional support remains somewhat "hidden" in the doc. Here's where you can almost always find it: Look in the "Class Support" section of a function's reference page. For example, here's what the Class Support section for watershed says:

"A can be a numeric or logical array of any dimension, and it must be nonsparse."

Whenever you see that "of any dimension" phrase, you'll know that the function offers multidimensional support.

PS. Don't ask me to send you the code for foobar. I lost it.

Get the MATLAB code

Published with MATLAB® 7.6

When we work on improving the speed of image processing functions in MATLAB, we are naturally interested in getting timing measurements. And we know customers are often timing their own code, typically by using the MATLAB functions tic and toc.

It turns out to be surprisingly difficult to get reliable, repeatable timing measurements on modern computers. This is true for any program, not just MATLAB.

I've written timing code many different times over the years, using a variety of methods. Along the way I've learned some important tips and pitfalls from people like Cleve Moler and Bill McKeeman. Late last year I was working on YABP (Yet Another Benchmark Program), and I wondered if I could encapsulate at least some of those tips into a timer utility function, one that anyone could use without thinking too hard or making too many choices. The result was the function timeit.m, which you can find on the MATLAB Central File Exchange.

timeit takes one input argument and returns one output argument:

   s = timeit(f)

f is a function handle, and s is the time, in seconds, required to call that function handle. timeit calls f with no input arguments. To use timeit successfully, it's very helpful to know how to make anonymous function handles. Here are a couple of examples:

How much time does it take to compute sum(A.' .* B, 1), where A is 12000-by-400 and B is 400-by-12000?

A = rand(12000, 400);
B = rand(400, 12000);
f = @() sum(A.' .* B, 1);
timeit(f)

ans =

    0.1125

How much time does it take to dilate the text.png image with a 25-by-25 all-ones structuring element?

bw = imread('text.png');
se = strel(ones(25, 25));
g = @() imdilate(bw, se);
timeit(g)

ans =

    0.0037

The function timeit is designed to handle three key timing issues automatically:

Also, since the calls to f all occur within a function, the large variability associated with timing code from the command prompt is eliminated.

The function is not perfect, and other people might make different choices than I did. For example, you could argue that I should be using the minimum time instead of the median time for the time trials. I chose median because in my experience I have found that it takes fewer repetitions to get a repeatable result by using the median. Also, timeit will call f at least 13 times, which is possibly more than necessary for some longer-running functions.

If you are interested in issues and techniques related to performance measurement, particularly in MATLAB, then I highly recommend that you read Bill's white paper on performance measurement.

If you have ideas about improving timeit, I'd be very interested to hear them.

Get the MATLAB code

Published with MATLAB® 7.5

Earlier this year I wrote about logical indexing and linear indexing.

Now I want to finish this little indexing series with an example of a technique I call neighbor indexing. I've written about neighbor indexing before, but in this post I'll show a more complete example of its application to image processing.

Many image processing operations can be classified as recursive neighborhood operations. (Thanks to Chris for suggesting this term to me.) 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, and neighbor indexing is a good way to implement this type of operation in MATLAB. If you implement your own image processing algorithms in MATLAB, you'll probably find neighbor indexing useful.

I'll start by reviewing the basics of neighbor indexing, and then I'll show how to use it to implement flood filling.

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

The linear indices can be used to extract the pixel values of each of the three seed locations:

A(idx)

ans =

     5     8    20

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! Subtracting 1 from the linear indexing of a pixel on the top row will just give you the pixel on the bottom row of the previous column.

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

Pad the image with ones. This step will prevent the flood filling operation from reaching the border of the image.

bw = padarray(bw, [1 1], 1);

Initialize a linear index array containing the "active pixels." This is the set of pixels to be changed from black to white. Remember that we have to add 1 to the original row and column coordinates because of the padding step above.

active_pixels = sub2ind(size(bw), r+1, c+1);
bw_out = bw;
M = size(bw, 1);
neighbor_offsets = [-1, M, 1, -M];

Our main processing loop will iterate until there are no more active pixels. In each loop iteration, we set the active pixels to 1 in the output image. Then we get all the neighbors of all the active pixels, and we make those neighbors the new active pixels.

That's the basic logic, but there are two special steps you should keep in mind if you use this technique yourself. First, we need to remove from the active pixel list the pixels that are already 1. If we don't do this, the flood fill operation will not stop at the white boundaries.

Second, we need to "dedup" the active pixel list at each iteration. That is, when you find all the neighbors of all the active pixels, there will be duplicates in the list. If you don't prune out the duplicates at each iteration, the number of duplicates can quickly grow to be quite large, causing the algorithm to run slowly and to use a lot of memory.

Here's the main processing loop:

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);
end

Remember to remove the extra padding:

bw_out = bw_out(2:end-1, 2:end-1);

And here's the final result:

imshow(bw_out)

Get the MATLAB code

Published with MATLAB® 7.5

Last week I posted an introduction to logical indexing. This week I want to continue with a brief discussion of linear indexing and its connection to image processing in MATLAB.

Lets start with a small matrix. Breaking with tradition, I'll use a Hilbert matrix instead of a magic square:

A = hilb(5)

A =

    1.0000    0.5000    0.3333    0.2500    0.2000
    0.5000    0.3333    0.2500    0.2000    0.1667
    0.3333    0.2500    0.2000    0.1667    0.1429
    0.2500    0.2000    0.1667    0.1429    0.1250
    0.2000    0.1667    0.1429    0.1250    0.1111

POP QUIZ: What does the MATLAB expression A(17) mean?

A(17)

ans =

    0.2000

MATLAB interprets a single subscript as an index into a vector containing all the values of A in column order. So A(17) is the same as A(2, 4).

A(2, 4)

ans =

    0.2000

This is called linear indexing.

So what's the connection to image processing? Suppose A is our image, and some previous computation has told us that we are interested in the pixel values at these row and column coordinates:

row = [2 1 5];
col = [1 3 4];

Can we extract the desired pixel values in a single expression? Let's try:

A(row, col)

ans =

    0.5000    0.2500    0.2000
    1.0000    0.3333    0.2500
    0.2000    0.1429    0.1250

Well, we got nine values out, not three, so that certainly isn't what we are looking for. MATLAB computes A(row, col) as the submatrix of A formed by the intersections of the specified rows and columns.

Instead of indexing into A using row and column subscripts, we need to index using a single subscript.

Here's how to convert from row-column subscripts into linear indices:

M = size(A, 1);
idx = M * (col - 1) + row

idx =

     2    11    20

The expression A(idx) pulls out just the three values we want:

A(idx)

ans =

    0.5000    0.3333    0.1250

Or we can assign to just those elements:

A(idx) = Inf

A =

    1.0000    0.5000       Inf    0.2500    0.2000
       Inf    0.3333    0.2500    0.2000    0.1667
    0.3333    0.2500    0.2000    0.1667    0.1429
    0.2500    0.2000    0.1667    0.1429    0.1250
    0.2000    0.1667    0.1429       Inf    0.1111

The functions sub2ind ("subscript" to "index") and ind2sub convert back and forth between row-column subscripts and linear indices.

idx = sub2ind(size(A), row, col)

idx =

     2    11    20

[r, c] = ind2sub(size(A), idx)

r =

     2     1     5


c =

     1     3     4

For a couple of detailed image processing examples illustrating linear indexing, see these previous posts:

Get the MATLAB code

Published with MATLAB® 7.5

Blog reader André recently asked me how to count the number of times each color appears in an image. One short solution involves optional output arguments from the function unique, as well as the relatively new function accumarray.

Let's look first at the function unique. Probably most people use it with a single output argument, like this:

v = [1 1 3 5 5 2 3 2 2]

v =

     1     1     3     5     5     2     3     2     2

u = unique(v)

u =

     1     2     3     5

unique can also return two additional outputs:

[u, m, n] = unique(v);
m

m =

     2     9     7     5

n

n =

     1     1     3     4     4     2     3     2     2

So what do m and n mean? Well, m tells you where you can find the unique values in the original vector. Specifically, m(1) is the index of the last occurrence of the value u(1) in the original vector v. Similarly, m(3) is the index of the last occurrent of the value u(3) in v.

fprintf('m(3) = %d; u(3) = %d, v(m(3)) = %d\n', m(3), u(3), ...
    v(m(3)))

m(3) = 7; u(3) = 3, v(m(3)) = 3

The output argument n gives you indices you can use to reconstruct the original vector from the vector of unique values. In orders, v equals u(n).

u(n)

ans =

     1     1     3     5     5     2     3     2     2

If you count how many times the value 2 appears in the vector n, that will tell you how many times the value u(2) appears in the vector v. That brings us to the second function I mentioned, accumarray. Here's the pertinent sentence from the function description in the doc:

A = accumarray(subs, val) creates an array A by accumulating elements of the vector val using the subscript in subs.

Passing n as subs and 1 as val to accumarray gives us a way to count how many times each value of n appears.

counts = accumarray(n(:), 1)

counts =

     2
     3
     2
     2

counts(1) is the number of times that the value u(1) appears in v. counts(2) is the number of times that the value u(2) appears in v, and so on.

for k = 1:numel(counts)
    fprintf('The value %d appears %d times in v.\n', ...
        u(k), counts(k));
end

The value 1 appears 2 times in v.
The value 2 appears 3 times in v.
The value 3 appears 2 times in v.
The value 5 appears 2 times in v.

So unique and accumarray help to count the number of occurrences of scalars. How do we adapt this technique to count colors? Specifically, how do we count how many times each unique color occurs in a M-by-N-by-3 truecolor image?

First, we reshape the image so that the columns correspond to three color channels.

rgb = imread('peppers.png');
imshow(rgb)

size(rgb)

ans =

   384   512     3

rgb_columns = reshape(rgb, [], 3);
size(column_data)

ans =

      196608           3

Then we use our procedure above with unique and accumarray, except that the we add the 'rows' option to unique.

[unique_colors, m, n] = unique(rgb_columns, 'rows');
color_counts = accumarray(n, 1);

fprintf('There are %d unique colors in the image.\n', ...
    size(unique_colors, 1))

There are 99059 unique colors in the image.

Finally, let's find out what color occurs most frequently.

[max_count, idx] = max(color_counts);

fprintf('The color [%d %d %d] occurs %d times.\n', ...
    unique_colors(idx, 1), unique_colors(idx, 2), ...
    unique_colors(idx, 3), max_count)

The color [254 254 254] occurs 344 times.

And where does this common color occur?

bw = n == idx;
bw = reshape(bw, size(rgb, 1), size(rgb, 2));
imshow(bw)

Get the MATLAB code

Published with MATLAB® 7.5