This is the third post about using lookup tables to process binary images. The first two posts were:

The most fun I've had with binary image lookup tables is playing around with Conway's Game of Life, the most famous example of a cellular automaton. Martin Gardner popularized the Game of Life in his Mathematical Games column in Scientific American:

Of course, there are MANY web sites devoted to this game.

In Conway's Game of Life, each cell in a rectangular grid is either "alive" or "dead" at generation N. Going from one generation to the next, cells "die" or are "born" according to the following rules (as stated in Gardner's first article):

1. Survivals. Every counter [cell] with two or three neighboring counters survives for the next generation.

2. Deaths. Each counter with four or more neighbors dies (is removed) from overpopulation. Every counter with one neighbor or none dies from isolation.

3. Births. Each empty cell adjacent to exactly three neighbors - no more, no fewer - is a birth cell. A counter is placed on it at the next move.

Gardner talks about "counters" because, at the time, computer simulation was simply not available to most readers of Scientific American. In his 1970 article, Gardner explains how to play: "To play life you must have a fairly large checkerboard and a plentiful supply of flat counters of two colors. (Small checkers or poker chips do nicely.) An Oriental "go" board can be used if you can find flat counters that are small enough to fit within its cells. (Go stones are unusable because they are not flat.) It is possible to work with pencil and graph paper but it is much easier, particularly for beginners, to use counters and a board."

OK, back to the present time, when we can use a computer instead of a checkerboard and poker chips. Conway's rules are based on a 3-by-3 neighborhood and so can be implemented using a lookup table. Here's a function which applies Conway's rules to a neighborhood:

type conway

function out = conway(nhood)

live = nhood(2,2) == 1;
P = sum(nhood(:)) - nhood(2,2);  % number of live neighbors

out = (live && ((2 <= P) && (P <= 3))) || ...
   (~live && (P == 3));

Using conway and makelut we can construct a lookup table:

conway_lut = makelut(@conway, 3);

Here's a very small "blinker" pattern. The pattern repeats itself with a period of two generations.

Generation 1:

bw = zeros(5, 5);
bw(3, 2:4) = 1

bw =

     0     0     0     0     0
     0     0     0     0     0
     0     1     1     1     0
     0     0     0     0     0
     0     0     0     0     0

Generation 2:

bw = applylut(bw, conway_lut)

bw =

     0     0     0     0     0
     0     0     1     0     0
     0     0     1     0     0
     0     0     1     0     0
     0     0     0     0     0

Generation 3:

bw = applylut(bw, conway_lut)

bw =

     0     0     0     0     0
     0     0     0     0     0
     0     1     1     1     0
     0     0     0     0     0
     0     0     0     0     0

According to the Wikipedia article on Conway's Game of Life, the first pattern discovered to grow without bound was the "Gosper glider gun." A "glider" is a small configuration of cells that periodically recreates itself one step over. With sufficient imagination, you can see that it "glides." Here is the glider gun:

url = 'http://blogs.mathworks.com/images/steve/2008/gosper_glider_gun.png';
bw = imread(url);
h = imshow(bw, 'InitialMagnification', 'fit')

h =

  349.1982

This pattern will periodically produce new gliders and send them shooting across space.

This code animates the Game of Life board for the Gosper glider gun.

for k = 1:300
    bw = applylut(bw, conway_lut);
    set(h, 'CData', bw);
    pause(0.05);
end

Here's an 18-second video clip showing the glider gun in action.

My last post (I think!) in this series will be about a lookup table optimization we introduced in a recent release of the Image Processing Toolbox. Look for it soon.

Get the MATLAB code

Published with MATLAB® 7.6

This is the second in a short series on using lookup tables for binary image neighborhood operations. See the first post for the basic idea.

The Image Processing Toolbox has two related functions:

Since you have to have a lookup table in order to use applylut, let's start with makelut. It's syntax is:

 lut = makelut(f, n)

f is a function handle, and n is either 2 or 3, depending on whether your lookup table is to be used with 2-by-2 or 3-by-3 neighborhoods.

For 3-by-3 neighborhoods, makelut will call your function handle f 512 times. In each call, makelut will pass as an argument to f one of the 512 possible 3-by-3 binary neighborhoods. f simply has to return the desired output value for a given neighborhood.

For example, suppose we want the output pixel to be 1 if and only if the corresponding input neighborhood looks like this:

  1 0 0
  0 1 0
  0 0 1

The you would call makelut as follows:

pattern = [1 0 0; 0 1 0; 0 0 1];
lut = makelut(@(neighborhood) isequal(neighborhood, pattern), 3);
size(lut)

ans =

   512     1

Let's try that on a small sample image.

bw = eye(5) | fliplr(eye(5))

bw =

     1     0     0     0     1
     0     1     0     1     0
     0     0     1     0     0
     0     1     0     1     0
     1     0     0     0     1

bw2 = applylut(bw, lut)

bw2 =

     0     0     0     0     0
     0     1     0     0     0
     0     0     0     0     0
     0     0     0     1     0
     0     0     0     0     0

So we see that only two input pixels, (2,2) and (4,4), have neighborhoods that exactly match the pattern.

Here's another example: Find all "endpoint" pixels. I'll define an endpoint pixel as a foreground pixel with exactly one foreground neighbor.

endpoint_fcn = @(nhood) (nhood(2,2) ~= 0) && (sum(nhood(:)) == 2);
endpoint_lut = makelut(endpoint_fcn, 3);

Here's the result on sample image from above:

bw3 = applylut(bw, endpoint_lut)

bw3 =

     1     0     0     0     1
     0     0     0     0     0
     0     0     0     0     0
     0     0     0     0     0
     1     0     0     0     1

Only the four corner pixels are endpoints.

In my next post on this topic, I'll use makelut and applylut to play Conway's Game of Life. Look for it in a week or so.

Get the MATLAB code

Published with MATLAB® 7.6

Today I'm introducing a short series about using lookup tables to implement neighborhood operations in binary images. This is a fast implementation technique that works for small neighborhoods, typically 2-by-2 or 3-by-3.

My tentative plan is to do three posts on this algorithm topic:

1. Explain the basic idea (today)

2. Show some examples using Image Processing Toolbox functions

3. Explain an algorithm optimization we recently put in the toolbox.

Let's talk 3-by-3 binary neighborhoods. Here's one:

bw = [0 0 1; 1 1 0; 1 0 0]

bw =

     0     0     1
     1     1     0
     1     0     0

Here's another:

bw = [1 0 0; 1 0 0; 0 1 0]

bw =

     1     0     0
     1     0     0
     0     1     0

OK, that might get old real fast.

QUESTION: How many different 3-by-3 binary neighborhoods are there?

Each pixel in the neighborhood can take on only two values, and there are nine pixels in the neighborhood, so the answer is:

2^9

ans =

   512

That's not so many different possibilities, although it would certainly be tedious to write them all out by hand! Is there an easy way to generate them?

I'm sure there are lots of reasonable methods. One way, used by the Image Processing Toolbox function makelut, is to use the 9-bit binary representation of the integers from 0 to 511.

label = 5;
bin = dec2bin(label, 9)

bin =

000000101

bin is a 9-character string. A logical comparison and a reshape turns it into a binary neighborhood:

bw = reshape(bin == '1', 3, 3)

bw =

     0     0     1
     0     0     0
     0     0     1

Just repeat this for the other 511 integers in the range [0,511] and you'll have all possible 3-by-3 neighborhoods. Compute the output of your operation for each one and save the result in a table.

Then use this procedure to implement your operation:

For each input pixel:
  Extract the 3-by-3 neighborhood
  Compute the table index corresponding to the neighborhood
  Insert the table value at that index into the output image

Computing the table index corresponding to a given neighborhood can be done by using bin2dec. Note that bin2dec takes a string, so we have to do a little work first to convert the binary neighborhood to a string of '0' and '1' characters.

str = repmat('0', 1, 9);
str(bw(:)') = '1';
bin2dec(str)

ans =

     5

And that's the heart of the lookup table algorithm for binary image neighborhood processes.

In my next post on this topic, I'll introduce the Image Processing Toolbox functions makelut and applylut and show examples illustrating their use.

Get the MATLAB code

Published with MATLAB® 7.6

I often find myself writing small functions that help visualize certain image processing algorithms. For example, my last three posts on bwlabel included identical snippets of code that performed the following steps on a label matrix:

1. Display a binary image using two light shades of gray.

2. Use regionprops to find the location of each labeled object.

3. Superimpose object text labels over each labeled object.

The visualization code, although not complicated, was longer than the algorithm code I was trying to explain. That can obscure the key points of the discussion. It would be nice if there were a single visualization function to call.

In the past, we haven't put algorithm visualization functions into the Image Processing Toolbox. They don't really increase what you can do with the toolbox, and there are many different visualization techniques that might be useful for any given algorithm. Even within a single basic technique, such as my bwlabel visualization, there are many possible variations.

Still, algorithm visualizations can be very useful. We've been talking recently about including "example" visualizations in the toolbox. My idea is that each example visualization should have simple code and a simple syntax. We shouldn't succumb to the temptation to overdesign these examples by supporting lots of options and variations, or by overoptimizing them. The primary purpose of the examples would be to provide inspiration and a starting point for our users to create their own visualizations.

This week I captured the visualization code from my recent posts into a function called VISLABELS, which I uploaded to the MATLAB Central File Exchange. Here's how to use it:

help vislabels

 VISLABELS Visualize labels of connected components
    VISLABELS is used to visualize the output of BWLABEL.
 
    VISLABELS(L), where L is a label matrix returned by BWLABEL,
    displays each object's label number on top of the object itself.
 
    Note: VISLABELS requires the Image Processing Toolbox.
 
    Example
    -------
        bw = imread('text.png');
        L = bwlabel(bw);
        vislabels(L)
        axis([1 70 1 70])

And here's what the example does:

bw = imread('text.png');
L = bwlabel(bw);
vislabels(L)
axis([1 70 1 70])

Do you have your own favorite ways of visualizing algorithms? I encourage you to share them on the File Exchange.

Get the MATLAB code

Published with MATLAB® 7.6

Several questions I've seen about bwlabel are about finding the correspondences between object labels in two images. In other words, if a particular pixel is in the foreground in two binary images, is that pixel labeled the same in both images, or is it different? In general, for foreground pixels in the same locations, how do object labels in one binary image match up to object labels in the second?

I'll show one method using the same small binary image that I've been using in my other recent bwlabel posts.

url = 'http://blogs.mathworks.com/images/steve/186/scanned_page.png';
bw = imread(url);
bw = ~bw(1107:1194, 17:135);
imshow(bw, 'InitialMagnification', 'fit')

Compute and display the labeled objects.

L1 = bwlabel(bw);

I = im2uint8(bw);
I(~bw) = 200;
I(bw) = 240;
s = regionprops(L1, 'extrema');
imshow(I, 'InitialMagnification', 'fit')
hold on
for k = 1:numel(s)
   e = s(k).Extrema;
   text(e(1,1), e(1,2), sprintf('%d', k));
end
hold off

Now let's try a morphological closing operation on this image.

bw2 = imclose(bw, ones(3,3));
imshow(bw2, 'InitialMagnification', 'fit')

Compute and display the labeled objects in the second image.

L2 = bwlabel(bw2);

I2 = im2uint8(bw2);
I2(~bw2) = 200;
I2(bw2) = 240;
s2 = regionprops(L2, 'extrema');
imshow(I2, 'InitialMagnification', 'fit')
hold on
for k = 1:numel(s2)
   e2 = s2(k).Extrema;
   text(e2(1,1), e2(1,2), sprintf('%d', k));
end
hold off

Compare the two sets of labels visually. You can see that the broken "f" character on the bottom line was labeled with 6 and 7 in the first image. In the second image, the broken pieces have been merged and labeled with 4.

On the first line of the first image, the "s" and "t" characters are distinct and are labeled with 2 and 4, respectively. In the second image, they have been merged into one object that was labeled 2.

Let's see how to compute all of the correspondences between the object labels in the two images. First, compute the set of pixels that belong to the foreground in both images.

overlap = bw & bw2;

Next, compute the corresponding label pairs by using logical indexing into the two label matrices.

pairs = [L1(overlap), L2(overlap)];

Eliminate the duplicates in the list of pairs.

pairs = unique(pairs, 'rows');

Let's look at a few of the pairs to see how to interpret their meaning:

pairs(1:5, :)

ans =

     1     1
     2     2
     3     3
     4     2
     5     1

Pixels belonging to objects 1 and 5 in the first image belong to a single object, labeled 1, in the second image. Similarly, pixels belonging to objects 2 and 4 in the first image belong to a single object, labeled 2, in the second image. Pixels labeled with a 3 in the first image are also labeled with a 3 in the second image.

Depending on your application, you might want to form an adjacency graph, represented as a sparse matrix.

S = sparse(pairs(:,1), pairs(:,2), 1);
spy(S)

One way to use the graph is to find all the object labels in image 1 that correspond to a given label in image 2. For example, here are all the first-image objects corresponding to the object labeled with 6 in the second image:

find(S(:, 6))

ans =

     9
    10

My other recent bwlabel posts were about search order and relabeling. You can also look the archive of all my posts about connected-component labeling.

Get the MATLAB code

Published with MATLAB® 7.6

The three most common questions I've been hearing about bwlabel, which is used for labeling connected components in a binary image, are about

Today I'll tackle relabeling.

In my previous post on search order, I showed how to postprocess the labeled objects to modify their order according to various criteria. Now I'll take that a step further and use the sorted output from regionprops to renumber the labels in the label matrix.

Blog reader Trung wanted to know if we could sort lexicographically by centroid, so I'll do that here.

url = 'http://blogs.mathworks.com/images/steve/186/scanned_page.png';
bw = imread(url);
bw = ~bw(1107:1194, 17:135);
imshow(bw, 'InitialMagnification', 'fit')

Here's a false-color view of the label matrix using label2rgb:

L = bwlabel(bw);
rgb = label2rgb(L, 'jet', [.95 .95 .95], 'shuffle');
imshow(rgb, 'InitialMagnification', 'fit')

Now let's use regionprops to get the centroids for each object. To facilitate the relabeling step, I'll get the 'PixelIdxList' for each object.

s = regionprops(L, {'Centroid', 'PixelIdxList'});

Next, sort lexicographically by the centroids, sorting first by the vertical coordinate.

centroids = cat(1, s.Centroid);
[sorted_centroids, sort_order] = sortrows(fliplr(centroids));
s2 = s(sort_order);

To visualize the sorted object order, we can display the object numbers on top of the objects like this:

% First, make an image with a light gray background instead
% of a black background, so that the numbers will be visible
% on top of it.
I = im2uint8(bw);
I(~bw) = 200;
I(bw) = 240;
imshow(I, 'InitialMagnification', 'fit')

% Now plot the number of each sorted object at the corresponding
% centroid:
hold on
for k = 1:numel(s2)
   centroid = s2(k).Centroid;
   text(centroid(1), centroid(2), sprintf('%d', k));
end
hold off

We sorted the output of regionprops, but we haven't touched the label matrix itself. We can relabel it according to the sort order by using linear indexing and the PixelIdxList of each object. An object's PixelIdxList is a vector of linear indices for the pixels belonging to the object. Here's the relabeling loop:

for k = 1:numel(s2)
   kth_object_idx_list = s2(k).PixelIdxList;
   L(kth_object_idx_list) = k;
end

In MATLAB 7.6 (R2008a), the publish feature lets you capture multiple graphics from within a loop. I'll use that new capability to show the location of the first few relabeled objects in L.

for k = 1:5

   imshow(L == k)

end

Soon I'll write another post that shows how how to match up labels in objects that overlap between two images.

Get the MATLAB code

Published with MATLAB® 7.6

Blog reader Ram asked a question last week that I hear fairly often: When you apply a spatial transformation to an image, how can you see the x-y coordinates when you display the image?

The answer has three steps:

1. Ask imtransform for the spatial coordinates of the output image

2. Give imshow the spatial coordinates when you display the output image

3. Turn on the axis box and tick labels

Let's walk through these steps using an affine transformation.

T = maketform('affine',[.5 0 0; .5 1.5 0; 100 200 1]);
I = imread('cameraman.tif');

The first output from imtransform is the transformed image. There are two more optional arguments, though: xdata and ydata. xdata is a two-element vector specifying the x-coordinates of the first and last image columns. Similarly, ydata specifies the y-coordinates of the first and last image rows.

[I2, xdata, ydata] = imtransform(I, T);
xdata

xdata =

   101   356

ydata

ydata =

  201.5000  584.5000

Step 2 is to pass the coordinate information to imshow. You can do that by using the 'XData' and 'YData' parameters, like this:

imshow(I2, 'XData', xdata, 'YData', ydata)

But we still can't see the coordinates! That's because imshow turns the axis display off. That brings us to step 3: Turn on the axis box and tick labels:

axis on

After you do these steps, you might also want to call impixelinfo, which turns on a display of pixel location and value in the figure window. The display updates as you move the mouse.

Get the MATLAB code

Published with MATLAB® 7.6