Steve on Image Processing

It's hard for me to believe, but I started writing this blog just over five years ago. I'm going to start occasionally recapping old blog posts, especially those with content that is still useful today and that more recent readers may not have seen.

In January 2006 I showed how to generate useful test images with some simple mathematical computations.

I also highlighted the MATLAB "publish" feature. This capability lets you easily publish simple MATLAB scripts as HTML, Word, PDF, LaTeX, and other types of documents. Command Window output and generated figures are captured automatically into the output. I use "publish" to write almost all of my blog posts.

I had just read Steve Squyres' book Roving Mars: Spirit, Opportunity, and the Exploration of the Red Planet, and I summarized its description of DIMES, the Descent Image Motion Estimation System. I didn't discover until later how MathWorks products had been used in the design of this and other Mars Exploration Rover systems. (See the June 2004 issue of MATLAB News and Notes.)

I started a long-running series of posts on the spatial (or geometric) transformation of images.

And I started another series of posts to explain all the different ways in which array values can get turned into image pixel colors on the screen. In 2006 we edited this series into a MATLAB Digest article called "How MATLAB Represents Pixel Colors."

In February 2006 I continued my two series on spatial transformations and pixel colors. (Later on I learned not to start two series at the same time.)

I described how developers often don't understand how long it takes for our users to become familiar with new features in the product. Partially this is because the "developer calendar" is shifted quite a lot from the user's experience. By the time the earliest downloaders of R2011a had installed that release, we were wrapping up feature work for R2011b and were well into planning for R2012a.

Finally, I exposed the limits of my artistic ability with the post called "Tracing George." (George is a person that all Georgia Tech students and alums know fondly.) This post demonstrated thresholding, morphological thinning, and boundary tracing.

Get the MATLAB code

Published with MATLAB® 7.12

A few months ago, MATLAB user Julia asked how to "fill 4-connected pixels exclusively" on MATLAB Answers. Here's a more precise statement of the problem: if a zero-valued pixel in a binary image has one-valued north, east, south, and west neighbors, then change that zero-valued pixel to a one-valued pixel.

I suggested bwhitmiss in my answer. The function bwhitmiss is a very useful tool for searching for specific pixel neighborhood configurations in a binary image. Here's how it works.

First, create one structuring element (as a matrix of zeros and ones) representing neighborhood pixels you want to "hit." For example, if you want the north, east, south, and west neighbors in a 3-by-3 neighborhood to be one-valued, use this structuring element:

se1 = [0 1 0
       1 0 1
       0 1 0]

se1 =

     0     1     0
     1     0     1
     0     1     0

Now create a second structuring element representing neighborhood pixels you want to "miss." For example, if you want the center pixel of a 3-by-3 neighborhood to be zero-valued, use this structuring element:

se2 = [0 0 0
       0 1 0
       0 0 0]

se2 =

     0     0     0
     0     1     0
     0     0     0

Then you'd pass your binary image and both of these structuring elements to bwhitmiss. Let's construct a small test image to experiment with.

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

bw =

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

pixel_matches = bwhitmiss(bw,se1,se2)

pixel_matches =

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

You can see that bwhitmiss found two places where a zero-valued pixel had one-valued north, east, south, and west neighbors. To modify the original image by turning these zero-valued pixels into one-valued pixels, use an elementwise OR operator:

bw2 = bw | pixel_matches

bw2 =

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

Sometimes you'll see references to "don't care" neighbors in discussions about the hit-miss operator. These are neighborhood pixels that can be either zero or one without affecting the match. In terms of the inputs to bwhitmiss, the "don't care" neighbors are the ones where neither se1 nor se1 is equal to one.

dont_care_neighbors = ~(se1 | se2)

dont_care_neighbors =

     1     0     1
     0     0     0
     1     0     1

Here's one more example to finish with. Let's use bwhitmiss to identify isolated pixels. These are one-valued pixels that are completely surrounded by zero-pixels.

se1 = [0 0 0
       0 1 0
       0 0 0];

se2 = [1 1 1
       1 0 1
       1 1 1];

bw3 = bwhitmiss(bw,se1,se2)

bw3 =

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

There was only one isolated pixel in our little sample image.

Get the MATLAB code

Published with MATLAB® 7.12

Blog reader Mutlu recently asked me this question:

"I have two images - one contains class labels (one of 1 or 2) and one containing image segment (region) unique IDs. I am trying to use regionprops (or something like regionprops) to figure out most frequently occurring label from the label image per region and generate a new image where each unique region contains this most frequently occurring label. Is there a way to do this?"

The answer is yes. I thought the techniques might be applicable to different kinds of problems, so I'm posting the description here for everyone.

I don't have a specific example from the Mutlu's own application, so I made up a sample image to illustrate.

url = 'http://blogs.mathworks.com/images/steve/2011/regions-classes-example.png';
bw = imread(url);
imshow(bw)

Suppose this image is divided up into three horizontal regions, as defined by this label matrix:

regions_url = 'http://blogs.mathworks.com/images/steve/2011/regions.png';
regions = imread(regions_url);
imshow(label2rgb(regions, 'jet', [.5 .5 .5]))

Further, suppose that the shapes in the original image are divided into three classes - the stars, the triangles, and the circles - as defined by this label matrix:

classes_url = 'http://blogs.mathworks.com/images/steve/2011/classes.png';
classes = imread(classes_url);
imshow(label2rgb(classes, 'jet', [.5 .5 .5]))

I can imagine different ways toward the solution. Here's the rough procedure I'm going to follow:

1. Determine the connected components (objects) in the image, including a list of pixel locations for each one.

2. Determine the class of each connected component.

3. Determine the region of each connected component.

4. For each region, determine the mode of the connected component classes contained in that region.

5. Construct the output image containing only the most frequently occurring object classes in each region.

Step 1: Determine the connected components (objects) in the image, including a list of pixel locations for each one. (I'm also going to compute the centroids for each object, but that's only because I'm going to use those for visualization steps along the way.)

cc = bwconncomp(bw);
s = regionprops(cc, 'PixelIdxList', 'Centroid');

Step 2: Determine the class of each connected component. I'll do this by using the pixel index list for each connected component to index into the classes label matrix. I'll tack on the class number to the existing structure array returned by regionprops; this is a pretty useful technique to remember.

for k = 1:cc.NumObjects
    s(k).ClassNumber = classes(s(k).PixelIdxList(1));
end

Let's plot the class numbers on top of the original image to make sure we did it right. We should have 1s for the stars, 2s for the triangles, and 3s for the circles.

imshow(bw)
hold on
for k = 1:numel(s)
    x = s(k).Centroid(1);
    y = s(k).Centroid(2);
    text(x, y, sprintf('%d', s(k).ClassNumber), 'Color', 'r', ...
        'FontWeight', 'bold');
end
hold off
title('Class number of each object')

Step 3: Determine the region number associated with each connected component. We have to think a little harder about this one, because there's a possibility that a connected component could lie in more than one region. I suggest that we find the region number of every pixel of each connected component and then use the mode function as a tiebreaker for components that lie in multiple regions.

for k = 1:cc.NumObjects
    s(k).RegionNumber = mode(single(regions(s(k).PixelIdxList)));
end

Once again, let's plot the region numbers as a sanity check.

imshow(bw)
hold on
for k = 1:numel(s)
    x = s(k).Centroid(1);
    y = s(k).Centroid(2);
    text(x, y, sprintf('%d', s(k).RegionNumber), 'Color', 'r', ...
        'FontWeight', 'bold');
end
hold off
title('Region number of each object')

Step 4: For each region, determine the most frequently occurring object class in that region. Again, we'll make use of the mode function.

num_regions = max(regions(:));
most_frequent_class = zeros(1, num_regions);
region_vector = [s.RegionNumber];
class_vector = [s.ClassNumber];
for k = 1:num_regions
    most_frequent_class(k) = mode(single(class_vector(region_vector == k)));
end

most_frequent_class

most_frequent_class =

     1     3     2

Step 5: Construct the output image containing only the most frequently occurring object classes in each region.

bw2 = false(size(bw));
for k = 1:numel(s)
    if most_frequent_class(s(k).RegionNumber) == s(k).ClassNumber
        % The class number of this object matches the most frequent class
        % number of objects in this region.
        bw2(s(k).PixelIdxList) = 1;
    end
end

imshow(bw2)
title('Final result')

You can see that only stars are in region 1, only circles are in region 2, and only triangles are in region 3.

Here's the complete computation in one block of code:

cc = bwconncomp(bw);
s = regionprops(cc, 'PixelIdxList', 'Centroid');

for k = 1:cc.NumObjects
    s(k).ClassNumber = classes(s(k).PixelIdxList(1));
end

for k = 1:cc.NumObjects
    s(k).RegionNumber = mode(single(regions(s(k).PixelIdxList)));
end

num_regions = max(regions(:));
most_frequent_class = zeros(1, num_regions);
region_vector = [s.RegionNumber];
class_vector = [s.ClassNumber];
for k = 1:num_regions
    most_frequent_class(k) = mode(single(class_vector(region_vector == k)));
end

bw2 = false(size(bw));
for k = 1:numel(s)
    if most_frequent_class(s(k).RegionNumber) == s(k).ClassNumber
        % The class number of this object matches the most frequent class
        % number of objects in this region.
        bw2(s(k).PixelIdxList) = 1;
    end
end

The entire computation takes about 11 ms on my laptop.

There were a lot of useful techniques demonstrated here, especially:

Get the MATLAB code

Published with MATLAB® 7.12

R2011a shipped back in April, but I haven't said too much about it, yet. I mentioned the new geotiffwrite function in Mapping Toolbox (18-Apr-2011 post), and I discussed the improved performance of certain kinds of automatic array growth in MATLAB (16-May-2011 and 20-May-2011 posts).

For new MATLAB releases, I usually like to highlight new features that are the most interesting to me personally, and that's what I want to do today. For additional information and perspectives on the release, I encourage you to check out the MATLAB R2011a release notes, Loren on the Art of MATLAB, and Mike on the MATLAB Desktop.

When I scan the MATLAB R2011a release notes, the items of most interest to me come from the MATLAB Math, Audio/Video, and Image & Scientific Formats teams. (I don't mean to overlook the work of other teams! These just happen to be some of the areas I know best.)

The MATLAB Math team has overhauled the interface for controlling the MATLAB random number generator. This was in response to persistent feedback that the previous interface was hard to understand and difficult to use. (I got to be an observer on some of the usability test sessions.) See the documentation for the new function rng for details. I'll tell Loren to get Peter to make a new guest blog appearance about rng.

The Math team always likes to make a bunch of stuff go faster, and this release was no different. Newly optimized functions include matrix transpose, many elementwise functions operating on single-precision inputs, several linear algebra functions, long-vector and large-matrix convolution with conv and conv2, and indexed assignment into sparse arrays.

The Audio/Video team added Motion JPEG support to VideoWriter.

The Image & Scientific Formats team added a suite of new functions that make it easier to work with NetCDF and HDF5 files. You can also now read in a subset of the pixels in FITS file.

I will probably discuss the new NetCDF and HDF5 functions in more detail in a future post.

Finally, I thought I would remind you that MathWorks now provides documentation for earlier releases on its web site. (You have to be a licensed user to access documentation for the earlier releases.) Start at mathworks.com, click on "Support," then "Product Documentation," and then "View documentation for other releases."

Get the MATLAB code

Published with MATLAB® 7.12

Earlier this week I wrote about the improvements in automatic array growth in MATLAB R2011a. The posted comments made me realize that I should give a few more details so that you can understand all the circumstances under which this improvement might benefit you.

To illustrate the improvement with a slightly different code fragment than the for-loop I used the other day:

  y = [];
  while not_done_yet()
      if some_condition()
          y(end+1) = next_value();
      end
  end

We don't know in advance how big y is going to be, so it's harder to preallocate space for it. It would really be nice if I could just write the above code and let MATLAB handle all the boring details automatically efficiently.

Well, prior to R2011a, MATLAB would all handle all the boring details of array growth automatically but not efficiently. Now it also handles array growth efficiently.

Here are some more details you might find helpful.

Also, for numeric array types, certain kinds of growth by concatenation are supported. Specifically, the case where a numeric array is concatenated with itself is now optimized:

   x = [];
   for k = 1:n
       x = [x k];
   end

One blog reader wondered if this optimization means that there might be situations where some very large vectors or arrays might be using twice as memory as necessary.

The answer is no, not really. MATLAB uses a smarter heuristic than simply doubling the allocated memory space whenever more is needed, so for large arrays the worst-case memory "overallocation" is much less than a factor of two. I don't intend to get into further details here because (a) I don't know them, and (b) I expect that we will continue to tune the heuristic and other aspects of automatic array growth with future releases.

I'll conclude by repeating my recommendation from earlier this week. If it's trivial to precalculate the final size of a growing array before entering a loop, then by all means use preallocation. It's the fastest way to go. But if it's not trivial, and if the pattern of growth matches one of the cases described above, then you can consider just letting MATLAB handle the array growth for you without worrying about bad performance scaling.

Get the MATLAB code

Published with MATLAB® 7.12

One of the most significant MATLAB improvements in the R2011a release is in how certain kinds of array growth operations are handled. This is a kind of "under-the-hood" change that I think is exciting and deserves some publicity.

For those of you who are MATLAB power users and understand preallocation and why it is necessary, I'll tell you the story's ending first: MATLAB R2011a performs much better than previous MATLAB versions when performing repeated array growth in a loop.

For everyone else, let's start with some not-so-hypothetical tech support questions:

The tech support engineer will start digging into the question to see why the unhappy MATLAB user has reached his conclusion:

"Well, I ran the profiler on my function flibbergibble, which calls foobar:

"You can see that foobar takes about 0.35 ms per call (70.93 / 200000), but when I call foobar 200000 times from the command line in a loop it only takes 8.5e-7 seconds per call:

 tic, for k = 1:200000, foobar(k); end, toc
 Elapsed time is 0.173287 seconds.

"That's more than 400 times faster! Why is foobar so slow in a function?"

Experienced MATLAB users know that the performance problem has little to do with foobar(k) on line 4 of the profile report and much to do with y(k) = on the same line.

Let's examine why. The first time through the loop, the variable y hasn't been assigned to yet, so MATLAB creates it through the assignment to y(1). The second time through the loop assigns to y(2), which doesn't exist yet. So MATLAB allocates enough new memory for two elements and then copies y(1) and assigns y(2) into the new space. The third time through the loop assigns to y(3), which doesn't exist yet. So MATLAB allocates enough new memory for three elements and then copies y(1) and y(2) and assigns y(3) into the new space.

See the pattern? On the k-th time through the loop, MATLAB has to make a copy of k elements into newly allocated memory. The time required to copy k elements is proportional to k. The time required to execute the entire loop is therefore proportional to the sum of the integers from 1 to n, which is n*(n+1)/2.

Bottom line: memory copying during the assignment statement causes the time required to execute the loop to be proportional to n^2.

Let's see some data. I timed flibbergibble for various values of n using MATLAB R2010b. Here's the data:

n = [100 200 500 1000 2000 10000 20000 40000 50000];
t = [1.2e-4 2.7e-4 8.6e-4 2.4e-3 2.8e-3 5.8e-2 .24 1.32 2.4];

plot(n,t)
title('n versus t, MATLAB R2010b')

That looks pretty parabola-like. I checked it out with cftool from the Curve Fitting Toolbox and found that a rough fit to the data is the curve 1.0e-9 * n^2. What if we tried n = 1000000? From our curve we could estimate that it would take about 17 minutes:

n = 1e6;
t_minutes = 1e-9 * n^2 / 60

t_minutes =

   16.6667

This is what triggers the tech support call where the user thinks that the function foobar has caused MATLAB to hang. MATLAB isn't hung, it's just working really hard reallocating all that memory.

The solution is to preallocate space for y before beginning the loop, like this:

type flibbergibble2

function y = flibbergibble2(n)

y = zeros(1, n);
for k = 1:n
    y(k) = foobar(k);
end

The code above does not have to reallocate any memory each time through the loop. In MATLAB R2010b on my computer it runs in less than a second:

 tic, flibbergibble2(1000000); toc
 Elapsed time is 0.737451 seconds.

The big change in MATLAB R2011a is that MATLAB uses a better strategy for reallocating memory as the array grows when it detects this growth pattern. As a result, the time required to execute the loop grows only linearly with n.

Here's how it takes for n = 1000000 in MATLAB R2011a on my computer:

 tic, flibbergibble(1000000); toc
 Elapsed time is 1.105256 seconds.

That's in the neighborhood of a thousand times faster than R2010b. You can see that the preallocation strategy (in flibbergibble2) is still a bit faster, but the performance penalty for not preallocating has been substantially reduced.

If the final size of your output array is trivial to calculate, then I would say go ahead and preallocate. It's just one extra code line and it's a little faster. But if the final size is complicated to calculate, then maybe you don't have to work as hard you did before to optimize your code. If you are growing a vector, or if the array growth is in the final dimension of a matrix or multidimensional array, you can now seriously consider just letting MATLAB do its own "automatic array growth" thing.

Get the MATLAB code

Published with MATLAB® 7.12

I just saw a very cool demo that I thought you might like. Teja Muppirala, an application engineer in our Japan office, created a tool that uses MATLAB, Image Processing Toolbox, Image Acquisition Toolbox, and a webcam to solve printed Sudoku puzzles. Here you can see Teja holding up an issue of MATLAB News and Notes in front of his webcam:

And here you can see a MATLAB figure showing that the puzzle has been read and solved:

Teja makes nice use of Image Processing Toolbox functions for block processing, binary morphology, region measurement, and geometric transforms. And he's got the whole thing working for streaming input, not just for a static image capture.

Take a look at the video Teja made. You can find all the code for the demo on the MATLAB Central File Exchange.

Nice work, Teja!

Get the MATLAB code

Published with MATLAB® 7.12

Some of you may have already noticed a new product in R2011a: Computer Vision System Toolbox. You may also have noticed that Video and Image Processing Blockset disappeared! These are not unrelated, as we took the blockset, added new computer vision algorithms, and changed the name to Computer Vision System Toolbox. Here’s a list of what’s new:

Calling this product a toolbox also allows us to clarify and highlight the MATLAB capabilities in the product that we’ve had since R2010a. Some of these algorithms overlap with Image Processing Toolbox, but provide support for C code generation and fixed-point modeling. Others are unique to Computer Vision System Toolbox, including these:

So, now we have one product for use in both MATLAB and Simulink that supports the design and simulation of computer vision and video processing systems. It contains MATLAB functions, MATLAB System objects, and Simulink blocks. You can learn more about these capabilities by looking at the documentation for Computer Vision System Toolbox.

Get the MATLAB code

Published with MATLAB® 7.12

This is my first post using R2011a, which shipped last week. (MathWorks ships product-line updates twice per year, usually in March and September.) There's a lot of interesting material to write about for this new release.

Because my blog readers have been asking for it for a long time, the first thing I want to mention is the new Mapping Toolbox function geotiffwrite.

Here are a couple of examples from the documentation. The first example reads a layer from a NASA WMS (Web Map Service) server and writes it out to a GeoTIFF file.

nasa = wmsfind('nasa', 'SearchField', 'serverurl');
layerName = 'bluemarbleng';
layer = nasa.refine(layerName,  'SearchField', 'layername', ...
  'MatchType', 'exact');
[A, R] = wmsread(layer(1));
filename = [layerName '.tif'];
geotiffwrite(filename, A, R)
worldmap world
geoshow(filename)

The next example reads an image from an existing GeoTIFF file, and then it writes out the first 1024 columns and last 1024 rows to a new GeoTIFF file.

[A, R] = geotiffread('boston.tif');
row = [size(A,1)-1024+1 size(A,1)];
col = [1 1024];
subImage = A(row(1):row(2), col(1):col(2), :);
xi = col + [-.5 .5];
yi = row + [-.5 .5];
[xlim, ylim] = intrinsicToWorld(R, xi, yi);
subR = R;
subR.RasterSize = size(subImage);
subR.XLimWorld = sort(xlim);
subR.YLimWorld = sort(ylim);
info = geotiffinfo('boston.tif');
filename = 'boston_subimage.tif';
geotiffwrite(filename, subImage, subR,  ...
    'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
mapshow(filename);

If you'd like to know more about the capabilities of this new function, I encourage you to take a look at the extensive new example online.

Stay tuned for more about R2011a...

Get the MATLAB code

Published with MATLAB® 7.12