This is the first post in a short series on index techniques that are particularly useful for image processing in MATLAB. I'll start with logical indexing today. Later I'll cover linear indexing, and then a technique I like to call neighbor indexing.

Every MATLAB user is familiar with ordinary matrix indexing notation.

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

A(2,3)

ans =

     7

A(2,3) extracts the 2nd row, 3rd column of the matrix A. You can extract more than one row and column at the same time:

A(2:4, 3:5)

ans =

     7    14    16
    13    20    22
    19    21     3

When an indexing expression appears on the left-hand side of the equals sign, it's assignment:

A(5,5) = 100

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   100

About every 13.6 days, someone asks this question on comp.soft-sys.matlab:

 How do I replace all the NaNs in my matrix B with 0s?

This is generally followed 4.8 minutes later with this reply from one of the newsgroup regulars:

 B(isnan(B)) = 0;

For example:

B = rand(3,3);
B(2, 2:3) = NaN

B =

    0.5470    0.1890    0.3685
    0.2963       NaN       NaN
    0.7447    0.1835    0.7802

Replace the NaNs with zeros:

B(isnan(B)) = 0

B =

    0.5470    0.1890    0.3685
    0.2963         0         0
    0.7447    0.1835    0.7802

The expression

 B(isnan(B))

is an example of logical indexing. Logical indexing is a compact and expressive notation that's very useful for many image processing operations.

Let's talk about the basic rules of logical indexing, and then we'll reexamine the expression B(isnan(B)).

If C and D are matrices, then C(D) is a logical indexing expression if C and D are the same size, and D is a logical matrix.

"Logical" is one of the builtin types, or classes, of MATLAB matrices. Relational operators, such as == or >, produce logical matrices.

C = hilb(4)

C =

    1.0000    0.5000    0.3333    0.2500
    0.5000    0.3333    0.2500    0.2000
    0.3333    0.2500    0.2000    0.1667
    0.2500    0.2000    0.1667    0.1429

D = C > 0.4

D =

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

whos D

  Name      Size            Bytes  Class      Attributes

  D         4x4                16  logical              

You can see from the output of whos that the class of the variable D is logical. The logical indexing expression C(D) extracts all the values of C corresponding to nonzero values of D and returns them as a column vector.

C(D)

ans =

    1.0000
    0.5000
    0.5000

Now we know enough to break down the B(isnan(B)) example to see how it works.

B = rand(3,3);
B(2, 2:3) = NaN;

nan_locations = isnan(B)

nan_locations =

     0     0     0
     0     1     1
     0     0     0

whos nan_locations

  Name               Size            Bytes  Class      Attributes

  nan_locations      3x3                 9  logical              

B(nan_locations)

ans =

   NaN
   NaN

B(nan_locations) = 0

B =

    0.0811    0.4868    0.3063
    0.9294         0         0
    0.7757    0.4468    0.5108

Functions in the Image Processing Toolbox, as well as the MATLAB functions imread and imwrite, follow the convention that logical matrices are treated as binary (black and white) images. For example, when you read a 1-bit image file using imread, it returns a logical matrix:

bw = imread('text.png');
whos bw

  Name        Size             Bytes  Class      Attributes

  bw        256x256            65536  logical              

This convention, together with logical indexing, makes it very convenient and expressive to use binary images as pixel masks for extracting or operating on sets of pixels.

Here's an example showing how to use logical indexing to compute the histogram of a subset of image pixels. Specifically, given a gray-scale image and a binary segmentation, compute the histogram of just the foreground pixels in the image.

Here's our original image:

I = imread('rice.png');
imshow(I)

Here's a segmentation result (computed and saved earlier), represented as a binary image:

url = 'http://blogs.mathworks.com/images/steve/192/rice_bw.png';
bw = imread(url);
imshow(bw)

Now use the segmentation result as a logical index into the original image to extract the foreground pixel values.

foreground_pixels = I(bw);
whos foreground_pixels

  Name                       Size            Bytes  Class    Attributes

  foreground_pixels      17597x1             17597  uint8              

Finally, compute the histogram of the foreground pixels.

imhist(foreground_pixels)

Or use logical indexing with the complement of the segmentation result to compute the histogram of the background pixels.

imhist(I(~bw))

Get the MATLAB code

Published with MATLAB® 7.5

Toward the end of last year, it occurred to me that I had accumulated enough material on this blog to construct a "master class" on image processing in MATLAB. The response to my suggestion was "Sounds good - here's a list of places we want you to go!"

So in a perverse way, it's the blog's fault that I've been slow posting recently. I've been frantically preparing a new image processing seminar, which I'll be presenting at the University of Illinois at Urbana-Champaign and at MIT later this week. (Please, no more snow storms for a few days!)

Have you ever heard of the linkaxes function? We introduced it to MATLAB several releases ago, sometime around version 7.0. It lets you plot into two different axes objects and then keep those axes automatically synchronized with each other. For example, as you zoom and pan in one axes, the other axes zooms and pans the same way. This is really useful for comparing input and output images.

Here's an example using linkaxes to examine closely the output of the edge function.

I = imread('rice.png');
imshow(I)

bw = edge(I, 'canny');
imshow(bw)

Let's use subplot to create two axes in the same figure, and then we'll link them together.

% Call subplot with an output argument to capture the axes
% handles, because we'll need the handles in the call to
% linkaxes.
ax1 = subplot(1, 2, 1);
imshow(I)

ax2 = subplot(1, 2, 2);
imshow(bw)

Now pass linkaxes a vector containing both axes handles.

linkaxes([ax1 ax2]);

If you then pan and zoom using the mouse in one image, the other image follows along. The linking also works if you change the axes limits from the command line, like this:

axis([50 90 200 240])

There's your quick tip for the week. I hope to pick up the blogging pace again starting next week.

Get the MATLAB code

Published with MATLAB® 7.5

Have you ever used the distance transform?

For a binary image, the distance transform is the distance from every pixel to the nearest foreground (nonzero) pixel. (Sometimes you'll see it defined the other way around. It doesn't really matter that much; you just need to pay attention to whatever convention is being used and complement the image as needed.)

Last week I posted a method for cleaning up scanned text. The method distinguished between small dots that were far away from or were close to pixels belonging to large characters. In that post I used dilation and logical operators to identify small dots that were far away from characters. Later it occurred to me that one could also use the distance transform for this purpose.

Here again are the first few steps. Use bwareaopen to remove the small dots, and use a logical operator to identify the pixels that got removed.

url = 'http://blogs.mathworks.com/images/steve/186/scanned_page.png';
page = imread(url);
bw = page(735:1280, 11:511);
imshow(bw)

bw2 = imcomplement(bw);
bw3 = bwareaopen(bw2, 8);
removed = xor(bw2, bw3);

Next, use the distance transform to identify all pixels that are within a certain distance from foreground pixels in the image bw3.

D = bwdist(bw3);
within_hailing_distance = D <= 10;
imshow(within_hailing_distance)

Which removed pixels do we want to put back?

put_back_pixels = removed & within_hailing_distance;

Use a logical OR to put the pixels back.

gotta_go_the_holiday_party_is_about_to_start = bw3 | put_back_pixels;
imshow(~gotta_go_the_holiday_party_is_about_to_start)

I suspect that not many people know about the "extra" feature that's tucked inside bwdist. Not only can bwdist compute the distance transform for you, but it can also compute a result sometimes called the "feature transform." For each pixel, the feature transform tells you which foreground pixel is nearest. You get the feature transform simply by using a second output argument when you call bwdist.

Do you use this capability? Or do you think you might have a use for it? Please let me know. I'd love to hear about it.

Get the MATLAB code

Published with MATLAB® 7.5

Earlier this year I exchanged e-mail with blog reader Craig Doolittle. Craig was writing MATLAB scripts to clean up scanned pages from old manuscripts. One of the samples he sent me was a page from "Fragmentation of Service Projectiles," N.F. Mott, J. H. Wilkinson, and T.H. Wise, Ministry of Supply, Armament Research Department, Theoretical Research Report No. 37/44, December 1944.

The image is too big to show at full resolution in this blog, so here's a thumbnail view.

url = 'http://blogs.mathworks.com/images/steve/186/scanned_page.png';
page = imread(url);
thumbnail = imresize(im2uint8(page), 'OutputSize', [256 NaN]);
imshow(thumbnail)

Craig wanted suggestions on how to clean up isolated "noise" dots without removing small dots that are part of characters. Let's look closely at a cropped portion of the page.

bw = page(735:1280, 11:511);
imshow(bw)

We could start by using bwareaopen to remove small dots. For For example:

bw2 = imcomplement(bw);
bw3 = bwareaopen(bw2, 8);
imshow(imcomplement(bw3))

Unfortunately, this approach has removed portions of some of the characters. Here's a method using bwlabel and regionprops to highlight the pixels that were removed.

removed = xor(bw2, bw3);
L = bwlabel(removed);
s = regionprops(L, 'Centroid');
centroids = cat(1, s.Centroid);
imshow(bw)
hold on
plot(centroids(:,1), centroids(:,2), 'ro')
hold off

You can see that some of the removed dots were noise, while others were parts of the characters "e", "i", "m", etc.

My suggestion to Craig was to restore removed dots that are "close" to the characters remaining after bwareaopen. We can do this using dilation and some logical operators.

bw4 = imdilate(bw3, strel('disk', 5));
imshow(bw4)

Now do a logical AND of the dilated characters with the pixels removed by bwareaopen. These are the pixels we are going to put back.

overlaps = bw4 & removed;
imshow(overlaps)

Use a logical OR to restore the removed pixels.

bwout = imcomplement(bw3 | overlaps);
imshow(bwout)

I also suggested using morphological recontruction to get all the pixels connected to the overlapping pixels found above. It doesn't seem to be really necessary here, though, so I'm going to save this technique for a future blog post, using a better example.

I'm sure there are lots of different ways to approach this text clean-up problem. Does anyone have suggestions for other approaches?

Thanks for letting me use your example, Craig.

Get the MATLAB code

Published with MATLAB® 7.5

The notion of neighbor connectivity is discussed in most image processing textbooks. Specifically, what is the set of neighbors of a pixel? For example, a commonly-used neighborhood connectivity is 4-connected, where each pixel has four neighbors. The neighborhood looks like this:

   1
 1 1 1
   1

Another common connectivity is 8-connected, where each pixel has eight neighbors. The neighborhood looks like this:

 1 1 1
 1 1 1
 1 1 1

Many Image Processing Toolbox functions depend on a definition of connectivity. Most of these functions allow you to specify the desired connectivity in a very general, flexible way.

Consider bwperim. This function computes perimeter pixels by finding the foreground pixels that are connected to background pixels. Some definition of connectivity is therefore required. Let's look at an example.

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

Find the foreground pixels that are 4-connected to the background.

perim4 = bwperim(bw, 4)

perim4 =

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

Now find the foreground pixels that are 8-connected to the background.

perim8 = bwperim(bw, 8)

perim8 =

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

In addition to specifying connectivity as "4" or "8", you can also use a 3-by-3 matrix of 0s and 1s to specify other kinds of connectivity. The 3-by-3 connectivity matrix below says that a pixel has only two neighbors: The pixel above and the pixel below.

conn2 = [ 0 1 0
          0 1 0
          0 1 0 ];

perim2 = bwperim(bw, conn2)

perim2 =

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

Note that the (4,2) and (4,6) pixels are not considered to be perimeter pixels with this connectivity. That's because their upper and lower neighbors are both part of the foreground, not the background.

Sometimes it's useful to define 6-connectivity. There are two different flavors, both of which can be expressed using the 3-by-3 matrix form:

conn6a = [ 1 1 0
           1 1 1
           0 1 1 ];

conn6b = [ 0 1 1
           1 1 1
           1 1 0 ];

For processing three-dimensional arrays, you specify 6-connectivity (face-adjacent pixels), 18-connectivity (face- or edge-adjacent pixels), or 26-connectivity (face-, edge-, or vertex-adjacent pixels). But you can specify your own flavor of connectivity by using a 3-by-3-by-3 matrix. Here's an example.

bw_3d = false(5,5,3);
bw_3d(2:4,2:4,:) = true;

% Label three-dimensional connected components.
L1 = bwlabeln(bw_3d)

L1(:,:,1) =

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


L1(:,:,2) =

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


L1(:,:,3) =

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

There's only one three-dimensional connected component. But that's using the default connectivity. Let's define a three-dimensional connectivity such that each pixel has neighbors only in the same plane.

conn_3d = cat(3, zeros(3,3), ones(3,3), zeros(3,3))

conn_3d(:,:,1) =

     0     0     0
     0     0     0
     0     0     0


conn_3d(:,:,2) =

     1     1     1
     1     1     1
     1     1     1


conn_3d(:,:,3) =

     0     0     0
     0     0     0
     0     0     0

L2 = bwlabeln(bw_3d, conn_3d)

L2(:,:,1) =

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


L2(:,:,2) =

     0     0     0     0     0
     0     2     2     2     0
     0     2     2     2     0
     0     2     2     2     0
     0     0     0     0     0


L2(:,:,3) =

     0     0     0     0     0
     0     3     3     3     0
     0     3     3     3     0
     0     3     3     3     0
     0     0     0     0     0

Now there are three connected components, one in each plane of the array.

The Image Processing Toolbox functions bwlabel, bweuler, bwboundaries, and bwtraceboundary can accept either 4 or 8 connectivity.

The toolbox functions bwareaopen, bwlabeln, bwperim, bwulterode, imclearborder, imextendedmax, imextendedmin, imhmax, imhmin, imimposemin, imreconstruct, imregionalmax, imregionalmin, and watershed can all accept general 3-by-3 (or higher dimensional) connectivity matrices.

Have you used this general connectivity capability before? If so, I'd like to hear about it.

Get the MATLAB code

Published with MATLAB® 7.5