Steve on Image Processing

In my previous post on dilation algorithms I discussed structuring element decomposition. We looked at the decomposition of structuring elements with rectangular domains. Today I want to follow up on the idea by looking at decompositions of other structuring element shapes.

The strel function offers a variety of flat structuring element shapes:

(The strel function also offers nonflat shapes, and shapes with arbitrary domains.)

Let's look at the 'diamond' shape. The doc says SE = strel('diamond', R) creates a flat, diamond-shaped structuring element, where R specifies the distance from the structuring element origin to the points of the diamond." For example:

se1 = strel('diamond', 5)

 
se1 =
 
Flat STREL object containing 61 neighbors.
Decomposition: 4 STREL objects containing a total of 17 neighbors

Neighborhood:
     0     0     0     0     0     1     0     0     0     0     0
     0     0     0     0     1     1     1     0     0     0     0
     0     0     0     1     1     1     1     1     0     0     0
     0     0     1     1     1     1     1     1     1     0     0
     0     1     1     1     1     1     1     1     1     1     0
     1     1     1     1     1     1     1     1     1     1     1
     0     1     1     1     1     1     1     1     1     1     0
     0     0     1     1     1     1     1     1     1     0     0
     0     0     0     1     1     1     1     1     0     0     0
     0     0     0     0     1     1     1     0     0     0     0
     0     0     0     0     0     1     0     0     0     0     0

 

The display says the structuring element has 61 neighbors in its domain, but the four decomposed structuring elements have a total of only 17 neighbors.

What shapes makes up the decomposition? We can use getsequence and getnhood to find out.

se1_seq = getsequence(se1);
for k = 1:numel(se1_seq)

    getnhood(se1_seq(k))

ans =

     0     1     0
     1     1     1
     0     1     0

ans =

     0     1     0
     1     0     1
     0     1     0

ans =

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

ans =

     0     1     0
     1     0     1
     0     1     0

end

This decomposition comes from Rein van den Boomgard and Richard van Balen, "Methods for Fast Morphological Image Transforms Using Bitmapped Binary Images," Computer Vision, Graphics, and Image Processing: Models and Image Processing, vol. 54, no. 3, May 1992, pp. 252-254. (That journal name sure is a mouthful!)

Now let's look at the 'disk' structuring element. By default, strel('disk',R) actually gives you an approximation to a disk shape, where the approximation is based on a technique called radial decomposition using periodic lines. The reference sources are:

se2 = strel('disk', 5)

 
se2 =
 
Flat STREL object containing 69 neighbors.
Decomposition: 6 STREL objects containing a total of 18 neighbors

Neighborhood:
     0     0     1     1     1     1     1     0     0
     0     1     1     1     1     1     1     1     0
     1     1     1     1     1     1     1     1     1
     1     1     1     1     1     1     1     1     1
     1     1     1     1     1     1     1     1     1
     1     1     1     1     1     1     1     1     1
     1     1     1     1     1     1     1     1     1
     0     1     1     1     1     1     1     1     0
     0     0     1     1     1     1     1     0     0

 

Similar to the 'diamond' shape above, the 'disk' shape is decomposed into a sequence of smaller structuring elements:

se2_seq = getsequence(se2);
for k = 1:numel(se2_seq)

    getnhood(se2_seq(k))

ans =

     1
     1
     1

ans =

     1     0     0
     0     1     0
     0     0     1

ans =

     1     1     1

ans =

     0     0     1
     0     1     0
     1     0     0

ans =

     1     1     1

ans =

     1
     1
     1

end

I mentioned above that strel('disk',R) gives you an approximation to a disk. Let's see what that approximation looks like with a larger radius:

se3 = strel('disk', 21);
imshow(getnhood(se3), 'InitialMagnification', 'fit')

Well, that's definitely an approximation. It looks a lot more like an octagon than a circle. The syntax strel('disk',R,N) gives you more control over the approximation. N can be 0, 4, 6, or 8. The values 4, 6, and 8 indicate different levels of approximation. The default value is 4. With N equal to 8, the result looks closer to a circle, but dilation is a little slower because the decomposition isn't as efficient.

se4 = strel('disk', 21, 8);
imshow(getnhood(se4), 'InitialMagnification', 'fit')

Setting N = 0 asks strel to construct a disk shape with no decomposition approximation.

se5 = strel('disk', 21, 0);
imshow(getnhood(se5), 'InitialMagnification', 'fit')

Using getsequence we can see that there is no decomposition for se5:

numel(getsequence(se5))

ans =

     1

Next time we'll do some timing experiments with the decomposed structuring elements.

Get the MATLAB code

Published with MATLAB® 7.6

This is the second post in my series on algorithm concepts behind the implementation of dilation and erosion in the Image Processing Toolbox. Today I want to talk about structuring element decomposition.

Some useful structuring elements can be represented as the dilation of two or more smaller structuring elements. Dilation implementations can exploit this decomposition of structuring elements into smaller structuring elements by taking advantage of the associativity property of dilation.

That is, suppose image g is the dilation of image f with structuring element b:

where the star operator means dilation. And further suppose that b is the dilation of two other structuring elements:

Then the original image dilation equation can be rewritten as:

In English, the equation above says we can dilate f with structuring element b by:

1. Dilating f first with structuring element b1.

2. Dilating the result of step 1 with structuring element b2.

Let's see how this might help make dilation faster. One of the most commonly-used structuring element shapes is the square. For example, here's a simple 5-by-5 array of 1s that can define a structuring element domain (or neighborhood):

b = ones(5, 5)

b =

     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1

It turns out that a square (or rectangular) structuring element can be decomposed into two structuring elements: a horizontal line, and a vertical line.

b1 = ones(1, 5)

b1 =

     1     1     1     1     1

b2 = ones(5, 1)

b2 =

     1
     1
     1
     1
     1

You can check that the dilation of b1 and b2 equals b by using the 'full' option of imdilate:

imdilate(b1, b2, 'full')

ans =

     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1

How does this help us? Well, as I discussed last time, computing dilation directly from its defining equation requires Q comparisons per pixel, where Q is the number of elements in the structuring element domain. The 5-by-5 square structuring element has 25 elements in its domain. But structuring elements b1 and b2 have only five elements each. So you could dilate by one and then again by the other and require only 40% of the comparisons.

Let's look at how structuring element decompositions are computed and represented in the Image Processing Toolbox. The function strel makes a structuring element object. You give strel the name of a shape, such as 'rectangle', plus any parameters required for the shape. When you use strel to make a rectangular structuring element, the resulting object displays information about itself, including the decomposition:

se = strel('rectangle', [5 5])

 
se =
 
Flat STREL object containing 25 neighbors.
Decomposition: 2 STREL objects containing a total of 10 neighbors

Neighborhood:
     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1

 

You can call the function getsequence (actually a method of the strel class) to get the structuring elements in the decomposition:

decomp = getsequence(se)

 
decomp =
 
2x1 array of STREL objects
 

You use ordinary MATLAB array indexing to get each of the structuring elements from the decomposition:

decomp(1)

 
ans =
 
Flat STREL object containing 5 neighbors.

Neighborhood:
     1
     1
     1
     1
     1

 

decomp(2)

 
ans =
 
Flat STREL object containing 5 neighbors.

Neighborhood:
     1     1     1     1     1

 

Ordinarily, of course, strel and imdilate together handle all of these details for you when you perform image dilation.

Next time I plan to look at some of the other shapes offered by strel and see how they get decomposed.

Get the MATLAB code

Published with MATLAB® 7.6

Post updated on September 4 based on feedback from Cris Luengo's comment.

Today I'm starting a new series covering some of the concepts underlying algorithms for performing morphological dilation (and erosion, which is very similar). Most of the time, when people talk about image dilation, they mean the form of dilation that is a local maximum operation on the neighbors of each pixel. For example, here's how to compute the local maximum, for each image pixel, with that pixel and its eight neighbors:

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

se = ones(3,3)

se =

     1     1     1
     1     1     1
     1     1     1

B = imdilate(A, se)

B =

    24    24    24    16    16
    24    24    24    22    22
    23    23    21    22    22
    18    25    25    25    22
    18    25    25    25    21

The set of neighborhood pixels involved in the max operator forms a shape called the structuring element. The structuring element doesn't have to be square, as above. Here is dilation with a short, 3-pixel line:

se = [1 1 1];
B2 = imdilate(A, se)

B2 =

    24    24    24    15    15
    23    23    14    16    16
     6    13    20    22    22
    12    19    21    21    21
    18    25    25    25     9

Let's consider how long it should take to compute dilation. Suppose an image has P pixels, and the structuring element contains Q neighbors. Computing the output for a particular pixel involves both Q memory reads and Q comparisons. Therefore we might expect the computation time to be proportional to Q.

We can try to verify that behavior using square structuring elements of varying sizes. (I'll be using my timeit function, which you can download from the MATLAB Central File Exchange.)

% Read in a small sample image and replicate it to make a 1k-by-1k
% test image.
I = repmat(imread('rice.png'), 4, 4);

% Measure dilation times using a set of square structuring elements
% ranging in size from 5-by-5 to 97-by-97.
w = 5:4:97;

% The number of neighbors, Q, in each structuring element will be
% w^2.  We are guessing that the running time will be proportional
% to Q.
Q = w.^2;

% Initialize vector containing recorded times.
times = [];

for wk = w
    % Not that many MATLAB users are familiar with the for-loop
    % construct above.  I'll send a MATLAB t-shirt to the first
    % person who adds a comment to this post with a brief, clear,
    % and correct explanation.

    se = strel(ones(wk, wk));
    f = @() imdilate(I, se);
    times(end+1) = timeit(f);
end

plot(Q, times)
ylim([0 1.5*max(times)])

That plot looks almost constant, with maybe just a hint of an increase with Q.

So is imdilate somehow computing dilation faster than expected?

That's the subject of this series. We'll start digging in deeper next time.

Get the MATLAB code

Published with MATLAB® 7.6

Transparent graphics objects can be used effectively to visualize image processing concepts. Two particularly useful techniques are:

Today I'll show how to highlight image regions with patches. For this example I'll use the 'Extrema' measurement returned by regionprops. The extrema for a given object are eight points: the left-most pixel on the bottom, the right-most pixel on the bottom, the top-most pixel on the right, the bottom-most pixel on the right, and so on.

I'll start with the rice image, segmenting it using techniques I've shown before.

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

Even out the illumination with the tophat operator, threshold, and then clean up the thresholded image a bit.

I2 = imtophat(I, ones(15, 15));
bw = im2bw(I2, graythresh(I2));
bw2 = bwareaopen(bw, 5);
bw3 = imclearborder(bw2);
imshow(bw3)

Label the binary objects and compute the extrema.

L = bwlabel(bw3);
s = regionprops(L, 'Extrema');

Each object has 8 extrema points associated with it.

s(1).Extrema

ans =

   11.5000   86.5000
   12.5000   86.5000
   34.5000  100.5000
   34.5000  102.5000
   33.5000  103.5000
   27.5000  103.5000
    9.5000   89.5000
    9.5000   87.5000

We can superimpose the extrema-bounded shapes on top of the original rice image by using patch objects.

imshow(I)

hold on
for k = 1:numel(s)
    x = s(k).Extrema(:,1);
    y = s(k).Extrema(:,2);
    patch(x, y, 'g')
end
hold off

The above visualization is pretty clear. If you zoom in on some of the odd, larger shapes, though, you can't really tell what's going on.

axis([120 200 1 75])

We solve that by displaying the patches transparently.

imshow(I)

hold on
for k = 1:numel(s)
    x = s(k).Extrema(:,1);
    y = s(k).Extrema(:,2);
    patch(x, y, 'g', 'FaceAlpha', 0.3)
end
hold off

Now if you zoom in on the same region, we can see exactly what caused the unusual region.

axis([120 200 1 75])

Two of the rice grains were touching.

Next time I'll show a couple of techniques for visualizing one image transparently superimposed on another.

Get the MATLAB code

Published with MATLAB® 7.6

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