Gray Scale Morphology
To visualize the gray scale morphology operators, it is useful to think as a gray scale image modeled as a Digital Elevation Model (DEM). The pixels values represents the altitude of the terrain. In the illustration below the same image is displayed as gray scale on the left and as a surface on the right.
f=imread('ggray.tif');
mmshow(f);
mmsurfshow(f);
|
In the case of real images, the noise and the higher resolution make the visualization by surface plot very difficulty, unless the image is strongly smoothed. An alternative surface visualization is by using mmsurf, which generates a top-view of the surface with a lateral illumination. This display gives the appearance of an engraved picture.
f=imread('teclas.tif');
mmshow(f);
fs=mmsurf(f);
mmshow(fs);
|
Extending the Binary Operators to Gray Scale
There are two main models to extend the binary operators to work with gray scale images: by the thresholding decomposition, and by the subgraph of a function. These two models are useful to understand the properties and the geometrical interpretation of gray scale operators, not necessarily to achieve efficient implementation.
Threshold Decomposition
A gray scale image f(x) can be decomposed in many binary images (cross sections) by thresholding it at each gray scale level. A cross section at level t is given by the set of all pixels greater or equal t.
Any gray scale image can be uniquely reconstructed from its cross sections.
An important property of the threshold decomposition is known as the stacking property.
As a consequence of this stacking property, the threshold decomposition can be illustrated by the MATLAB plot with the option 'stack'. In the example below, a small 1-D gray scale signal f is first plotted as a bar graph. The signal is decomposed by its cross sections and stacked together using the option 'stack' of bar. Each color represents a cross section (binary image).
f=[3 0 5 4 5 6 3 1 4]'; bar(f); % F1 = f>=1; F2 = f>=2; F3 = f>=3; F4 = f>=4; F5 = f>=5; F6 = f>=6; bar([F1, F2, F3, F4, F5, F6],'stacked'); colormap(jet);
|
Stack Operators
Note that to reconstruct a gray scale image from its cross sections, the stacking property is crucil. Any binary operator that preserves the stacking property when applied to the cross sections can be extended to gray scale. This class of operators is know as stack operators. Note that any binary operator that is increasing preserves the stacking property and hence can be extended to gray scale using the thresholding decomposition model.
Since the union, intersection, dilation and erosion are increasing operators, they can be readly extended to gray scale. The example below illustrate the gray scale flat dilation using the threshold decomposition. Note that mmdil accepts both binary and gray scale images.
f=[3 0 5 4 5 6 3 1 4]'; bar(f); % F1 = uint8(f>=1); F2 = uint8(f>=2); F3 = uint8(f>=3); F4 = uint8(f>=4); F5 = uint8(f>=5); F6 = uint8(f>=6); bar(double([F1, F2, F3, F4, F5, F6]),'stacked'); colormap(jet); % bar(double([mmdil(F1), mmdil(F2), mmdil(F3), mmdil(F4), ... mmdil(F5), mmdil(F6)]),'stacked'); % bar(double(mmdil(uint8(f))));
|
Flat gray scale dilation
The thresholding decomposition allow us to think the stack operators as working in all possible cross sections of the image. This means that if a sequence of thresholding followed by a dilation is applied, it may be better to apply the dilation first and then the thresholding. The gray scale dilation stores the dilations of all possible thresholdings.
In the example below, the dilation is applied to three cross sections using a real 2-D gray scale image. First the original image is displayed both in gray scale and as surface model.
f=imread('teclas.tif');
mmshow(f);
fs=mmsurf(f);
mmshow(fs);
|
Three cross sections, at 30, 100, and 200 are overlayed on the surface view.
F30=mmcmp(f,'>=',uint8(30)); mmshow(mmsurf(f),F30); F100=mmcmp(f,'>=', uint8(100)); mmshow(mmsurf(f),F100); F200=mmcmp(f,'>=', uint8(200)); mmshow(mmsurf(f),F200);
|
The gray scale image is dilated by an Euclidean disk of radius 4.
B=mmsedisk(4); g=mmdil(f,B); mmshow(g); mmshow(mmsurf(g));
|
The three cross sections of the dilated image are overlayed on the surface view of the dilated image. These cross sections are the dilated versions of the original cross sections.
G30=mmcmp(g,'>=', uint8(30)); mmshow(mmsurf(g),G30); G100=mmcmp(g,'>=', uint8(100)); mmshow(mmsurf(g),G100); G200=mmcmp(g,'>=', uint8(200)); mmshow(mmsurf(g),G200);
|
Extending Gray Scale Operators Using Top and Subgraph
A gray scale image f(x) is a mapping of a subset of 
(called domain of f) into an interval of nonnegative integers K. A gray scale image can be model as sets using its subgraph. The subgraph of a function f, SG(f) is the set of points (x,t) below the function.
It is possible to recover the function from its subgraph using the top or top surface. The top of a set is given by:
Note that to convert a gray scale function to a set it is necessary to add an additional dimension. An 1-D signal is mapped into a 2-D binary image. In the same sense, a 2-D gray scale image is mapped into a 3-D binary volume. This model is useful for geometrical interpretation of gray scale operators, but may not be appropriated for efficient implementation.
In the example below, a subgraph of an 1-D gray scale signal is displayed. Two auxiliary functions mmsg and mmtop was created to implement the subgraph of a function and top of a set. A subgraph of a small 1-D signal is displayed using the MATLAB bar plot and also explicitly converted to a set (binary image) and converted back to a function.
f=[3 0 5 4 5 6 3 1 4]'; bar(f,1,'k'); A = mmsg(f,8); mmshow(A); fr = mmtop(A); fr' ans = 3 0 5 4 5 6 3 1 4 whos A fr Name Size Bytes Class A 8x9 72 uint8 array (logical) fr 9x1 9 uint8 array Grand total is 81 elements using 81 bytes
|
Gray scale dilation
The gray scale dilation is the top of the binary dilation of the subgraph of the input gray scale function.
In the example below, an ECG signal is processed. Its subgraph is created and displayed as an image. The SE has 10 pixels to the left and 30 pixels to the right (B). The dilation of the subgraph is displayed in white superposed by the subgraph of the original signal in red. Next, the top of the dilated subgraph is compared with the gray scale dilation to confirm that the model is correctly implemented. Finally, the original and dilated signals are visualized using the plot command. Note that you can think as the SE as stamping all the points of the subgraph, as done with the geometrical interpretation of the binary dilation. The SE can be visualized in the dilated signal at the peak of the signal.
f=imread('heart_uint8.tif');
Bimg=logical(uint8(zeros(1,61)));
Bimg(20:61)=1;
B=mmimg2se(Bimg);
SGf=mmsg(f');
mmshow(SGf);
SGdil=mmdil(SGf,B);
mmshow(SGdil, SGf);
gt = mmtop(SGdil)';
g=mmdil(f,B);
mmis(g,'==',gt)
ans =
1
x=1:size(f,2);
plot(x,f,'-r', x,g,'-b');
|
Non-flat Structuring Element
An advantage of the use of top and subgraph is its ability to work with structuring elements that are not sets (binary), but structuring function. A structuring function is a function mapping a domain given by a structuring element to integers. A structring function is also called non-flat structuring element.
Below is an example of a non-flat structuring element, it is 1-D semi-sphere. As this structuring function is 1-D, we can visualize it by its subgraph.
b = mmsedisk(20,'1D','EUCLIDEAN','NON-FLAT'); bf = mmseshow(b); SGb = mmsg(bf,20); mmshow(SGb);
|
The same experiment of processing the ECG signal using a flat SE is repeated below, but using the semi-spherical non-flat SE. First, a 2D SE is created by concatenating a zero semi-spherical below the subgraph. This is necessary as the origin coordinates of the image being converted to a SE is always at its center.
Bimg = [SGb; logical(uint8(zeros(19,41)))]; B = mmimg2se(Bimg); mmshow(mmseshow(B)); mmshow(mmseshow(B,'EXPAND')); SGg=mmdil(SGf,B); mmshow(SGg,SGf); g=mmdil(f,b); gt = mmtop(SGg)'; g=mmdil(f,b); mmis(g,'==',gt) ans = 1 plot(x,f,'-r', x,g,'-b');
|
The equation of the gray scale dilation by a non-flat SE is the following:
To understand why the addition used above is different than the arithmetic addition (see mmdil documentation), consider the following example, where part of the signal has value 0. The dilation using the subgraph is repeated here. In order to get the gray-scale dilation to reproduce this behavor at the values 0, the result of the addition must be zero if the image value is zero, no matter the value of the structuring function.
f1=f; f1(170:250)=0; SGf1=mmsg(f1'); mmshow(SGf1); SGg1=mmdil(SGf1,B); mmshow(SGg1,SGf1);
|
The visualization of 2D non-flat SE as an uint8 image is not simple. The solution adopted here is to visualize the response to an impulse of a gray scale image with value 1. In the example below, a city-block non-flat disk of radius 2 is created. To interpret its visualization by mmseshow, the zeros values does not belong to the domain, so their don't enter in the operation. It is equivalent if their values were minus infinity. The other values must be subtracted by one.
b=mmsedisk(2,'2D','CITY-BLOCK','NON-FLAT');
mmseshow(b)
ans =
0 0 1 0 0
0 1 2 1 0
1 2 3 2 1
0 1 2 1 0
0 0 1 0 0
|
Exercises
- Give an example of a binary operator that does not build a stack operator.
- Give the geometrical interpretation of the following equation of the flat gray scale dilation.
- Build a MATLAB script to illustrate the geometrical interpretation of the gray scale dilation, using the same technique applied in the binary dilation.
- How can you test if a SE is flat or non-flat?
- Create a function to display the non-flat SE with the MATLAB
-Inf value when the points do not belong to the domain of the SE and the correct values for the rest of points of the structuring function.
-Inf value when the points do not belong to the domain of the SE and the correct values for the rest of points of the structuring function.