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.

  • murat: Hi Steve, I have an rgb image of a kind of cream and it contains some small black particles (black dots). In...
  • Steve: Ernest—Look at setting the FaceColor property. The code for setting that is shown on the page you asked...
  • Ernest Miller: Hi Steve, Understood. However, can you explain how to change the colors? Thanks, Ernest
  • Jan: Hi Steve Very useful code, yet what if I parts of my rotated+translated object are outside the original...
  • Steve: MoHDa—It might be possible. You’ll need to use one of the options that produces closed edge...
  • MoHDa: I have one question about the ROIPOLY: I have an image with stripes, I use the “edge” command for...
  • Steve: Shahn—My November 17, 2006 post shows you how to do it.
  • Steve: Kay-Uwe—Thanks for following up. I am planning to make it easier to use test directories in a package....
  • shahn: Hello Steve Instead of superimposing a star on the image to show the centroide. How would you superimpose a...
  • Kay-Uwe: Having TestSuite.fromPackag e() would be nice to have, but so far using simple “test” subdirs...

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