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
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.
Omar—Nice work.
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
Jim—Thanks for adding your comment showing how the syntax works for matrices.
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)
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.
Interesting,
this is a good article to help with algorithms…
keep up the good work
Thanks for writing about it
Hi Steve,
Thanks for the tips and explanations.
Is there any way to go back to the original image from the dilated image?
Say,
I have dilated as:
Now, I want to go back to the original data from the p_dilated data. How can I do it?
Any help would be appreciated.
Thanks.
Rashed—No, in general you can’t recover the original image from a dilated image.
hi Steve
i want to dilate the image which contains the Y of YCbCr to propagate the high values and erode for the low ones. please tell me hw could i perorm it.
im=imread('image058.jpg'); YCbCr= rgb2ycbcr(im); y=YCbCr(:,:,1); figure;imshow(y);how could i perfom dilation nd erosion on this.
regards
Hira—Use imdilate and imerode.
what sort of structuring element should i use???
after applying dilation and erosion i have to divide both images “dy”(dilated image) and “ey”( eroded image). dy/ey gives error. please guide how to divide the images.
regards
Hira—The choice of structuring element depends on what effect you are trying to accomplish. You haven’t said. For element-wise division use the ./ operator.
hi
actually i want to make a MapLum = dil(y)/ero(y) as doen in this paper “Face and Eyes Detection to Improve Natural Human-Computer Dialogue”.
for the MapLum will give high values around the eyes. i m trying to implement it but stucked in this MapLum as i m not getting the desired results as shown in the paper. after performing the dilation and erosion when i divide the images it gives me a totally black image. i dont know where is the problem. may be its the matter of the structuring element.
please help me. kindly take a look of the paper and guide me how i could perform it.the link for the paper is:
http://www.eurasip.org/Proceedings/Ext/ISCCSP2006/defevent/…/cr1070.pdf
Hira—I don’t generally have time to consult on individual image analysis or algorithm development projects, or to help interpret research papers. I suggest that you take a look at the numerical pixel values at each step. Is the output of the dilation what you expect? And the output of the erosion? And then check the output of the elementwise division. If that is not what you expect, then it may be a data type or scaling issue. Check out the MATLAB documentation on the how arithmetic operators work on integer data types, and check out the Image Processing Toolbox documentation on the default dynamic ranges for different image types and data types. Finally, your formula dil(y)/ero(y) is missing any mention of structuring element, so you might want to review the paper again to try to determine what structuring element they are using.
thanks
Hi Steve,
How can we apply morphological operators (esp. Dilation) to 3D triangulated meshes similar to :
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.36.9118&rep=rep1&type=pdf
Thanks
Danish—The Image Processing Toolbox supports morphology on arrays (rectlinear grids) only. You’ll need to develop your own implementations for other geometries.
Hi Steve, how can i get dilation using element of
hemispheric at selected scale.(sory for bad english)
Anton—Take a look at the ‘ball’ option of strel.