Steve on Image Processing

September 3rd, 2008

Dilation algorithms—Introduction

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

7 Responses to “Dilation algorithms—Introduction”

  1. Omar replied on :

    for wk = w

    w is a 1-by-n vector that contains n indices. wk will iterate n-times w and use each value in w as a current index. E.g.,

    wk = [1, 3, 5, 7]

    wk will have values (1,3,5,7) at each of the four iterations.

  2. Steve replied on :

    Omar—Nice work.

  3. Jim replied on :

    for i = A
    …statements…
    end;

    Description:
    The …statements… are executed (as MATLAB commands) for i equal sequentially to each column in A. If A is 3D or greater, the …statements… are executed for i equal sequentially to the vector of first index values for each combination of the other indices of A.

    Example:
    A = [1 2 3; 4 5 6];
    for i = A
    this_is_i = i
    end;

    Result:
    this_is_i =

    1.00
    4.00

    this_is_i =

    2.00
    5.00

    this_is_i =

    3.00
    6.00

  4. Steve replied on :

    Jim—Thanks for adding your comment showing how the syntax works for matrices.

  5. Cris Luengo replied on :

    Not to spoil your upcoming bog entry too much, but if you scale the first graph (times vs Q) by setting the lower ylim to 0, you’ll see it is pretty close to a horizontal line. There’s only a 10% difference from low to high.

    However, isn’t this supposed to be an O(1) algorithm? Where’s the increase coming from? Actually, only yesterday I saw some paper with a timing plot, where operations with a linear SE actually had lower times for longer SE. That one seriously surprises to me! (Urbach and Wilkinson, IEEE TIP 17(1):1-8, fig 9)

  6. Steve replied on :

    Cris—You’re right, I should have caught the plot scaling issue. I wasn’t actually trying hide the point. I updated the post to make it more clear. Thanks. Also, thanks for the TrIP pointer; I’ll take a look.

  7. software developer replied on :

    Interesting,

    this is a good article to help with algorithms…

    keep up the good work

    Thanks for writing about it

Leave a Reply

Wrap code fragments inside <pre> tags, like this:

<pre class="code">
a = magic(3);
sum(a)
</pre>

If you have a "<" character in your code, either follow it with a space or replace it with "&lt;" (including the semicolon).


Steve Eddins manages the Image & Geospatial development team at The MathWorks and coauthored Digital Image Processing Using MATLAB. He writes here about image processing concepts, algorithm implementations, and MATLAB.

  • Sana: hi steve, could you explain to me how i would be able to use the dir function, to do a loop through a directory...
  • Nishtha: Sir, I have preprocessed the image in following steps: [1] adaptive histogram equalization [2] thresholding...
  • Kristof: I also strongly support the idea. I have just recently bumped into the problem that im2single was not...
  • Steve: David—I’ m glad you found it useful!
  • David Lalejini: I found your example very useful for finding connected nodes in a large set of input pairs. I start...
  • tommy: Dear Steve, I have a question,please if you are kind to help me regarding the accumulator array dimensions of...
  • Steve: Abc—I don’t know how to distinguish the faces. You might try posting your question in the MATLAB...
  • Manju: well if we have a few ovals within each other like in a cell how do we measure the distance from the center...
  • Steve: Manju—What do you mean? How is each region defined?
  • Manju: if we have 2-3 regions within each other how do we measure the regions of each one?

These postings are the author's and don't necessarily represent the opinions of The MathWorks.