# Orientation-independent empirical mode decomposition for images based on unconstrained optimization

Marcelo Colominas, Anne Humeau-Heurtier, Gastòn Schlotthauer

<figure><img src="image_1.png" /></figure>

## To cite this version:

Marcelo Colominas, Anne Humeau-Heurtier, Gastòn Schlotthauer. Orientation-independent empirical mode decomposition for images based on unconstrained optimization. IEEE Transactions on Image Processing, 2016, 25 (5), pp.2288-2297. 10.1109/TIP.2016.2541959. hal-01388824

HAL Id: hal-01388824

https://hal.science/hal-01388824

Submitted on 22 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Orientation-Independent Empirical Mode Decomposition for Images Based on Unconstrained Optimization

Marcelo A. Colominas, Anne Humeau-Heurtier, and Gastón Schlotthauer, Member, IEEE

**Abstract**—This paper introduces a 2D extension of the empirical mode decomposition (EMD), through a novel approach based on unconstrained optimization. EMD is a fully data-driven method that locally separates, in a completely data-driven and unsupervised manner, signals into fast and slow oscillations. The present proposal implements the method in a very simple and fast way, and it is compared with the state-of-the-art methods evidencing the advantages of being computationally efficient, orientation-independent, and leads to better performances for the decomposition of amplitude modulated-frequency modulated (AM-FM) images. The resulting genuine 2D method is successfully tested on artificial AM-FM images and its capabilities are illustrated on a biomedical example. The proposed framework leaves room for an nD extension (n > 2).

**Index Terms**—Empirical mode decomposition, unconstrained optimization, data-driven, non-stationary image.

## I. INTRODUCTION

The separation of non-stationary images into locally, but not globally, frequency disjoint components is often necessary for applications such as texture classification [1]–[3], face recognition [4]–[6], and digital watermarking [7]. The main goal is to extract local quantities, such as amplitudes, phases, and frequencies. In all cases, the images that are processed can be seen as a sum of a usually small number of amplitude modulated-frequency modulated (AM-FM) components plus a final trend:

$$X(x, y) = \sum_{k=1}^{K} \alpha_k(x, y) \cos \phi_k(x, y) + A_K(x, y), \quad (1)$$

M. A. Colominas is with the Laboratorio de Señales y Dinámicas no Lineales, Facultad de Ingeniería, Universidad Nacional de Entre Ríos: Ruta Prov. 11 Km. 10, CP: 3100. Oro Verde, Entre Ríos, Argentina, and also with the Consejo Nacional de Investigaciones Científicas y Técnicas: Godoy Cruz 2290 (C1425FQB) CABA, Argentina (e-mail: macolominas@bioingenieria.edu.ar).

A. Humeau-Heurtier is with the Laboratoire Angevin de Recherche en Ingénierie des Systèmes, Université d'Angers, Angers 49000, France (e-mail: anne.humeau@univ-angers.fr).

G. Schlotthauer is with the Laboratorio de Señales y Dinámicas no Lineales, Facultad de Ingeniería, Universidad Nacional de Entre Ríos, Argentina, also with the Consejo Nacional de Investigaciones Científicas y Técnicas, Argentina, also with the Centro de Investigaciones y Transferencia de Entre Ríos (CITER) Av. Monseñor Tavella 1450 Sede UNER Alimentos CP 3200. Concordia, Entre Ríos, Argentina (e-mail: gschlotthauer@conicet.gov.ar).

where every local detail $D_k(x, y) = \alpha_k(x, y) \cos \phi_k(x, y)$ is riding on a local mean $A_k = \sum_{i=k+1}^{K} D_i + A_K$.

Empirical Mode Decomposition (EMD) [8] provides an approach to the form of Eq. (1) for 1D-signals. It locally separates, in a completely data-driven and unsupervised manner, signals into fast and slow oscillations. A post-processing via Hilbert Transform (HT) provides an estimation of the local frequencies and amplitudes. The global procedure of EMD plus HT is the so-called Hilbert Huang Transform (HHT). In the EMD jargon the local details are called Intrinsic Mode Functions (IMFs) or simply modes: functions symmetrically oscillating around zero, albeit with modulation both in amplitude and frequency. The final trend, usually a monotonic function, is called the residue.

The adaptivity and flexibility of EMD encouraged several researchers to adapt this method for 2D-data. Nunes et al. [9] developed a bidimensional EMD where they used morphological operators to detect regional maxima and radial basis function for surface interpolation. Damerval et al. [10] introduced a fast bidimensional EMD where a piecewise cubic interpolation of a Delaunay triangulation is used to find the envelopes. Linderhed [11] introduced an image EMD where the envelopes are found through thin plate splines. Wu et al. [12] proposed a multi-dimensional approach where they decomposed rows and columns independently and combined the results based on a minimal-scale combination strategy. However, all of these algorithms rely on an iterative procedure called sifting process which lacks of convergence guarantees. To avoid this iterative process, several optimization-based strategies have recently been developed for 1D-signals [13], [14]. In [15], the authors introduced an unconstrained optimization approach to EMD which is simpler and faster than the other proposals. The main advantages of the latter method compared to others are the following: the explicit computation of envelopes to find the local mean is not needed; there is no use of explicit spline interpolations; the proposed method provides an analytical solution. Moreover, the computational cost is similar to the one of EMD and only one parameter has to be set.

Nevertheless, such an approach for bidimensional data has not been proposed yet. We therefore herein introduce
an orientation-independent EMD for images that is based on an unconstrained optimization problem. We show that, compared to other 2D-EMD-based methods, our framework has the advantages of being computationally efficient, orientation-independent, and leads to better performances for the decomposition of AM-FM images.

The paper is organized as follows. In Section II, we recall the unconstrained optimization approach to EMD for 1D-signals. Section III introduces the framework that we propose for the decomposition of bidimensional data. Experimental results for synthetic and experimental data are presented and discussed in Section IV. We end with a Conclusion.

## II. UNCONSTRAINED OPTIMIZATION APPROACH TO EMD

The separation of a signal in detail (mode) plus approximation (trend or local mean) can be done in several ways. Traditional EMD [8] uses the so-called sifting process, which heavily relies on cubic spline interpolations and does not have a guaranteed convergence. On the other hand, the optimization approach to EMD introduced in [15] solves the following convex optimization problem to find the local mean $a \in \mathbb{R}^N$ of a given signal $x \in \mathbb{R}^N$:

$$ (P) \min_{a} \|P_{x}(x-a)\|_{2}^{2} + \lambda \|La\|_{2}^{2}, \quad (2) $$

where matrix $P_x \in \mathbb{R}^{N \times N}$ is an operator, depending only on the locations of local extrema of $x$, which models the penalization imposed on the mode $d = x - a$ at each extrema, matrix $L \in \mathbb{R}^{N \times N}$ is a second-order difference matrix, $a$ and $x$ are considered column vectors, and $|| \cdot ||_2$ stands for the $\ell_2$-norm. This problem has a unique solution:

$$ a^* = (P_x^T P_x + \lambda L^T L)^{-1} P_x^T P_x x. \quad (3) $$

The matrix $P_x$ is a key element in the problem statement. It has as many non-zero rows as local extrema of the signal $x$. Let $t_l$, $l = 1, \dots, L$, with $t_l < \dots < t_l < t_{l+1} < \dots < t_L$, be the locations of the local extrema of $x$, which act as estimates of the local extrema of the mode $d$. The matrix evaluates the following function in every local extrema:

$$ \langle p_l, d \rangle = d(t_l) + \frac{d(t_{l+1})(t_l - t_{l-1})}{t_{l+1} - t_{l-1}} + \frac{d(t_{l-1})(t_{l+1} - t_l)}{t_{l+1} - t_{l-1}}, \quad (4) $$

where the vector $p_l \in \mathbb{R}^N$ is the $t_l$-th row of matrix $P_x$, and it has non-zeros elements only on the $t_{l-1}$, $t_l$, $t_{l+1}$-th positions:

$$ p_l = \left[ \dots \frac{t_{l+1} - t_l}{t_{l+1} - t_{l-1}} \dots 1 \dots \frac{t_l - t_{l-1}}{t_{l+1} - t_{l-1}} \dots \right]. \quad (5) $$

The goal of Eq. (4) is to compare the extremum $d(t_l)$ with its mirror point on the would-be other envelope which is locally defined thanks to $d(t_{l-1})$ and $d(t_{l+1})$ [16], thus favoring the symmetry of the mode. It can be clearly seen that if $t_l$ is a local minimum (resp. maximum), then both $t_{l-1}$ and $t_{l+1}$ are local maximum (resp. minimum). In this approach the IMF conditions (envelope symmetry) are evaluated only on the signal local extrema and not throughout the whole time span. Equations (3) and (5) allow to perform the decomposition in a deflationary scheme, which is outlined in Algorithm 1 [15].

## Algorithm 1 UOA-EMD for 1D Data (Signals)

1. Assign $a_0 = x \in \mathbb{R}^N$ (signal).

2. Define the second-order difference matrix

$$ L \in \mathbb{R}^{N \times N} = \begin{bmatrix} -2 & 1 & & & \\ 1 & -2 & 1 & & \\ & \ddots & \ddots & \ddots & \\ & & 1 & -2 & 1 \\ & & & 1 & -2 \end{bmatrix}. $$

3. **for** $k = 1, 2, \dots, \text{do}$

4. Construct the matrix $P_{a_{k-1}}$, according to (5), by using the local extrema from $a_{k-1}$.

5. Calculate the local mean $a_k$ for a $\lambda_k > 0$

$$ a_k = (P_{a_{k-1}}^T P_{a_{k-1}} + \lambda_k L^T L)^{-1} (P_{a_{k-1}}^T P_{a_{k-1}}) a_{k-1}. $$

6. Calculate the mode $d_k = a_{k-1} - a_k$.

7. **end for**

## III. PROPOSED FRAMEWORK

### A. Bidimensional Approach

Let us consider now a function of two real variables $X = X(x, y)$, with $X \in \mathbb{R}^{N_1 \times N_2}$. The countable set of isolated local extrema is

$$ \mathcal{L} = \left\{ z_l = (x_l, y_l) / X(z_l) > X(z) \lor X(z_l) < X(z), \quad 0 < |z - z_l| < \delta_l \right\}. \quad (6) $$

this definition of local extremum may vary depending on whether a strict inequality is required or not. Our desire is to find an approximation $A \in \mathbb{R}^{N_1 \times N_2}$ such that the penalization on every local extremum of the desired mode $D = X - A$ is minimized. To do that, we must apply Eq. (4) in a 2D manner. If $z_l$ is a local maximum (resp. minimum), let us define $z_i$, $i = 1, 2, 3$, as the three closest non-collinear local minima (resp. maxima) [17]. Then, the penalization on $z_l$ is

$$ \langle P_l, D \rangle = D(z_l) + c_1 D(z_1) + c_2 D(z_2) + c_3 D(z_3), \quad (7) $$

where $\langle \cdot, \cdot \rangle$ stands here for the inner product between matrices, $P_l \in \mathbb{R}^{N_1 \times N_2}$ is a matrix with $P_l(z_i) = c_i$, $i = 1, 2, 3$, and the coefficients $c_i$ are those necessary to compare the function $D(z_l)$ with the plane intersecting $D(z_i)$:

$$ c_2 = \frac{(x_l - x_1)(y_3 - y_1) - (x_3 - x_1)(y_l - y_1)}{(x_2 - x_1)(y_3 - y_1) - (x_3 - x_1)(y_2 - y_1)} \quad (8a) $$

$$ c_3 = \frac{(x_2 - x_1)(y_l - y_1) - (x_l - x_1)(y_2 - y_1)}{(x_2 - x_1)(y_3 - y_1) - (x_3 - x_1)(y_2 - y_1)} \quad (8b) $$

$$ c_1 = 1 - c_2 - c_3. \quad (8c) $$

To favour smooth solutions, we must use a roughening matrix. We choose a discrete approximation to the Laplacian operator $L_{i,j} \in \mathbb{R}^{N_1 \times N_2}$ which consists in a matrix of the same size as $X$ with all zeros except for the kernel

$$ K = \begin{bmatrix} 0.5 & 1 & 0.5 \\ 1 & -6 & 1 \\ 0.5 & 1 & 0.5 \end{bmatrix} \quad (9) $$
centered at the point i, j. When i, j is at the border of the matrix, we simply cropped the kernel, and modify the central value so that the sum of the kernel entries is zero.

With Eq. (7) and the Laplacian operator, we can formulate the following optimization problem:

$$ (\text{P1-2D}) \quad \min_A \sum_{l/z_l \in \mathcal{L}} \langle P_l, (X-A) \rangle^2 + \lambda \sum_{i,j} \langle L_{i,j}, A \rangle^2, \quad (10) $$

in which the first term models the penalization in every local extrema and the second term favours smooth solutions. However, problem P1-2D can be difficult to solve. We must modify it. The inner product between matrices can be computed as $\langle A, B \rangle = \text{Tr}(A^T B) = (\text{vec}(A))^T (\text{vec}(B))$, where $\text{vec}(\cdot)$ is the classical vectorization operation, in which the columns of a matrix are vertically concatenated to form a vector; and $\text{Tr}(\cdot)$ stands for the trace of a matrix. Let us define the vectors $\hat{a} \in \mathbb{R}^{N_1 N_2} = \text{vec}(A)$ and $\hat{x} \in \mathbb{R}^{N_1 N_2} = \text{vec}(X)$, the matrix $\hat{P}_X \in \mathbb{R}^{N_1 N_2 \times N_1 N_2}$ which contains $(\text{vec}(P_l))^T$ as rows and $\hat{L} \in \mathbb{R}^{N_1 N_2 \times N_1 N_2}$ which contains $(\text{vec}(L_{i,j}))^T$ as rows. The following optimization problem

$$ (\text{P-2D}) \quad \min_{\hat{a}} \|\hat{P}_X(\hat{x} - \hat{a})\|_2^2 + \lambda \|\hat{L}\hat{a}\|_2^2, \quad (11) $$

is equivalent to (P1-2D) and has a unique minimizer at

$$ \hat{a}^* = (\hat{P}_X^T \hat{P}_X + \lambda \hat{L}^T \hat{L})^{-1} \hat{P}_X^T \hat{P}_X \hat{x}, \quad (12) $$

from which the desired local mean $A$ can be constructed.

Problem (P-2D) arises as a natural 2D extension of (P), and paves the road for nD extensions, simply by defining the penalization tensor for every local extrema, and a roughening tensor on nD. It also constitutes a true genuine-2D approach, as a contrast to pseudo-2D approaches, because it does not separately constraints rows, columns and diagonals. The result is an orientation-independent image processing tool.

## B. Implementation Issues

The sizes of the matrices involved in problem (P-2D) might look overwhelming at first sight, and therefore the inversion in Eq. (12) might appear as very expensive in terms of computational cost. Sparsity must be taken into account to significantly reduce the cost. Matrices $P_l$ have at most 4 non-zero elements, and therefore $\hat{P}_X$ has a proportion of $4|\mathcal{L}|/N_1^2 N_2^2$ non-zero elements, with $|\mathcal{L}|$ being the cardinality of the set $\mathcal{L}$, and $|\mathcal{L}| \ll N_1 N_2$. The same applies to $L_{i,j}$, which has at most 9 non-zero elements, and therefore $\hat{L}$ has a proportion of $9/N_1 N_2$ of non-zero elements. Efficient routines to invert sparse matrices make our proposal a very low-cost method when compared with other bidimensional EMD implementations. From Eqs. (7)-(9) and (12), we construct a deflationary scheme to decompose an image $X$, as shown in Algorithm 2. We will call our proposal UOA-EMD, for unconstrained optimization approach.

# IV. EXPERIMENTS AND RESULTS

## A. 2D Dirac Function

Genuine 2D approaches must be orientation-independent. In order to test this property on UOA-EMD, we present as

## Algorithm 2 UOA-EMD for 2D Data (Images)

1: Given $X \in \mathbb{R}^{N_1 \times N_2}$ (image), assign $A_0 = X$ and $\hat{a}_0 = \text{vec}(X) \in \mathbb{R}^{N_1 N_2}$.

2: Define the sparse roughening matrix $\hat{L} \in \mathbb{R}^{N_1 N_2 \times N_1 N_2}$

3: **for** $k = 1, 2, \dots, \text{do}$

4: Determine the set $\mathcal{L}_k$ of local extrema of $A_k$.

5: Define the sparse matrix $\hat{P}_{A_0} \in \mathbb{R}^{N_1 N_2 \times N_1 N_2}$.

6: **for** $l = 1, 2, \dots, |\mathcal{L}_k|$ **do**

7: Define the sparse matrix $P_l \in \mathbb{R}^{N_1 \times N_2}$ and the coefficients $P(z_l) = 1$ and $P(z_i) = c_i$ as in (8).

8: $\hat{P}_{A_0}(y_l(N_1 - 1) + x_l, \cdot) \leftarrow \text{vec}(P_l)^T$.

9: **end for**

10: Calculate the local mean $\hat{a}_k$ for a $\lambda_k > 0$.

$$ \hat{a}_k = (\hat{P}_{A_k}^T \hat{P}_{A_k} + \lambda_k \hat{L}^T \hat{L})^{-1} \hat{P}_{A_k}^T \hat{P}_{A_k} \hat{a}_{k-1}. $$

11: Construct $A_k$ from $\hat{a}_k$ and calculate the mode $D_k = A_{k-1} - A_k$.

12: **end for**

Notations: we used the notation $M(n; \cdot)$ to denote the n-th row of matrix $M$.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Dirac function decomposition. Average results of one hundred decompositions of noisy 2D Dirac functions.</figcaption></figure>

a first experiment the decomposition of one hundred noisy copies of a 2D Dirac function and the averaging of the final results, in an Ensemble EMD fashion [18]. Given $\Delta \in \mathbb{R}^{64 \times 64}$, with $\Delta(32, 32) = 1$, we decomposed one hundred realizations of $\Delta^{(i)} = \Delta + \beta W^{(i)}$, with $W^{(i)} \sim \mathcal{N}(0, 1)$ and $\beta = 0.001$, into two modes and one residue.

The results are shown in Fig. 1. The latter figure confirms the orientation-independence of UOA-EMD.
<figure><img src="image_4.png" /><figcaption>Fig. 2. Results for Experiment 1. Best results for the four analyzed methods when decomposing an image mixing components described in Eqs. (13).</figcaption></figure>

## B. AM-FM Images

In this subsection we present four different experiments decomposing AM-FM images, measuring both performance in component retrieval and computational time as a function of the only parameter $\lambda$. Moreover, we compare our results with those of the state of the art methods: Image EMD (IEMD) [11], Prox-EMD [17] in its pseudo-2D implementation, and TV-G texture-geometry decomposition [19].

IEMD uses sifting process to extract every mode, iterating until $|\bar{e}(x, y)| < \varepsilon$, $\forall (x, y)$, where $\bar{e}(x, y)$ stands for the mean envelope. The toolbox available on http://aquador.vovve.net/IEMD/ requires $0 < \varepsilon \le 1$. As $\varepsilon$ becomes smaller, more sifting iterations are performed.

Prox-EMD achieves a decomposition by solving a convex non-smooth criterion which has texture and geometry terms. The toolbox is available on http://perso.ens-lyon.fr/nelly.pustelnik/logicielE.html.

TV-G texture-geometry decomposition provides a separation where the texture component is related with finer scales and thus can be considered as detail and the geometry component contents coarser scales and is related with the trend. The toolbox is available on http://www-rohan.sdsu.edu/~jegilles.

Both Prox-EMD and TV-G texture-geometry decomposition contain two parameters to be set: $v$ for the texture and $\rho$ for the geometry. Based on [17], we use $v = \{1, 5, 10, 50\}$ and $1 \le \rho \le 20$.

For the four methods (UOA-EMD, IEMD, Prox-EMD, and TV-G texture-geometry decomposition), we decompose the given image $X = A+D$ into two components: $X \mapsto [\hat{A}, \hat{D}]$.

For UOA-EMD and IEMD, we have $X = \hat{A} + \hat{D}$; for Prox-EMD and TV-G we have $X \approx \hat{A} + \hat{D}$. Therefore, to measure recovering capabilities of UOA-EMD and IEMD we can use either the error for $A$ or the error for $D$, since $||A - \hat{A}||_F = ||D - \hat{D}||_F$, where $|| \cdot ||_F$ stands for the Frobenius norm. For Prox-EMD and TV-G we use $(||D - \hat{D}||_F + ||A - \hat{A}||_F)/2$. For all cases, we also show the decomposition results for the parameters that minimize the error, with the images being scaled to use the full colormap.

1) Experiment 1: This experiment consists in decomposing $X = X_1 + X_2$, with

$$X_1 = \cos(6\pi(x+y)^2) + \cos(18\pi(x-y)) \\ X_2 = \cos(1.5\pi(x+y)^2), \qquad (13)$$

with $-1 \le x \le 1$ and $-1 \le y \le 1$. The image has a size of 128 x 128 pixels. The results from the four methods are shown in Fig. 2. UOA-EMD and IEMD achieve very similar results, showing good recovering capabilities. On the other hand, Prox-EMD and TV-G show rather poor performances. Figure 3 presents both the errors when estimating the components, and the computational times. UOA-EMD gives the lowest error for two thirds of the considered ranges, and a constant computational time. IEMD evidences three different regimes: when the threshold value $\varepsilon$ goes from 0 to 1, three different number of sifting iterations are used. Moreover, this first experiment shows that Prox-EMD and TV-G are unable to successfully recover AM-FM components for the image processed. Furthermore, the tuning of their two main parameters can be tricky.
<figure><img src="image_5.png" /><figcaption>Fig. 3. Errors and computational times for Experiment 1 (see text for the details).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Results for Experiment 2. Best results for the four analyzed methods when decomposing an image mixing components described Eqs. (14).</figcaption></figure>

2) Experiment 2: For this experiment we consider the decomposition of $X = X_1 + X_2$, with

$$X_1 = \begin{cases} \cos(6\pi(x+y)^2) + \cos(18\pi(x-y)), & \text{if } 0 \le y \le 1 \\ \cos(2\pi(x+y)^2), & \text{if } -1 \le y < 0 \end{cases} \\ X_2 = \begin{cases} \cos(2\pi(x+y)^2), & \text{if } 0 \le y \le 1 \\ \cos(0.5\pi(x+y)^2), & \text{if } -1 \le y < 0 \end{cases} \quad (14)$$

The bottom half of $X_1$ is equal to the top half of $X_2$. Since $X_2$ is trend to $X_1$, the same spatial frequencies are trend in one part of the image and detail in the other. A simple linear filtering is unable to separate the components. Results obtained with the four methods are shown in Fig. 4. UOA-EMD presents the best results, among others for the bottom half of the residue. The errors for Prox-EMD are close

to those of UOA-EMD (Fig. 5). However, when analyzing the residue, Prox-EMD seems to favor piecewise constant solutions, which are not suited to catch AM-FM characteristics. The best performances for the computational times are obtained with UOA-EMD.

3) Experiment 3: In this experiment we introduce amplitude modulations:

$$X_1 = \begin{cases} \cos(24\pi(x + 0.1y)), & \text{if } 0 \le y \le 1 \\ e^{-3(x^2+y^2)}(\cos(6\pi(x+y)) + \cos(6\pi(x-y))), & \text{if } -1 \le y < 0 \end{cases} \\ X_2 = \begin{cases} e^{-3(x^2+y^2)}(\cos(6\pi(x+y)) + \cos(6\pi(x-y))), & \text{if } 0 \le y \le 1 \\ e^{-3(x^2+y^2)}, & \text{if } -1 \le y < 0 \end{cases} \quad (15)$$
<figure><img src="image_7.png" /><figcaption>Fig. 5. Errors and computational times for Experiment 2 (see text for the details).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Results for Experiment 3. Best results for the four analyzed methods when decomposing an image mixing those of Eqs. (15).</figcaption></figure>

As before, the bottom half of $X_1$ equals the top half of $X_2$. The decomposition given by the four methods are presented in Fig. 6. The errors and computational times are presented in Fig. 7. For this third experiment, results from UOA-EMD and IEMD are similar. The errors are of the same order of magnitude, except for $\epsilon = 1$ for IEMD where the error is significantly higher. Although errors for Prox-EMD are low for small values of geometry parameter $\rho$, the results deserve special attention. The Prox-EMD residue presents isolated squares of almost constant value, which are not present in the original component. The first mode of Prox-EMD presents in its bottom half a pattern of vertical, horizontal and diagonal lines, which do not appear in $X_1$.

This effect is a clear consequence of a pseudo-2D approach, in which rows, columns and diagonal are separately constrained. The top half of TV-G first mode is almost equal to that of the original image $X$. Moreover, computational time remains almost constant for UOA-EMD, being thirteen times smaller than that of IEMD.

4) Experiment 4: As fourth and final experiment for this subsection, we masked the images from Experiment 3 with

$$M = \begin{cases} 1, & \text{if } x^2 + y^2 < \frac{16}{25} \\ 0, & \text{otherwise,} \end{cases} \quad (16)$$
<figure><img src="image_9.png" /><figcaption>Fig. 7. Errors and computational times for Experiment 3.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Experiment 4 Original image X</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Results for Experiment 4. Best results for the four analyzed methods when decomposing an image mixing those of Eqs. (15) masked with Eq. (16).</figcaption></figure>

in order to test the behaviour of the methods on an image with a piecewise constant zone.

The results are presented in Fig. 8. Same features as in the previous example appears. The errors (see Fig. 9) for IEMD increase, while the performance of UOA-EMD remains almost unchanged.

### C. Real Data

Real standard images are widely used by image processors in order to test their performances. We herein present the decomposition of a 512 × 512 version of the image “Barbara” via the four methods. In every case, the image

was decomposed into two modes plus residue. We used $\lambda_1 = \lambda_2 = 1$ for UOA-EMD, $\varepsilon = 0.5$ and $K = 2$ for IEMD, and $\nu_1 = \nu_2 = 1$, $\rho_1 = \rho_2 = 10$ and $K = 2$ for both Prox-EMD and TV-G. The decomposition needed around 5 minutes for UOA-EMD, five and a half hours for IEMD, 17 minutes for Prox-EMD and 50 seconds for TV-G. The results of the decompositions for the four methods are shown in Fig. 10.

The first mode of UOA-EMD captures the highest frequencies. The stripes of the woman's clothes and the tablecloth are perfectly conserved, along with the edges of the rest of the image. The second mode of UOA-EMD catches lower
<figure><img src="image_12.png" /><figcaption>Fig. 9. Errors and computational times for Experiment 4.</figcaption></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 10. Decompositions of Barbara. Results of the decompositions *via* the four methods.</figcaption></figure>

frequencies but conserves important features such as the folds of the woman's pants, of the cloth around her head and of the tablecloth. A portion of the original image is pointed out. The same zone in modes first and second is zoomed. In the zoomed zone of the first mode the stripes of the fabric can be appreciated, with no traces of the dark areas. The stripes of the fabrics are absent from this zone in the second mode, and only the folds can be seen. All the contrast of the original image is present in the residue.

On the other hand, the other three analyzed methods seem unable to capture such features. IEMD first mode includes not

only the edges and fabric stripes but also coarser oscillations such as the folds of the fabrics. The second mode is too blurry for fine structures to be identified. It also contains most part of the image contrast, leaving the residue with almost no relevant information. Prox-EMD retrieves a first mode where some edges are hard to detect, such as those of the legs of the table. The second mode and the residue are difficult to analyze. While the second mode resembles a watercolour version of the image with no fine details, the residue is too blurry and no structure can be easily identified (we must remember that Prox-EMD does not provide an exact reconstruction).
<figure><img src="image_15.png" /><figcaption>Fig. 11. Decomposition of a magnetic resonance image with UOA-EMD, from fine to coarse scales.</figcaption></figure>

Finally, TV-G first mode is the most similar, except for the bright, to the original image among first modes. The legs of the table are significantly darker than the floor, proving more than just edge detection; the same can be said about the woman's chest. Second mode and residue provide little information about the original image (as with Prox-EMD, TV-G reconstruction is not exact).

## D. Biomedical Application

We herein present the decomposition of a 256 × 256 medical image. The latter is extracted from the website of Harvard Medical School for Radiology at http://www.med.harvard.edu/AANLIB/cases/case1/mr1-tl4/029.html [12]. This image, a T2-weighted magnetic resonance image (MRI), is shown in Fig. 11. It has been recorded in a woman suffering from a glioma. The decomposition of the data into three components obtained with UOA-EMD is shown in Fig. 11, where $λ_k = 1, 1 ≤ k ≤ 3$ was used. We observe that the three components represent each a separate scale range: the first component corresponds to the highest frequencies (the finest texture); as the number of decomposition grows, the largest scales are discovered. The trend is observed on the residue. The first mode clearly reveals the contours of the lesion. These results show that UOA-EMD can also be useful for the medical field.

## V. CONCLUSION

Local and auto-adaptive decompositions such as EMD have gained an increasing attention both for signals (1D data) and images (2D data). For the process of images, most EMD-based decompositions rely on an iterative procedure (sifting process) that has the drawback of lacking convergence guarantees. The recent optimization-based approaches to EMD, on the other hand, makes no use of such process, paving the road for mathematical foundations of EMD.

In this work we proposed a 2D extension of an unconstrained optimization approach to EMD (UOA-EMD), and compared its performance to three other decompositions methods: IEMD, Prox-EMD, and TV-G. Our results on synthetic and real data show that the proposed framework overcomes the existing algorithms: it is orientation-independent, shows a low temporal complexity, and gives better performances in decomposing AM-FM images.

For the artificial AM-FM images we confirmed the locality of the method. Applying the here proposed technique to a well-known standard image ("Barbara") the capabilities of

our method, when different spatial frequencies are related to different structures, were illustrated. Finally, we showed a potential application of our method decomposing a biomedical image (MRI data).

The simplicity of the proposed framework leaves room for an nD extension ($n > 2$).

## REFERENCES

[1] R. Maani, S. Kalra, and Y.-H. Yang, "Rotation invariant local frequency descriptors for texture classification," *IEEE Trans. Image Process.*, vol. 22, no. 6, pp. 2409-2419, Jun. 2013.

[2] Z. Guo, L. Zhang, and D. Zhang, "Rotation invariant texture classification using LBP variance (LBPV) with global matching," *Pattern Recognit.*, vol. 43, no. 3, pp. 706-719, 2010.

[3] J. Bruna and S. Mallat, "Classification with scattering operators," in *Proc. IEEE Conf. Comput. Vis. Pattern Recognit.*, Jun. 2010, pp. 1561-1566.

[4] Z. Lei, T. Ahonen, M. Pietikäinen, and S. Z. Li, "Local frequency descriptor for low-resolution face recognition," in *Proc. IEEE Int. Conf. Autom. Face Gesture Recognit. Workshops (FG)*, Mar. 2011, pp. 161-166.

[5] X. Tan and B. Triggs, "Enhanced local texture feature sets for face recognition under difficult lighting conditions," *IEEE Trans. Image Process.*, vol. 19, no. 6, pp. 1635-1650, Jun. 2010.

[6] T. Ahonen, E. Rahtu, V. Ojansivu, and J. Heikkilä, "Recognition of blurred faces using local phase quantization," in *Proc. 19th Int. Conf. Pattern Recognit. (ICPR)*, Dec. 2008, pp. 1-4.

[7] S. Stanković, I. Orović, and N. Zarić, "An application of multidimensional time-frequency analysis as a base for the unified watermarking approach," *IEEE Trans. Image Process.*, vol. 19, no. 3, pp. 736-745, Mar. 2010.

[8] N. E. Huang et al., "The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis," *Proc. Roy. Soc. London, A Math. Phys. Eng. Sci.*, vol. 454, no. 1971, pp. 903-995, 1998.

[9] J. C. Nunes, Y. Bouaoune, E. Delechelle, O. Niang, and P. Bunel, "Image analysis by bidimensional empirical mode decomposition," *Image Vis. Comput.*, vol. 21, no. 12, pp. 1019-1026, 2003.

[10] C. Damerval, S. Meignen, and V. Perrier, "A fast algorithm for bidimensional EMD," *IEEE Signal Process. Lett.*, vol. 12, no. 10, pp. 701-704, Oct. 2005.

[11] A. Linderhed, "Image empirical mode decomposition: A new tool for image processing," *Adv. Adapt. Data Anal.*, vol. 1, no. 2, pp. 265-294, 2009.

[12] Z. Wu, N. E. Huang, and X. Chen, "The multi-dimensional ensemble empirical mode decomposition method," *Adv. Adapt. Data Anal.*, vol. 1, no. 3, pp. 339-372, 2009.

[13] N. Pustelnik, P. Borgnat, and P. Flandrin, "Empirical mode decomposition revisited by multicomponent non-smooth convex optimization," *Signal Process.*, vol. 102, pp. 313-331, Sep. 2014.

[14] T. Oberlin, S. Meignen, and V. Perrier, "An alternative formulation for the empirical mode decomposition," *IEEE Trans. Signal Process.*, vol. 60, no. 5, pp. 2236-2246, May 2012.

[15] M. A. Colominas, G. Schlottthauer, and M. E. Torres, "An unconstrained optimization approach to empirical mode decomposition," *Digit. Signal Process.*, vol. 40, pp. 164-175, May 2015.

[16] N. Pustelnik, P. Borgnat, and P. Flandrin, "A multicomponent proximal algorithm for empirical mode decomposition," in *Proc. 20th Eur. Signal Process. Conf. (EUSIPCO)*, Aug. 2012, pp. 1880-1884.
[17] J. Schmitt, N. Pustelnik, P. Borgnat, P. Flandrin, and L. Condat, "2D Prony-Huang transform: A new tool for 2D spectral analysis," *IEEE Trans. Image Process.*, vol. 23, no. 12, pp. 5233-5248, Dec. 2014.

[18] Z. Wu and N. E. Huang, "Ensemble empirical mode decomposition: A noise-assisted data analysis method," *Adv. Adapt. Data Anal.*, vol. 1, no. 1, pp. 1-41, 2008.

[19] J. Gilles and S. Osher, "Bregman implementation of Meyer's G-norm for cartoon + textures decomposition," Dept. UCLA Comput. Appl. Math. Rep., UCLA, Los Angeles, CA, USA, Tech. Rep. CAM11-73, Nov. 2011.