Steve on Image Processing

November 26th, 2008

A structuring element decomposition “gotcha”

In my earlier posts about dilation and erosion algorithms, I wrote about exploiting structuring element decomposition for improving speed. Today I want to discuss a potential decomposition gotcha: If you don't implement decomposed dilation and erosion carefully, you can easily end up with wrong answers near the image boundaries.

Here's an example to illustrate what can go wrong, starting with a single binary image containing only one foreground pixel;

A = false(5,5);
A(5,3) = true
A =

     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

Here's a structuring element:

b = [1 0 0; 1 1 0; 1 1 1; 0 1 1; 0 0 1]
b =

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

Compute the dilation of A by b using imdilate:

A_b = imdilate(A, b)
A_b =

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

The structuring element b can be decomposed into two smaller structuring elements, b1 and b2.

b1 = [1; 1; 1]
b1 =

     1
     1
     1

b2 = [1 0 0; 0 1 0; 0 0 1]
b2 =

     1     0     0
     0     1     0
     0     0     1

Morphology theory (actually, just the distributive property you learned back in high school) suggests that we should be able to dilate A by b1, and then dilate the result by b2, to get the same result as dilating A by b. Let's try that.

A_b1_b2 = imdilate(imdilate(A, b1), b2)
A_b1_b2 =

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

Uh oh. Compare A_b and A_b1_b2 carefully to see they are not the same. So where did we go wrong? Roughly speaking, we failed to take into consideration the effect of dilation outside the image boundaries. If you take the entire image plane into account, dilation by b1 affects pixels outside the original image boundaries. The out-of-bounds pixels can then affect in-bounds results when we dilate by b2. But in the call imdilate(imdilate(A, b1), b2) above, each call to imdilate returns a result that's only the size of the input image; affected out-of-bounds pixels get cropped out.

So let's try that again, but this time we'll pre-pad the input image with 0s.

A_p = padarray(A, [1 1], 0, 'both')
A_p =

     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     0     0     1     0     0     0
     0     0     0     0     0     0     0

Now dilate with b1 and b2 in succession:

A_p_b1_b2 = imdilate(imdilate(A_p, b1), b2)
A_p_b1_b2 =

     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     1     1     0     0     0
     0     0     1     1     1     0     0
     0     0     0     1     1     0     0

And trim away the padded rows and columns:

A_p_b1_b2_cropped = A_p_b1_b2(2:end-1, 2:end-1)
A_p_b1_b2_cropped =

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

Success! That result equals A_b.

That means we have a trade-off to make. We can get significant speed improvements using structural element decomposition, but we can guarantee correct results at the image boundaries only if we use extra memory to make a padded copy of the input image. That's what imdilate and imerode do.

I have two more brief topics in mind for this series on dilation and erosion algorithms. I'll try to wrap things up in December.


Get the MATLAB code

Published with MATLAB® 7.7

8 Responses to “A structuring element decomposition “gotcha””

  1. Ken Eaton replied on :

    Steve,

    I’ve found your posts about dilation and structure elements very interesting, since I’ve been doing some things with it recently. I came across something strange that I wasn’t expecting, and I was hoping you might explain it for me…

    I had some operations were I wanted to dilate points with a disk shape, and found that using IMFILTER instead of IMDILATE was actually much faster:

      a = zeros(9);
      a(5,5) = 1;
    
      % using imdilate:
    
      disk = strel('disk',5);
      tic;
      for iCount = 1:100,
        b_dilate = imdilate(a,disk);
      end
      toc;
    
      % using imfilter:
    
      disk = fspecial('disk',4);
      tic;
      for iCount = 1:100,
        b_filter = double(imfilter(a,disk,'replicate') > 0);
      end
      toc;
    
      % using a thresholded imfilter:
    
      disk = double(fspecial('disk',4) > 0);
      tic;
      for iCount = 1:100,
        b_thresh = imfilter(a,disk,'replicate');
      end
      toc;
    

    Using MATLAB version 7.1, I get the following from the above code, and all the results are equal:

    Elapsed time is 0.341824 seconds.
    Elapsed time is 0.089447 seconds.
    Elapsed time is 0.082622 seconds.

    Why is this happening? Is this the same for newer MATLAB versions?

    Thanks,
    Ken

  2. Oliver Woodford replied on :

    Hi Steve

    On the topic of structural elements, why can they be used for IMDILATE, but not for BWLABEL? Currently I need to use them in the latter (see http://www.mathworks.com/matlabcentral/newsreader/view_thread/240944), but I can’t see how to achieve this using the IPT. Is there a way?

    Many thanks,
    Oliver

  3. Steve replied on :

    Oliver—bwlabeln supports labeling operations with an arbitrary connectivity definition. See the doc regarding the CONN input argument. You can use bwlabeln on two-dimensional inputs.

  4. Oliver Woodford replied on :

    Thanks, Steve. However, the CONN matrix doesn’t allow arbitrary strel neighbourhoods. In particular, the neighbourhoods are limited to 3×3 (for 2d) and must be symmetric. For the problem I have (with a 9×9 neighbourhood), BWLABELN can’t be used.

  5. Steve replied on :

    Oliver—I didn’t say bwlabel allows arbitrary strels, just arbitrary connectivities. I would be interested to hear more about your application, if you care to share it. For example, when would it make sense to say that the two 1-valued pixels below are connected?

    1 0 0 0 0 0 0 0 1
    

    Anyway, some of the material in my series on connected component labeling algorithms might be useful to you.

  6. Steve replied on :

    Ken—I get similar results. Clearly, we have a lot of room to improve our dilation implementation! We’re looking into it.

  7. Alex replied on :

    Hi.. I cannot understand what is written in matlab help. Could you give me a brief description on how to select the appropriate structuring element parameter? wheter is disc or rectangle or wad.

  8. Steve replied on :

    Alex—The appropriate structuring element shape depends on your application and data set, just like choosing between a lowpass filter and a highpass filter. You might want to get a good textbook on image processing.

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.