#### Fuzzy Morphological Polynomial Image Representation

EURASIP Journal on Advances in Signal Processing
Hindawi Publishing Corporation
Fuzzy Morphological Polynomial Image Representation
Chin-Pan Huang 1
Luis F. Chaparro (EURASIP Member) 0
Srdjan Stankovic
0 Department of Electrical and Computer Engineering, University of Pittsburgh , Pittsburgh, PA 15261 , USA
1 Department of Computer and Communication Engineering, Ming Chuan University , Taoyuan 333 , Taiwan
A novel signal representation using fuzzy mathematical morphology is developed. We take advantage of the optimum fuzzy fitting and the efficient implementation of morphological operators to extract geometric information from signals. The new representation provides results analogous to those given by the polynomial transform. Geometrical decomposition of a signal is achieved by windowing and applying sequentially fuzzy morphological opening with structuring functions. The resulting representation is made to resemble an orthogonal expansion by constraining the results of opening to equate adapted structuring functions. Properties of the geometric decomposition are considered and used to calculate the adaptation parameters. Our procedure provides an efficient and flexible representation which can be efficiently implemented in parallel. The application of the representation is illustrated in data compression and fractal dimension estimation temporal signals and images.
1. Introduction
Signal representation is an area of great interest in the
signal and image processing. Many representation techniques
currently available are well developed and offer
satisfactory performance in many applications [
1–6
]. However, in
general, drawbacks of these techniques include intensive
computation, sequential implementation, and disregard of
geometrical information present in signals. In this paper, we
propose fuzzy mathematical morphology [7] to represent
one- and two-dimensional signals. Fuzzy morphological
operators, similar to morphological operator, are nonlinear
but well suited for efficient implementation in parallel.
Furthermore, they allow to extract geometrical information
in signals by appropriate transformations.
Among recently introduced representation techniques,
Martens [
5, 6
] proposes a linear combination of polynomial
to represent signals. Although the method achieves high
performance in data compression, it has high computational
complexity and a sequential implementation. Pitas and
Venetsanopaulos [
8–10
] propose a morphological signal
decomposition method to decompose a signal into a set of
morphologically simple function. Song and Delp [11] use
multiple structuring functions instead of a single function
to enhance the performance of morphological filters. In
developing the Fuzzy Morphological Polynomial (FMP)
representation [
12
], we take the advantage of polynomial
transform idea [
5, 6
], the morphological decomposition
recursive procedures [
8–10
] and using multiple structuring
functions [11] overcoming some of the problems mentioned
before.
Binary morphology, as developed in [
13
], is based on
the concept of fitting a structuring element to the signal. Its
extension to multilevel morphology was achieved by treating
the space underneath the signal (“umbra”) as a binary signal.
Recently, Shinha and Dougherty [
7
] propose an alternate
mathematical morphology based on fuzzy set theory. The
morphological operations are modeled on a “fuzzy” notion
of fitting, the umbra concept is not required and as such
binary morphology becomes a special case. The fuzzy fitting
yields one or zero for crisp fitting, and between zero and
unity for a partial fit. The closer to unity, the higher the
degree of fit.
The rest of the paper is organized as follows. In Section 2,
we briefly review fuzzy mathematical morphology. Then in
Section 3, we develop the one-dimensional FMP
representation based on a recursive geometric decomposition for a
given signal membership function. The properties of our
FMP representation are investigated to develop algorithms
to compute the adaptive parameters. In Section 4, we extend
our algorithm to two-dimensions. We use the tensor product
of the two one-dimensional functions as fuzzy
structuring functions. The order of the two-dimensional fuzzy
structuring functions is explored and a rational solution
is recommended. In Section 5, we present our experiment
results of applying our algorithm to data compression and
fractal dimension estimation for one- and two-dimensional
signals, demonstrating the advantages of our algorithm. The
conclusion is in Section 6. Some of the results in this paper
were presented before in [
12
].
2. Fuzzy Mathematical Morphology
Recently, Shinha and Dougherty [
7
] proposed to consider
fuzzy set theory [
14
] instead of the classical set theory
to develop mathematical morphology. They have in fact,
obtained a new approach that considers simultaneously
binary and multilevel morphology. The concept of “umbra”
is no longer needed to develop the multilevel case.
Morphological operations are then developed on the “fuzzy” fitting
so that for crisp sets the fitting still remains characterized
as either 0 or 1, but for fuzzy or noncrisp sets it is possible
to have a fitting characterized by a value between 0 and 1.
The closer to unity, the better the fitting of the structuring
element. As in the classical morphology, fuzzy morphology
[
7
] also consists in transforming a fuzzy set into another.
Such a transformation is performed by means of a fuzzy
structuring set containing the desired geometric structure.
If we let X be the universe of discourse and x be its generic
element, the difference between crisp and fuzzy sets is the
characteristic function of a crisp set C which is defined as μC :
X → {0, 1} while the membership function μF : X → [0, 1]
of a fuzzy set F is defined so that μF (x) denotes the degree to
which x belongs to the set F. Among the different operations
on fuzzy sets [
15
], the following are important operations to
be used later.
(a) Complement operation:
(b) Translation of a fuzzy set F by a vector v ∈ X:
(c) Reflection of a set F:
(d) Bold union of two sets F and G:
(e) Bold intersection F∇G:
μFΔG(x) = min 1, μF (x) + μG(x) .
μF∇G(x) = max 0, μF (x) + μG(x) − 1 .
The degree of fitting of a set A into a set B is measured
by an inclusion grade operator
I(A, B) = xin∈fX μAcΔB(x)
= 1 + min 0, inf x∈X μB(x) − μA (x) ,
where Δ is the bold union operator. According to
the above index, the degree of subsethood of two
crisp sets A, B is either 0 or 1, while for fuzzy sets
C and D I(C, D) ∈ [0, 1]. Moreover, if C ⊆ D then
I(C, D) = 1 and in general 0 ≤ I(C, D) ≤ 1.
Using such an index [
7
] has shown that the erosion
operation can be defined, and from it the dilation,
opening, and closing operators are obtained. In fact,
if f (n) is a multilevel and k(n) is a structuring
element with supports F and K and membership
function μ f (n) and μk(n), then we have
Erosion:
μ f Θk(n) = I T (k; n), f
Dilation:
μ f ⊕k(n) = μ( f cΘ−k)c (n)
Opening: Closing:
μFc (x) = 1 − μF (x).
μT (F;v)(x) = μF (x − v).
μ−F (x) = μF (−x).
(6)
(7)
(8)
(9)
(10)
(11)
(12)
(13)
= mi∈iKn min 1, 1 − μk(i) + μ f (n + i) ,
= mi∈aKx max 0, μk(i) + μ f (n − i) − 1 ,
μ f ◦k(n) = μ( f Θk)⊕k(n),
μ f •k(n) = μ( f ⊕k)Θk(n),
2.1. Fuzzification. To apply the above fuzzy morphological
operators, the multilevel signal must be converted to its
membership function. Thus, in case of an image, the
membership function will determine the degree of belonging
to categories ranging over different intensities. The
membership function μ f thus maps the image intensity range R into
[0, 1],
(1)
(2)
(3)
(4)
(5)
S-function
⎧
⎪⎪⎪⎪ 2
μS f ; a, b, r = ⎪⎨⎪
f − a
r − a
2
,
⎪
⎪⎪⎪⎪⎪⎩ 1 − 2
f − r
r − a
a ≤ f ≤ b,
2
, b ≤ f ≤ r,
μ f : R −→ [0, 1]
f x, y −→ [0, 1],
μL f ; a, r = rf −− aa ,
according to some criteria. Some of them are the Linear (L)
function
and π-function [
16, 17
]. In signal representation, the
fuzzification needed to be single valued and as such we will
not consider the π-function fuzzification. We propose a new
fuzzification method, called Z-function fuzzification, which
has the inverse effect of the S-function fuzzification. The
Zfunction is
a ≤ f ≤ b,
, b ≤ f ≤ r,
where μ is the membership function and f , a, b, and r are
the grey level, and the minimum, middle, and maximum
value of the image grey level. Notice that these three
functions are all single valued and monotonically increasing
in the analysis interval (minimum to maximum of the
signal). The fuzzification functions for image is shown in
Figure 1. Given that the fuzzification functions chosen are
single valued, the defuzzification process is easily achieved by
the inverse mapping.
The application of these fuzzification techniques in signal
analysis varies. For instance, the S-function enhances the
image contrast about the b level, and as such it could be used
for edge enhancement. On the other hand, the Z-function
decreases the contrast or blurs the image about the b level.
The linear fuzzification does not alter the contrast, it simply
normalizes the image to a range of [0, 1].
To illustrate the application of S-function in edge
detection, consider the fuzzy morphological gradient (FMG)
μd(n) = μ f ⊕k(n) − μ f Θk(n)
(15)
obtained as an extension of the classical morphological
gradient first proposed by Serra [
18
] and Evans and Liu
[
19
]. Applying the FMG to the real signal of an image in
Figure 2(a), we obtain the lower figure where the peaks
correspond to the edges. To enhance these edges, we
consider the fuzzification of this result to achieve considerable
enhancement of the edges as shown in Figure 2(b).
3. Fuzzy Morphological Polynomial (FMP)
Representation
The FMP representation is analogous to the morphological
polynomial transform [
20
] and the orthogonal polynomial
representation [
5
]. Using fuzzy morphological opening we
obtain a representation similar to a polynomial
representation by means of a geometrical decomposition of the
signal. One of the difficulties encountered in the process
was the selection of the structuring functions, which can
be either arbitrary or derived from the signal. In our case,
we get them from a complete set of ordered real-valued
orthogonal polynomials in 0 ≤ n ≤ N − 1. In the
examples, we use the discrete Legendre orthogonal (DLO)
polynomials [
21
]. It should be noted that the corresponding
membership functions are not necessarily orthogonal. Let
μ f (n) be the membership function of a given signal, and
μki (n) : 0 ≤ i ≤ N − 1} be one-dimensional fuzzy structuring
functions, such that 0 ≤ μki (n) ≤ 1. Let {ai} be adaptive
parameters used to make the fuzzy structuring function
fit μ f (n) closely. To consider all possibilities, the fuzzy
structuring functions {μki (n)} are derived from a shifted and
normalized set of orthogonal polynomials {μgi (n)} and its
complementary functions {μgic (n)}. Figure 3 illustrates the
shifted and normalized functions {μgi (n)} when we consider
the discrete Legendre orthogonal polynomials for N = 5.
The geometric decomposition of the given membership
function μ f (n) is obtained recursively as follows.
(14)
(i) Windowing with W (n):
(16)
(17)
(18)
(19)
(20)
1
0.9
0.8
0.7
0.6
(ii) Adaptive recursive approximation of μz0 (n):
μzv0 (n) = μ f (n) × W (n − vN ).
v v
μzi+1 (n) = μzvi (n) − μzi◦aiki (n),
where i = 0, 1, . . . , N − 1 relates to the structuring
functions aiμki (n), v = 0, 1, 2, . . . refers to the window
W , and ai are in [0, 1] adaptive parameter. Each
window is processed similarly.
The term a0μk0 (n) is very important in the above
decomposition as it provides a coarse approximation to the
signal membership function while {aiμki (n), i > 0} gives the
fine information of μzv0 (n). Applying (17) recursively we have
μzv0 (n) =
v
μzi◦aiki (n) + μzvN (n),
where the last term corresponds to the residual or the part of
the signal that cannot be well represented with N function
μki (n). We will show in the next section that the above
representation can be considerably simplified by choosing
values of {ai} such that
v
μzi◦aiki (n) = aiμki (n)
to convert (18) into
N−1
i=0
N−1
i=0
μzv0 (n) =
aiμki (n) + μzvN (n),
0.8
0.7
0.9
0.8
0.7
representation analogous to a polynomial representation of
the windowed signal.
Properties. The following propositions will give insight on
how the FMP representation works and how to develop
formulas to calculate the {ai} coefficients. Here, we work
on a frame signal only, and thus the superscript v can be
omitted. Also, we use min , max to stand for min0≤ ≤N−1
and max0≤ ≤N−1. In the following propositions, we assume
μz(n), μk(n) are both defined on 0 ≤ n ≤ N − 1, and
{ai} ∈ [0, 1] then
Proposition 1. The opening
v
μzi◦aiki (n) = max 0, aiμki (n) + μc − 1 ,
(21)
where μc = 1 + min{0, min [μzi ( ) − aiμki ( )]}.
n
(b)
This proposition provides a simplification of the
nonlinear opening and yield μc, an index of the degree of fitting of
the structuring function in the signal membership function.
Proposition 2. There exists an optimum {ai} ∈ [0, 1]
(denoted as ai∗) such that μzvi◦ai∗ki (n) = ai∗μki , 0 ≤ n ≤ N − 1,
if and only if the following optimum condition is satisfied
The value ai∗ is calculated as
min μzi ( ) − ai∗μki ( ) = 0.
ai∗ =
min
0≤l≤N−1
μki (l) =/ 0
μzi ( )
μki ( )
.
Let J(ai) = min [μzi ( ) − aiμki ( )] be the fitting cost function,
and refer to
J ai∗
= 0
as the optimum condition. There are two direct corollaries
that follow from this proposition.
Corollary 1. If the optimum condition is met, then, there are
the following.
(i) If ai > ai∗, then J(ai) < 0, 0 ≤ n ≤ N − 1,
(ii) If ai < ai∗, then J(ai) > 0, 0 ≤ n ≤ N − 1,
(iii) ai∗ = minn[μz0 (n)].
(iv) 0 ≤ ai∗ ≤ max[μzi (n)] ≤ 1, 0 ≤ i ≤ N − 1.
Corollary 2. If the optimum condition is met for μz0c (n), then
a0∗c = max μz0 (n) .
n
Proposition 2 gives a solution to the optimal fitting
problem. In the next section, we develop a geometrically
intuitive solution.
(22)
(23)
(24)
(25)
Proposition 3. For μzi (n), i = 0, 1, . . . , N − 1 in (17), the
following relation holds
(i) 0 ≤ μzi+1 (n) ≤ μzi (n) ≤ 1,
This proposition shows that the residual at each
decomposition step satisfies the membership function conditions,
decreases, and at least in one point is zero.
Proposition 4. If the optimum condition is met for ai and ai ,
then the only value of ai is zero to satisfy the following equations
μzi◦aiki (n) = aiμki (n), μzi ◦ai ki (n) = ai μki (n), for 0 ≤ i ≤ N −1,
where μzi (n)
μzi (n) − aiμki (n).
This proposition establishes that once geometrical
features are decomposed from an input membership function,
opening with previously used fuzzy structuring function
gives zero.
3.1. Implementation. The representation of μ f (n) according
to (20) requires the calculation of the adaptive coefficient
{ai} and choosing the appropriate structuring functions.
There are two possible methods to find the coefficients {ai},
the first one is based on a geometrically intuitive argument,
while the second uses Proposition 2 given before. Consider a
given window, and for simplicity let us not indicate it in the
equation.
3.1.1. Calculation of {ai}
Iterative Method. Let aiq be the qth iteration of ai (ai∗ as
the optimum) when attempting to optimize the fitting cost
function
J(ai) = 0≤m≤iNn−1 μzi ( ) − ai μki ( ) ,
with respect to ai, the optimum value ai∗ is obtained when
the cost is zero, and such that μzi◦ai∗ki (n) = ai∗μki (n).
According to (26) and Corollary 1, the following
algorithm can be used to find the optimum value of ai:
ai0 = max μzi (n) ,
J aiq = 0≤m≤iNn−1 μzi ( ) − aiq μki ( ) ,
aiq+1
= aiq + J aiq ,
where q ≥ 0. That J(aiq) converges toward zero when q
increases can be established. If we assume that aiq > ai∗, then
J(aiq) ≤ 0 (Corollary 1), and, thus
q+1
J ai
min
= 0≤ ≤N−1
min
≥ 0≤ ≤N−1
q q
μzi ( ) − ai + J ai
μzi ( ) − aiq μki ( ) = J aiq .
Direct Method. As in Proposition 2, we can compute
ai∗ using (23). This method gives the same ai∗ as the iterative
method, but in a faster way.
(29)
(30)
(31)
(32)
(33)
(34)
3.1.2. Choosing Structuring Functions. Given shifted and
normalized orthogonal polynomials {μgi (n)} and their
complements {μgic (n)}, we need to determine which of these two
should be used as the structuring functions {μki (n)} for the
representation. This needs to be done due to the positivity
condition on the adaptation coefficients {ai}. We decide this
by comparing the reconstruction error corresponding to the
coefficients attached to each of these structuring functions.
Let the coefficients ai and ai be the optimum values
for μgi (n) and μgic (n), respectively. We want to choose the
optimum value which gives us the smaller reconstruction
error. Let the reconstruction error membership function
corresponding to ai and ai be, respectively,
μe1 (n) = μzi (n) − aiμgi (n),
μe2 (n) = μzi (n) − aiμgic (n).
It can be easily shown that
where αi is found to be
In that case, we then let αi = ai and μki (n) = 1 − μgi (n).
Otherwise, we choose αi = ai and μki (n) = μgi (n).
4. Two-Dimensional Fuzzy Morphological
Polynomial Representation
The FMP representation can be easily generalized to two
dimensions. Let μ f (m, n) be the given signal membership
function and, μk∅ (m, n), 0 ≤ ∅ ≤ MN − 1 be ordered
twodimensional fuzzy structuring functions, based on
orthogonal polynomials on 0 ≤ m ≤ M − 1, 0 ≤ n ≤ N − 1. The
geometrical decomposition algorithm becomes
μu,v(m, n) = μ f (m, n) × W (m − uM, n − vN ),
z0
μzu∅,v+1 (m, n) = μzu,v(m, n) − μzu∅,v◦a∅ μk∅ (m, n),
∅
where ∅ = 0, 1, . . . , MN −1 is an index related to the adaptive
tshtreubctloucrkinbgeifnugncctoionnsidae∅reμdk.∅T(hme,anb)o;vaenpdrouc,evduarree isinrdeipceeasteodf
until the residual is insignificant or ∅ = MN − 1. Each
αi
N−1
n=0
αi =
N−1
n=0
μgi (n) =
μgic (n),
N
n=0 μgi (n) − 1
N−1
ai ≤ aiαi,
N−1
n=0
N−1
n=0
(notice α0 = 0 and 1 − μk0 (n) = 0, so i ≥ 1). If
then we have that
μe1 (n) ≥
μe2 (n).
block is decomposed similarly. As in (18)–(20), our
twodimensional FMP representation for a frame signal is
MN−1
∅=0
μu,v(m, n) =
z0
a μu,v (m, n) + μzuM,vN (m, n).
∅ k∅
(35)
Notice that the ∅, abbreviation of ∅(i, j), is an
ordering function used for the two-dimensional structuring
functions. The properties of one-dimensioned FMP can
be extended to two dimensions easily, thus, we omit the
derivation here. As in the one dimension, the optimum
condition in the two dimensions is
min μz∅ (s, t) − a∗∅μk∅ (s, t) = 0,
s,t
(36)
where a∗∅ is an optimum value to satisfy this equation.
It is understood from the previous section that the
one-dimensional algorithm can be extended to
twodimensional provided the generation and ordering of the
two-dimensional structuring functions are determined.
Separable and nonseparable bivariate orthogonal polynomials
may be used to generate two-dimensional structuring
functions. Consider the separable structuring functions
μki,j (m, n) = μki (m)μkj (n)
(37)
obtained from the one-dimensional structuring functions.
Figure 4 shows an example of two-dimensional separable
structuring functions with size of 5 × 5.
An inherent problem in two dimensions is the ordering
of the structuring functions, which in the one-dimensional
case occurs naturally. Our approach first investigates the
structuring function properties and establishs the possible
guides to the order of the two-dimensional structuring
functions and then comes out a procedure to get the solution.
Proposition 5. If one multiplies two normalized one variate
structuring functions, derived from the discrete Legendre
orthogonal polynomials, μki (m) of size M and μkj (n) of size N
as two-dimensional structuring function, that is, μki,j (m, n) =
μki (m)μkj (n) with dimension of M × N , then one can have the
following properties:
(i) μk0,j (m, n) ≥ μki,j (m, n) for all i, j, m, n.
(ii) μki,0 (m, n) ≥ μki,j (m, n) for all i, j, m, n.
(iii) μk0,0 (m, n) ≥ μki,j (m, n) for all i, j, m, n.
This proposition gives properties of the separable
two-dimensional structuring functions based on
onedimensional DLO polynomials.
Proposition 6. If optimum condition is met for a∗∅(i,j) and
a∗∅(s,t) and μk∅(i,j) (m, n) ≤ μk∅(s,t) (m, n) for a pair (i, j), (s, t)
and for all m, n, then the only value of a∗∅(s,t) is zero to satisfy
the following equations:
μz∅(i,j)a∗∅(i,j) μk∅(i,j) (m, n) = a∗∅(i,j) μk∅(i,j) (m, n),
μz∅(s,t)a∗∅(s,t) μk∅(s,t) (m, n) = a∗∅(s,t) μk∅(s,t) (m, n),
(38)
that 0 ≤ μz∅(s,t) (m, n) ≤ 1.
where μz∅(s,t) (m, n) μz∅(i,j) (m, n) − a∗∅(i,j) μk∅(i,j) (m, n) −
(m, n), with 0 ≤ (m, n) ≤ 1, for all m, n and arbitrary such
This proposition gives guides of the structuring function
order. If the structuring function has higher amplitude values
than the other one at every point, then the structuring
function must have a lower order otherwise the decomposition
gets zero.
By observation, the factor in addition to the amplitude
which affect the natural order of the one-dimensional
structuring function is its complexity. Based on the properties and
the observation, we come up with our rationale solution.
We first define the amplitude and complexity index to
quantitatively measure the structuring function amplitude
and complexity characteristic.
Definition 1. An amplitude index (AI) of a geometrical
structuring function μki,j (m, n) of size M × N is defined by
the summation of amplitude at every pixel as
N−1 M−1
Definition 2. A complexity index (CI) of a geometrical
structuring function μki,j (m, n) of size M × N is defined by
the distance between the adjacent pixel in the horizontal,
vertical, and diagonal directions as
CI i, j =
μki,j (m, n) − μki,j (m, n + 1)
2
2
μki,j (m, n) − μki,j (m + 1, n)
μki,j (m, n) − μki,j (m + 1, n + 1)
2
μki,j (m, n) − μki,j (m − 1, n + 1) .
2
(40)
M−1 N−1
The effects of the amplitude and complexity index of
the structuring function to order are that the smaller the
amplitude index, the higher the order of the structuring
function, and the greater the complexity index, the higher
the order of that structuring function. We combine these
definitions and their effects and come out, the structuring
index (SI).
Definition 3. A structuring index (SI) of a structuring
function μki,j (m, n) with amplitude index AI(i, j) and complexity
index CI(i, j) is defined as
SI i, j = 0.5 N AI i, j −1 + N CI i, j
,
(41)
where 0.5 is average factor and N (·) is a normalization
operator over the structuring functions.
i
Notice that the value of SI is in [0, 1].
Based on SI(i, j), the order of the fuzzy structuring
functions is
We can get the order of the two-dimensional fuzzy
structuring functions from one variate DLO polynomials. As an
example, this will order 5 × 5 two-dimensional structuring
functions as Table 1. Our ordering method, considering both
the amplitude and intrinsic geometrical complexity of the
(42)
5. Applications
structuring function, is more reasonable than the commonly
used close-neighbor ordering method [
22
], considering only
the index of the structuring functions.
In this section, we show how the FMP representation is
applied to data compression and fractal dimension
estimation. We compare the data compression results with those
using the discrete cosine transform (DCT) method [
23
], and
the fractal dimension estimation results with those using
morphological covering (MC) method [
24
] and differential
box-counting (DBC) method [
25
].
L
4
3
5
3
5.1. Data Compression. The application of the FMP
representation for data compression is shown in Figure 5. The
signal f (n) is fuzzified by F and then processed by the FMP
decomposition to get the adaptive coefficients. The signal
membership function is reconstructed and f (n) is recovered
by the defuzzifier D . The block diagram for two dimensional
signals is similar. The pepper image with size of 512 × 512
shown in Figure 6 is used as a test image.
In the first example, the one-dimensional FMP algorithm
is used to process the test image horizontally. We consider
different fuzzification methods, window lengths. The fuzzy
structure functions are obtained from DLO polynomials as
shown before. The performance of our representation is
evaluated by the “peak-to-peak” signal to noise ratio (SNR
dB) and the entropy-based compression ratio (ECR). The
entropy-based compression ratio (ECR) is defined as
total l.c.of bits of compressed signal
total l.c.of bits of original signal
=
iN=−01 Mili ,
MT lT
(43)
where N is the number of subblock signals, Mi is the number
of samples of the subblock i, li is the bits/sample required to
code subblock i, lT is the bits/sample required for the original
signal, MT is the total number of samples of the original
signal. The average bits/sample li required to code a subblock
signal is defined by entropy as
G−1
j=0
li = −
p j log2 p j ,
(44)
where p j is a probability of a sample with amplitude j, G
is the greatest amplitude of the signal. In Table 2,
signalto-noise ratio (SNR dB) values for different fuzzification
methods and window sizes are shown. In Table 3, the
entropy-based compression ratio (ECR) for different
fuzzification methods and window sizes are shown. These results
show that our one-dimensional algorithm has a high data
compression when using L or Z fuzzification methods. The
results also indicate that by using the Z fuzzification we
achieve a higher compression ratio with good SNR than
those results with L fuzzification.
In the second example, we apply our two-dimensional
FMP algorithm to process block by block the test image. The
fuzzy structuring functions are generated by multiplying two
one-dimensional structuring functions derived from DLO
3
polynomials. The order is determined by the structuring
index method discussed before. In Table 4, we show the
signal-to-noise ratio (SNR dB) for different fuzzification
methods and window sizes. In Table 5, the entropy-based
compression ratio (ECR) for different fuzzification methods
and window sizes is shown. Those results indicate that our
two-dimensional algorithm has a higher performance when
using L and Z fuzzification methods. The results also indicate
that the Z fuzzification achieves a higher compression ratio
with good SNR than the L fuzzification. To illustrate the
results, we show in Figures 7(a) and 7(c) and Figures 7(b)
and 7(d) the FMP component and the corresponding error
images, when using window size of 1 × 3 and 3 × 3, for
Lfunction fuzzification.
Since there is no published papers, as we know, using
the fuzzy morphology approach for data compression, we
want to show the high performance of our algorithm by
comparing with those obtained by commonly used high
performance method such as discrete cosine transform
(DCT) [
23
]. In order to have a fair comparison, we use
the coefficients with high energy of the DCT and FMP
such that the “peak-to-peak” SNR of the reconstructed
f (n)
F
μzv1 (n)
image which is closely used in [
26, 27
], then we compare
their data compression ratio at each different window size.
The comparisons of SNR (dB) and ECR are shown in
Tables 6 and 7 for one and two dimension, respectively.
Those results indicate that our FMP representation using Z
fuzzification achieves the higher data compression ratio that
of DCT. We also provide, as an example, the one- and
twodimensional reconstructed images in Figure 8 using the FMP
representation with Z fuzzification and the DCT method,
with window size of 3 and 3 × 3 cases, for visual quality
comparison purpose.
Figures 8(a) and 8(c) and Figures 8(b) and 8(d) show
the reconstructed images for one and two dimensions,
respectively. Figures 8(a) and 8(b) show the reconstructed
images using FMP with Z fuzzification. Figures 8(c) and
8(d) show the the reconstructed images using DCT. Those
results indicate that the visual quality of using our FMP
representation with Z fuzzification and DCT method is
similar. However, our FMP representation with Z fuzzification
achieves higher data compression ratio (fewer bits/pixel).
For computation complexity comparison, Table 8 shows the
required operations including addition/subtraction,
multiplication/division, and minimum/maximum and look-up
table/check-up for FMP and DCT. Assume a window size
of N is used. Notice that the multiplication/division and
addition/subtraction operations contribute more
computation complexity than the other operations. We can see clearly
that the FMP requires less number of operations than that of
the DCT.
5.2. Fractal Dimension Estimation Using FMP. To estimate
the fractal dimension (FD) of one- or two-dimensional
signals, we obtain the FMP representation of signal frames
of increasing dimension. This representation is used to find
an approximate cover of the windowed signal. It will be
shown that when using the discrete Legendre polynomials as
structuring functions, the covers can be obtained recursively
providing a better FD estimation at each recursion. Our
procedure can be easily extended to two dimensions.
5.2.1. One-Dimensional Estimation. Let μ f (n), 0 ≤ n ≤
M − 1, be the membership function of a signal f (n). For
each of the signal frames, we will use {μki (n)}, 0 ≤ n ≤
[0, N − 1], as structuring functions based on the discrete
Legendre polynomials (see Figure 3 when N = 5). The
FMP algorithm provides the adaptation parameters {ai}, for
each frame membership function μz(n). Let then the support
length of μ f (n) be S = M − 1 (if we know the sampling
period Ts then S = (M − 1)Ts) which is divided into an
integer number of windows of increasing length r = N − 1
(or r = (N − 1)Ts when Ts is known).
For each of the frames, we will attempt to come up with
a cover that encloses the signal as tightly as possible. The
length of the cover can be calculated recursively from the
FMP adaptive coefficients and the values of the structuring
functions. For a frame v with corresponding length r, we
have that (See Figure 9 when r = 4)
r
0v(r) = S ,
r
j=1
iv(r) =
μki j − μki j − 1
2 aiv 2 +
iv−1(r) 2
r
,
(45)
where 0v(r) is the geometric length corresponding to the
constant FMP decomposition, and i = 1, 2, . . . , I, I ≤
r corresponds to the ith FMP geometric decomposition.
j = 1, 2, . . . , r corresponds to the point of the structuring
function. The length calculation is done in each of the r
segments in which the window is divided (See Figure 9). The
height of the the cover dv(r) ∈ [0, 1] can be found exactly in
the case where μzv(n) is either a constant or has a great deal
of variation. In the first case, dv(r) = 0, and in the second,
dv(r) = 1. However, it should be noted that in these two
cases, we will have aiv = 0, i > 0 that is only a constant
approximation is possible. In cases different from the above
ones, we cannot calculate dv(r) exactly although a good
estimate of it can be obtained as the difference between the
maximum and minimum of the residual μzvI+1 (n) (see (18))
obtained after the Ith decomposition. Thus, the diagonal of
the cover in frame v
obtained by doing the the linear fitting piecewise, obtaining
a better estimate for linear fitting for small rs.
(2) If μ f (n) is so smooth that the FMP gives aiv = 0,
0 < i ≤ I, for every v, then we have that dv(r) ≈ 0, and
v(r)/ 0v(r) = 1 and, therefore, Q(r) = S/r for any value of r.
The estimation of the slope of log(Q(r)) versus log(1/r) for
two different values of r gives that the fractal dimension D is
one.
(3) On the other hand, if μ f (n) varies widely everywhere,
then the FMP gives aiv = 0, 0 ≤ i ≤ I, for every v, but due to
the large variation dv(r) ≈ 1, in which case Iv(r) = 0v(r). We
will then get that
v(r)
0v(r) =
1 +
S2 S
r2 ≈ r
(48)
Using least-square fitting in the log(Q(r)) versus log(1/r)
graph for various values of r, the slope of the line will
correspond to an estimate of the FD.
Remarks. (1) In order for the above discrete algorithm to
work properly, we need to choose r so that S/r is an integer
for every chosen value of r. Results vary for different choices
of r due to the global linear fitting. A better estimate might be
after substituting dv(r) and 0v(r), and using the fact that
S/r 1 (i.e., we divide S into several windows) which will
give us Q(r) = S2/r2, so that the slope estimation for two
different values will give a fractal dimension close to two.
(4) According to 2 and 3 above mentioned we have that,
except for very smooth and very rough signal, the FD will be
a value between 1 and 2.
D = 1.8
D = 1.5
D = 1.2
from which estimation of the FD can be done as before.
To test our algorithm, artificial fractal signals with known
FD are generated using the weistrass cosine function (WCF)
method [
24, 28
]. Figure 10 shows some WCF generated
signals, and when we apply our procedure to them. The
results are shown in Table 9. For comparison purpose, the
results of using MC method [24] are also shown in the table.
5.2.2. Two-Dimensional Estimation. Different from the
onedimensional case, the estimation of the FD of a
twodimensional signal can only be done using the constant and
linear term of the FMP approximation. Using just the first
three terms of the approximation, we have for an r × r block
of a square image of dimension S × S that get the following
recursive formulas:
A0u,v(r) =
r 2
S
where {aiu,v} are coefficients found using the FMP. Notice
that these equations are similar to those in the
onedimensional case, except that in this case the block is
not subdivided as we did with the window in the
onedimensional case. A cube covering the signal can be thought
of having a base with area Au,v(r) and an approximate height
I
of du,v(r) equal to the difference between the maximum and
the minimum of μzu,v(n). The diagonal plane of this cube has
an area equal to
Au,v(r) =
AIu,v(r) 2 + AIu,v(r)(du,v(r))2,
a formulation analogous to that in the one-dimensional case,
when we use only two components
Similar comments as those made in the one-dimensional case
can be made here. As before, S/r must be an integer and the
FD varies for different choices of r. Likewise a very smooth
signal has an FD close to 2, and a very rough signal has an FD
close to 3.
Finally if in the FMP approximation we only use the a0
component, the above algorithm simplifies to
A0u,v(r) =
r 2
S
,
Au,v(r) =
A0u,v(r) 2 + A0u,v(r)(du,v(r))2,
(53)
from which we can find an estimate of the FD of the given
two-dimensional signal.
We then apply our algorithm to estimate the FD of the
Brodatz texture images [
29
] shown in Figure 11. The code of
texture image is same as in the Brodatz album. The results
are shown in Table 10. This example shows the applicability
of our algorithm to estimate FD to images. We also show the
results of using DBC method [
25
] for comparison purpose.
The results of using these two methods are similar.
6. Conclusion
A novel signal representation using fuzzy morphological
approach has been proposed in this paper. Using the fuzzy
morphological operators and a set of structuring functions,
a decomposition and reconstruction procedure, similar to
polynomial transform [
5, 6
], is developed for the fuzzified
signals. Through using the fuzzy morphological approach, a
signal can be efficiently represented with several additional
advantages, such as lower computation complexity and
flexible in using the different fuzzification methods to extract
signal geometrical information for a better signal
representation. Furthermore, our representation can be
implemented very fast by parallel. We successfully use the fuzzy
mathematical morphology [7] approach to extend the work
D03
D24
D25
D03
2.55
2.60
D04
D28
D68
D05
D33
D84
of the Pitas and Venetsanopoulos [
8–10
] and of Song and
Delp [
11
] on morphological signal representation. We have
applied our representation to data compression and fractal
dimension estimation for one- and two-dimensional signals,
the experimental results have shown the high performance
in data compression and applicability in estimating fractal
dimension as compared with those using the DCT [
23
], MC
[
24
], and DBC [
25
] methods.
Appendix
Proof of Proposition 1. By definition
μzi◦aiki = μ(ziΘaiki)⊕aiki .
(A.1)
According to fuzzy morphological erosion
μziΘaiki (n) = m,ni+n min 1, 1 − aiμki ( ) + μzi (n + )
= min min 1, 1 − aiμki ( ) + μzi ( )
= μcδ(n).
Now
μc⊕aiki (n) = m,na−x max 0, aiμki ( ) + μcδ(n − ) − 1
= max 0, aiμki (n) + μc − 1 .
Thus, μzi◦aiki (n) = max[0, aiμki (n) + μc − 1].
D09
D54
D92
D92
Proof of Proposition 2. From Proposition 1,
μzi◦aiki (n) =
⎩ 0,
⎧⎨ aiμki (n) + μc − 1, aiμki (n) + μc − 1 ≥ 0,
One want μzi◦aiki (n) = aiμki (n). Since aiμki (n) ≥ 0, for all n,
then a sufficient condition for this to happen is to set
μc = 1. From Proposition 1, one can have μc = 1 +
min{0, min [μzi ( ) − aiμki ( )]}. That μc = 1 means that
min{0, min [μzi ( ) − aiμki ( )]} = 0, which implies that
min μzi ( ) − aiμki ( ) ≥ 0.
(A.5)
The optimum (maximum) value of ai (denoted as ai∗)
satisfying (A.5) is reached when the equality holds.
(1) Existence of ai∗.
By contradiction, assume for all ai ∈ [0, 1] such that
min [μzi ( ) − aiμki ( )] =/ 0. There are two cases:
(a) min [μzi ( ) − aiμki ( )] > 0. ⇒ μzi ( ) >
aiμki ( ) for all , which is not true, in particular
when ai = 1 and μki ( ) = 1 for some .
(b) when min [μzi ( ) − aiμki ( )] < 0 is similarly
shown.
(2) Uniqueness of ai∗.
min [μzi ( ) − aiμki ( )] = 0 ⇒
If μki ( ) = 0, for some ∈ [0, N − 1], then (A.6) is
always satisfied, so one needs to consider only when
μki ( ) > 0, for which ai ≤ μzi ( )/μki ( ), Therefore,
μzi ( ) ≥ aiμki ( ),
∀ .
min
ai∗ = 0≤ ≤N−1
μki ( ) =/ 0
μzi ( )
.
(A.6)
(A.7)
Proof of Corollary 1. (i) It follows from Proposition 2. If ai >
ai∗, then
μzi ( ) − aiμki ( ) < μzi ( ) − ai∗μki ( ),
(A.8)
taking minimum over , we get J(ai) < J(ai∗) = 0.
(ii) If ai < ai∗ can be similarly shown.
(iii) By Proposition 2 and μk0 (n) = 1.
(iv) By the antiextensive property of opening and
optimum condition μzi◦ai∗ki (n) = ai∗μki (n) ≤ μzi (n).
By Proposition 1μzi◦ai∗ki (n) ≥ 0 ⇒ 0 ≤ ai∗μki (n) ≤
μzi (n). Thus, 0 ≤ ai∗maxn[μki (n)] ≤ maxn[μzi (n)] ≤ 1, and
0 ≤ ai∗ ≤ maxn[μzi (n)] ≤ 1.
Proof of Corollary 2. By Proposition 2 and μk0 (n) = 1. a0∗ =
minn[μz0c (n)]. Thus, a0∗c = maxn[μz0 (n)].
Proof of Proposition 3. (i) According to the antiextensive
property of opening μzi◦aiki (n) ≤ μzi (n) then μzi+1 (n) =
μzi (n) − μzi◦aiki (n) ≥ 0, so μzi+1 (n) ≤ μzi (n).
(ii) By (17) and Proposition 2, we get
min μzi+1 (n) = min μzi (n) − μzi◦aiki (n) = 0.
n n
(A.9)
Proof of Proposition 4. This can be proven easily using (22)
and (23).
Proof of Proposition 5. By definition,
μki,j (m, n) = μki (m)μkj (n).
We have
μk0,j (m, n) = μk0 (m)μkj (n)
≥ μki (m)μkj (n),
∀m, n, i, j.
Notice that μk0 ( )≥ μki ( ), for all i, .
(ii) and (iii) are similarly shown.
Proof of Proposition 6. By optimum condition,
μz∅(s,t)a∗∅(s,t) μk∅(s,t) (m, n) = a∗∅(s,t) μk∅(s,t) (m, n)
aφ∗(s,t)
=
=
≤
≤
min
p,q
μkφ(s,t) (p,q) =/ 0
min
p,q
μkφ(s,t) (p,q) =/ 0
min
p,q
μkφ(s,t) (p,q) =/ 0
min
p,q
μkφ(i,j) (p,q) =/ 0
= aφ∗(i,j) − aφ∗(i,j)
Therefore, aφ∗(s,t)=0.
μzφ(s,t) p, q
μkφ(s,t) p, q
≥ 0
μzφ(i,j) p, q − aφ∗(i,j)μkφ(i,j) p, q − ∈ p, q
μkφ(s,t) p, q
μzφ(i,j) p, q − aφ∗(i,j)μkφ(i,j) p, q
μkφ(s,t) p, q
μzφ(i,j) p, q − aφ∗(i,j)μkφ(i,j) p, q
μkφ(i,j) p, q
(A.10)
(A.11)
(A.12)
[3] S. G. Mallat , “ Theory for multiresolution signal decomposition: the wavelet representation , ” IEEE Transactions on Pattern Analysis and Machine Intelligence , vol. 11 , no. 7 , pp. 674 - 693 , 1989 .
[4] S. G. Mallat , “ Multifrequency channel decompositions of images and wavelet models , ” IEEE Transactions on Acoustics, Speech, and Signal Processing , vol. 37 , no. 12 , pp. 2091 - 2110 , 1989 .
[5] J. Martens , “ The Hermite transform-theory,” IEEE Transactions on Acoustics, Speech, and Signal Processing , vol. 38 , no. 9 , pp. 1595 - 1606 , 1990 .
[6] J. Martens , “ The Hermite transform-applications,” IEEE Transactions on Acoustics, Speech, and Signal Processing , vol. 38 , no. 9 , pp. 1607 - 1618 , 1990 .
[7] D. Shinha and E. R. Dougherty , “Fuzzy mathematical morphology, ” Journel of Visual Communication and Image Processing , pp. 286 - 302 , 1992 .
[8] I. Pitas and A. N. Venetsanopoulos , “Morphological shape decomposition, ” IEEE Transactions on Pattern Analysis and Machine Intelligence , vol. 12 , no. 1 , pp. 38 - 45 , 1990 .
[9] I. Pitas , “Morphological signal decomposition,” in Proceedings of the International Conference on Acoustics, Speech, and Signal Processing (ASSP '90) , pp. 2169 - 2172 , Albuquerque, NM , USA, April 1990 .
[10] I. Pitas and A. N. Venetsanopoulos , “Morphological shape representation,” in Proceedings of the International Conference on Acoustics, Speech, and Signal Processing (ICASSP '91) , pp. 2381 - 2384 , Toronto, Canada, May 1991 .
[11] J. Song and E. J. Delp , “ The analysis of morphological filters with multiple structuring elements,” Computer Vision , Graphics and Image Processing , vol. 50 , no. 3 , pp. 308 - 328 , 1990 .
[12] C. Huang and L. F. Chaparro , “ Signal representation using fuzzy morphology,” in Proceedings of the 3rd International Symposium on Uncertainty Modeling and Analysis and Annual Conference of the North American Fuzzy Information Processing Society , (ISUMA-NAFIPS '95) , pp. 607 - 612 , IEEE, September 1995 .
[13] G. Matheron , Random Sets and Integral Geometry , Wiely, New York, NY, USA, 1973 .
[14] L. A. Zadeh , “Fuzzy sets, ” Information and Control , vol. 8 , no. 3 , pp. 338 - 353 , 1965 .
[15] D. Dubois and H. Prade , Fuzzy Sets and Systems Theory and Applications , Academic Press, New York, NY, USA, 1980 .
[16] H. Bandemer and W. Nather , Fuzzy Data Analysis , Kluwer Academic Publishers, Dordrecht, The Netherlands, 1992 .
[17] L. A. Zadeh , K. S. Fu , K. Tanaka , and M. Shimura , Fuzzy Sets and Their Applications to Cognitive and Decision Processes , Academic Press, New York, NY, USA, 1975 .
[18] J. Serra , Image Analysis and Mathematical Morphology , Academic Press, New York, NY, USA, 1982 .
[19] A. N. Evans and X. U. Liu , “ A morphological gradient approach to color edge detection , ” IEEE Transactions on Image Processing , vol. 15 , no. 6 , pp. 1454 - 1463 , 2006 .
[20] H. Cha and L. F. Chaparro , “ A morphological polynomial transform ,” in Proceedings of the International Conference on Acoustics, Speech, and Signal Processing (ICASSP '93) , vol. 5 , pp. 173 - 176 , April 1993 .
[21] C. P. Neuman and D. I. Schonbach , “ Discrete (legendre) orthogonal polynomials-a survey,” International Journal for Numerical Methods in Engineering , vol. 8 , pp. 743 - 770 , 1974 .
[22] L. F. Chaparro and M. Boudaoud , “ Two-dimensional linear prediction covariance method and its recursive solution , ” IEEE Transactions on Systems, Man and Cybernetics , vol. 17 , no. 4 , pp. 617 - 621 , 1987 .
[23] R. Gonzales and R. E. Woods , Digital Image Processing, Prentice-Hall, Englewood Cliffs, NJ, USA, 3rd edition, 2008 .
[24] P. Maragos and F. Sun , “ Measuring the fractal dimension of signals: morphological covers and iterative optimization , ” IEEE Transactions on Signal Processing , vol. 41 , no. 1 , pp. 108 - 121 , 1993 .
[25] N. Sarkar and B. B. Chauduri , “ Efficient differential boxcounting approach to compute fractal dimension of image,” IEEE Transactions on Systems, Man and Cybernetics , vol. 24 , no. 1 , pp. 115 - 120 , 1994 .
[26] S.-H. Jung , S. K. Mitra , and D. Mukherjee , “ Subband DCT: definition, analysis, and applications , ” IEEE Transactions on Circuits and Systems for Video Technology , vol. 6 , no. 3 , pp. 273 - 286 , 1996 .
[27] B. Kosko , Neural Network and Fuzzy System , Prentice-Hall, Englewood Cliffs, NJ, USA, 1991 .
[28] B. B. Mandelbrot , The Fractal Geometry of Nature, W. H. Freeman , San Francisco, Calif, USA, 1982 .
[29] P. Brodatz , Textures, Dover, New York, NY, USA, 1966 .