A MATLAB user recently asked in the MATLAB newsgroup how to fill "small" holes in a binary image. The function imfill can be used to fill all holes, but this user only wanted to fill holes having an area smaller than some threshold.

That's an interesting question. It can be done using a combination of imfill, bwareaopen, and MATLAB logical operators. Here's how.

Step 1: Fill all holes using imfill:

original = imread('circbw.tif');
imshow(original)

filled = imfill(original, 'holes');
imshow(filled)
title('All holes filled')

Step 2: Identify the hole pixels using logical operators:

holes = filled & ~original;
imshow(holes)
title('Hole pixels identified')

Step 3: Use bwareaopen on the holes image to eliminate small holes:

bigholes = bwareaopen(holes, 200);
imshow(bigholes)
title('Only the big holes')

Step 4: Use logical operators to identify small holes:

smallholes = holes & ~bigholes;
imshow(smallholes)
title('Only the small holes')

Step 5: Use a logical operator to fill in the small holes in the original image:

new = original | smallholes;
imshow(new)
title('Small holes filled')

All done!

Get the MATLAB code

Published with MATLAB® 7.6

I've previously blogged about image restoration methods based on inverse filtering and Wiener filtering. Wiener filtering gives really impressive results compared to inverse filtering. However, Wiener filtering assumes that the image is modeled as a random process for which 2nd-order statistics are known, along with the noise variance. This viewpoint may not be attractive to some users, and the required information may not always be available.

Regularized restoration provides similar results to Wiener filtering but is justified by a very different viewpoint. Also, less prior information is required to apply regularized restoration. (Note: The Image Processing Toolbox function that implements regularized restoration is called deconvreg.)

As before, I will assume a shift-invariant blur model with independent additive noise, which is given by

y(m,n) = h(m,n)*x(m,n) + u(m,n)

where * is 2-D convolution, h(m,n) is the point-spread function (PSF), x(m,n) is the original image, and u(m,n) is noise.

Regularization trades off two desirable goals -- 1) the closeness of the model fit and 2) the closeness of the model behavior to something that would be expected in the absence of specific knowledge of the model parameters or data. Remember that inverse filtering minimizes:

which fits the model to the data as closely as possible.

As a reminder, let's look at the horrible results that this seemingly logical procedure can yield:

% form PSF as a disk of radius 3 pixels
h = fspecial('disk',4);
% read image and convert to double for FFT
cam = im2double(imread('cameraman.tif'));
hf = psf2otf(h,size(cam));
cam_blur = real(ifft2(hf.*fft2(cam)));
imshow(cam_blur)
title('blurred image')
xlabel('(original image courtesy Massachusetts Institute of Technology)')

sigma_u = 10^(-40/20)*abs(1-0);
cam_blur_noise = cam_blur + sigma_u*randn(size(cam_blur));
imshow(cam_blur_noise)
title('noisy image')

cam_inv = real(ifft2(fft2(cam_blur_noise)./hf));
imshow(cam_inv)
title('inverse filter restoration')

The problem is that the solution fits the noise as well as the underlying data. The inverse filter divides by some very small values where noise in the numerator is relatively large compared to the attenuated signal. As a result, the solution becomes far too rough. Regularization adds another term to the minimization criterion to force the image to be somewhat smooth:

The second term

is a filtered version of the estimated image that we expect to be small. In image processing, we usually expect a high level of local correlation or smoothness in the image. So the filter is chosen as a highpass filter so that we minimize the high-frequency content or "roughness" in the solution. The parameter

controls the degree of the roughness penalty. The larger it is, the smoother the restored image will be.

We can rewrite the minimization criterion in the DFT domain as

Taking a derivative with respect to the image DFT coefficients and setting the result to zero yields the regularized restoration in the DFT domain:

The key to the performance of this filter is the extra term in the denominator. If this extra term is zero, the filter reverts to an inverse filter. However, for a non-zero term, the denominator is prevented from growing too small, especially at higher frequencies where the attenuated signal is likely to be quite small. Thus, noise amplification resulting from division by very small numbers is significantly reduced.

The regularization filter is often chosen to be a discrete Laplacian:

The regularization parameter alpha can be chosen either by trial-and-error or by one of several different statistical methods. Some of these methods require no additional prior information.

lf = fft2([0 1 0; 1 -4 1; 0 1 0],256,256);
alpha = 0.01;
cam_inv = real(ifft2(fft2(cam_blur_noise)./hf));
cam_reg = deconvreg(cam_blur_noise,h,0,alpha);

imshow(cam_reg)
title('regularized filter restoration')

This filter can be understood as an approximation of a Wiener filter. If we choose

then the result approximates the Wiener filter expression.

Let's explore the role of

alpha = 0.0005;
cam_reg = deconvreg(cam_blur_noise,h,0,alpha);
imshow(cam_reg)
title('\alpha = 0.0005')

alpha = 0.01;
cam_reg = deconvreg(cam_blur_noise,h,0,alpha);
imshow(cam_reg)
title('\alpha = 0.01')

alpha = 0.2;
cam_reg = deconvreg(cam_blur_noise,h,0,alpha);
imshow(cam_reg)
title('\alpha = 0.2')

It is evident from the images that a smaller alpha results in a noisier but sharper image while larger alpha results in a cleaner but blurrier image.

While the choice of the regularization filter does have an influence on the results, even a non-optimal choice can yield good results. If we use an impulse (no filtering) instead of the high-pass discrete Laplacian, we get the following results:

alpha = 0.01;
regop = 1;
cam_reg = deconvreg(cam_blur_noise,h,0,alpha,regop);
imshow(cam_reg)
title('l(m,n) = no filtering')

alpha = 0.01;
lf = fft2([0 1 0; 1 -4 1; 0 1 0],256,256);
cam_reg = deconvreg(cam_blur_noise,h,0,alpha);
imshow(cam_reg)
title('l(m,n) = discrete Laplacian')

The results are pretty insensitive to the type of filter used to penalize roughness. Even no filtering (penalizing only the sum of squared pixel values) yields results that are nearly as good as the discrete Laplacian.

Regularization can be a useful tool when statistical information is unavailable. Furthermore, this framework can be extended to adapt to image edges, spatially varying noise, and other challenges. But those are topics for other blogs...

- by Stan Reeves, Department of Electrical and Computer Engineering, Auburn University

Get the MATLAB code

Published with MATLAB® 7.6

Today I want to show you a morphological operation called "opening by reconstruction."

The normal morphological opening is an erosion followed by a dilation. The erosion "shrinks" an image according to the shape of the structuring element, removing objects that are smaller than the shape. Then the dilation step "regrows" the remaining objects by the same shape.

Here's an example using a fragment of text from the book Digital Image Processing Using MATLAB.

url = 'http://blogs.mathworks.com/images/steve/2008/book_text.png';
text = imread(url);
bw = text(1:500, 1:500);
imshow(bw)

Suppose we want to identify characters containing a tall vertical segment. We can do this by opening with a vertical structuring element.

Erode first:

se = strel(ones(51, 1));
bw2 = imerode(bw, se);
imshow(bw2)

Then dilate:

bw3 = imdilate(bw2, se);
imshow(bw3)

Or you can do the opening in a single step by calling imopen:

bw3 = imopen(bw, se);
imshow(bw3)

The dilation step in the opening operation restored the vertical strokes, but the other strokes of the characters are missing. How can we get the entire characters containing vertical strokes?

The answer is to use morphological reconstruction. For binary images, reconstruction starts from a set of starting pixels (or "seed" pixels) and then grows in flood-fill fashion to include complete connected components.

To get ready to use reconstruction, first define a "marker" image. This is the image containing the starting or seed locations. For our text example, the marker image will the output of the erosion.

marker = imerode(bw, se);
imshow(marker)

Next, define mask image. The flood-filling will be constrained to spread only to foreground pixels in the mask image. We can use the original text image as our reconstruction mask.

mask = bw;

Finally, call imreconstruct to perform the operation.

characters = imreconstruct(marker, mask);
imshow(characters)

Performing morphological reconstruction, using the eroded image as the marker and the original image as the mask, is called "opening by reconstruction."

Do you have other uses for morphological reconstruction in your own applications? Tell us about it: Click on the "Comment" link below.

Get the MATLAB code

Published with MATLAB® 7.6

The local standard deviation operator is often used as a measure of "busy-ness" throughout the image. For each pixel, the standard deviation of that pixel's neighbors is computed:

sqrt( sum ( (x - xbar).^2 ) )

(Scale factors omitted.)

Mathematically, the summation can be rewritten as:

sum(x.^2) - sum(x).^2/N

Both of these local sums can be computed by using imfilter with an all-ones filter.

I = im2double(imread('cameraman.tif'));
imshow(I)

Original image credit: Massachusetts Institute of Technology

To compute the sum(x.^2) term, we square (elementwise) the input to imfilter.

h = ones(5,5);
term1 = imfilter(I.^2, h);
imshow(term1, [])

Notice the dark band around the edge. This is because imfilter zero-pads by default. We might want to use the 'symmetric' option instead.

term1 = imfilter(I.^2, h, 'symmetric');
imshow(term1, [])

To compute the sum(x).^2 term, we square the output of imfilter.

term2 = imfilter(I, h, 'symmetric').^2 / numel(h);
imshow(term2, [])

Then we subtract the second term from the first and take the square root.

local_std = sqrt(term1 - term2);  % scale factor omitted
imshow(local_std, [])

Cautionary notes

local_std = sqrt(max(term1 - term2, 0));

Recent releases of the Image Processing Toolbox include the function stdfilt, which does all this work for you.

The local geometric mean filter multiplies together all the pixel values in the neighborhood and then takes the N-th root, where N is the number of pixels in the neighborhood. The geometric mean filter is said to be slightly better than the arithmetic mean at preserving image detail.

Use the old logarithm trick to express the geometric mean in terms of a summation. Then imfilter can be used to compute the neighborhood summation, like this.

geo_mean = imfilter(log(I), h, 'replicate');
geo_mean = exp(geo_mean);
geo_mean = geo_mean .^ (1/numel(h));

imshow(geo_mean, [])

Does anyone have other creative uses of imfilter to share?

Copyright 2008 The MathWorks, Inc.
Get the MATLAB code

Published with MATLAB® 7.6

Today I want to demonstrate a useful technique to produce a false-color visualization of different sets of binary image objects. Here's the sample image that we'll use:

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

Let's look at two sets of objects: Those that touch the image border, and those that do not.

notouch = imclearborder(bw);
imshow(notouch)

The objects that do touch the border can be computed using logical operators:

touch = bw & ~notouch;
imshow(touch)

Here's one way to turn these two sets of objects into a single, false-color, indexed image. First, initialize the index matrix:

X = zeros(size(bw), 'uint8'); % Did you know about
                              % this way to call zeros?

Now assign 1 to the elements of X corresponding to the border-touching set, and assign 2 to the elements corresponding to the interior set.

X(touch) = 1;  % Logical indexing!
X(notouch) = 2;

Now we just need to pick some colors for the color map. I'll make the background white:

map(1,:) = [1 1 1];

Make the touching objects be purple-ish.

map(2,:) = [0.7 0.3 0.8];

And use a green shade for the removed objects.

map(3,:) = [0.4 0.8 0.7];

Now we can display the resulting indexed image.

imshow(X,map)

Get the MATLAB code

Published with MATLAB® 7.6

All of the toolbox sample images are stored in the same subdirectory of the MATLAB installation directory. We can use matlabroot and fullfile to find the complete path.

matlabroot

ans =

C:\Program Files\MATLAB\R2008a

demos_folder = fullfile(matlabroot, 'toolbox', 'images', 'imdemos')

demos_folder =

C:\Program Files\MATLAB\R2008a\toolbox\images\imdemos

The purpose of fullfile is to produce a valid path string on all MATLAB platforms. That way, you don't have to worry about whether your code will be running on a platform that uses "/" or "\" as the path separator.

Did you know you can automatically determine all of the file formats supported by imread and imwrite? The secret is to use a little-known function called imformats. This function returns a struct array containing information about format descriptions, recognized file extensions, and functions to be used for reading and writing each format.

f = imformats

f = 

1x16 struct array with fields:
    ext
    isa
    info
    read
    write
    alpha
    description

The number of elements of the structure array is the number of recognized image file formats.

num_formats = numel(f)

num_formats =

    16

Each element of the structure contains information about a particular format. This information is used by imread, imwrite, and imfinfo.

f(1)

ans = 

            ext: {'bmp'}
            isa: @isbmp
           info: @imbmpinfo
           read: @readbmp
          write: @writebmp
          alpha: 0
    description: 'Windows Bitmap (BMP)'

The ext field of the struct array is a cell array containing the recognized filename extensions for each format.

f(1).ext

ans = 

    'bmp'

The reason ext is a cell array is that multiple extensions are recognized for certain formats. For example:

f(3).ext

ans = 

    'fts'    'fits'

f(7).ext

ans = 

    'jpg'    'jpeg'

Use the comma-separated list syntax to create a single list of all recognized extensions.

extensions = [f.ext]

extensions = 

  Columns 1 through 7

    'bmp'    'cur'    'fts'    'fits'    'gif'    'hdf'    'ico'

  Columns 8 through 14

    'jpg'    'jpeg'    'pbm'    'pcx'    'pgm'    'png'    'pnm'

  Columns 15 through 19

    'ppm'    'ras'    'tif'    'tiff'    'xwd'

Now build up a list of image files in the Image Processing Toolbox demos folder. Use dir in a loop to do directory listings of *.jpg, *.tif, *.bmp, etc.

% Initialize an empty struct array.
files = struct([]);
for k = 1:numel(extensions)
    files_k = dir(fullfile(demos_folder, sprintf('*.%s', extensions{k})));
    files = [files; files_k];
end

files

files = 

64x1 struct array with fields:
    name
    date
    bytes
    isdir
    datenum

files(1)

ans = 

       name: 'solarspectra.fts'
       date: '13-Oct-2002 08:48:02'
      bytes: 66240
      isdir: 0
    datenum: 7.3150e+005

files(2)

ans = 

       name: 'football.jpg'
       date: '01-Mar-2001 09:52:38'
      bytes: 27130
      isdir: 0
    datenum: 7.3091e+005

We might want to sort the image files alphabetically. Do this by:

1. Make a cell array of file names.

2. Sort the cell array of file names and save the second output from sort. (The second output is called the permutation vector of the sort.)

3. Rearrange the structure array using the permutation vector.

Step 1:

filenames = {files.name};

Step 2:

[filenames_sorted, idx] = sort(filenames);

Step 3:

files = files(idx);

Now we construct HTML code for a page that displays all the sample image files in a bulleted list.

html = cell(0,0);
html{end+1} = '<html>';
html{end+1} = '<head>';
html{end+1} = '<title>Image Processing Toolbox Demo Images</title>';
html{end+1} = '</head>';
html{end+1} = '<body>';
html{end+1} = '<h2>Image Processing Toolbox Demo Images</h2>';
html{end+1} = '<ul>';

for k = 1:numel(files)
    html{end+1} = '<li>';
    html{end+1} = files(k).name;
    html{end+1} = '</li>';
end

html{end+1} = '</ul>';
html{end+1} = '</body>';
html{end+1} = '</html>';

We could save the generated HTML into a file. But we can also display the HTML directly in the MATLAB Web Browser without saving it. To do that, we need to:

1. Construct a single character array that contains all of the HTML.

2. Prepend 'text://' to the HTML code.

3. Pass the resulting string to the web function.

html_str = sprintf('%s\n', html{:});
web(['text://', html_str])

Here's a screen shot showing the results in the MATLAB Web Browser:

Get the MATLAB code

Published with MATLAB® 7.6

I thought I would finish my discussion of applylut and makelut by describing a performance optimization we implemented for version 6 (R2007b) of the Image Processing Toolbox.

A couple of years ago, a developer on our team wanted to know the details of the 'thin' option of bwmorph. This syntax "thins" connected strokes in a binary image. For example:

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

bwt = bwmorph(bw, 'thin');
imshow(bwt)

With the syntax bwmorph(bw, 'thin', N), you can tell bwmorph to repeat the thinning operation N times. If you use bwmorph(bw, 'thin', Inf), then bwmorph will repeat the operation until the image stops changing.

bwt = bwmorph(bw, 'thin', Inf);
imshow(bwt)

So, in this long-ago conversation about bwmorph, I explained that many of the operations supported by bwmorph are implemented using applylut. About halfway through explaining applylut, I had a light bulb moment. It dawned on me that the thinning operation of bwmorph never changes a 0-valued pixel to a 1! So for any 0-valued input pixel, it is kind of pointless to examine the surrounding 3-by-3 neighborhood, which is what applylut does. One could just immediately determine that the corresponding output pixel should be 0.

There are other operations ('thicken', for example) that never change a 1-valued pixel to a 0.

After some further discussion, we realized that applylut could examine any lookup table to determine whether it had one of these properties. So that's what we did. Starting in version 6 of the toolbox, applylut examines the input lookup table to see if it always transforms a 0-valued pixel to 0, or a 1-valued pixel to 1. If so, special-case code is used that skips the neighborhood lookup step for pixels that are 0 (or 1).

Here's the time required on my Windows laptop to apply the thinning operation to text.png. (You can find the function timeit on the MATLAB Central File Exchange.)

f = @() bwmorph(bw, 'thin', Inf);
timeit(f)

ans =

    0.0044

Running with an older version of the Image Processing Toolbox, the same operation takes about 104 milliseconds. That's about a factor of 23 reduction in execution time!

[Note added June 16, 2008:] The applylut optimization was not the only change made to speed up the 'thin' option of bwmorph. We also realized that work being done using several lookup tables could be combined in a single lookup table.

These optimizations went into the product in R2007b, which shipped in September 2007. Your mileage will certainly vary depending on your particular data. The speedup you see depends on how many 0-valued (or 1-valued) pixels the image contains.

For your reference, here are my previous posts on applylut and makelut:

Get the MATLAB code

Published with MATLAB® 7.6