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



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
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
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.
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.
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?
Anyway, some of the material in my series on connected component labeling algorithms might be useful to you.
Ken—I get similar results. Clearly, we have a lot of room to improve our dilation implementation! We’re looking into it.
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.
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.