Paul Desombre

Haïfa Farès

Yves Louët

Performance Comparison of Digital Modulations in the Presence of Gaussian Phase Noise in the Sub-THz Context

Keywords: Sub-THz communications, Phase Noise, Waveform

Driven by the need of more available frequency bands, an increasing attention has been paid to the use of sub-THz bands. In such high frequency bands, the signal is expected to be impaired by phase noise, I/Q mismatch, power amplifier non linearities, etc. It has been shown that at such frequencies, with possibly large bandwidth (∼ 50 GHz), Gaussian phase noise models accurately the phase noise. This paper studies the effect of Gaussian phase noise on several waveforms. Here, five performance metrics are used to compare them: Bit Error Rate (BER), Peak-to-Average Power Ratio (PAPR), Adjacent Channel Power Ratio (ACPR), spectral efficiency (SE) and Error Vector Magnitude (EVM). The aim of this comparison is to give recommandations on which waveform to use in this context, depending on the impact of the phase noise on the performance of such transmissions.

I. INTRODUCTION

Considering the latest report of the European Communications Office (ECO) Frequency Information System (FIS) [START_REF]European Communications Office Frequency Information System[END_REF], frequencies over 71GHz have yet to be allocated to wideband data communications. In the last decades, the research community payed more attention to the sub-THz bands (200 GHz -1 THz) [START_REF] Nagatsuma | Advances in terahertz communications accelerated by photonics technologies[END_REF]. However, with higher frequencies comes several impairments, such as the Phase Noise (PN), which is the subject of this paper.

Phase noise is the main distorsion experienced by the signal in the sub-THz bands, which makes its recovery impossible. Therefore, its impact has been widely studied [START_REF] Bicaïs | Phase Noise Model Selection for Sub-THz Communications[END_REF]- [START_REF] Hillger | Toward Mobile Integrated Electronic Systems at THz Frequencies[END_REF]. Indeed, in [START_REF] Bicaïs | Phase Noise Model Selection for Sub-THz Communications[END_REF], Bicais et al. compared two models of PN: correlated Wiener plus Gaussian PN, and uncorrelated Gaussian PN (GPN). They showed that, for broadband communications, such as sub-THz communications, uncorrelated GPN is an accurate model. An other study of PN [START_REF] Khanzadi | Calculation of the Performance of Communication Systems from Measured Oscillator Phase Noise[END_REF] gives a model for PN simulations, with a thourough statistical study with the goal of being more accurate than other models, at the cost of complexity. The authors from [START_REF] Bicaïs | On the Optimum Demodulation in the Presence of Gaussian Phase Noise[END_REF], [START_REF] Bicaïs | Design of Digital Communications for Strong Phase Noise Channels[END_REF] gave a design of digital communications for strong PN channels: their contributions showed that, by using polar metric in the detector, the performance in terms of Bit Error Rate (BER) is improved. An other digital modulation framework has been proposed in [START_REF] Bicaïs | An Algebraic Framework for Digital Envelope Modulation[END_REF], where PN robustness has been considered as the performance metric to optimize.

In the review of the THz communications advances drawn in [START_REF] Nagatsuma | Advances in terahertz communications accelerated by photonics technologies[END_REF], the author explained that publications have been numerous since early 2000. Communication systems using all-electronics (such as [START_REF] Hillger | Toward Mobile Integrated Electronic Systems at THz Frequencies[END_REF]) or photonics are using carrier frequencies in the range 100 GHz -700 GHz, though the most studied frequency bands lies between 200 GHz and 400 GHz. From this review, researchers have reached over 100 Gbps. These systems used OOK, QPSK and 16-QAM modulation schemes.

The aim of this paper is to compare multiple waveforms in the presence of GPN using the following five performance metrics: Bit Error Rate (BER), Peak to Average Power Ratio (PAPR), Adjacent Channel Power Ratio (ACPR), Spectral Efficiency (SE) and Error Vector Magnitude (EVM). Some well-known Phase-Amplitude Modulations (PAM) modulations were used: BPSK, QPSK, π/4-DQPSK, 8-PSK, 16-APSK, 16-QAM, 256-QAM and OOK. The OOK scheme was taken from the framework presented in [START_REF] Bicaïs | An Algebraic Framework for Digital Envelope Modulation[END_REF]. Alongside widely used modulation schemes, Continuous Phase Modulation (CPM) have been considered into the comparison study due to their high energy efficiency. Another promising modulation scheme was used: Parallel Sequence Spread Spectrum (PSSS) [START_REF] Krishnegowda | Towards 100 Gbps Wireless Communication in THz Band with PSSS Modulation: A Promising Hardware in the Loop Experiment[END_REF] combined with BPSK and QPSK modulations. The goal of this comparison is to show the performance of different modulation schemes impaired by PN, with PN modelled as GPN. To the best of our knowledge, no comparison of such scale has been published so far.

The remainder of this paper is organized as follows. Section II briefly describes the channel and GPN models. Section III provides simulations scheme and results. Section IV illustrates the main outcomes of this comparison study. Finally, conclusions are drawn in Section V.

II. SYSTEM MODEL

In this work, we considered a SISO (Single-Input, Single-Output) system with an Additive White Gaussian Noise (AWGN) channel including PN at both sides (transmitter and receiver). The received signal is defined as

r = e jφpn,r (s • e jφpn,t + n), (1)
with r the received signal after the receiver's oscillator, n the circular, zero-mean, complex white Gaussian noise with variance σ 2 n , s the transmitted signal, and Φ pn,t and Φ pn,r the PN at the transmitter and the receiver side respectively. The transmitted signal is composed of symbols which take their values from the aforementioned modulations. Each of the two instances of PN follow an independant zero-mean Gaussian process φ pn ∼ N (0, σ 2 pn ). In our simulations, we used several modulation schemes: BPSK, QPSK, π/4-DQPSK, 8-PSK, 16-APSK, 16-QAM, 256-QAM, OOK, GMSK, PSSS BP SK and PSSS QP SK . To compute the PSSS, we need a spreading sequence (maximum length sequence), we used the first 2 8 -length sequence taken from [START_REF] Davies | System Identification for self-adaptive control[END_REF]. All of our simulations were baseband, filtered by a Square-Root-Raised-Cosine (SRRC) filter with a roll-off factor of α = 0.25, except for the OOK modulation (Non-Return-to-Zero filter, NRZ) and GMSK (Gaussian filter, with parameters L= 4 and BT= 0.25). The variance of the GPN was: σ 2 pn = 10 -4 for low PN, σ 2 pn = 10 -2 for medium PN and σ 2 pn = 10 -1 for high PN. No Forward Error Correction (FEC) nor equalization algorithms were used in our simulation setup.

The BER simulations used Monte-Carlo (MC) scheme, with E b /N 0 = [0, 1, ..., 30] dB, and frames of 2 15 bits. To process the PAPR, we calculated the Complementary Cumulative Distribution Function (CCDF). The ACPR was obtained by looking at the normalized frequency of the first channel offset. In a normalized frequency setup, with the first channel offset taken at the normalized frequency value of F of f set = 0.16. Spectral efficiency was investigated for a BER of 10 -5 . And EVM was processed at E b /N 0 = [0, 5, 10, 15, 30] dB. As EVM for the GMSK is not trivial to obtain, we chose to calculate its Mean Quadratic Error (MQE), which is given by

MQE = 1 N N i=1 (φ s,t (i) -φ s,r (i)) 2 , (2)
with Φ s,t and Φ s,r are the phases for the ith sample at the transmitter and the receiver side, respectively, for the i th sample. The EVM for the other modulations were obtained thanks to the MATLAB comm.EVM object, with reference signal normalization.

In the following section, results are illustrated for medium and high level of PN for which the impact on performance metrics could be observed.

III. SIMULATIONS

A. Bit error rate

Figure 1 shows the BER as a function of the E b /N 0 (dB) for the simulated systems. All of the 11 modulations are shown, without and with PN (medium level) to enable comparison. As expected, medium PN renders 256-QAM unrecoverable. As highlighted in [START_REF] Bicaïs | On the Optimum Demodulation in the Presence of Gaussian Phase Noise[END_REF], the effect of PN on the signal's BER is that it reaches an error floor. With medium PN, the error floor is reached at 10 -2 for 256-QAM and a little over 10 -8 for 16-APSK. Low-order modulation waveforms have their respective error floor lower than 10 -8 .

In a system with high PN, the BER curves are shown on Figure 2. Most of the modulations exhibit an error floor. At Regarding only the BER performance, the modulations we should use in PN-impaired systems such as sub-THz communications should be low-order modulations. The PSSSmodulated BPSK and QPSK are remarkably robust to PN.

B. Peak to Average Power Ratio

If we look at the PAPRs represented by their CCDFs in Figure 3, the two PSSS are of the same order as Multi-Carrier waveforms. As shown in [START_REF] Dore | Above-90 GHz Spectrum and Single-Carrier Waveform as Enablers for Efficient Tbit/s Wireless Communications[END_REF], in sub-THz systems, it is more preferable to use Single-Carrier (SC) systems. As expected, the GMSK is not illustrated in the CCDF Figure , as it has a null-PAPR. Regarding only the PAPR of the waveforms, the PSSS should be excluded of the possible waveforms for sub-THz communications, whereas the other simulated waveforms could all be used.

C. Adjacent Channel Power Ratio

As the ACPR is processed on the Power Spectral Density (PSD) of the waveforms, it is expected that all the simulated waveforms using the same SRRC filter have exactly the same performance (Fig. 4). As the Gaussian filter used in the GMSK has slightly higher rejection than the SRRC in the vicinity of the main lobe, it is also not surprizing that the GMSK has better performance. Furthermore, as for the OOK, the use of an NRZ filter is causing spikes in its PSD, this explains its higher ACPR.

It should be noted that the PSD of a White Gaussian noise is constant. As the level of the PN PSD is higher than the level of the signal PSD in the rejection band, it explains the high increase of ACPR with the increase of PN.

Looking only at the ACPR, GMSK has the best behavior for these parameters. It is not excluded that, with a better filter, the other modulations would have better performance.

D. Spectral efficiency

Without PN (Fig. 5), SE is well known and should serve as a comparator for PN-impaired waveforms.

As seen earlier, the 256-QAM is excluded from the list of considered waveforms. Its BER does not reach the 10 -5 threshold within the considered E b /N 0 interval for medium PN (Fig. 6) and high PN (Fig. 7).

As for the other waveforms, with medium PN (Fig. 6), the E b /N 0 level is higher than without PN, appart from the lower order modulations. This highlights the robustness to PN of low-order modulations. As expected, with higher PN (Fig. 7), the same trend is observed, except for QPSK, which has as an E b /N 0 of 15.33 dB for a BER of 10 -5 .

Then, if we tune on for an E b /N 0 threshold of 16 dB, which is the upper bound of E b /N 0 for high PN, only the lower order modulations should be considered for a highly PN-impaired system.

E. Error vector magnitude

At high E b /N 0 , the EVM is expected to be near 0%, which can be observed in Figure 8. Without PN, for high order modulations, the EVM is indeed near 0%. With low order modulations, the EVM is higher. The higher the E b /N 0 , the higher the increase of EVM in systems impaired by PN.

It can be said the same about MQE for GMSK modulation. However, the level of MQE is significantly lower for GMSK than the level of EVM for any low order modulation. This makes GMSK a good candidate for sub-THz communications.

Considering the EVM and MQE, GMSK is far ahead of the other waveforms.

IV. RESULTS AND DISCUSSION Now, instead of considering each criterion alone we consider them all together. To compare the waveforms, we graded objectively each of them on a scale from 0 to 10, with 0 the worst candidate for sub-THz communications and 10 the best candidate (Fig. 10).

Given the results of the grading, in a system with little to no PN, depending on the channel conditions and other impairments, any modulation could be used. However, in a system impaired with medium to high PN, we advise to use waveforms with a low-order modulation, such as the QPSK or CPM such as GMSK.

QPSK has a relatively good robustness to PN, which is surely one of the reasons it is used in existing systems in the literature [START_REF] Nagatsuma | Advances in terahertz communications accelerated by photonics technologies[END_REF]. QPSK also has objectively good EVM and PAPR performance. The ACPR depends on the shaping filter, which can be optimized. SE is low compared to higher order modulations, but with large bandwidth, even low SE can enable high throughput.

On the other hand, GMSK has better performance compared to the QPSK in presence of high PN. It also has the best PAPR possible (since it has a constant complex envelope). Its MQE is excellent, and it has the best ACPR among the studied waveforms. Nevertheless, its complexity depends on the number of states in the Viterbi algorithm needed to perform Maximum Likelihood detection. Therefore, its complexity is obviously higher than that of the QPSK scheme, based on a simple separator. And yet this GMSK remains of great interest thanks to its high energy efficiency. Indeed, Since GMSK has a constant complex envelope, it would enable the use of Power Amplifiers (PA) at saturation, which is not the case with the other waveforms.

V. CONCLUSION In this paper, we presented simulated sub-THz communications impaired by phase noise. We assumed that phase noise is modelled as a Gaussian process. We provided simulations results and compared 11 waveforms using five performance metrics (BER, PAPR, ACPR, SE and EVM). Eventually we discussed the results and then we recommanded to use GMSK or QPSK waveforms in medium to high phase noise impaired system.

Fig. 1 .

 1 Fig. 1. BER curves for no PN and medium PN.

Fig. 2 .

 2 Fig. 2. BER curves for no PN and high PN.

Fig. 3 .

 3 Fig. 3. PAPR CCDF of the simulated waveforms.

Fig. 4 .

 4 Fig. 4. ACPR of the simulated waveforms.

Fig. 5 .

 5 Fig. 5. SE performance without PN and BER = 10 -5 .

Fig. 6 .

 6 Fig. 6. SE performance with medium PN and BER = 10 -5 .

Fig. 8 .

 8 Fig. 8. EVM for E b /No = 15 dB.

Fig. 9 .

 9 Fig. 9. MQE for GMSK modulation.

Fig. 10 .

 10 Fig. 10. Comparison of each waveform with (a) no phase noise, (b) medium phase noise, (c) high phase noise.

ACKNOWLEDGMENT

The research leading to these results received funding from the French LabEx CominLabs and French administrative region Bretagne within the frame of the project WASSAP.


Olivier Lezoray
email: olivier.lezoray@unicaen.fr

Vinh Thong
email: vinhthong.ta@unicaen.fr

Abderrahim Elmoataz
email: abderrahim.elmoataz-billah@unicaen.fr

Nonlocal

Nonlocal graph regularization for image colorization

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Colorization is the process of adding color to monochrome images. It is usually made by hand by color experts but this process is tedious and very timeconsuming. In recent years, several methods have been proposed for colorization [START_REF] Levin | Colorization using optimization[END_REF][START_REF] Yatziv | Fast image and video colorization using chrominance blending[END_REF] that less require intensive manual efforts. These techniques colorize the image based on the user's input color scribbles and are mainly based on a diffusion process. However, most of these diffusion processes only use local pixel interactions that cannot properly describe complex structures expressed by nonlocal interactions. In this paper, we propose to use nonlocal graph regularization for image colorization. Regularization by variational methods has shown its effectiveness for many applications. Since the advent of the nonlocal means filters [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF][START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF], the use of nonlocal interaction to capture complex structures of the data has received a lot of attention. The nonlocal approach has also shown to be very effective and more flexible in the regularization process. We have recently proposed a nonlocal discrete p-Laplacian regularization framework for the processing of images and manifolds represented by weighted graphs of the arbitrary topologies [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing[END_REF]. The proposed methodology enables local and nonlocal regularization by using appropriated graphs topologies and edge weights. In this work, we consider the image colorization problem within this nonlocal graph-based framework.

Nonlocal graph regularization

In this Section, basic definitions on graphs are recalled, nonlocal discrete differential operators and a family of weighted p-Laplace operators are introduced. The p-Laplace regularization on weighted graphs problem is presented and the associated filter family.

Preliminary definitions

A weighted graph G = (V, E, w) consists in a finite set V = {v 1 , . . . , v N } of N vertices and a finite set E ⊂ V × V of weighted edges. We assume G to be undirected, with no self-loops and no multiple edges. Let (u, v) be the edge of E that connects vertices u and v of V . Its weight, denoted by w(u, v), represents the similarity between its vertices. Similarities are usually computed by using a positive symmetric function w :

V × V → R + satisfying w(u, v) = 0 if (u, v)
is not an edge of E. The notation u ∼ v is also used to denote two adjacent vertices. Let H(V ) be the Hilbert space of real-valued functions defined on the vertices of a graph. A function f :

V → R of H(V ) assigns a real value f (v) to each vertex v ∈ V . The function space H(V ) is endowed with the usual inner product f, h H(V ) = v∈V f (v)h(v)
, where f, h : V → R. Similarly, one can define H(E), the space of real-valued functions defined on the edges of G.

Weighted gradient and divergence operators

Let G = (V, E, w) be a weighted graph, and let f :

V → R be a function of H(V ). The difference operator of f , noted d : H(V ) → H(E), is defined on an edge (u, v) ∈ E by: (df )(u, v) = w(u, v)(f (v) -f (u)).
The directional derivative (or edge derivative) of f , at a vertex u ∈ V , along an edge e = (u, v), is defined as:

∂ v f (u) = (df )(u, v).
The weighted gradient operator of a function f ∈ H(V ), at a vertex u ∈ V , is the vector operator defined by:

∇ w f (u) = [∂ v f (u) : v ∼ u] T .
The L 2 -norm of this vector represents the local variation of the function f at a vertex of the graph. It is de-

fined by |∇ w f (u)| = v∼u w(u, v)(f (v) -f (u)) 2 .
The local variation is a semi-norm that measures the regularity of a function around a vertex of the graph. The adjoint of the difference operator, noted d * :

H(E) → H(V ), is a linear operator defined by df, H H(E) = f, d * H H(V ) , for all f ∈ H(V ) and all H ∈ H(E).
Using the definitions of the inner products in H(V ) and H(E), and of the difference operator, we obtain the expression of

d * at vertex u by (d * H)(u) = v∼u w(u, v)(H(v, u) -H(u, v)).
The divergence operator, defined by -d * , measures the net outflow of a function of H(E), at each vertex of the graph.

A family of weighted p-Laplace operators

Let p ∈ (0, +∞) be a real number. The weighted p-Laplace operator of a function f ∈ H(V ), noted

∆ p w : H(V ) → H(V ), is defined by ∆ p w f (u) = 1 2 d * (|∇ w f (u)| p-2 df (u, v)).
The p-Laplace operator of f ∈ H(V ), at a vertex u ∈ V , can be computed by:

∆ p w f (u) = 1 2 v∼u γ f w (u, v)(f (u) -f (v)), (1)
with

γ f w (u, v) = w(u, v)(| w f (v)| p-2 + | w f (u)| p-2
). The p-Laplace operator is nonlinear, with the exception of p = 2. When p = 2, it corresponds to the combinatorial graph Laplacian. When p = 1, it corresponds to the weighted curvature of the function f on the graph.

p-Laplace regularization on weighted graphs

In this Section, one considers a general function f 0 : V ⊂ R n → R m defined on graphs of the arbitrary topologies and we want to regularize this function. The regularization of such a function corresponds to an optimization problem which can be formalized by the minimization of a weighted sum of two energy terms:

min f ∈H(V ) E p w (f, f 0 , λ) = R p w (f ) + λ 2 f -f 0 2 .
(2) with R p w (f ) = 1 p u∈V |∇ w f (u)| p a regularization functional. The first term in (2) is the smoothness term or regularizer, meanwhile the second is the fitting term. The parameter λ ≥ 0 is a fidelity parameter, called the Lagrange multiplier, which specifies the trade-off between the two competing terms. Both terms of the energy E p w are strictly convex functions of f . In particular, by standard arguments in convex analysis, the problem (2) has a unique solution, for p = 1 and p = 2 , which satisfies

∂E p w (f,f 0 ,λ) ∂f (u)
= 0, ∀u ∈ V . This is rewritten as

∂R p w (f ) ∂f (u) + λ(f (u) -f 0 (u)) = 0, ∀u ∈ V, (3)
that is equivalent to:

2∆ p w f (u) + λ(f (u) -f 0 (u)) = 0, ∀u ∈ V, (4)
since

∂R p w (f ) ∂f (u) = 2∆ p w f (u), ∀u ∈ V .
Substituting the expression of the p-Laplace operator in (3), we obtain:

λ + v∼u γ f w (u, v) f (u)- v∼u γ f w (u, v)f (v) = λf 0 (u).
(5) We propose to use the linearized Gauss-Jacobi iterative method to solve the system (5). Let t be an iteration step, and let f (t) be the solution of (5) at the step t. Then, the method is given by the following algorithm ∀u ∈ V :

       f (0) (u) = f 0 (u) f (t+1) (u) = λf 0 (u) + v∼u γ f (t) w (u, v)f (t) (v) λ + v∼u γ f (t)
w (u, v) .

(6) It describes a family of discrete diffusion processes, which is parameterized by the structure of the graph (topology and weight function), the parameter p, and the parameter λ. Also, the stopping time can be given a priori, or can be determined by a stopping criterion. To get the convergence of the process, a classical stopping criterion is f (t+1) -f (t) < τ , where τ → 0 is a small fixed constant. At each iteration of the algorithm [START_REF] Yatziv | Fast image and video colorization using chrominance blending[END_REF], the new value f (t+1) (u) depends on two quantities: the original value f 0 (u), and a weighted average of the filtered values of f (t) in a neighborhood of u. This shows that the proposed filter, obtained by iterating [START_REF] Yatziv | Fast image and video colorization using chrominance blending[END_REF], is a low-pass filter which can be adapted to many graph structures and weight functions.

Construction of graphs

The minimization problem (2), and the discrete diffusion processes [START_REF] Yatziv | Fast image and video colorization using chrominance blending[END_REF], can be used to regularize any function defined on a finite set V of discrete data. This is achieved by constructing a weighted graph G = (V, E, w), and by considering the function to be regularized as a function f 0 : V ⊂ R n → R m , defined on the vertices of G. A typical graph in image processing is the 8-adjacency graph where vertices are associated to pixels and edges correspond to pixel adjacency relationships: this corresponds to a classical local processing. By changing the graph topology and the edge weights, we naturally obtain an expression of nonlocal processing methods which is embedded in the graph structure. Indeed, once edges are added between vertices, they are considered as direct neighbors and the processing is local over the graph. Figure 1 resumes all these elements. A classical 8-adjacency grid graph is considered and expresses only local interactions through local edges. An example of nonlocal interaction is depicted for one vertex by adding nonlocal edges between vertices not spatially connected (for a 5 × 5 neighborhood in Figure 1). Let F (f 0 , v) ∈ R q denote a feature vector associated to each vertex v ∈ V . This feature vector associated to vertices can be the initial function value:

F (f 0 , v) = f 0 (v) (for a local configura- tion) or a vector F (f 0 , v) = [f 0 (u) : u ∈ B v,s ] T
(for a nonlocal configuration). For this latter case, F (f 0 , v) is a patch where B v,s denotes a bounding box of size (2s + 1) × (2s + 1) centered at v (s = 1 in Figure 1). Weights are computed according to a measure of similarity w(u, v) = exp -F (f 0 ,u)-F (f 0 ,v) 2 σ 2 . To have a parameterless weight function, σ has to be computed locally at each vertex of the graph [START_REF] Lezoray | Parameterless discrete regularization on graphs for color image filtering[END_REF].

A local edge A nonlocal edge F (f 0 , v) = f 0 (v) F (f 0 , v) = [f 0 (u) : u ∈ Bv,s] T Vertices to compare

Application to image colorization

We now explain how to perform image colorization with the proposed framework. From a gray level image f l : V → R, a user provides an image of color scribbles f s : V s ⊂ V → R 3 that defines a mapping from vertices to a vector of RGB color channels:

f s (v) = [f s 1 (v), f s 2 (v), f s 3 (v)] T where f s i : V → R is the i-th component of f s (v).
From these functions, one computes f c : V → R 3 that defines a mapping from the vertices to a vector of chrominances:

     f c (v) = f s 1 (v) f l (v) , f s 2 (v) f l (v) , f s 3 (v) f l (v) T , ∀v ∈ V s . f c (v) = [0, 0, 0] T , ∀v / ∈ V s . (7)
We then consider the regularization of function f c (v) by applying algorithm [START_REF] Yatziv | Fast image and video colorization using chrominance blending[END_REF] where γ f image:

w(u, v) = exp -F (f l ,u)-F (f l ,v) 2 σ 2 l , with σ l estimated from f l . At convergence, final colors are ob- tained by f l (v) × f c (t) 1 (v), f c (t) 2 (v), f c (t) 3 (v)
T , ∀v ∈ V . For all the experiments λ = 0.01, p = 1 and the data-fitting term is measured only on vertices v ∈ V s . This is needed to ensure that the original scribbles do not change too much of color. Figures 2 and3 present colorization results on two gray level images obtained from color images given a set of color scribbles as input. Local (8-adjacency grid graph with F (f l , v) = f l (v)) and nonlocal (99-adjacency grid graph that corresponds to add nonlocal edges in a 11 × 11 neighborhood with F (f l , v) as a 5 × 5 patch) are presented. The benefits of nonlocal processing are evident for the boy (Figure 2): the eyes and several areas of the bib are not properly colored and have diffused over straight edges. On the opposite, nonlocal colorization has successfully colored these areas thanks to its ability to discover similar textures and fine details. For the girl (Figure 3), even with precise locations of color scribbles, local results are less good than nonlocal ones: too much diffusion around the hairs and the shirt (see the zoomed areas in Figure 3). To have a quantitative evaluation between local and nonlocal colorization, PSNR measures between the original color image (before its conversion into a gray level one) and the colored one are provided and show once again the benefits of nonlocal colorization.

Conclusion

In this paper, we first presented a framework that enables local or nonlocal regularization by using appropriated graph topologies and edge weights. We have shown how the proposed framework can be used to perform nonlocal image colorization. Moreover, results were provided that demonstrate the benefits of nonlocal image colorization compared to local image colorization.

Figure 2 .

 2 Figure 2. Colorization results in local or nonlocal configurations for a boy image.

Figure 1 .

 1 Figure 1. Graph topology and feature vector illustration.

Figure 3 .

 3 Figure 3. Colorization results in local or nonlocal configurations for a girl image.

* This work was supported under a research grant of the ANR Foundation (ANR-06-MDCA-008-01/FOGRIMMI)


Olivier Lezoray
email: olivier.lezoray@unicaen.fr

Vinh Thong
email: vinhthong.ta@unicaen.fr

Abderrahim Elmoataz
email: abderrahim.elmoataz-billah@unicaen.fr

Impulse Noise Removal by Spectral Clustering and Regularization on Graphs *

In this paper we present a method for impulse noise removal that makes use of spectral clustering and graph regularization. The image is modeled as a graph and local spectral analysis is performed to identify noisy and noise free pixels. On the set of noise free pixels, a topology adapted graph regularization is performed. Experimental results show the benefits of the proposed approach regarding the standard VMF when noise proportion is high.

Introduction

Images are often corrupted by impulse noise due to noisy sensors or channel transmission errors. There are many works on the restoration of images corrupted by impulse noise. The most popular approach for removing impulse noise is the Vector Median Filter (VMF) [START_REF] Astola | Vector median filters[END_REF] because of its good denoising power and computational efficiency. Many improvements of the VMF have been proposed to avoid the modification of pixels not affected by impulse noise [START_REF] Smolka | Fast detection and impulsive noise removal in color images[END_REF]. Various decision-based filters have been proposed where noisy pixels are first identified and then replaced by using the VMF [START_REF] Lukac | A taxonomy of color image filtering and enhancement solutions[END_REF]. However, such approaches still have the same drawbacks of the VMF: blurring and low performances when noise ratio is high. In the spirit of [START_REF] Chan | Salt-and-pepper noise removal by median-type noise detectors and detailpreserving regularization[END_REF], we propose an iterative approach that, at each iteration, identifies corrupted pixels by spectral clustering on the filtering window and performs graph regularization by topology adaptation to take only into account uncorrupted pixels.

Preliminaries on graphs

We provide some basic definitions on graph theory. A graph G is a couple G = (V, E) where V is a finite set of vertices and E is a set of edges included in a subset of V × V . Two vertices u and v in a graph are adjacent if the edge (u, v) exists in E. u ∼ v denotes the set of vertices u connected to the vertex v via the edges (u, v) ∈ E. In the rest of this paper, we consider only simple graphs for which maximum one edge can link two vertices. These simple graphs are always assumed to be connected and undirected. A graph, as defined above, is said to be weighted if it is associated with a weight function w :

E → R + satisfying w(u, v) > 0 if (u, v) ∈ E, w(u, v) = 0 if (u, v) / ∈ E.

The volume of a set of vertices A can be computed by vol

(A) = vi∈A d(v i ) with d(v i ) = u∼vi w(u, v i ) the degree of a vertex.

Noise detection by spectral clustering

The transition probability from a vertex v i to a vertex v j is provided by p(v i , v j ) = w(vi,vj ) d(vi) . The associated transition matrix P is then defined by P = D -1 S where D is the degree matrix and S a similarity matrix associated to the graph G. The eigenvectors v of P are obtained by solving P v = λv. The associated eigenvalues of P are

λ 1 = 1 ≥ λ 2 ≥ • • • λ n ≥ -1 and one has P = yV y i=1 λ i v i v T i .
The eigenvector v 2 of P associated to λ 2 is known as the FIEDLER vector of P and provides most of geometrical information. Indeed, according to the sign of the FIEDLER vector, the graph can be partitioned into two sets. This type of analysis is known as spectral clustering [START_REF] Luxburg | A tutorial on spectral clustering[END_REF]. Spectral clustering can be used to partition a filtering window into two sets to identify noisy pixels. Given an arbitrary graph G = (V, E) that models an image where vertices correspond to pixels, the neighborhood set N (G, u) of a vertex u is defined as N (G, u) = {v ∈ V : (v, u) ∈ E} ∪ {u}. N (G, u) provides the set of vertices in a filter window centered on u. Then, we define the complete graph associated to a filtering window centered on u by

F (G, u) = (N (G, u), N (G, u) × N (G, u)). A spectral clustering of graph F (G, u
) is performed and according to the sign of the FIEDLER vector v 2 , the set of vertices of

F (G, u) is partitioned into two disjoint sets N 1 (G, u) = {v ∈ N (G, u) : v 2 (v) ≥ 0} and N 2 (G, u) = {v ∈ N (G, u) : v 2 (v) < 0}.
This can be considered as a natural nonparametric impulse noise detector since we use a complete graph and noise free pixel tend to form a single cluster while noisy pixels are outliers. From these sets, we then want to obtain two sets N c (G, u) and N uc (G, u) that respectively denote noisy and noise free candidates for a filtering window centered on u in a graph G that models an image. To that aim, we apply the following rules:

• If N 1 (G, u) or N 2 (G, u) has a cardinality lower or equal to 1, it is considered as being N c (G, u).
This corresponds to evident outliers.

• If vol(N1(G,u)) vol(N (G,u)) ≤ vol(N2(G,u)) vol(N (G,u)) , N c (G, u) = N 1 (G, u) and N c (G, u) = N 2 (G, u) otherwise.
N uc (G, u) is then obtained as the set N 1 (G, u) or N 2 (G, u) that was not retained as N c (G, u). These rules first treat the cases of one or less noisy pixels. Otherwise, the noisy set is determined as the less coherent set of vertices.

Topology adapted graph regularization

We consider a general function f 0 : V ⊂ R n → R m defined on graphs of the arbitrary topologies and we want to regularize this function. The regularization of such a function corresponds to an optimization problem which can be formalized by the minimization of a weighted sum of two energy terms:

1 2 u∈V |∇ w f (u)| 2 + λ f (u) -f 0 (u) 2 2 (1)
With the following definition of the weighted gradient

|∇ w f (u)| = v∼u w(u, v)(f (v) -f (u)) 2
, this problem (1) has a unique solution and can be solved with the following algorithm ∀u ∈ V [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing[END_REF]:

     f (0) (u) = f 0 (u) f (t+1) (u) = λf 0 (u) + v∼u w(u, v)f (t) (v) λ + v∼u w(u, v) . ( 2
)
Such a graph regularization enables local and nonlocal [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] regularization by using appropriated graphs topologies and edge weights [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing[END_REF]. However, it efficient for Gaussian noise removal but not for impulse noise. Therefore, we propose to combine graph regularization with spectral clustering to obtain a general algorithm of impulse noise removal. The key point is that the topology of the graph changes along the iterations and this is performed by retaining only noise free pixels (the N uc (G, u) set) on a given neighborhood centered on a vertex u. The major problem we have to face with is the fact that the vertex u does not necessarily belong to N uc (G, u) since it can be a corrupted pixel and Algorithm (2) is not directly applicable. Once we get the set N uc (G, u), an adapted neighborhood graph G(u) is constructed with N uc (G, u) as vertices and with {(v, v med ) : v ∈ N uc (G, u), v = v med } as edges. v med denotes the vertex that correspond to the vector median of the set N uc (G, u). First part of Figure 2 illustrates these different steps on a toy filter window. This corresponds to consider the following energy for one vertex u:

8 > > < > > : 1 2 y∇ G(u) w f (v med )y 2 + λ 2 f (u) -f 0 (u) 2 2 if u ∈ N uc (G, u) 1 2 y∇ G(u) w f (v med )y 2 if u / ∈ N uc (G, u) (3)
where |∇ G(u) w f (v med )| denotes the weighted graph gradient operator on the adapted neighborhood graph G(u). One has to note that the data fitting term is taken into account only for noise free vertices. Finally we obtain the following algorithm, ∀u ∈ V :

f (t+1) (u) = χ(u ∈ G(u))λf 0 (u) + P (v,v med )∈G(u) w(v, v med )f (t) (v med ) χ(u ∈ G(u))λ + P (v,v med )∈G(u) w(v, v med ) (4)
where χ : V → {0, 1} is the indicator function.

Experimental results

To test our algorithm, we have considered a part of the Lena image that has been corrupted by impulse noise expressed as [START_REF] Lukac | Adaptive vector median filtering[END_REF]:

x i,j = v with probability p v o i,j with probability 1 -p v
where i, j characterize the sample position, o i,j is the original sample, x i,j represents the sample from the noisy image, p v is a corruption probability and v = (v R , v G , v B ) is a noise vector of intensity random values. For the experiments, the considered degree of the impulse noise corruption p v has ranged from 5% to 90%. To evaluate the achieved results, objective criteria as Peak Signal-to-Noise Ratio (PSNR) and Normalized Color Difference (NCD) have been used [START_REF] Lukac | Adaptive vector median filtering[END_REF]. Figure 1 presents these results for the standard VMF, the proposed approach with 8-adjacency grid graphs (called RSVMF for Regularized Spectral VMF) and with a 24adjacency grid graph (called NLRSVMF for Nonlocal Regularized Spectral VMF). Second part of Figure 2 presents visual results of denoising. Weights are com- puted according to a measure of similarity w(u, v) = exp -F (f 0 ,u)-F (f 0 ,v) 2 σ 2

. F (f 0 , u) ∈ R q denotes a feature vector associated to each vertex u ∈ V . The feature vector associated to vertices F (f 0 , v) can be the initial function value: F (f 0 , v) = f 0 (v) (classical local processing: RSVMF) or a vector F (f 0 , v) = [f 0 (u) : u ∈ B v,s ] T (nonlocal processing: NLRSVMF). For this latter case F (f 0 , v) is a patch where B v,s denotes a bounding box of size (2s + 1) × (2s + 1) centered at v (3×3 in all our experiments). For all the experiments we set λ = 0.05. Since the proposed filters are iterative, we provide results at convergence. As depicted in Figure 1 and in the second part of Figure 2, our approach outperforms the VMF once the noise proportion is higher than 20%. Lower performances are obtained for lower noise proportions since our approach tends to sharpen the image while reducing the VMF blurry effect. For very high proportions of noise, the nonlocal version of the proposed approach is very effective. Moreover, nonlocal results show its remarkable robustness as noise proportion increases. Our approach is therefore very effective for removing very high proportions of impulse noise in local or nonlocal configurations.

Conclusion

In this paper, a new algorithm for impulse noise removal has been proposed. The image is modeled as a graph and a regularization process is performed. The latter has the property to locally adapt the topology of the graph under study by a spectral analysis to retain noise free pixels. Experimental results show that the new filter outperforms the VMF when high proportions of impulse noise should be eliminated.

Figure 1 .

 1 Figure 1. PSNR and NCD according to the percentage of impulse noise.

Figure 2 .

 2 Figure 2. First part: spectral clustering and adapted graph neighborhood construction. Second part: impulse noise removal results. First line presents corrupted images for different amounts of noise. Next lines present filtering results with the standard VMF, the Regularized Spectral VMF (RSVMF) on a 8-adjacency grid graph with F (f 0 , v) = f 0 (v) (local processing) and the Nonlocal Regularized Spectral VMF (NLRSVMF) on a 24-adjacency grid graph with F (f 0 , v) = [f 0 (u) : u ∈ B v,2 ] T (nonlocal processing).

* This work was supported under a research grant of the ANR Foundation (ANR-06-MDCA-008-01/FOGRIMMI)