I've written often here about various computational and visualization techniques involving labeling connected components in binary images. Sometimes I use the function label2rgb to convert a label matrix into a color image with a different color assigned to each label.

Here's an example.

bw = imread('http://blogs.mathworks.com/images/steve/2012/rice-bw.png');
imshow(bw)

Now compute the connected components and the corresponding label matrix.

cc = bwconncomp(bw);
L = labelmatrix(cc);

In the label matrix, each foreground object in the original binary image is assigned a unique positive integer. Here, for instance, is how to display the tenth object.

imshow(L == 10)

Use the function label2rgb to "colorize" the label matrix.

rgb = label2rgb(L);
imshow(rgb)

That's a nice effect, but because of the way the colors are assigned, object near each other tend to have very similar colors. It might be better sometimes to assign the colors differently. That's what the 'shuffle' argument to label2rgb is for.

Here's the full syntax including 'shuffle':

rgb = label2rgb(L,map,zerocolor,'shuffle');

zerocolor is a three-element vector specifying what color is used for the background pixels.

Let's try it with the jet colormap and a light gray background color.

rgb = label2rgb(L,'jet',[.7 .7 .7],'shuffle');
imshow(rgb)

Now it's easier to see where two objects might be actually touching and so receive the same label. Let's zoom in closer to see:

xlim([128 185]);
ylim([5 62]);

I hope you find this useful!

Get the MATLAB code

Published with MATLAB® 7.14

On my list of potential blog topics today I saw just this cryptic item labeled dftmtx. Hmm, the MATLAB dftmtx function. But have I written about this function before? I better double-check by searching the old blog postings:

Oh, I forgot about Loren's post! She showed how the discrete Fourier transform, which MATLAB users normally compute by calling fft, can also be computed via a matrix multiply.

x = rand(1000,1);
X = fft(x);

T = dftmtx(1000);
X2 = T*x;

max(abs(X2(:) - X(:)))

ans =

   3.9790e-13

The difference is just floating-point round-off error.

My old post talked about the value of having an independent computation method available when you are testing your algorithm.

So today let's do something a little different. Let's compare the performance of computing the discrete Fourier transform using the DFT matrix versus using the fast Fourier transform.

To help with the timing, I'm going to use a function I wrote called timeit that you can download from the MATLAB Central File Exchange.

clear
n = 100:50:3000;
for k = 1:length(n)
    nk = n(k);
    x = rand(nk,1);
    T = dftmtx(nk);

    f = @() T*x;
    g = @() fft(x);

    times_f(k) = timeit(f);
    times_g(k) = timeit(g);
end

plot(n,times_f,n,times_g)
legend({'Using T*x', 'Using fft(x)'})

That's a pretty dramatic difference. The blue curve, showing the computation time using T*x, is an n^2 curve. Compared to that, the green curve, showing the computation time using fft(x), is so low that you can hardly see it.

Let's expand the y axis.

ylim([0 0.0005])

The lower green curve is an n*log(n) curve. The dramatic difference between that and the n^2 curve is why everyone got so excited when the fast Fourier transform algorithm was invented (or re-invented) a few decades ago.

Get the MATLAB code

Published with MATLAB® 7.14

http://www.youtube.com/watch?v=ywWBy6J5gz8

I said there was some room for improvement because I thought the last 20 seconds or so of the algorithm could be optimized away. However, math developer Bobby informed me that those last few dance steps are necessary to make sure that the right thing happens when the input contains NaNs. (Thanks also to Bobby for providing the title of today's post.)

Last year I posted a few times about the "Jähne test pattern," something that readers then told me was called a zone plate. This is a radial sinusoid pattern with low frequency in the middle of the image and high frequency toward the edges:

In my 02-Jul-2012 post, Filtering fun, I used the zone plate image to illustrate the effects of different kinds of frequency-selective filtering.

Lowpass:

Bandpass:

Directional:

Last week I submitted a function called imzoneplate to the MATLAB Central File Exchange.

I resisted the urge to give imzoneplate a lot of optional parameters. Instead there's only one, the size of the image.

imshow(imzoneplate(301))

Here's a smaller zone plate image and its horizontal cross-section.

I = imzoneplate(151);
imshow(I)

plot(I(76,:))

Let's zoom in at the edge of the plot.

xlim([140 151])

You can see that, with the parameters I've chosen for imzoneplate, the sinusoidal period decreases to about 3 samples at the edge of the image.

Be sure to download imzoneplate and give it a try. Post a comment here if you have a particular application for it.

Get the MATLAB code

Published with MATLAB® 7.14

Start by clicking on Support.

Then click on Bug Reports.

To find out which bugs have been fixed in R2012a in the Image Processing Toolbox, for example, perform the following search:

Here's the result:

This icon:

means that there is a way to work around the problem if you are using an earlier release.

I want to mention the iradon issue in particular and acknowledge Bob Alvarez, who called our attention to it with a blog comment last July. (However, you should generally report product bugs directly to MathWorks technical support rather than posting them as comments on this or other MATLAB Central blogs. My usual practice is to delete comments that aren't directly relevant to the topic of the blog post.)

We are, of course, continually trying to improve the quality of our products and therefore reduce the need for published bug reports.

Some of my fellow MATLAB Central bloggers have already talked about the MATLAB R2012a release. Usually, I like to scan the release notes and pick out just a few things that particularly interest me. I encourage you to take a look at the release notes to see what else is there that might interest you.

I was particularly happy to see that the web browser system built into MATLAB (used for various purposes, including the help system and the publish feature) was completely overhauled so that it looks a lot better on 64-bit Windows systems.

Speaking of the publish feature, have you heard of it? It's a great way to turn a MATLAB script (or function) into a web page, or a PDF, or a Word document. I use it to create almost all of my blog posts. You can find a lot of examples online by searching for the phrase "Published with MATLAB", such as Calculus with MATLAB. Anyway, you can include syntax-highlighted code in your published output. For example:

bw = imread('text.png');
s = regionprops(bw, {'Area', 'Centroid'});
imshow(bw)
hold on
for k = 1:numel(s)
  plot(s(k).Centroid(1), s(k).Centroid(2), 'b*');
end
hold off

There are several enhancements in the math area that interest me. For example, there are new numerical integration functions (integral, integral2, integral3) that have more capabilities than the functions they replace (quad and its cousins). Try an integration over an unbounded interval using integral:

fun = @(x) exp(-x.^2).*log(x).^2;
integral(fun,0,inf)

ans =

    1.9475

With quad that didn't work so well:

quad(fun,0,inf)

Warning: Infinite or Not-a-Number function value encountered. 

ans =

   NaN

The griddata function now supports natural neighbor interpolation.

And the unique function now has an additional option to preserve the original order of the data. I really like this one. The output of unique is normally sorted. Although you can postprocess the output to restore the original data order, it's tricky and I always puzzled over it a bit. With the new optional argument, it's really easy.

data = [3 2 9 4 8 1 1 10 1 2];
unique(data)

ans =

     1     2     3     4     8     9    10

You can see that the duplicates were removed and the result is sorted.

unique(data,'stable')

ans =

     3     2     9     4     8     1    10

Now the original data order has been preserved.

There are several nice enhancements in the file import/export area, including:

Finally, if you are really into object-oriented programming using MATLAB, check out the new ability to control the set of allowed subclasses, as well as the ability to control access to class members.

Get the MATLAB code

Published with MATLAB® 7.14

Image Processing Toolbox Enhancements

Version 8.0 includes these enhancements.

Intensity-Based Image Registration

The new imregister function lets you automatically align two images using intensity values, even when the images were created by two different devices (multimodal). With intensity-based registration, you do not need to specify control points.

You use the new imregconfig function to create the optimizer and the metric that imregister uses to specify the desired registration parameters.

Two New Functions to Visually Compare Images

The toolbox includes two new functions for visually comparing images: imshowpair and imfuse.

imshowpair creates a composite of two images and displays them in a figure.

imfuse creates a composite of two images and returns a third image that is a numeric matrix containing a fused version of the original images.

Circle Detection Using the Circular Hough Transform

The new imfindcircles function uses the Hough transform to find circular elements in grayscale, RGB, or binary images. To view the circles that have been detected, overlaid on the original image, use the viscircles function.

Performance Improvements

The performance of the imlincomb function has improved by taking advantage of multicore capabilities.

New and Updated Demos

The toolbox includes these new and updated demos.

• Registering Multimodal MRI Images (ipexknee)
• Finding the Rotation and Scale of a Distorted Image (ipexrotate)
• Measuring the Radius of a Roll of Tape (ipexradius) – Updated to use the new imfindcircles function

Computer Vision System Toolbox

First, I want to note that, as R2012a, the Computer Vision System Toolbox no longer requires the DSP System Toolbox or the Signal Processing Toolbox.

Version 5.0 includes these enhancements.

Audio Output Sampling Mode Added to the From Multimedia File Block

The From Multimedia File block now includes a new parameter, which allows you to select frame- or sample-based audio output. If you do not have a DSP System Toolbox license and you set this parameter for frame-based processing, your model will return an error. The Computer Vision System Toolbox software uses only sample-based processing.

New Viola-Jones Cascade Object Detector

The vision.CascadeObjectDetector System object uses the Viola-Jones algorithm to detect objects in an image. This detector includes Haar-like features and a cascade of classifiers. The cascade object detector is pretrained to detect faces, noses and other objects.

New MSER Feature Detector

The detectMSERFeatures function detects maximally stable extremal regions (MSER) features in a grayscale image. You can use the MSERRegions object, returned by the function, to manipulate and plot MSER features.

New CAMShift Histogram-Based Tracker

The vision.HistogramBasedTracker System object uses the continuously adaptive mean shift (CAMShift) algorithm for tracking objects. It uses the histogram of pixel values to identify the object.

New Integral Image Computation and Box Filtering

The integralKernel object with the integralImage and integralFilter functions use integral images for filtering an image with box filters. The speed of the filtering operation is independent of the filter size, making it ideally suited for fast analysis of images at different scales.

New Demo to Detect and Track a Face

This release provides a new demo, Face Detection and Tracking. This example shows you how to develop a simple face tracking system by detecting a face, identifying its facial features, and tracking it.

Improved MATLAB Compiler Support

MATLAB Compiler™ now supports detectSURFFeatures and disparity functions.

Code Generation Support

The vision.HistogramBasedTracker and vision.CornerDetector System objects now support code generation. See About MATLAB Coder for more information about code generation.

Image Acquisition Toolbox Enhancements

Version 4.2 includes these enhancements.

GenTL Support for GigE Vision, USB, and FireWire Devices

A new adaptor type, GenTL, is now available for GigE Vision, USB, and FireWire cameras.

Windows 64 Support on DCAM

The DCAM adaptor support is now extended to 64-bit Windows®.

VideoDevice System Object

The Image Acquisition Toolbox™ introduces the VideoDevice System Object, which allows single-frame image acquisition and code generation from MATLAB.

You use the imaq.VideoDevice function to create the System Object. It supports the same adaptors and hardware that the videoinput object supports; however, it has different functions and properties associated with it. For example, the System Object uses the step function to acquire single frames.

For more information on using the System Object, use this command in MATLAB:

help imaq.VideoDevice

From Video Device Block Enhancements

Several enhancements have been added to the From Video Device block for doing image acquisition in Simulink. Options for setting color space and Bayer Sensor Alignment are now included in the block properties.

New Hardware Support – National Instruments

Support has been added for additional National Instruments® hardware.

• PCI-1433
• PXI-1435

Test Suite for Third-Party Adaptor Developers and Camera Vendors

As part of the Image Acquisition Toolbox Adaptor Kit, we now offer a test suite for third-party adaptor developers and camera vendors to test adaptors and hardware against the toolbox.

DCAM Adaptor Improvement

Support has been added for strobe outputs to the DCAM adaptor on Windows.

Linuxvideo Adaptor Improvement

Improved support has been added for V4L1 cameras with the linuxvideo adaptor when using the V4L1 compatibility library provided by the libv4l project.

This week I want to highlight a new article on mathworks.com called "Applying Modern PDE Techniques to Digital Image Restoration," by Carola-Bibiane Schönlieb of the University of Cambridge. Dr. Schönlieb describes PDE-based inpainting, or image interpolation, used to reconstruct missing parts of images.

Her inpainting work began on a project to restore some 14th-century Viennese frescoes that were rediscovered in 1979:

She also describes other applications, such as filling in obscured portions of roads in satellite imagery.

The article describes inpainting using Fourier's heat equation as well as more sophisticated methods that can better handle discontinous image features such as edges. Here's a result using the heat equation. (Click here for the full size image.)

Check out the article if you're interested. You might also want to look at roifill, which implements a simple form of PDE-based inpainting.

Get the MATLAB code

Published with MATLAB® 7.14

 fid = fopen(filename, 'rb');
 assert(fid > 0, 'Could not open file.');

And . . .

 assert(isfield(meta, 'sizes') && ...
        isfield(meta, 'dimension') && ...
        isfield(meta, 'encoding') && ...
        isfield(meta, 'endian'), ...
        'Missing required metadata fields.')

onCleanup - Why worry about trying to remember to clean up resources? Let the onCleanup class take care of it for you. Construct one of these objects by giving it an anonymous function that closes your file handle when the object goes out of scope—whether from an error or at the end of the function.

 cleaner = onCleanup(@() fclose(fid));

regexp - Use MATLAB's regular expression engine to handle complicated text parsing for you.

 theLine = fgetl(fid);

 % "fieldname:= value" or "fieldname: value" or "fieldname:value"
 parsedLine = regexp(theLine, ':=?\s*', 'split', 'once');

Dynamic structure field indexing - If you have a string that's a legal MATLAB identifier, there's no need to write complicated logic just to use it as a field name in a structure. Simply use the .(string) construct.

 field = lower(parsedLine{1});
 value = parsedLine{2};

 field(isspace(field)) = '';  % Remove embedded spaces.
 meta(1).(field) = value;

Using temporary files to decompress data - The NRRD format supports storing the image data as raw bytes, human readable text, or GZIP-compressed byte streams. When a file contains compressed or encapsulated data and MATLAB has a file reader capable of handling that, it's easiest just to write the data to a temporary file and use the supported reader. Consider the readData() subfunction that recursively handles three different kinds of encoding:

 function data = readData(fidIn, meta, datatype)

 switch (meta.encoding)
  case {'raw'}

   data = fread(fidIn, inf, [datatype '=>' datatype]);

  case {'gzip', 'gz'}

   tmpBase = tempname();
   tmpFile = [tmpBase '.gz'];
   fidTmp = fopen(tmpFile, 'wb');
   assert(fidTmp > 3, ...
      'Could not open temporary file for GZIP decompression')

   tmp = fread(fidIn, inf, 'uint8=>uint8');
   fwrite(fidTmp, tmp, 'uint8');
   fclose(fidTmp);

   gunzip(tmpFile)

   fidTmp = fopen(tmpBase, 'rb');
   cleaner = onCleanup(@() fclose(fidTmp));

   meta.encoding = 'raw';
   data = readData(fidTmp, meta, datatype);

  case {'txt', 'text', 'ascii'}

   data = fscanf(fidIn, '%f');
   data = cast(data, datatype);

  otherwise
   assert(false, 'Unsupported encoding')
 end

swapbytes - Like many formats, NRRD supports little-endian and big-endian byte ordering. The swapbytes function makes it dead simple to change endianness, and the computer function helps you determine whether swapping is necessary. Here's the pattern, which uses the "endian" metadata value read from the .nrrd file:

 function data = adjustEndian(data, meta)

 [~,~,endian] = computer();

 needToSwap = (isequal(endian, 'B') && ...
                isequal(lower(meta.endian), 'little')) || ...
              (isequal(endian, 'L') && ...
               isequal(lower(meta.endian), 'big'));

 if (needToSwap)
     data = swapbytes(data);
 end

Happy coding!

- Jeff Mather

Get the MATLAB code

Published with MATLAB® 7.13

The MATLAB function imwrite writes image data to a variety of different formats including: reads image data from a variety of formats, including JPEG, JPEG 2000, TIFF, PNG, and GIF (as well as a few older formats no longer commonly used).

Probably most people use form where you just provide the image array and a filename.

I'll read in one of the sample images in Image Processing Toolbox to experiment with.

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

imwrite(rgb,'mypeppers.tif')

In the above call to imwrite, the function decides which image format to use based on the extension (.tif) in the filename.

Sometimes, though, you might want to use a filename extension for a different purpose. In that case, you call pass a third argument to imwrite to let it know which image format to use.

imwrite(rgb, 'test_image.study_013', 'tif')

Most image file formats support some form of compression. That means that the image takes up less space on disk than you might think based on a simple calculation. For example, the rgb array above takes up about 576 kilobytes in memory.

bytes_in_memory_peppers = numel(rgb)

bytes_in_memory_peppers =

      589824

But it takes less space on disk when you save it to an image format that uses compression. Let's try writing to a PNG file:

imwrite(rgb,'mypeppers.png')

info = imfinfo('mypeppers.png');
bytes_on_disk_peppers_png = info.FileSize

bytes_on_disk_peppers_png =

      287589

compression_ratio = bytes_in_memory_peppers / bytes_on_disk_peppers_png

compression_ratio =

    2.0509

The PNG format uses lossless compression. This means that you can recover the original image data exactly.

peppers2 = imread('mypeppers.png');
isequal(rgb, peppers2)

ans =

     1

Lossless compression methods save a lot more space when the image has large constant areas or contains certain kinds of repetitive, predictable patterns. For example, this image is a screen shot of a MATLAB plot:

url = 'http://blogs.mathworks.com/images/steve/2012/plot_screen_shot.png';
rgb_plot = imread(url);
bytes_in_memory_plot = numel(rgb_plot)

bytes_in_memory_plot =

      336000

info = imfinfo(url);
bytes_on_disk_plot = info.FileSize

bytes_on_disk_plot =

        3284

compression_ratio = bytes_in_memory_plot / bytes_on_disk_plot

compression_ratio =

  102.3143

That's a pretty big compression ratio!

Other image file formats, such as JPEG, use lossy compression. (There is a lossless form of JPEG, but it is obscure and rarely used.) With lossy compression, properties of the human visual system are exploited to reduce the disk size even further, but at a cost --- the original pixels are not exactly recoverable. Let's write out the peppers image using JPEG compression.

imwrite(rgb,'peppers.jpg');
info = imfinfo('peppers.jpg');
bytes_on_disk_peppers_jpg = info.FileSize

bytes_on_disk_peppers_jpg =

       23509

compression_ratio_peppers_jpg = bytes_in_memory_peppers / ...
    bytes_on_disk_peppers_jpg

compression_ratio_peppers_jpg =

   25.0893

That's a much bigger compression ratio than we got with the PNG format above. The cost is in reduced image quality. JPEG does a pretty good job of reducing space without adversely affecting image quality:

imshow('peppers.jpg')

But it isn't perfect, as you can see by zooming in.

subplot(1,2,1)
imshow('peppers.png')
limits = [232 276 215 248];
axis(limits)
title('Original')

subplot(1,2,2)
imshow('peppers.jpg')
axis(limits)
title('JPEG compressed')

You can see the compression artifacts around the stem.

Compression artifacts tend to be more visible on images such as the plot screen shot above.

imwrite(rgb_plot,'plot_screen_shot.jpg')
rgb_plot_compressed = imread('plot_screen_shot.jpg');

subplot(2,2,1)
imshow(rgb_plot)
title('Original image')
limits = [80 105 210 230];
axis(limits);

subplot(2,2,2)
imshow(rgb_plot_compressed)
title('JPEG compressed image')
axis(limits);

subplot(2,2,3)
imshow(rgb_plot)
limits = [345 380 240 265];
axis(limits);

subplot(2,2,4)
imshow(rgb_plot_compressed)
axis(limits);

That's why I usually use PNG for screen shots of many MATLAB graphics. For more complex images, on the other hand, I usually use JPEG.

For more information, see Section 2.4 of Digital Image Processing Using MATLAB.

See also the reference pages for imwrite, as well as the section Reading and Writing Image Data in the Image Processing Toolbox User's Guide.

Get the MATLAB code

Published with MATLAB® 7.13