# The 2D Spectral Intrinsic Decomposition Method Applied to Image Analysis

Authors: Samba Sidibe; Oumar Niang; Abdoulaye Thioune; Abdoul-Dalibou Abdou; Ndeye Fatou Ngom
Journal: Journal of Electrical and Computer Engineering
 (2017)

---

## Abstract

We propose a new method for autoadaptive image decomposition and recomposition based on the two-dimensional version of the Spectral Intrinsic Decomposition (SID). We introduce a faster diffusivity function for the computation of the mean envelope operator which provides the components of the SID algorithm for any signal. The 2D version of SID algorithm is implemented and applied to some very known images test. We extracted relevant components and obtained promising results in images analysis applications.


---

## Body

## 1. Introduction

The need of components extraction and reconstruction in signal and image processing in time frequency analysis is very strong for many fields of application. Notorious methods that have been proposed include Fourier technics, wavelet decomposition, and Empirical Mode Decomposition. While Fourier transform is localized in frequency, wavelets are localized in both time and frequency; EMD is autoadaptive. EMD decomposes a signal in AM-FM components called Intrinsic Mode Functions (IMF) and a residue. This nonlinear and nonstationary decomposition works on 1D signals [1] and 2D signals such as images [2, 3]. The EMD algorithm is based on a procedure called sifting process which iteratively uses the upper and lower envelopes to extract IMFs. To create a mathematical model to compute envelopes directly, an envelope operator has been proposed in [4] and from this operator a new decomposition method called Spectral Intrinsic Decomposition, SID, was proposed in [5].The SID method allows decomposing any signal into a superposition of Spectral Proper Mode Functions (SPMFs) [5]. This method has been presented in a 1D version and depends on an operator interpolating the characteristics points of the signal to be decomposed. In this paper, the two-dimensional version of the Spectral Intrinsic Decomposition for images analysis is introduced. An algorithm for a faster spectral decomposition is proposed and illustrated with some images. We first recall the SID principle in one dimension and propose a faster method to determine the signal characteristics points in Section 2. Than an algorithm of the two-dimensional SID is presented in Section 3. Applications on grayscale images are depicted in Section 3.1.


## 2. The Spectral Intrinsic Decomposition Method

The Spectral Intrinsic Decomposition Method decomposes any signal into a combination of eigenvectors of a Partial Differential Equation (PDE) interpolation operator as presented in [5, 6].


### 2.1. The PDEs System Interpolator

For a given signals0, the upper (s+) and lower (s-) envelope are the asymptotic solution of the following PDEs system:(1)∂s±x,t∂t=-g±x,tα∂2s±x,t∂2x+1-α∂4s±x,t∂4x,s±x,0=s0.α is the tension parameter which ranges from 0 to 1. g+ and g- are diffusivity functions for upper and lower envelope which are equal to zero at characteristic points of s0 and range from zero to one. g± based on Maximum Curvature Points (MCP) can be computed as follows:(2)g±x=19sgn∂3s0x∂3x±sgn∂2s0x∂2x+12,where sgn denotes the sign function.Equation (1) is resolved numerically in its discrete implicit unconditionally stable scheme as follows:(3)Sk+1=Sk+ΔtASk+1,S0=S0,where Δt is the time step, Sk+1=s(x,kΔt) is signal value at step k+1, and A is a matrix formed with finite difference approximation coefficients of second- and fourth-order differential operators (resp., D2 and D4), as A=GαD2-1-αD4, with G being the diagonal matrix constructed with discrete version of stopping function values g(x), exactly, G(i,i)=g(i).So the explicit form leads to the following numerical resolution:(4)Sk+1=I-ΔtA-1Sk,S0=s0,k≥0with I being the identity matrix. Finally (1) can be decomposed into a linear system from implicit numerical scheme (4) by(5)Sk+1=L-1Sk,S0=s0,k≥0.L is given by L=I-ΔtA.The operator matrix,L, has real-valued eigenvalues that are always greater than or equal to 1. Then, eigenvalues, λn, of L-1 are always smaller than or equal to 1  (0<λn≤1); see [4].


### 2.2. On the Asymptotic Solution

Iterative scheme (5) can be rewritten in terms of initial solution s0 as(6)Sk=L-1ks0,k≥1.After convergence (see [7]), the asymptotic solution, S∞, is given by(7)S∞=L-1∞s0.LetV be a matrix of L-1’s sequence of eigenvectors (Vn) and D a diagonal matrix having L-1’s sequence of eigenvalues (λn), at the diagonal. So we have the following decomposition L-1=VDV-1. It is easy to see that(8)L-1k=VDV-1k=VDkV-1.So, the asymptotic solution in (7) is given by(9)S∞=VD∞V-1S0.The asymptotic eigenvalue matrixD∞ is a diagonal matrix with eigenvalues λn∞=1 only at loci where matrix G is zeroed, and λn∞=0, where g[n]>0. Finally, the asymptotic solution of the PDE interpolator system is a linear combination of fixed vector point of upper and lower envelope operators.


### 2.3. A Faster Stopping Function for Discrete Signal

For image processing we will consider region boundaries as characteristic points. The characteristic points of the upper envelope will be the local maximums and the limits of the regions where the value of the gray level of the pixel is equal to or greater than the gray level of all the pixels in their neighborhood represented, for example, by a rectangular window.We define the diffusion functiongM- for lower envelope to be equal to 1 everywhere except in characteristic points of the lower envelope where it will be equal to 0.Similarly the characteristic points of the lower envelope are local minimums and region boundaries where the pixel value is equal to or less than the gray level of all pixels in their neighbors. We define the diffusion functiongM+ for upper envelope to be equal to 1 everywhere except in characteristics points of the upper envelope where it will be equal to 0.The diffusivity function called stopping functiongM± is calculated by using morphological dilation and erosion operations [8].Letb=-1,0,1 be a structured element; the grayscale dilation of s by b at x is given by(10)s⊕bx=maxa∈bsx+a.The grayscale erosion of s by b at x is given by(11)s⊖bx=mina∈bsx+a.s
⊕
b is equal to s at local maxima and when s is locally constant; s⊖b is equal to s at local minima and when s is locally constant.g
M
+ is zeroed at points which are invariant to morphological dilation and variant to morphological erosion; gM+ is equal to 1 for any other point.Similarly,gM- is zeroed at points which are invariant to morphological erosion and variant to morphological dilation; gM- is equal to 1 for any other point.Let(12)A=x∣s⊕bx=sx,B=x∣s⊖bx=sx.We have(13)gM+x=0,ifx∈A∖B,1,otherwise,gM-x=0,ifx∈B∖A,1,otherwiseg
M functions are faster than g to compute and give satisfying results for the computation of envelope operators for real images (Table 1).Table 1
Image
Width
Height
Time forg±
Time forgM±
Figure4(a)
300
300
0,4684466
0.320695
Figure3(a)
400
266
0,826974833
0,2877035
Figure2(a)
400
272
0,807968
0,284108Figure1 shows an example of calculating the envelope using gM, gM+ for upper envelope and gM- for lower envelope.Figure 1
The effectiveness of envelope computation usinggM instead of g.


### 2.4. The Spectral Intrinsic Decomposition

In the following,E denotes either the upper or lower envelope operator. The upper and lower envelope of the signal are calculated with the eigenvectors associated with eigenvalue λ=1. Hence, S∞ in formula (9) is a linear combination of all 1-eigenvectors of the envelope operator E weighted by the signal amplitude. Instead of focusing only on the envelope calculus, we now consider all the set of eigenvalues of the envelope operator E.The Spectral Intrinsic Decomposition procedure is defined as the calculus of all the SPMFs for a given signal.E=L-1. The same procedure can be performed with the lower envelope. The eigendecomposition of E gives [VE,LE]=eig(E), where VE=[V1,…,Vsize(s0)] and LE=[L1,…,Lsize(s0)] (with the possibility of zeros to complete the size of the vector) are, respectively, the set of eigenvectors and the set of eigenvalues of E. The coefficient reconstruction of s0 is given by C=LEVE-1s0⊥, with s0⊥ being the transpose of s0. Hence s0 is computed by the formula s0=VC (Algorithm 2). The Spectral Intrinsic Decomposition of s0 is given as follows:(14)s0=∑k=1NVkCk,N=sizes0.This decomposition is intrinsic and depends only on the position of the characteristic points ofs0 that define the diffusivity function in the interpolation operator. We notice that the SID with lower envelope works like the SID with the upper envelope and has the same reconstruction ability.


### 2.5. Advantage and Disadvantage of SID

The SID is adaptive and depends on the position of the characteristics points of the signal. It is autoadaptive and works for nonlinear and nonstationary signal. SID can decompose an IMF and can be used to separate mixing mode [9] in EMD.However, the main disadvantage of the proposed SID algorithm is the computation time when the size of signal is a large. This is due to matrix inversion in the algorithm. Thus a faster algorithm is proposed in the next section.


## 2.1. The PDEs System Interpolator

For a given signals0, the upper (s+) and lower (s-) envelope are the asymptotic solution of the following PDEs system:(1)∂s±x,t∂t=-g±x,tα∂2s±x,t∂2x+1-α∂4s±x,t∂4x,s±x,0=s0.α is the tension parameter which ranges from 0 to 1. g+ and g- are diffusivity functions for upper and lower envelope which are equal to zero at characteristic points of s0 and range from zero to one. g± based on Maximum Curvature Points (MCP) can be computed as follows:(2)g±x=19sgn∂3s0x∂3x±sgn∂2s0x∂2x+12,where sgn denotes the sign function.Equation (1) is resolved numerically in its discrete implicit unconditionally stable scheme as follows:(3)Sk+1=Sk+ΔtASk+1,S0=S0,where Δt is the time step, Sk+1=s(x,kΔt) is signal value at step k+1, and A is a matrix formed with finite difference approximation coefficients of second- and fourth-order differential operators (resp., D2 and D4), as A=GαD2-1-αD4, with G being the diagonal matrix constructed with discrete version of stopping function values g(x), exactly, G(i,i)=g(i).So the explicit form leads to the following numerical resolution:(4)Sk+1=I-ΔtA-1Sk,S0=s0,k≥0with I being the identity matrix. Finally (1) can be decomposed into a linear system from implicit numerical scheme (4) by(5)Sk+1=L-1Sk,S0=s0,k≥0.L is given by L=I-ΔtA.The operator matrix,L, has real-valued eigenvalues that are always greater than or equal to 1. Then, eigenvalues, λn, of L-1 are always smaller than or equal to 1  (0<λn≤1); see [4].


## 2.2. On the Asymptotic Solution

Iterative scheme (5) can be rewritten in terms of initial solution s0 as(6)Sk=L-1ks0,k≥1.After convergence (see [7]), the asymptotic solution, S∞, is given by(7)S∞=L-1∞s0.LetV be a matrix of L-1’s sequence of eigenvectors (Vn) and D a diagonal matrix having L-1’s sequence of eigenvalues (λn), at the diagonal. So we have the following decomposition L-1=VDV-1. It is easy to see that(8)L-1k=VDV-1k=VDkV-1.So, the asymptotic solution in (7) is given by(9)S∞=VD∞V-1S0.The asymptotic eigenvalue matrixD∞ is a diagonal matrix with eigenvalues λn∞=1 only at loci where matrix G is zeroed, and λn∞=0, where g[n]>0. Finally, the asymptotic solution of the PDE interpolator system is a linear combination of fixed vector point of upper and lower envelope operators.


## 2.3. A Faster Stopping Function for Discrete Signal

For image processing we will consider region boundaries as characteristic points. The characteristic points of the upper envelope will be the local maximums and the limits of the regions where the value of the gray level of the pixel is equal to or greater than the gray level of all the pixels in their neighborhood represented, for example, by a rectangular window.We define the diffusion functiongM- for lower envelope to be equal to 1 everywhere except in characteristic points of the lower envelope where it will be equal to 0.Similarly the characteristic points of the lower envelope are local minimums and region boundaries where the pixel value is equal to or less than the gray level of all pixels in their neighbors. We define the diffusion functiongM+ for upper envelope to be equal to 1 everywhere except in characteristics points of the upper envelope where it will be equal to 0.The diffusivity function called stopping functiongM± is calculated by using morphological dilation and erosion operations [8].Letb=-1,0,1 be a structured element; the grayscale dilation of s by b at x is given by(10)s⊕bx=maxa∈bsx+a.The grayscale erosion of s by b at x is given by(11)s⊖bx=mina∈bsx+a.s
⊕
b is equal to s at local maxima and when s is locally constant; s⊖b is equal to s at local minima and when s is locally constant.g
M
+ is zeroed at points which are invariant to morphological dilation and variant to morphological erosion; gM+ is equal to 1 for any other point.Similarly,gM- is zeroed at points which are invariant to morphological erosion and variant to morphological dilation; gM- is equal to 1 for any other point.Let(12)A=x∣s⊕bx=sx,B=x∣s⊖bx=sx.We have(13)gM+x=0,ifx∈A∖B,1,otherwise,gM-x=0,ifx∈B∖A,1,otherwiseg
M functions are faster than g to compute and give satisfying results for the computation of envelope operators for real images (Table 1).Table 1
Image
Width
Height
Time forg±
Time forgM±
Figure4(a)
300
300
0,4684466
0.320695
Figure3(a)
400
266
0,826974833
0,2877035
Figure2(a)
400
272
0,807968
0,284108Figure1 shows an example of calculating the envelope using gM, gM+ for upper envelope and gM- for lower envelope.Figure 1
The effectiveness of envelope computation usinggM instead of g.


## 2.4. The Spectral Intrinsic Decomposition

In the following,E denotes either the upper or lower envelope operator. The upper and lower envelope of the signal are calculated with the eigenvectors associated with eigenvalue λ=1. Hence, S∞ in formula (9) is a linear combination of all 1-eigenvectors of the envelope operator E weighted by the signal amplitude. Instead of focusing only on the envelope calculus, we now consider all the set of eigenvalues of the envelope operator E.The Spectral Intrinsic Decomposition procedure is defined as the calculus of all the SPMFs for a given signal.E=L-1. The same procedure can be performed with the lower envelope. The eigendecomposition of E gives [VE,LE]=eig(E), where VE=[V1,…,Vsize(s0)] and LE=[L1,…,Lsize(s0)] (with the possibility of zeros to complete the size of the vector) are, respectively, the set of eigenvectors and the set of eigenvalues of E. The coefficient reconstruction of s0 is given by C=LEVE-1s0⊥, with s0⊥ being the transpose of s0. Hence s0 is computed by the formula s0=VC (Algorithm 2). The Spectral Intrinsic Decomposition of s0 is given as follows:(14)s0=∑k=1NVkCk,N=sizes0.This decomposition is intrinsic and depends only on the position of the characteristic points ofs0 that define the diffusivity function in the interpolation operator. We notice that the SID with lower envelope works like the SID with the upper envelope and has the same reconstruction ability.


## 2.5. Advantage and Disadvantage of SID

The SID is adaptive and depends on the position of the characteristics points of the signal. It is autoadaptive and works for nonlinear and nonstationary signal. SID can decompose an IMF and can be used to separate mixing mode [9] in EMD.However, the main disadvantage of the proposed SID algorithm is the computation time when the size of signal is a large. This is due to matrix inversion in the algorithm. Thus a faster algorithm is proposed in the next section.


## 3. The 2D Spectral Intrinsic Decomposition Algorithm

In this section we present, in Algorithm1, the 2D SID procedure and implement it for images decomposition-recomposition.Algorithm 1:Pseudocode for decomposition and recomposition.
(1) M← transforme image to data
(2) NL ← number of lines of M
(3) NC ← number of columns of M
(4) MR← NL by NC zeroed matrix
(5) for  i←1, NL  do
(6)  si← line i of M
(7)  compute gi+ from si
(8)  compute Ei
(9)  decompose Ei, [VEi,LEi]←eig(E)
(10)   compute Ci, Ci=LEiVEi-1si′
(11)    for  j←1, NC  do
(12)      MR(i,j)←Vi(j)∗Ci(j)
(13)   end for
(14)  end for
(15) recompose image from MRAlgorithm 2:Pseudocode to extract SPMFs between two eigenvalues.
(1) [λmin,λmax]← eigenvalues interval to extract
(2) M← transforms image to data
(3) NL ← number of lines of M
(4) NC ← number of columns of M
(5) MR← NL by NC zeroed matrix
(6) for  i←1, NL  do
(7)  Si← line i of M
(8)  compute gi+ from si
(9)  compute Ei
(10)  decompose Ei, [VEi,LEi]←eig(E)
(11)   LRi←LEi
(12)   for  j←1, NC  do
(13)     if  LRi(j)≤λminorLRi(j)≥λmax  then
(14)       LRi(j)←0
(15)     end if
(16)   end for
(17)   compute CRi, CRi←LRiVEi-1Si′
(18)   for  j←1, NC  do
(19)     MR(i,j)←MR(i,j)+VEi(j)∗CRi(j)
(20)   end for
(21) end for


### 3.1. The Algorithm of Image Decomposition and Recomposition by SID

Let us consider an image as represented by a matrixM(Line (1)). Each row i of the matrix M(Line (5)) can be seen as a one-dimensional signal si. Thus, we can apply the spectral decomposition of si(Lines (6) to (10)) as described in [4].Each line can be recomposed to build a matrixMR(Lines (11) to (13)) and finally reconstruct the image from MR(Line (15)). Image decomposition and recomposition are described in Algorithm 1. We can also make decomposition along columns and catch more properties of the image to be analyzed.The elementary components of a spectral decomposition are the modulation of eigenvectors by their coefficients as can be seen in (14). It is clear that these components depend on eigenvalues of the envelope operator.The range of smallest eigenvalues catches higher frequencies contents of the reconstructed signal with the smallest modulated amplitude.So we can associate components which have similar frequency and amplitude by summing the components that have same eigenvalues; this method works well for signal in one dimension. For images, it is possible numerically to have missing eigenvalues in many lines. Hence, to avoid this drawback, elementary components which have the same eigenvalues will be associated with the same belonging to a specific range of eigenvalues.


### 3.2. Application to Image Components Extraction

In the following, 2D SID is applied to very known images test to demonstrate the ability of this pectoral intrinsic decomposition for images analysis, particularly in components extraction. In Figures4(a), 3(a), and 2(a), we have a representation of high frequency components on Figures 4(b), 3(b), and 2(b); we note that they are more sensitive to small variations of the intensity of the image than low frequencies components in Figures 4(c), 3(c), and 2(c).Figure 2
Decomposition (a), component extraction and representation of high frequency component (λ∈[0,0.1]) (b), medium frequency component (λ∈[0.5,0.7]) (c), high and medium frequency component (λ<1) (d), low frequency component (λ=1) (e), and recomposition of all components (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 3
Decomposition (a), component extraction and representation of high frequency component (λ∈[0,0.1]) (b), medium frequency component (λ∈[0.5,0.7]) (c), high and medium frequency component (λ<1) (d), low frequency component (λ=1) (e), and recomposition of all components (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 4
Decomposition (a), component extraction and representation of high frequency component (λ∈[0,0.1]) (b), medium frequency component (λ∈[0.5,0.7]) (c), high and medium frequency component (λ<1) (d), low frequency component (λ=1) (e), and recomposition of all components (f).
(a)
(b)
(c)
(d)
(e)
(f)


## 3.1. The Algorithm of Image Decomposition and Recomposition by SID

Let us consider an image as represented by a matrixM(Line (1)). Each row i of the matrix M(Line (5)) can be seen as a one-dimensional signal si. Thus, we can apply the spectral decomposition of si(Lines (6) to (10)) as described in [4].Each line can be recomposed to build a matrixMR(Lines (11) to (13)) and finally reconstruct the image from MR(Line (15)). Image decomposition and recomposition are described in Algorithm 1. We can also make decomposition along columns and catch more properties of the image to be analyzed.The elementary components of a spectral decomposition are the modulation of eigenvectors by their coefficients as can be seen in (14). It is clear that these components depend on eigenvalues of the envelope operator.The range of smallest eigenvalues catches higher frequencies contents of the reconstructed signal with the smallest modulated amplitude.So we can associate components which have similar frequency and amplitude by summing the components that have same eigenvalues; this method works well for signal in one dimension. For images, it is possible numerically to have missing eigenvalues in many lines. Hence, to avoid this drawback, elementary components which have the same eigenvalues will be associated with the same belonging to a specific range of eigenvalues.


## 3.2. Application to Image Components Extraction

In the following, 2D SID is applied to very known images test to demonstrate the ability of this pectoral intrinsic decomposition for images analysis, particularly in components extraction. In Figures4(a), 3(a), and 2(a), we have a representation of high frequency components on Figures 4(b), 3(b), and 2(b); we note that they are more sensitive to small variations of the intensity of the image than low frequencies components in Figures 4(c), 3(c), and 2(c).Figure 2
Decomposition (a), component extraction and representation of high frequency component (λ∈[0,0.1]) (b), medium frequency component (λ∈[0.5,0.7]) (c), high and medium frequency component (λ<1) (d), low frequency component (λ=1) (e), and recomposition of all components (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 3
Decomposition (a), component extraction and representation of high frequency component (λ∈[0,0.1]) (b), medium frequency component (λ∈[0.5,0.7]) (c), high and medium frequency component (λ<1) (d), low frequency component (λ=1) (e), and recomposition of all components (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 4
Decomposition (a), component extraction and representation of high frequency component (λ∈[0,0.1]) (b), medium frequency component (λ∈[0.5,0.7]) (c), high and medium frequency component (λ<1) (d), low frequency component (λ=1) (e), and recomposition of all components (f).
(a)
(b)
(c)
(d)
(e)
(f)


## 4. Conclusion

In this paper, we have presented a new method for autoadaptive image representation called two-dimensional version of the Spectral Intrinsic Decomposition. A new faster diffusivity function in the computation of the mean envelope operator is also provided. In future works, we will investigate how to use the Spectral Proper Mode Functions (SPMFs) to do signals classification or to treat other aspects of image processing, like edge detection, segmentation, and so on.



---

*Source: 7948571-2017-12-20.xml*