Marjaneh Taghavi
email: mtaghavi@kish.sharif.edu

Mohammad Bagher-Shamsollahi

M B Shamsollahi

L Senhadji
email: lotfi.senhadji@univ-rennes1.fr

Keywords: Ensemble Empirical Mode Decomposition, denoising, Mutual Information, Instantaneous Half Period, Interval Thresholding, ECG

HAL is

I. INTRODUCTION

LECTROCARDIOGRAM (ECG) signal is a recording of the electric potential produced by rhythmical cardiac activities. ECG is widely used in clinical applications where the accuracy of information from ECG is very important. In [1]- [START_REF] Thakor | Application of adaptive filtering to ECG analysis: Noise cancellation and arrythmia detection[END_REF], numerous methods have been proposed to denoise ECG signal. Empirical Mode Decomposition (EMD) is a typical tool for denoising signals, since EMD is very suitable for nonlinear and nonstationary signals. EMD was initially introduced by Huang [1], and immediately used in biomedical engineering applications [START_REF] Huang | Engineering analysis of biological variables: an example of blood pressure over 1 day[END_REF], [START_REF] Liang | Artifact reduction in electrogastrogram based on the empirical model decomposition method[END_REF]. EMD decomposes a signal into several Intrinsic Mode Functions (IMFs). The EMD denoising is based on the fact that noise components are often found in the first several IMFs. Thus, denoised signal is obtained by detecting noisy IMFs and removing or filtering them and then the signal is reconstructed by filtered IMFs.

EMD has some limitations in decomposing signal which has the close amplitude-frequency ranges to each other. This problem is termed 'mode mixing', and it can be happened when different modes of oscillation are appeared in a single IMF. We use Ensemble Empirical Mode Decomposition (EEMD) from [START_REF] Wu | Ensemble empirical Mode decomposition: a noise-assisted data analysis method[END_REF] to solve this problem. In EEMD, an ensemble mean is taken over a number of IMFs extracted from EMD applied to some new noisy data obtained by adding a white noise series each time to the original data. Numerous methods such as soft and hard thresholding are proposed to filter noisy IMFs [START_REF] Boudraa | EMD-based signal noise reduction[END_REF]- [START_REF] She | De-noisng of ECG based on EMD improved-thresholding and mathematical morphology operation[END_REF]. In these methods it is assumed that all IMFs are noisy and thresholding operations are applied to them. Hence, some useful information is lost especially in higher order IMFs which contain more signal information than noise.

In this paper, we combine two EMD based interval thresholding, mutual information technique and correlation coefficient test for ECG denoising. EEMD is applied to the noisy signal to obtain IMFs and then mutual information is used to classify IMFs into high frequency (HF) and low frequency (LF) groups [START_REF] Omitaomu | Empirical Mode Decomposition Technique With Conditional Mutual Information for Denoising Operational Sensor[END_REF]. We consider HF group as a strong noisy IMFs. Then we use Block Energy (BE) model [START_REF] Zhang | An Empirical Mode Decomposition Signal Denoising Method Based on Block Energy[END_REF] to filter HF group. Some of IMFs in the LF group may have dominant frequency component of the original signal. Therefore, we use correlation coefficient test to detect noise free IMFs [START_REF] Siwal | A new approach of denoising the regular and chaotic signals using Empirical Mode Decomposition: comparison and application[END_REF] and then we filter the noisy IMFs of LF group by Instantaneous Half Period model (IHP) model [START_REF] Fang | Stress wave signal denoising using ensemble empirical mode decomposition and an instantaneous half period model[END_REF]. Finally, the denoised signals are reconstructed by adding filtered IMFs from LF and HF groups to the remaining nonfiltered IMFs which are in the LF group. Therefore in this method the useful information of signal in each IMF can be preserved. In order to demonstrate the performance of our proposed method, we compare its performance with IHP model and BE model [START_REF] Zhang | An Empirical Mode Decomposition Signal Denoising Method Based on Block Energy[END_REF] in terms of improvement in signal to noise ratio (SNRimp) and Mean Square Error (MSE). For this comparison, we have used different type of ECG signals from the MIT-BIH database.

The outline of the paper is as follows. In section II, a brief review of theoretical EMD, EEMD, BE model, IHP model and mutual information are presented. Furthermore, our proposed method is explained. Section III illustrates the results of the three methods and their performance will be compared. Finally, conclusion is made in section IV.

II. METHODOLOGY

A. Description of EMD and EEMD

EMD is able to decompose a given signal into sum of IMFs by iterative sifting process. An IMF should satisfy the following two conditions: (1) the number of extrema and the number of zero crossings are either equal or differ at most by one; and (2) the mean of its upper and lower envelopes equals zero [ The above procedure to extract the IMF is called the sifting process. 4. Repeat step 1 to 3 until the resulting signal, , has two conditions of IMF, then the first IMF, , is obtained. There are some criteria for terminating the sifting process. A commonly used one is the Standard Deviation (SD):

where k is the number of sifting iteration. The value of SD is usually in the range [0.2, 0.3].

Treat the residue signal

as a new signal and apply the sifting process again; 6. Repeat step 1 to 5 to obtain IMFs until the residue of signal, step 5, is constant. The signal is expressed as:

where I is the number of IMFs and is the final residue.

There are a number of considerations used in implementing EMD such as number of sifting, stop criteria and envelop interpolation [1], [START_REF] Rilling | On empirical mode decomposition and its algorithms[END_REF]. As an improved version of EMD, EEMD was originally developed to reduce the mode-mixing effect [START_REF] Wu | Ensemble empirical Mode decomposition: a noise-assisted data analysis method[END_REF]. In this method white noise is added to the original signal for many times with different series. Then average IMFs are obtained by calculating mean of IMFs which are obtained by applying EMD on each trial, the final IMFs are:

where I is the number of IMFs and M is the number of trials. The added noises are different in each trial and their effect can be cancelled out on average, under the condition that M is high enough [START_REF] Huang | Engineering analysis of biological variables: an example of blood pressure over 1 day[END_REF].

B. EMD Based denoising using interval thresholding

In this section, we study two methods to detect noisy parts of IMFs and then filter them to obtain denoised IMFs.

1) Instantaneous Half Period (IHP)

In [START_REF] Fang | Stress wave signal denoising using ensemble empirical mode decomposition and an instantaneous half period model[END_REF] a criterion is proposed to detect the noisy parts of IMFs based on the time interval between two zero-crossings of each IMF. In this method, the noisy parts are detected if the time interval between two zero-crossings is smaller than predefined threshold. Then, the noisy parts are removed and the signal is reconstructed by remaining parts of IMFs [START_REF] Fang | Stress wave signal denoising using ensemble empirical mode decomposition and an instantaneous half period model[END_REF].

Steps of this method are as follows: 1. Find zero crossings of each IMF; 2. Obtain time interval between two adjacent zero-crossings of each IMF;

3. Detect noisy parts of each IMF as follows and replace these part with zero:

4. Denoised signal can be obtained as:

where is the denoised signal, is the denoised IMF,

, is the number of IMFs. is zero-crossing of IMF.

is the time when zero crossing of IMF emerges. is the time interval between two zero-crossings of IMF, . In this method, finding the value of threshold is very important. The optimum threshold is obtained based on the maximum frequency of signal [START_REF] Fang | Stress wave signal denoising using ensemble empirical mode decomposition and an instantaneous half period model[END_REF]:

where denotes the maximum frequency of the signal and is the constant coefficient used to determine the frequency range of the retained oscillation [START_REF] Fang | Stress wave signal denoising using ensemble empirical mode decomposition and an instantaneous half period model[END_REF].

2) Block Energy Model In [START_REF] Zhang | An Empirical Mode Decomposition Signal Denoising Method Based on Block Energy[END_REF], an interval thresholding method is proposed to detect noisy parts of IMFs based on the energy of samples which are between two zero-crossings of each IMF. The threshold value is determined for each IMF. The waveform between the two adjacent zero crossings of IMF is retained when the energy of samples which are between the two adjacent zero crossings is larger than the predefined threshold, whereas the waveform is set to zero. Steps of interval thresholding are as follows: 1. Find zero crossings of each IMF; 2. Determine the samples which are between the two adjacent zero crossing of IMF:

where is the sample points set between two adjacent zero crossing of IMF and is the number of sample points. where denotes the filtered IMF, is the energy of sample point set between and zero crossings of IMF, Ti is the threshold value for IMF. The denoised signal is obtained by using question [START_REF] Liang | Artifact reduction in electrogastrogram based on the empirical model decomposition method[END_REF].

In [START_REF] Flandrin | EMD equivalent filter banks, from interpetation to applications[END_REF] a denoising method based on the IMF energies is proposed. The IMF energy should be an exponentially [START_REF] Liang | Artifact reduction in electrogastrogram based on the empirical model decomposition method[END_REF] (

(5) [START_REF] Wu | Ensemble empirical Mode decomposition: a noise-assisted data analysis method[END_REF] (1) decreasing function of the IMF index [START_REF] Flandrin | EMD equivalent filter banks, from interpetation to applications[END_REF]. Therefore the first IMF has the highest amount of energy. As in this paper we consider white Gaussian noise, the noise-only IMF energies for can be approximated as a function of the first IMF energy [START_REF] Flandrin | EMD equivalent filter banks, from interpetation to applications[END_REF] according to where E1 is the energy of the first IMF, is a scaling factor and is estimated as 2 [START_REF] Flandrin | EMD equivalent filter banks, from interpetation to applications[END_REF], is experimentally estimated in [START_REF] Flandrin | EMD equivalent filter banks, from interpetation to applications[END_REF] and H is the Hurst exponent as for specific EMD implementation [START_REF] Flandrin | EMD equivalent filter banks, from interpetation to applications[END_REF]. The energy equation of the first IMF is given:

The optimal value of threshold in this method is suggested in [START_REF] Zhang | An Empirical Mode Decomposition Signal Denoising Method Based on Block Energy[END_REF]:

where Ti is the determined threshold value for i th IMF, c is a constant experimentally found to take values 1 to 0.7, Vi is the noise variance of i th IMF which is calculated in [START_REF] Omitaomu | Empirical Mode Decomposition Technique With Conditional Mutual Information for Denoising Operational Sensor[END_REF] for i>1 and N is the number of IMF samples. We assume that most of the noise energy is in the first IMF [START_REF] Tsolis | Signal Denoising Using Empirical Mode Decomposition and Higher Order Statistics[END_REF], so we use the absolute median deviation of the first IMF for V1 as in [START_REF] Boudraa | Denoising via empirical mode decomposition[END_REF].

C. IMF classification based on Mutual information

Mutual information quantifies the relationship between two random variables [START_REF] Cover | Elements of information theory[END_REF]. Mutual information has been used in different applications such as feature selection [START_REF] Fleuret | Fast binary feature selection with conditional mutual information[END_REF] and signal denoising [START_REF] Omitaomu | Empirical Mode Decomposition Technique With Conditional Mutual Information for Denoising Operational Sensor[END_REF]. The mutual information is used in EMD based denoising for separating IMFs into high frequency and low frequency groups. Steps of classification of IMFs based on mutual information are as follows [START_REF] Omitaomu | Empirical Mode Decomposition Technique With Conditional Mutual Information for Denoising Operational Sensor[END_REF] where is the original signal, and are the IMF in the HF and LF groups, I is the number of IMFs and residue is concidered as an IMF in the LF group.

D. Proposed method

In this paper, we propose a denoising method based on the combination of EMD based interval threshoding, mutual information and correlation coefficient test. EEMD is applied to the noisy ECG signal to decompose it into several IMFs. Then we use mutual information to classify IMFs into high frequency and low frequency groups as explained in section C. We use BE model to filter HF group as in section B. In this interval thresholding, a determined threshold value for each IMF is defined in order to gain more accuracy and performance. Moreover it is possible to guess the samples are either noise or signal, based on the energy of samples. Whereas based on the absolute amplitude of the i th IMF, it is impossible to know whether samples are noise or not [START_REF] Zhang | An Empirical Mode Decomposition Signal Denoising Method Based on Block Energy[END_REF]. The IMFs in LF group contain low frequencies of the original signal or noise or both of them. Therefore, we use correlation coefficient test (c.c) to detect IMFs with dominant frequency of the original signal [START_REF] Siwal | A new approach of denoising the regular and chaotic signals using Empirical Mode Decomposition: comparison and application[END_REF].

where is the correlation coefficient between IMFs ( ) and noisy signal in frequency domain, and are their total energies, and are their DFTs and is the number of data points.

The correlation coefficient ranges from values +1 to -1. If the absolute value of correlation coefficient of IMFs in the LF group is greater than 0.2, it indicates that the IMFs have dominant frequency component of the original signal [START_REF] Siwal | A new approach of denoising the regular and chaotic signals using Empirical Mode Decomposition: comparison and application[END_REF]. The remaining IMFs contain noise component and then should be filtered. The frequency range of signal is often lower than the frequency of noise. Therefore the time interval between two adjacent zero crossings of signal dominated oscillation is longer than noise dominated oscillation [START_REF] Fang | Stress wave signal denoising using ensemble empirical mode decomposition and an instantaneous half period model[END_REF]. So we use IHP model to filter noisy LF IMFs as explained in section B. Finally, the denoised signal is obtained by adding filtered IMFs from two interval thresholding output and the remaining non-filtered IMFs in LF group. The flowchart of the proposed algorithm is shown in Fig. 1.

III. RESULTS

Our method is applied to the several standard data sets from PhysioBank [START_REF] Physiobank | physiologic signal archives for biomedical research[END_REF], including MIT-BIH Arrhythmia database (DB1) [22], MITBIH normal sinus rhythm database (DB2) [START_REF]The MIT-BIH Normal Sinus Rhythm Database[END_REF] and the MIT-BIH noise stress test database (DB3) [START_REF]The MIT-BIH Noise Stress Test Database[END_REF]. Ten records are arbitrarily chosen from DB1. Each signal contains 60 seconds recording and is digitized at 360 Hz. The first 60 seconds of all 18 records of DB2 are used too. White Gaussian noises are added artificially to the ECG signals to make -5dB, 0dB, 5dB, input SNRs. The performance of the proposed method is evaluated by (

improvement of signal to noise ratio and mean square error (MSE). For each input SNR, 100 times runs are performed to get an averaged and MSE value. The can be represented as follows:

where and are formulated as follows:

where is the clean ECG signal, is the noisy signal, is the denoised signal and is the noise. We have also used MSE criterion for evaluating quality of the information which is preserved in the denoised signal, it can be represented as follows:

where L is the length of the signal. Fig. 2 shows an example of a clean signal, noisy signal and denoised signal using our proposed method. The tested signal in Fig. 2 is the record 100 from DB1 with an additive Gaussian noise that makes a 0 dB input SNR. We apply EEMD on this record to obtain IMFs. The total number of IMFs is 15 considering residue as well. Then we use mutual information to classify IMFs into HF and LF group. The IPI index is 7 and the index is 9. So the sum of 8 first components (IMF1 to IMF8) gives the final HF group and the sum of remaining IMFs gives the final LF group. Then we calculate the correlation coefficient to detect the IMFs in the LF group which should be filtered. For this signal, the IMFs 9, 12, 14, 15 are obtained. So the remaining components in LF group (IMF10, IMF11 and IMF13) need to be filtered. We use BE model to filter IMFs in the HF group with a determined threshold for each IMFs and IHP model is applied to IMF10, IMF11 and IMF13. Finally we add filtered IMFs and the IMFs with index 9, 12, 14, 15 to obtain denoised signal as shown in Fig. 2.

Table I shows the values for output SNR for different values of input SNR of the ECG records of DB1 for three methods: our proposed method, EEMD-based IHP model (IHP) [START_REF] Fang | Stress wave signal denoising using ensemble empirical mode decomposition and an instantaneous half period model[END_REF] and EEMD-based BE model [START_REF] Zhang | An Empirical Mode Decomposition Signal Denoising Method Based on Block Energy[END_REF]. We can compare their performance in SNR improvement from [START_REF] Boudraa | Denoising via empirical mode decomposition[END_REF].

As can be seen in Table I, in all three methods, by increasing the input SNR, ECG signal becomes cleaner and therefore the SNR improvement is lower. From Table I, it is clear that our proposed method has better performance than the other two methods for all inputs SNR (i.e.-5, 0, 5, dB) and all ECG records.

Table II shows the average values of MSE for ten ECG records from DB1 as in Table I. We can see the denoised signal which is obtained by our method has a less error than two other methods. We plot the output SNR versus different input SNRs for all 18 records of DB2 in Fig. 3. The IHP and BE model assume that the noise in the signal is spread over all IMFs and the threshold value is determined for each IMF. Therefore, some useful information may be lost, especially when input SNR is increased. We can see the BE model has larger output SNR than IHP model. Therefore, according to the [START_REF] Boudraa | Denoising via empirical mode decomposition[END_REF], the BE model has a better SNR improvement. It is clear that BE model has more accuracy in detecting noisy parts of IMFs than the IHP model. So, we use BE model to filter IMFs in the HF group which contain more noise components and we use IHP model to filter noisy IMFs in the LF group to preserve signal information. In our proposed method all IMFs are filtered with minimum loss of useful information, and maximum noise removal is achieved in comparison to the single method of IHP model or BE model. This is obtained by using mutual information, correlation coefficient test, and two different interval thresholding approaches. So, this Model has more accuracy in recognizing and filtering noisy IMFs. The performance of our proposed method is evaluated by real noise. Fig. 4 shows denoising results for different input ECG records with an additive EMG noise. The input SNR of records 118e24 and 118e18 from DB3 is 6dB and for record 116 from dB1 is 10dB. The output SNRs which are obtained by applying our method is 16.7dB for 118e24, 20.4dB for 116 and 17.4dB for 118e18. As it can be seen, the denoised signal follows the clean ECG morphology.

IV. CONCLUSION AND DISCUSSION

This paper presents a method based on EEMD to reduce noise from corrupted ECG signals. In this method, EEMD is used to resolve mod mixing problem of EMD. Moreover, two interval thresholding approaches in combination with mutual information and correlation coefficient test are used to preserve signal information. Our method provides a larger SNR improvement, especially for fewer input SNRs, where the original signal is lost in noise. The denoised signal obtained by our method follows better the original signal and has a less MSE than the methods based on IHP or BE model. Thus, by using our method, the noise components will be removed from the signal with minimum loss of useful information in all input SNRs.

[1] N. E. Huang

3 .

 3 Determine the energy of each sample point set: 4. Determine threshold value for each IMF by using equation (13); 5. Detect noisy parts of each IMF and filter it as follows:

Fig. 1 .

 1 Fig. 1. The flowchart of the proposed algorithm for performing denoising of a given signal.

Fig. 2 .

 2 Fig. 2. Clean ECG, Noisy ECG with 5 dB SNR and ECG denoised by the proposed algorithm

Fig. 3 .

 3 Fig. 3. Output SNR versus different input SNRs for the first 60 seconds of all 18 records of DB2.

Fig. 4 .

 4 Fig. 4. Typical filtering results of proposed method for different input ECGs with an additive real EMG noise. Dash-dot is the clean signal, solid line is the denoised signal and dash line is the nosiy signal. (a) Record 118e24 from dB3 (corresponding to an input SNR of 6dB). (b) Record 116 from dB1 (corresponding to an input SNR of 10 dB). (c) Record 118e18 from DB3 (corresponding to an input SNR of 6 dB).

TABLE I OUTPUT

 I SNR FOR PROPOSED METHODS, IHP MODEL AND BE MODEL TABLE II PERFORMANCE EVALUATINO OF PROPOSED MTHODS, IHP MODEL AND BE MODEL BASED ON MSE

  , Z. Shen, S.R. Long, M.C. Wu, H.H. Shih, Q. Zheng, N.-C. Yen, C.C. Tung, H.H. Liu, "The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis", Proc. R. Soc. London., vol. 454, pp. 903-995, 1998.

	Algorithm	-5	Input SNR (dB) 0	5
	Proposed	0.091	0.064	0.036
	BE	0.112	0.075	0.047
	IHP	0.147	0.092	0.085


Thierry Géraud

Yongchao Xu

Edwin Carlinet

Nicolas Boutry

Introducing the Dahu Pseudo-Distance

Keywords: discrete topology, minimal path, minimum barrier distance, set-valued maps

The minimum barrier (MB) distance is defined as the minimal interval of gray-level values in an image along a path between two points, where the image is considered as a vertex-valued graph. Yet this definition does not fit with the interpretation of an image as an elevation map, i.e. a somehow continuous landscape. In this paper, based on the discrete set-valued continuity setting, we present a new discrete definition for this distance, which is compatible with this interpretation, while being free from digital topology issues. Amazingly, we show that the proposed distance is related to the morphological tree of shapes, which in addition allows for a fast and exact computation of this distance. That contrasts with the classical definition of the MB distance, where its fast computation is only an approximation.

Introduction

Que la montagne de pixels est belle. Jean Serrat.

The minimum barrier (MB) distance [START_REF] Strand | The minimum barrier distance[END_REF] is defined as the minimal interval of gray-level values in an image along a path between two points. For instance, the image of Fig. 1(a) can be seen as a graph depicted in 1(b), and the path between the red points which is depicted in blue is minimal: the sequence of values is 1, 0, 0, 0, 2 so the interval is [0, 2] and the distance is 2.

Basic Notions

Let us consider a discrete finite set X. The set of all subsets of X is denoted by P(X), so E ∈ P(X) means E ⊆ X. X is endowed with a discrete neighborhood relation; the set of neighbors of x ∈ X is denoted by N X (x) ⊂ X, and satisfies:

x ∈ N X (x) iff x = x and x ∈ N X (x ).

A discrete path π of X is a sequence π = π 1 , .., π i , .., π k , where k is the length of the sequence, such as π i ∈ X and π i+1 ∈ N X (π i ). A sequence x, .., x is called a path between x and x ; the set of all paths between x and x is denoted by Π(x, x ), and we use π(x, x ) to denote an element of Π(x, x ).

A connected component Γ of E ⊂ X is a subset of E which is path-connected (∀ x = x ∈ Γ, Π(x, x ) = ∅) and maximal (∀ x ∈ Γ, ∀ x ∈ E \Γ, Π(x, x ) = ∅). The set of the connected components of E is denoted by CC(E), and given x ∈ X, the connected component of E containing x is denoted by C(E, x) ∈ CC(E)-by convention, if x ∈ E then C(E, x) = ∅. Given a particular element x ∞ ∈ X, the cavity-fill-in operator (in 2D, hole-filling operator) is denoted by Sat and defined by: Sat

(E, x ∞ ) = X \ C(X \E, x ∞ ).
Throughout this paper, we will only consider functions having for domain a set X being discrete, finite, and path-connected.

Scalar vs. Set-Valued Maps

A gray-level image is typically a bounded scalar function such as X is a subset D ⊂ Z n with n ≥ 2, and its codomain is Z or R. It is the case of the image depicted in 1(a). Now, if we look at the surface of Fig. 1 A set-valued map u : X → P(Y ) is characterized by its graph, Gra(u) = {(x, y) ∈ X × Y ; y ∈ u(x)} [START_REF] Aubin | Set-Valued Analysis[END_REF]. One definition of the inverse by a set-valued map u of a subset M ⊂ Y is u (M ) = {x ∈ X; u(x) ⊂ M }, and it is called the core of M by u. Let us assume that X and Y are metric spaces. The "natural" extension of the continuity of single-valued functions to set-valued maps is characterized by the following property: u is continuous iff the core of any open subset of Y is an open set of X.

In [START_REF] Géraud | A quasi-linear algorithm to compute the tree of shapes of n-D images[END_REF], the authors have defined the notion of level sets for set-valued maps; at λ ∈ Y , the lower and upper level sets of u are respectively:

[u λ] = {x ∈ X; ∀ y ∈ u(x), y < λ} and [u λ] = {x ∈ X; ∀ y ∈ u(x), y > λ}.
(1) Such definitions are fundamental because they are some means to build some morphological tools on set-valued maps (we will use them later in this paper in Sec. 2.4 and 4.1).

A scalar function u : X → Y can be simply translated into a set-valued map [START_REF] Mazo | Digital imaging: A unified topological framework[END_REF]). Fig. 2. A cubical complex for space X (left) and a tree of shapes for u (right).

• u : X → I Y with ∀ x ∈ X, • u(x) = {u(x)}. Between set-valued maps of X → h h' ' h' h h' ' h' (a) Faces of H n (from
P(Y ) we can define the relation < -by:

u 1 < -u 2 ⇔ ∀ x ∈ X, u 1 (x) ⊆ u 2 (x)
. By abuse of notation, we will write u <u to state that • u <u, meaning that we have ∀ x ∈ X, u(x) ∈ u(x). An illustration can be seen later: the map u of Fig. 3(d) and the map u of Fig. 3(c) are such that u < -u. In the following, we will only consider interval-valued maps, that is, maps of X → I Y .

Cubical Complexes

From the sets H 1 0 = {{a}; a ∈ Z} and H 1 1 = {{a, a + 1}; a ∈ Z}, we can define

H 1 = H 1 0 ∪ H 1 1
and the set H n as the n-ary Cartesian power of H 1 . If an element h ⊂ Z n is the Cartesian product of d elements of H 1 1 and n-d elements of H 1 0 , we say that h is a d-face of H n and that d is the dimension of h; it is denoted dim(h). The set of all faces, H n , is called the nD space of cubical complexes. Figure 2(a) depicts a set of faces {h, h , h } ⊂ H 2 where h = {0}×{1}, h = {1}×{0, 1}, and h = {0, 1}×{0, 1}, the dimension of these faces being respectively 0, 1, and 2; they are depicted as subsets of Z 2 (left) and as geometrical objects (right).

The pair (H n , ⊆) forms a poset, from which can be derived a T0-Alexandroff topology on H n . The closure operator on subsets of H n is denoted by cl . The n-faces are the minimal open sets of H n , and the set of n-faces, denoted by H n 1 , is the n-Cartesian product of H 1 . The faces of H n can be arranged onto a grid, called Khalimsky's grid, and the inclusion between faces leads to a neighborhood relationship between them, depicted in orange in Figure 2(a) (right). This neighborhood relationship will be used in Sec. 3.2 in order to define paths in H n .

Tree of Shapes

Given a gray-level image u : X → Y and any scalar λ ∈ Y , the lower level sets are defined as [u < λ] = {x ∈ X; u(x) < λ}, and the upper level sets as [u ≥ λ] = {x ∈ X; u(x) ≥ λ}. Considering the connected components of these sets, and using the cavity-fill-in operator, the tree of shapes (ToS) of an image u is classically [START_REF] Monasse | Fast computation of a contrast-invariant image representation[END_REF] defined by:

S(u) = { Sat(Γ ); Γ ∈ CC([u < λ]) ∪ CC([u ≥ λ]) } λ . (2)
An image and its tree of shapes are depicted in Fig. 2. An element of S(u) is called a shape; it is a connected component of X with no cavity, and its boundary is a level line of u. Every shape corresponds to a node of the tree; for instance, in Fig. 2(b) (right), the sub-tree rooted at node "B" corresponds to the shape B ∪ D ∪ E. Keeping the level of every node-such as displayed in Fig. 2(b) (right)-allows to reconstruct the image from its tree. The tree of shapes of an image u is a morphological representation of u which makes it easier to deal with the image contents [START_REF] Caselles | Geometric Description of Images as Topographic Maps[END_REF]. Storing [START_REF] Carlinet | A comparative review of component tree computation algorithms[END_REF] and computing [START_REF] Géraud | A quasi-linear algorithm to compute the tree of shapes of n-D images[END_REF][START_REF] Crozet | A first parallel algorithm to compute the morphological tree of shapes of nD images[END_REF] the tree of shapes can be done very efficiently. Last, let us mention that the tree of shapes is a versatile tool to perform image filtering [START_REF] Xu | Connected filtering on tree-based shape-spaces[END_REF], and a very relevant structure to perform some computer vision tasks [START_REF] Caselles | Topographic maps and local contrast changes in natural images[END_REF][START_REF] Cao | A Theory of Shape Identification[END_REF][START_REF] Xu | Hierarchical segmentation using tree-based shape spaces[END_REF].

Let us now recall some results from [START_REF] Najman | Discrete set-valued continuity and interpolation[END_REF] about defining the tree of shapes on set-valued maps. Given an interval-valued map u : H n → I Y , using Eq. 1, its set of level sets is defined by

T(u) = λ CC([u λ]) ∪ CC([u λ]
). The set of connected components of X defined by:

S(u) = { Sat(Γ ); Γ ∈ T(u) } (3)
can be arranged into a tree for some very particular maps u. Every intervalvalued maps u does not have a tree of shapes. Yet, if all the level sets of u are "well-composed" [START_REF] Boutry | How to make nd functions well-composed in a self-dual way[END_REF], then the tree exists.

Minimum Barrier Distance

As defined in [START_REF] Strand | The minimum barrier distance[END_REF], the barrier of a path π in a gray-level image u is τ u (π) = max πi∈π u(π i ) -min πi∈π u(π i ), and the minimum barrier distance between x and x in u is:

d MB u (x, x ) = min π∈Π(x,x ) τ u (π). ( 4
)
It is actually a pseudo-distance since, ∀ u and ∀ x, x , x ∈ X, it verifies:

d MB u (x) ≥ 0 (non-negativity); d MB u (x, x) = 0 (identity); d MB u (x, x ) = d MB u (x , x) (symmetry); and d MB u (x, x ) ≤ d MB u (x, x ) + d MB u (x , x ) (subadditivity
). Yet it is not a distance because the positivity property (x = x ⇒ d MB u (x, x ) > 0) does not hold. Indeed we can have d MB u (x, x ) = 0 for some x = x , e.g., in flat zones. This pseudo-distance have been recently used in [START_REF] Zhang | Minimum barrier salient object detection at 80 FPS[END_REF] and [START_REF] Tu | Real-time salient object detection with a minimum spanning tree[END_REF] to perform salient object segmentation, and it allows to achieve some leading performance (in terms of efficiency and accuracy) as compared to state-of-the-art methods.

The Dahu Pseudo-Distance

To tackle the issue explained in the introduction, we need to map a discrete image into a continuous representation such as 1(c); yet for computations to remain tractable we want to stay in a discrete setting.

A Convenient Continuous Yet Discrete Representation

[2,4] [0,2] [2,4] {2} [0,3] [0,3] [2,4] 0 2 2 3 2 3 6 [0,3] [2,3] [2,6] [3,6] [0,2] {2} [0,3] [0,3] [3,6]
(a) u median . (c) u . the 2D case. In the general nD case, these parts are called critical configurations, and an image without any critical configuration is said "digitally well-composed" (DWC) [START_REF] Boutry | On making nD images well-composed by a self-dual local interpolation[END_REF]. Let us denote by u an interpolation of the scalar image u : Z n → Y , defined on the subdivided space Z 2 n , and taking its values either in Y or in another scalar space Y . A continuous yet discrete representation of u is now u :

(u : Z n → Y ) step 1 ----→ (u : Z 2 n → Y ) step 2 ----→ ( u : H 2 n → I Y ). (7)
Under the assumption of u being digitally well-composed, the interval-valued image u allows for representing u in a continuous yet discrete way, and without any topological inconsistency. We will say that such maps u are path-consistent. Now we can provide a new definition of the Dahu pseudo-distance:

∀ x, x ∈ Z n , â u (x, x ) = min u < -u d MB u (h x , h x ). ( 8
)
As compared with the initial (naive) definition given by Eq. 6, we have just replaced u by u . Several subdivision-based interpolations of nD scalar images are known to produce DWC images. It is the case of the interpolation with the max operator, used in [START_REF] Géraud | Self-duality and discrete topology: Links between the morphological tree of shapes and well-composed gray-level images[END_REF] for the quasi-linear computation of the tree of shapes of a scalar image, as defined by Eq. 2. In the 2D case-and not in nD with n > 2-the only self-dual local interpolation verifying strong invariant properties leading to DWC images is obtained with the median operator [START_REF] Géraud | Self-duality and discrete topology: Links between the morphological tree of shapes and well-composed gray-level images[END_REF]. It is depicted in Fig. 5(a) (after the step 2 u → u being applied); Fig. 5(b) shows that the level lines do not cross each other. For the general nD case, a self-dual non-local interpolation has been defined [START_REF] Boutry | How to make nd functions well-composed in a self-dual way[END_REF], which is illustrated in Fig. 5(c) (again, after applying the step 2 to get an interval-valued representation).

Computing the Dahu Pseudo-Distance

The tree of shapes of an image is also called topographic map. Since the Dahu distance is related to the notion of topography, the tree of shapes may be a good candidate tool to compute any Dahu distance â u (x, x ).

Tree of Shapes of Interval-Valued Maps

Since the scalar interpolations u are digitally well-composed, the tree of shapes of the interval-valued maps u , namely S( u ), as defined by Eq. 3, exists. One can easily understand that it is due to the fact that u is path-consistent: reminding that a level line is the boundary of a shape, since two level lines at different levels cannot cross each other, it implies that shapes are either disjoint or nested. Actually, the set S( u ) verifies the following properties: i) Every element of this set is a connected component of H n , having no cavity, and being a regular open set. ii) The set S( u ) is a self-dual decomposition of H n with respect to u. iii) The boundary of every component is an iso-level line of u, so:

∀ S ∈ S( u ), ∃ µ ∈ Y such that ∀ h ∈ cl (S) \S, µ ∈ u (h).
iv) A mapping that associates a level, denoted by µ(S), to every element S ∈ S( u ) makes this set being a representation of u . We will say that µ(S) is the level of the shape S (we will use this mapping in the next section).

When the Dahu Pseudo-Distance and the Tree of Shapes Meet

Intuitively it is easy to understand that the minimal path between two points of an image u, w.r.t. the Dahu pseudo-distance, corresponds to a path of nodes on the tree of shapes S( u ). If we look back at Fig. 2(b) (left), the path between the two points (x, x ) indicated by red bullets in u starts from region B, then goes through A and C, and finally ends in region F. When considering the tree of shapes, in Fig. 2(b) (right), this path is exactly the same, and it crosses the level lines, depicted by dashed lines both in the image and on the tree. Such a path is minimal because every path in Π(x, x ) should at least cross this same set of level lines to go from x to x . So the Dahu pseudo-distance corresponds to the level dynamics of this set of lines. In both the continuous version in Fig. 1(c) and its discrete representation in Fig. 3(b), we can directly read that the distance between the red points is 1 (the level lines are the dashed green lines).

Given a node t of a tree, let us denote by par(t) the parent node of t in the tree, and by lca(t, t ) the lowest common ancestor of the nodes t and t . Let us denote by t x the node of S( u ) corresponding to x ∈ X; this node corresponds to the smallest shape of S( u ) containing h x . We can concatenate the two sequences of nodes π x = t x , par(t x ), par 2 (t x ), ..., lca(t x , t x ) and π x = t x , par(t x ), par 2 (t x ), ..., lca(t x , t x ) , to form the new sequence:

π(t x , t x ) = π x π -1 x = t x , par(t x ), ..., lca(t x , t x ), ..., par(t x ), t x . (9)
Note that, if t x is an ancester of t x (or the reverse way), this writing is slightly different: π(t x , t x ) = t x , par(t x ), ..., t x . In the case of the example of Fig. 2(b), the lca of the nodes t B and t F is t A , and the sequence of nodes is t B , t A , t C , t F (depicted on the right), and it corresponds on the tree to the minimal path in the image space (depicted on the left).

To every shape S ∈ S( u ) corresponds a node, say t S , of the tree of shapes. Let us then write abusively µ(t S ) = µ(S), where µ(S) is the level of the shape as defined in Sec. 4.1. Eventually, the Dahu pseudo-distance can then be reexpressed by:

â u (x, x ) = max t ∈ π(tx,t x ) µ(t) - min t ∈ π(tx,t x ) µ(t), (10)
where the nodes considered and the mapping µ are related to the tree S( u ). Practically, to compute â u (x, x ) we just have to update the interval values (min and max of µ) starting respectively from the nodes of x and of x , up to their lca. For instance, on π x we iteratively compute: µ min ( t x , ..., par i+1 (t x ) ) = min( µ min ( t x , ..., par i (t x ) ), µ(par i+1 (t x )) ), and likewise for µ max . Eventually, â u (x, x ) is the length of the span of both intervals, that is: â u (x, x ) = max( µ max (π x ), µ max (π x ) ) -min( µ min (π x ), µ min (π x ) ). On a regular PC, running GNU Linux with a 2.67GHz CPU, computing one million distances between points randomly taken in the lena image takes about 1 second. Some code is available from: http://publications.lrde.epita.fr/geraud.17.ismm

Illustration on Document Image Segmentation

Figure 6 gives an illustration on a camera-based document image taken from the SmartDoc competition organized at ICDAR 2015 [START_REF] Burie | ICDAR 2015 competition on smartphone document capture and OCR (SmartDoc)[END_REF]. Let us first recall that the tree of shapes has already been used to such an application [START_REF] Carlinet | MToS: A tree of shapes for multivariate images[END_REF][START_REF] Xu | Hierarchical segmentation using tree-based shape spaces[END_REF].

From a scalar image u (Fig. 6(a)), we compute its tree of shapes S( u ) (Fig. 6(b)), which allows for computing very efficiently a distance map from the top left point (0, 0) (Fig. 6(c)) using the Dahu distance. One can see that the document is now well contrasted w.r.t. its environment. A saliency map corresponding to the Dahu distance is depicted in Fig. 6(d).

In addition to the references given throughout the paper, related either to the context of our work or to its theoretical background, several other works shall be mentioned. Actually the notion of path on tree nodes also appear in [START_REF] Carlinet | Morphological object picking based on the color tree of shapes[END_REF] to assign to each node a label and then perform some scribble-based segmentation, in [START_REF] Carlinet | MToS: A tree of shapes for multivariate images[END_REF] to count level lines and build a multi-variate tree of shapes, and in [START_REF] Dubrovina | Image editing using level set trees[END_REF] to compute some curvilinear variation and separate an object from its background. The notion of barrier used on the tree of shapes is also close to the one of shape saliency used in [START_REF] Xu | Hierarchical image simplification and segmentation based on Mumford-Shah-salient level line selection[END_REF]. Last, about distance transforms and minimal path approaches, the reader can refer to the recent survey [START_REF] Saha | Digital topology and geometry in medical imaging: A survey[END_REF].

Conclusion

We have presented a new pseudo-distance, the Dahu distance, which is a variant of the minimum barrier (MB) distance on a continuous-yet discreterepresentation of images. There are two major advantages of this distance as compared to the "classical" MB distance:

-It considers that the input discrete images are actually defined in a continuous world, that is, that their elevation map is a surface.

-Computing the Dahu distance for many couples of points (x, x ) is very efficient, thanks the pre-computation of a tree of shapes; furthermore, the distances â u (x, x ) obtained using this tree are exact. That contrasts with the efficient but approximate computation of the MB distance given in [START_REF] Tu | Real-time salient object detection with a minimum spanning tree[END_REF].

Some other practical interests of this new distance over its classical definition remain to be proved. Yet we believe that there are numerous perspectives for applications, which are made possible thanks to the fast tree-based computation.

Another perspective is to adapt this distance to color images; that can be done using the multi-variate tree of shapes (MToS) defined in [START_REF] Carlinet | MToS: A tree of shapes for multivariate images[END_REF].

Acknowledgments. The authors want to thank Guillaume Tochon (new evangelist with δ γ ε

) and Laurent W. Najman for their valuable comments on an erstwhile version of this paper, and the reviewers for their helpful remarks. This work has been conducted within the mobidem project, part of the "Systematic Paris-Region" and "Images & Network" Clusters. This project is partially funded by the French Government and its economic development agencies.

  as a graph... (c) or as a surface.

Fig. 1 .

 1 Fig. 1. Two different minimal barrier distances: (b) gives 2, while (c) gives 1.

  (c), its elevation cannot be defined by a scalar function. Yet it can be a mapping in R 2 → I R , where I Y denotes the space of the intervals of Y . Given two intervals a = [y a , y a ] and b = [y b , y b ], the interval [min(y a , y b ), max(y a , y b )] is denoted by span{a, b}. Now let us recall some definitions related to set-valued maps.

  An image u and its tree of shapes.

  This section explain how to represent a gray-value image u : D ⊂ Z n → Y by an interval-valued map u : D H ⊂ H n → I Y , and by a new tree of shapes S( u).

  0

Fig. 5 .

 5 Fig. 5. Interpolation-based interval-valued maps which are path-consistent.

  (a) Input image u. (b) Level lines of u (every 5 levels). (c) Distance map âu((0, 0), x).(d) Saliency map.

Fig. 6 .

 6 Fig. 6. Using the Dahu pseudo-distance on an actual application.

First, to every

A simple example of a transform of a gray-level image u, given in Fig. 3(a), into an interval-valued map u is given in Fig. 3(c). As a result of Eq. 5, we can see that u | H n 1 (the pinkish part) looks like the set-valued version • u of the original image u). The span computation is displayed by the gray and olive triangles, respectively for the dimensions d = 1 and d = 0. Eventually, only faces with a dimension less than n can be non-degenerated intervals (they are then displayed in orange).

A 3D version of the discrete interval-valued representation u is depicted in Fig. 3(b); one can easily see that it is a discrete equivalent of the continuous surface of Fig. 1(c). We will define the tree of shapes of such u maps in Sec. 4.1.

Definition of the Naive Dahu Pseudo-Distance

Relying on an interval-valued representation u :

we define the Dahu pseudo-distance by:

Let us just recall that the notation u < -u means that ∀ x, u(x) ∈ u(x), so the argument of the minimum operator is a scalar function u : means that we actually search for a minimal path in the cubical complex D H (not in D), with the classical definition of the minimum barrier distance d MB (that is, a distance defined on scalar functions), and considering all the possible scalar functions u that are "included" in the interval-valued map u. The definition of the Dahu pseudo-distance is thus combinatorial w.r.t. all the scalar images u "included" in u (see the parts highlighted in blue in Eq. 6). The Dahu pseudodistance can be interpreted as the best minimum barrier distance that we can have considering that the input function is continuous.

The Dahu function is a pseudo-distance since it verifies the non-negativity, identity, symmetry, and subadditivity properties. In addition, we have the property:

). An example is given in Fig. 3(a) for an input function u, in Fig. 3(c) for its "continuous" representation u, and in Fig. 3(d) for a scalar function u < -u. A minimal Dahu path between two original points of u is depicted in blue Fig. 3(d); it is actually the same path as the one depicted on the continuous 3D surface in Fig. 1(c).

Solving a Discrete Topology Issue

In 2D discrete images, we can have saddle points. An example is given in Fig. 4(a), from which we deduce the continuous representation (given in Fig. 4(b)) in order to compute some minimal paths. We have â naive u (x a , x a ) = 0 which is obtained with u a (Fig. 4(c)), and â naive u (x b , x b ) = 2 which is obtained with u b (Fig. 4(d)). An interpretation of this situation "in the continuous world" is the following. On one hand, having â naive u (x a , x a ) = 0 implies not only that the level in u a of the saddle point (the 0-face) is 0, but also that there is a 0-level line joining x a and x a in u. On the other hand, since we have u(x b ) = 4 and u(x b ) = 6, u b considers that the level 5 exists between x b and x b , which is in contradiction with the previous conclusion. We thus end up with some inconsistencies, if we consider two distinct discrete paths in u, and if we try to interpret these paths in a continuous way in u. The main issue, here in 2D, is that the representation u cannot prevent lines to cross at different levels, which is the case at saddle points. From this example, we can draw a negative conclusion: The representation u, as defined in Eq. 5, cannot give us a continuous interpretation of the set of all paths in u without topological flaws.

A solution to this problem is to add an intermediate step which computes an "interpolation" of u that has no pathological parts such as the saddle points as in


Lê Duy Huỳnh

Yongchao Xu

Thierry Géraud

Morphological Hierarchical Image Decomposition Based on Laplacian 0-Crossings

Keywords: morphological Laplace operator, well-composed images, tree of shapes, hierarchical decomposition, text detection

A method of text detection in natural images, to be turned into an effective embedded software on a mobile device, shall be both efficient and lightweight. We observed that a simple method based on the morphological Laplace operator can do the trick: we can construct in quasi-linear time a hierarchical image decomposition / simplification based on its 0-crossings, and search for some text in the resulting tree. Yet, for this decomposition to be sound, we need "0-crossings" to be Jordan curves, and to that aim, we rely on some discrete topology tools. Eventually, the hierarchical representation is the morphological tree of shapes of the Laplacian sign (ToSL). Moreover, we provide an algorithm with linear time complexity to compute this representation. We expect that the proposed hierarchical representation can be useful in some applications other than text detection.

Introduction

In [START_REF] Huỳnh | Morphology-based hierarchical representation with application to text segmentation in natural images[END_REF] we have proposed a method to perform robust text segmentation in natural images; see Fig. 1 for an illustration. This method is interesting for several reasons. 1. To select candidate regions for characters, we rely on the morphological Laplace operator. We can observe that this operator is relatively well robust to poor contrast and uneven illumination; the experimental results show that it outperforms the "more classical" component-based methods-namely the Stroke Width Transform (SWT), the Toggle Mapping Morphological Segmentation (TMMS), and the Maximally Stable Extrema Region (MSER) methods-see [START_REF] Huỳnh | Morphology-based hierarchical representation with application to text segmentation in natural images[END_REF] for details. 2. Based on the 0-crossings of the morphological Laplace operator, we compute a hierarchical representation of the image contents, so that regions form an inclusion tree. Such a structure is great because we have a strong simplification of the image contents, and because dealing with a simple tree structure allows for powerful decision taking when grouping regions / characters into text boxes. 3. This text segmentation method is suitable to real-time implementation on mobile devices. Indeed computing the hierarchical representation eventually translates to a simple labeling algorithm that gives both an inclusion tree and a label image; in addition, identifying and grouping text components is achieved by an easy tree-based processing.

In this present paper we focus on two points which are not detailed in [START_REF] Huỳnh | Morphology-based hierarchical representation with application to text segmentation in natural images[END_REF]: how we can rely on some discrete topology tools to get a sound definition of a hierarchical decomposition of an image into regions, and how we take benefit from some properties to get an efficient algorithm to compute such a representation.

Section 2 starts by giving an overview of the hierarchical representation and of the segmentation method. Section 3 recalls some theoretical background, and Section 4 details the algorithm that computes the hierarchical representation. Last we conclude in Section 5.

Method Overview

The method that we propose to represent an image by a hierarchical decomposition and to extract text in natural images is very simple: put shortly, text components are selected among the 0-crossing lines of the morphological Laplace operator of the gray-level input image. The method is composed of several steps, depicted in Fig. 2, and described and justified just below.

We start with a gray-level input image. If the primary data is a color image, we just take the luminance value of its pixels (we lose color information but it almost never negatively affects text retrieval). We then compute its morphological dilation δ and erosion ε to directly deduce two images. First, we obtain the morphological thick gradient ∇ = δ -ε ; it is used later to discard contours that are not enough contrasted. Second, we obtain the morphological Laplace operator ∆ = δ + ε -2 id. The text character boundaries are expected to belong to the 0-crossing contours of this non-linear operator; actually they are, and their localization is precise.

We virtually3 compute a particular interpolated image, ∆ wc , of the Laplacian image ∆ , having 4 times more pixels than the original. This resulting image is well-composed, meaning that the boundaries of every components of any threshold set are Jordan curves. As a consequence, the (boundaries of the) 0-crossings are simple closed curves: they cannot have the shape of a '8'. In addition, they are disjoint, and this set of curves can be organized in an inclusion tree.

Due to the fact that ∆ wc is well-composed, the regions delimited by the 0crossings can be labeled very efficiently (by the classical blob labeling algorithm), and their inclusion tree is built. In addition, many 0-crossings are discarded on the fly during the labeling, because they are not contrasted enough (we use ∇ ), or because they do no satisfy some geometrical criteria (when they are too small for instance). The result "inclusion tree + label image" is the hierarchical decomposition of the image contents into regions; it is depicted on the bottomright part of Fig. 2. We end up with two structures: the inclusion tree, encoded by a parenthood relationship between labels, and an image of labels, so that each pixel is assigned to a region. In addition, we have some additional information related to every regions, that are computed during the labeling: the region area, its bounding box, etc.

From the resulting hierarchical representation, we have derived an application described in [START_REF] Huỳnh | Morphology-based hierarchical representation with application to text segmentation in natural images[END_REF]: text detection. Indeed, we can group regions, that are separated components, together to form text boxes. For that, we only consider the bottom of this tree (the leaves and sometimes their parent): for each component, we search spatially in the label image what are their left and right components to be grouped into a text box. In this step, we highly take advantage of the tree structure: it allows very easily to discard many regions as non-text, and to determine if a leaf region is a character hole or a plain character.

3 Theoretical Background h. The set of all faces, H n , is called the nD space of cubical complexes. Fig. 3(a) depicts a set of faces {h, h , h } ⊂ H 2 where h = {0}×{1}, h = {1}×{0, 1}, and h = {0, 1}×{0, 1}; the dimension of those faces are respectively 0, 1, and 2. The set of faces of H n is arranged onto a grid, called the Khalimsky's grid [START_REF] Khalimsky | Computer graphics and connected topologies on finite ordered sets[END_REF], depicted in orange in Fig. 3(a) (right). Consequently, it means that we just can store values defined on such a grid in a simple matrix structure.

This grid is a way to 1. get a discrete topology, and 2. represent what lies in-between pixels, for instance contours. In Fig. 3(b), an image is depicted with integer values. The green and red regions are open sets corresponding to components having respectively positive and negative values. A 1D contour (thus with 0-faces and 1-faces only) is depicted in blue: it is precisely the discrete boundary of the non-positive region on the right.

Morphological Laplace Operator

In order to detect text in images, many methods first look for candidate regions for characters. A seminal method, based on contour detection, is to consider the 0-crossings of a discrete Laplace operator: given a gray-level image u, the contours of interest are given by ∆u = u xx + u yy = 0. This method is interesting for several reasons: 1. it is a very simple approach; 2. it provides closed contours; 3. labeling the components of the image having the same sign, resp. positive and negative; gives a segmentation; 4) it is self-dual, i.e., it processes dark objects and bright ones the same way.

The morphological Laplace operator has been defined in [START_REF] Van Vliet | An edge detection model based on non-linear laplace filtering[END_REF] by ∆ = (δ -id) -(id -ε ), relying on the elementary dilation (δ) and erosion (ε) morphological operators. Actually, the morphological non-linear version features a much higher fidelity to actual object contours than the linear version; in addition, salient contours hardly depends on the size of the structuring element.

The image depicted in Fig. 4(a) has been created by the authors of [START_REF] Blayvas | Efficient computation of adaptive threshold surfaces for image binarization[END_REF] to make classical binarization methods fail, due to the presence of uneven illumination. On Fig. 4(c), one can observe that the object boundaries belong to the 0-crossings of the morphological Laplace operator.

Tree of Shapes

The tree of shapes is a morphological self-dual representation of an image; see [START_REF] Géraud | A quasi-linear algorithm to compute the tree of shapes of n-D images[END_REF] for history, implementation, and references. This tree encodes the inclusion of the level sets, i.e., the connected components obtained by thresholding. An illustration is given on a very simple image by Fig. 5(a). The node B lists the pixels of the region B. The contour of this node is a level line, and the component, called shape, associated with this node is actually the sub-tree rooted at this node. The node B thus represents the shape given by B ∪ D ∪ E. In [START_REF] Carlinet | A comparative review of component tree computation algorithms[END_REF] the reader can find more details about how to store and process efficiently such a tree structure.

With X representing a 2D discrete space, given an image u : X → Z and any scalar λ ∈ Z, the lower and upper level sets are respectively defined by

[ u < λ ] = { x ∈ X | u(x) < λ } and [ u ≥ λ ] = { x ∈ X | u(x) ≥ λ }.
Using CC, the operator that takes a set and gives its set of connected components, and Sat, the cavity-fill-in operator, we can define:

S(u) = { Sat(Γ ); Γ ∈ CC([u < λ]) ∪ CC([u ≥ λ]) } λ ,
called the tree of shapes of u.

Well-Composed Sets and Images

A sub-class of sets defined on the cubical grid, called well-composed, has been proposed in [START_REF] Latecki | Well-composed sets[END_REF], where all connectivities are equivalent, thus avoiding many topological problems. A very easy characterization of well-composed sets is based on the notion of "critical configurations"; in the 2D case, a set is wellcomposed if the configurations and do not appear. The notion of wellcomposedness has been extended in [START_REF] Latecki | 3D well-composed pictures[END_REF] from sets to functions, i.e., gray-level images, and has been generalized in the nD case in [START_REF] Boutry | On making nD images well-composed by a self-dual local interpolation[END_REF]. A gray-level image u is well

-composed if any set [u ≥ λ] is well-composed. A straightforward character- D E B A C F O A O F B C D E > 2 2 > 2 > 2 > 0 2 1 > 4 > > > (a)
The tree of shapes of an image. max(v,w) . An image is not a priori well-composed. To get a well-composed image from a primary image, one can compute an interpolation of the primary image that is wellcomposed. Fig. 5(b) gives an example of an image which is not well-composed (left), but whose interpolation is well-composed (right). In [START_REF] Boutry | How to make nD functions well-composed in a self-dual way[END_REF] the authors have proposed a simple method to get a self-dual well-composed interpolation, that makes sense when considering the tree of shapes of an image [START_REF] Géraud | Self-duality and discrete topology[END_REF].

Putting Things Together

As explained in Section 2, from a gray-level input image u, we consider the Laplacian image once well-composed; both images are depicted respectively in Figs. 6(a) and 6(b). The tree of shapes of the Laplacian sign image (ToSL), S(sign( ∆ wc (u) )), contains nodes corresponding to positive, negative, and null regions; it is depicted in Fig. 6(d).

Yet the 0-crossings, corresponding to nodes at level 0 (depicted in white), can be "thick", that is, they can contain pixels (2-faces) instead of being only defined by 1D objects (set of 0-faces and 1-faces). In Fig. 6(c) for instance, we can see that some pixels of the 's' contour belong to the 0-crossing region. Since we want to obtain a partition of the image into "positive" and "negative" regions, we merge null regions with their parent regions; that way, instead of the tree of shapes depicted in Fig. 6(d), we consider the simplified one, depicted in Fig. 6(e). Eventually, the actual "0-crossings" we are looking at are the boundaries of the shapes of this final tree, so they are effectively 1D objects. It is illustrated in Fig. 5(b) (and also in Fig. 3(b)): the null region is merged with the negative one, so we consider the blue contour to be the 1D "0-crossing" separating regions having different signs.

In the next section, to compute the hierarchical representation, we do not consider that we have a cubical complex as the space structure: we just ignore that 0-faces and 1-faces exist. Though, from a theoretical point of view, the contours / 0-crossings expressed in terms of 0-faces and 1-faces really are Jordan curves.

Fast Computation of the Hierarchical Representation

A Particular Well-Composed Non-Local Interpolation

In the following, we do not need to consider the cubical complex, so we only deal with pixels; they are, for instance, the valued 2-faces in Fig. 5

(b) (right).

The hierarchical representation is computed from an interpolated image, ∆ wc (u), which is a very particular well-composed version of the Laplacian image ∆ (u). This particular interpolation takes its origin from the work in [START_REF] Géraud | A quasi-linear algorithm to compute the tree of shapes of n-D images[END_REF], and is detailed in [START_REF] Boutry | How to make nD functions well-composed in a self-dual way[END_REF]. Briefly put, it is a non-local interpolation driven by a propagation from the border of the image, which browses the nodes of its tree of shapes from the root to the leaves. The interpolated pixels are assigned with the current gray-level value, which evolves "continuously" during the process. This interpolation has two important features [START_REF] Géraud | Self-duality and discrete topology[END_REF]: it is related to the tree of shapes of the input image, so it actually follows the same scheme as a blob labeling algorithm where a blob would be a level set, and its topological behavior is deterministic.

There are two main consequences: this particular interpolation makes sense in the present context, since we want to label regions that are in an inclusion relationship, and we can optimize the computation of the hierarchical representation (the label image and the inclusion tree) by actually emulating the interpolation.

An example on the image of Laplacian sign given in Fig. 7(a) is depicted in Fig. 7(b) 4 . We can observe in Fig. 7(b) that we obtain the desired properties: the boundaries of the regions are Jordan curves, and the regions are in an unambiguous spatial inclusion relationship. The inclusion tree of the Laplacian sign image, called ToSL, is thus a hierarchical image decomposition.

Labeling Interpolated Sign of ∆ to Construct ToSL

Thanks to the fact that ∆ wc is well-composed, the regions delimited by the 0crossings can be labeled very efficiently (by the classical blob labeling algorithm), and their inclusion tree is built. This resulting inclusion tree is the tree of shapes of the sign of the Laplacian, a ternary-valued image-with pixels set to -1 (red), 0 (gray), or 1 (green) as depicted in Fig. 7. The whole labeling process is depicted by the black part in Algorithm 1. It computes a label image L that labels all components delimited by the 0-crossings of ∆ and a tree structure encoded in an array of parenthood parent (e.g., having parent(l 1 ) = l 2 means that the region with label l 1 is included in the region with label l 2 ). Q is a queue of pixels, border is an auxiliary image that marks the inner component contours as active q ← Q.pop(), N ← N4 24:

if border (q) = undef then 25:

N ← N8 /*optimized version*/ 26: for all n ∈ N(q) do 27:

if

L(n) = 0 and ∆ (p) × ∆ (n) ≥ 0 then 28: L(n) ← , Q.push(n) 29: else 30: border (n) ← a
follow contour (p) 31: P.init(), border (p) ← ā, Q.push(p) 32: while Q is not empty do 33:

q ← Q.pop(), P.update(q) 34:

for all n ∈ N(q) do /*N8 or N4*/ 35:

if border (n) = a then 36:

Q.push(n), border (n) ← ā 37: return P =0 a or inactive ā state, nlabels is the temporary number of labels, is the current label value, and N represents either 4-connectivity (N 4 ) or 8-connectivity (N 8 ).

The core of the algorithm is an alternate application of the following two routines (depicted on the right side of Algorithm 1). blob labeling is a classical queue-based "blob labeling" algorithm that labels the underlying connected component with current label value , and also marks the inner component contours as active state a. Note that actually we do not want regions representing null values in the final tree, so we merge nodes corresponding to 0-crossings with their parent (see line 27 and Fig. 6). follow contour is also a queue-based process which is similar with previous "blob labeling" algorithm, but applied on the underlying active contour of an unlabeled connected component. Instead of labeling the active contour, this routine browses it and marks it as inactive ā (see line 31 and 36), and collects a measure P (e.g., length) characterizing the contour. Note that both these two routines are very efficient thanks to the queue-based "blob labeling".

More precisely, for the complete algorithm depicted on the left side of Algorithm 1, we browse the pixels in raster scan order (main loop, line 5). when we reach an unlabeled pixel p, if this is the first pixel p 0 in the scan order (i.e., the top left pixel), we know that its label value is 1, and the label value of its parent (i.e., itself) is also 1 (line 8 to 10), because p 0 is in the root node thanks to the added external boundary described in previous section. For all other unlabeled pixel p = p 0 , we follow the contour of the unlabeled region, which is a hole in the label image, thanks to the border image. The follow contour routine computes on the fly a contour-based measure P characterizing the hole, such as the average of gradient's magnitude along the contour and the bounding box of the hole. If this contour-based measure P does not satisfy some given criterion, for instance, it is not contrasted enough, or it does not satisfy some geometrical criteria (when the hole is too small for instance), we discard this region by not creating a new label value for this unlabeled region. Let p -1 be the pixel just before p in the raster scan order (p -1 is guaranteed to be labeled), we assign the label value of p -1 to the underlying hole region, which means we merge it with its parent region. If the contour measure P satisfies the corresponding criterion, we create a new label value to label the hole region, and update the parenthood relationship of this new label value to L(p -1 ) (see line 13 to line 18). Then we proceed the routine blob labeling to label the connected set of pixels having the same Laplacian sign as p or being null and update the auxiliary border image. Note that since we use the 4-connectivity neighborhood (N 4 ) in the routine blob labeling to label regions and mark neighboring pixels having different sign of ∆ , we need to use the 8-connectivity neighborhood (N 8 ) to follow completely the active contour of an unlabeled region (see also the longest contour in the left image in Fig. 7(e)).

An example of the proposed algorithm on the interpolated image in Fig. 7(b) is depicted in Fig. 7(c-e). Different stages of the algorithm are depicted in Fig. 7(e). Note that the pixels having null Laplacian sign are grouped with the parent region. The two small regions inside cyan and respectively yellow region are also discarded, they are grouped with the parent region. The resulting "label image + tree" are depicted in Fig. 7(c) and 7(d).

Optimization of ToSL Construction

The well-composed interpolation described in Section 4.1 resolves all the topological issues at critical configurations with the cost of quadrupling the number of pixels. Yet, in practice, we do not need to apply this interpolation, which can be emulated efficiently without subdivising the image domain thanks to the particular non-local way of interpolation described in Section 4.1.

The optimized version of the proposed algorithm by emulating the interpolation is depicted in Algorithm 1 by the black and red parts. More precisely, the used particular non-local interpolation is based on the inclusion relationship of components. The interpolated values at critical configurations are given by the values of outer region. Besides, the proposed algorithm for such interpolated image labels regions from outside to inside. Consequently, it is equivalent to use 8-connectivity (N 8 ) when we use blob labeling algorithm to label pixels that are not yet activated thanks to the border image (see line 24 and line 25 in the scope of blob labeling on the right side of Algorithm 1). This makes the other two pixels of critical configurations having a different sign of Laplacian disjoint and marked as active contours. Consequently, when we proceed the follow contour routine in the following, we need to use 4-connectivity (N 4 ) (see line 34 in the scope of follow contour in Algorithm 1) to avoid connecting two disjoint unlabeled regions. An example is given in Fig. 7(f) which shows different stages of the proposed optimized algorithm without using interpolation. Note that the green pixels (resp. red pixels) of the critical configuration in the top (resp. bottom) red rectangle in Fig. 7(a) are considered as connected using 8-connectivity (N 8 ) when we proceed the routine blob labeling. The active contours (black pixels) in Fig. 7(f) are processed using 4-connectivity (N 4 ). Both the algorithm depicted in Fig. 7(e) and Fig. 7(f) result in the same "label image + tree" depicted in Fig. 7(c) and Fig. 7(d).

Application to Text Segmentation and Detection

As an application, the ToSL has been used for text detection and segmentation; the results have been further detailed in [START_REF] Huỳnh | Morphology-based hierarchical representation with application to text segmentation in natural images[END_REF]. Briefly put, we compute the ToSL, and then group nodes together to form candidates for text boxes.

During the follow contour process, we discard some regions w.r.t. their contour average gradient magnitude ∇ (to keep only contrasted regions), the contour length (to remove noise, i.e., too small regions), and the bounding box size B (to keep components that reasonably look like characters).

The ToSL is such that sibling nodes are objects sharing the same background in the original image. Consequently, the grouping process is performed efficiently because the tree structure helps to limit the sets of candidates that can be grouped together. We found that most of text components are leaves, and sometimes leaf parents (in the case of characters with holes). To group text components together into text boxes, starting from each leaf, we look in the image space, that is in the label image, for a sibling, and we control the grouping thanks to some geometric information.

Some results are depicted in Fig. 8; more details and a quantitative evaluation are given in [START_REF] Huỳnh | Morphology-based hierarchical representation with application to text segmentation in natural images[END_REF].

In this paper, we have presented a hierarchical decomposition of the image contents into regions. Although we use the very "classical" idea of relying on the 0-crossings of the Laplace operator image to obtain the objects of interest, our version is innovative for several reasons. First, we rely on the morphological Laplace operator, which performs well in the case of uneven illumination, which is a common defect in natural images. Second, we ensure that the "0-crossings" are really 1D objects. For that, we use a well-composed Laplacian image, we compute the tree of shapes of its sign, and we consider the 1D topological boundary of shapes. As a consequence, the positive and negative regions can be organized into a tree without topological ambiguity. Last we present a linear time complexity algorithm, which is a hardly more sophisticated blob labeling algorithm, to compute the hierarchical structure. We also provide a way, directly during the computation process, to ignore some regions if they are not relevant, and an optimization that mimics well-composedness and then avoids to duplicate pixels.

Fig. 1 .

 1 Fig. 1. A hierarchical image decomposition (center) leading to text detection (right).

Fig. 2 .

 2 Fig. 2. Overview of the proposed method to get a hierarchical image decomposition.

  Two regions and a separating contour.

Fig. 3 .

 3 Fig. 3. Representation of images on the Khalimsky's grid.

Fig. 4 .

 4 Fig. 4. Contour characterization by morphological operators; 4(b) is the morphological Laplacian ∆B = δB +εB -2id; 4(c) depicts on ∆B = 0 the gradient values ∇B = δB -εB (inverted) showing that the contours of actual objects are effectively salient.

  Level sets contours as Jordan curves.

Fig. 5 .

 5 Fig. 5. Topological representations and some topological considerations.

  (a) u. (b) ∆ wc (u).(c) ∆ wc (u) = 0. S(sign( ∆ wc (u) )).

  Final simplified S, called ToSL.

Fig. 6 .

 6 Fig.6. Tree of shapes of Laplacian sign (ToSL): positive and negative regions are respectively green and red nodes of the ToS, and null regions are white nodes.

  Different stages of the proposed algorithm based on interpolation.(f) Different stages of the proposed fast algorithm without using interpolation.

Fig. 7 .

 7 Fig. 7. An example of the proposed labeling algorithm. Black pixels: contours of components inside labeled ones. White pixels: pixels that are not yet labeled and are also not marked as inside contours at the current stage of the labeling process.

Fig. 8 .

 8 Fig. 8. Qualitative results using "ICDAR 2015 Robust Reading" database: input (top), label image (middle), text boxes (bottom).

  Algorithm 1 Computation of ToSL by labeling interpolated image (black part) and its fast version (adding red part) without using interpolation.

		labeling (∆ ,∇ )	blob labeling (p, )
	1: for all p do	21: L(p) ← , Q.push(p)
	2:	L(p) ← 0	22: while Q is not empty do
	3:	border (p) ← undef	23:
	4: nlabels ← 1
	5: for all p do
	6:	if L(p) = 0 then
	7:	continue
	8:	if p = p0 then
	9:	← 1
	10:	parent[ ] = 1
	11:	else
	12:	P ← follow contour(p)
	13:	if evaluate(P) then
	14:	nlabels ← nlabels + 1
	15:	← nlabels
	16:	parent[ ] ← L(p-1)
	17:	else
	18:	← L(p-1)
	19:	blob labeling (p, )
	20: return parent, L

† This work has been conducted in the context of the mobidem project, part of the "Systematic Paris-Region" and "Images & Network" Clusters (France). This project is partially funded by the French Gov. and its economic development agencies.

We will see later that we actually do not interpolate the Laplacian image, but proceed as if there were an interpolation. Practically, it means that we avoid the need of multiplying by

the number of pixels in the process.

Note that the two identical local configurations enclosed by red rectangles in Fig.7(a)do not lead to the same interpolation; this is due to the non-local interpolation process that depends of the outer region, which is different in the two cases: respectively negative for the top configuration, and positive for the bottom one.