Laurent Condat
email: laurent.condat@ensicaen.fr

GREYC Saleh Mosaddegh

JOINT DEMOSAICKING AND DENOISING BY TOTAL VARIATION MINIMIZATION

Keywords: Demosaicking, denoising, Bayer color filter array, frequency selection, spatio-spectral sampling

Joint demosaicking and denoising consists in reconstructing a color image from the noisy raw data output by the sensor of a digital camera. We adopt a variational formulation in which the reconstructed image has minimal total variation under the constraint of consistency with the available measurements. This way, the recovered color image has smooth chrominance but the sharp edges are maintained and the noise is transferred to the luminance channel. This channel is denoised subsequently.

INTRODUCTION

Color images are acquired in digital cameras by means of a single sensor on which the Bayer color filter array (CFA) is overlaid [1]. In order to reconstruct a full-color image from the raw data delivered by the sensor, an interpolation process called demosaicking is performed. There is an abundant literature on demosaicking and we direct the readers to the good survey by Menon [START_REF] Menon | Color image reconstruction for digital cameras[END_REF]. However, most demosaicking methods are developed under the unrealistic assumption of noise-free data. In the presence of noise, the performances of the algorithms degrade drastically, since their sophisticated nonlinear mechanisms are generally not robust to noise. Moreover, denoising after demosaicking is untractable, because demosaicking distorts the characteristics of the noise in a complex and hardly computable form. Thus, demosaicking and denoising have to be handled jointly. We refer to the introduction of [START_REF] Condat | A simple, fast and efficient approach to denoisaicking: Joint demosaicking and denoising[END_REF] for a survey of the relevant literature.

To formulate the problem, let us first introduce some notations. Boldface letters denote vectors, e.g. a = [a 1 , a 2 ] T ∈ C 2 with norm |a| = |a 1 | 2 + |a 2 | 2 . We define the color image u = (u[k]) k∈Z 2 as the ground-truth to be estimated. For every k, u

[k] = u R [k], u G [k], u B [k]
T is the color of the pixel of u at location k, in the canonical R, G, B (red, green, blue) basis. In this paper, we adopt an additive white Gaussian noise (AWGN) model; that is, we have at our disposal the noisy mosaicked image v such that

v[k] = u X[k] [k] + ε[k], ∀k ∈ Z 2 , ( 1
)
where X[k] ∈ {R, G, B} is the color of the filter in the Bayer pattern at location k (see Fig. 1a), ε[k] ∼ N (0, σ 2 ) for every k and σ 2 is the noise variance. Then, the problem is to reconstruct a color image d from v, which is a good estimate of u. In real conditions, the AWGN assumption is not met; real noise is more accurately modeled by the sum of a Gaussian and a Poissonian components [START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for singleimage raw-data[END_REF]. Moreover, the observed values are photon counts, which have to be tone mapped/gamma corrected. However, variance stabilization techniques can be efficiently employed [START_REF] Foi | Clipped noisy images: heteroskedastic modeling and practical denoising[END_REF], so that the problem can be recast in the AWGN context. The article is organized as follows. In sect. 2, we recall the properties of the joint demosaicking and denoising approach by frequency selection we proposed in [START_REF] Condat | A simple, fast and efficient approach to denoisaicking: Joint demosaicking and denoising[END_REF]. To improve upon this method, we formulate a new optimization problem using total variation and we propose an algorithm to solve it in sect. 3. The approach is validated by experiments in sect. 4.

SPATIO-SPECTRAL MODEL OF SAMPLING AND DEMOSAICKING BY FREQUENCY SELECTION

It is well known that the R, G, B components of natural images are strongly correlated [1]. That is why we define the components of luminance, green/magenta and red/blue chrominances of a color image a as a L = a, L , a G/M = a, C G/M , and a R/B = a, C R/B , respectively, using the orthonormal basis

L = 1 √ 3 [1, 1, 1] T , C G/M = 1 √ 6 [-1, 2, -1] T , C R/B = 1 √ 2 [1, 0, -1] T .
In first approximation, the luminance and chrominance channels of natural images in this basis are statistically decorrelated. A major contribution of Alleysson et al. [START_REF] Alleysson | Linear demosaicing inspired by the human visual system[END_REF] consisted in showing that the basis L, C G/M , C R/B is appropriate to characterize the Bayer CFA and that the mosaicked image is the sum of the modulated luminance and chrominance components of u:

v(ω) = 1 √ 3 ûL (ω) + 1 √ 24 ûG/M (ω) + √ 6 4 ûG/M (ω-[π, π] T ) + √ 2 4 ûR/B (ω-[0, π] T ) - √ 2 4 ûR/B (ω -[π, 0] T ) + ε(ω), (2)
for every ω ∈ R 2 , where the Fourier transform â(ω) of an image a is defined as

â(ω) = k∈Z 2 a[k]e -jω
T k and an image with finite support is implicitely extended to an infinite one by zero-padding.

This frequency analysis of the spatio-spectral sampling induced by the Bayer CFA, illustrated in Fig. 1b, sheds an interesting light on the problem: it aims at separating the three images u L , u R/B and u G/M from their noisy mixing in v. This is exactly what demosaicking by frequency selection does in the noise-free case, as proposed by Dubois [START_REF] Dubois | Frequency-domain methods for demosaicking of Bayer-sampled color images[END_REF]. We extended the method to the noisy case in our last paper [START_REF] Condat | A simple, fast and efficient approach to denoisaicking: Joint demosaicking and denoising[END_REF]. In short, the approach consists in estimating u R/B and u G/M by modulation and lowpass filtering. Then, the residual in v is known to be an estimate of 1 √ 3 u L + ε, from eqn. [START_REF] Menon | Color image reconstruction for digital cameras[END_REF]. Thus, in first approximation, the noise is completely contained in the luminance channel of the demosaicked image and this channel can be subsequently denoised using any method adapted to grayscale images and AWGN. This approach yields stateof-the-art results but, due to the estimation of the chrominance using very selective lowpass filters, the sharp color edges are blurred. So, there is still room for improvement over the method of [START_REF] Condat | A simple, fast and efficient approach to denoisaicking: Joint demosaicking and denoising[END_REF], which motivates this work.

In [START_REF] Condat | Le dématric ¸age par sélection spectrale revisité[END_REF], we have shown that, with minor modifications, in the noise-free case, demosaicking by frequency selection yields the solution of the following variational problem considered in [START_REF]A generic variational approach for demosaicking from an arbitrary color filter array[END_REF]:

d = argmin a µ ∇a L 2 ℓ2 + ∇a G/M 2 ℓ2 + ∇a R/B 2 ℓ2 (3) subject to a X[k] [k] = v[k], ∀k ∈ Z 2 ,
where we introduce the discrete gradient vector using finite differences as

∇a[k] = a[k] -a[k 1 -1, k 2 ], a[k] -a[k 1 , k 2 -1] T , (4)
for every k ∈ Z 2 , using Neumann boundary conditions; that is, a finite difference is set to zero when it involves a pixel value outside the image domain. We also introduce the discrete divergence operator div = -∇ * as

div a[k] = a 1 [k 1 +1, k 2 ]-a 1 [k]+a 2 [k 1 , k 2 +1]-a 2 [k]. (5)
The parameter µ in (3) plays a crucial role; it controls the balance between the smoothness of the luminance and of the chrominance in the reconstructed image. If µ is close to zero, then d will be close to a monochrome image, since all the high frequency energy of v will be assigned to d L . On the contrary, for µ = 1, the regularization functional is diagonal in the R, G, B basis and the process amounts to reconstructing the R, G, B channels independently by interpolation, a naive solution which yields bad results. Consequently, µ should be chosen relatively small in order to get a smooth hue [START_REF]A generic variational approach for demosaicking from an arbitrary color filter array[END_REF]. This way, the inter-correlations between color channels in natural images are automatically taken into account.

In the noisy case, by further reducing µ and keeping exact consistency with v, we reconstruct an image d with almost all the noise of v assigned to the luminance band d L . However, when decreasing µ, the high frequencies of chrominance are mistakenly assigned to the luminance channel. Thus, at sharp color transitions, the edges are over-smoothed and zipper artifacts appear in the luminance. This motivates a formulation of the problem using non-quadratic regularization, namely total variation, whose superiority in restoration problems in keeping sharp edges is well known.

A VARIATIONAL FORMULATION BASED ON TOTAL VARIATION MINIMIZATION

First, we remark that the regularization in (3) is invariant with respect to a rotation in every iso-luminance color plane; that is, there is no privileged color axis and the problem formulation is independent on the choice of the chrominance basis C G/M , C R/B . We have the equality

∇a G/M 2 ℓ2 + ∇a R/B 2 ℓ2 = ∇a C 2 ℓ2 ( 6
)
where a C = a G/M + j.a R/B and j is the complex square root of -1. In this work, we propose a new formulation of the total variation (TV) of a color image:

a TV = µ ∇a L ℓ1 + ∇a C ℓ1 (7)
for some parameter µ > 0, where the ℓ 1 norm of a vectorvalued image a is a ℓ1 = k∈Z 2 a[k] H a[k] and • H indicates the complex conjugate of the transpose of a vector. Then, we formulate demosaicking as the following optimization problem, for some parameter 0 < µ < 1 to choose:

d = argmin a a TV s. t. a X[k] [k] = v[k], ∀k ∈ Z 2 .
(8) Note the difference with usual variational formulations in which we trade the fit to the data and the smoothness of the solution. Here we keep the exact consistency with the data, so that the noise is not removed but is transfered to the luminance channel d L .

Minimizing a quadratic penalty like in (3) boils down to solving a linear system and a direct implementation can be designed [START_REF] Condat | A simple, fast and efficient approach to denoisaicking: Joint demosaicking and denoising[END_REF][START_REF] Condat | Le dématric ¸age par sélection spectrale revisité[END_REF]. With the non-quadratic TV, there is no direct way to obtain the solution and an iterative method has to be designed. Note that the TV is convex so that the problem ( 8) is well posed. However, the TV is not differentiable so that conventional smooth optimization techniques are not applicable. Actually, very few methods are available to minimize the TV under an affine constraint, see [START_REF] Alfonso | An augmented Lagrangian approach of the constrained optimization formulation of imaging inverse problems[END_REF] and references therein. Very recently, a breakthrough in the field has appeared under the form of new primal-dual methods proposed independently by several authors [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Esser | A general framework for a class of first order primal-dual algorithms for convex optimization in imaging science[END_REF][START_REF] Zhang | A unified primal-dual algorithm framework based on Bregman iteration[END_REF]. In this work, we apply the algorithm of [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] to our problem. This yields the following implementation:

Demosaicking Algorithm

1. Choose the initial estimate d (0) and the constant α > 0

2. n := 0; b (0) := d (0) ; β := 1/(8.01 α); ∀k ∈ Z 2 , ∀X ∈ {R, G, B}, a X (0) [k] := [0, 0] T 3. Repeat until stopping criterion is met 4. ∀X ∈ {R, G, B}, a X (n+1) := a X (n) + α∇b X (n) 5. ∀k ∈ Z 2 , a L (n+1) [k] := a L (n+1) [k] max(1, |a L (n+1) [k]|/µ) 6. ∀k ∈ Z 2 , a C (n+1) [k] := a C (n+1) [k] max(1, |a C (n+1) [k]|) 7. ∀X ∈ {R, G, B}, d X (n+1) := d X (n) + β div a X (n+1) 8. ∀k ∈ Z 2 , d X[k] (n+1) [k] := v[k] 9. b (n+1) := 2d (n+1) -d (n) 10. n := n + 1
It can be shown [11, Theorem 1] that this algorithm converges to the solution d of (8).

EXPERIMENTAL VALIDATION

In [START_REF] Condat | A simple, fast and efficient approach to denoisaicking: Joint demosaicking and denoising[END_REF], different methods of the literature were compared and the method proposed in [START_REF] Condat | A simple, fast and efficient approach to denoisaicking: Joint demosaicking and denoising[END_REF] was shown to be significantly better than every other method. Therefore, by lack of space, we only compare the proposed TV approach to the method of [START_REF] Condat | A simple, fast and efficient approach to denoisaicking: Joint demosaicking and denoising[END_REF]. The 24 color images of the classical Kodak test set 1 1 The vertical images were first rotated by 90 o counterclockwise to simulate an acquisition with a digital camera whose sensor is aligned horizontally. were mosaicked 2 with the Bayer CFA and corrupted with different noise levels. 50 iterations of the proposed algorithm were run, α was set to 0.1 and the initial guess d (0) was set as the result of the method of [START_REF] Condat | A simple, fast and efficient approach to denoisaicking: Joint demosaicking and denoising[END_REF]. The state-of-the-art BM3D denoising method [START_REF] Dabov | Image denoising by sparse 3D transform-domain collaborative filtering[END_REF] was used to denoise the luminance image d L . In Tab. 1, we report the CPSNR 3 between u and d, averaged over the 24 images. We also propose an extended variant of our approach in which the reconstructed image is mosaicked again and then demosaicked using the method of [START_REF] Zhang | Color demosaicking via directional linear minimum mean square-error estimation[END_REF]. Indeed, the re-mosaicked image is relatively free of noise and a classical demosaicking method can be applied to it. From the numerical results in Tab. 1, we see that the proposed approach is only slightly better than the method of [START_REF] Condat | A simple, fast and efficient approach to denoisaicking: Joint demosaicking and denoising[END_REF]. However, as shown in Fig. 2, the proposed TV-based method yields cleaner edges. The extended variant provides a significant gain for low noise levels, where the denoising process has a negligible effect and the nonlinear demosaicking methods can show their superiority. This gain vanishes as the noise level increases.

CONCLUSION

In this work, we proposed a new expression of the total variation of a color image and a new primal-dual algorithm to minimize it under a consistency constraint. Applied to the illposed inverse problem of demosaicking noisy data, the proposed variational approach yields even better results that the state-of-the-art method of [START_REF] Condat | A simple, fast and efficient approach to denoisaicking: Joint demosaicking and denoising[END_REF], with sharp color edges free of noticeable artifacts. Moreover, our method is generic and can be applied to images acquired with an arbitrary CFA [START_REF] Condat | A new color filter array with optimal properties for noiseless and noisy color image acquisition[END_REF]. A Matlab implementation of our approach is available online.

[1] B. K. Gunturk, J. Glotzbach, Y. Altunbasak, R. W. Schaffer, and R. M. Mersereau, "Demosaicking: Color filter array interpolation," IEEE Signal Processing Mag., vol. 22, no. 1, pp. 44-54, Jan. 2005. 2 The filter for the top-left pixel was green and its right neighbor was red. 3 The CPSNR between the color images u 1 and u 2 is equal to 10 log

Fig. 1 .

 1 Fig. 1. (a) The Bayer CFA and (b) schematic representation of the spectrum of a mosaicked image, with the spectrum of the luminance u L in the baseband and the replicas of the spectra of chrominance with hatched fill.

  10 (255 2 /MSE) with MSE = P X=R,G,B P k (u X 1 [k]u X 2 [k]) 2 /(3N.M ) for images of size N × M .A 20 pixel-wide band around the border of the images was ignored when computing the CPSNR, to get rid of boundary effects. The pixel values were rounded to integers in 0 . . . 255 before computing the CPSNR.

Fig. 2 .

 2 Fig. 2. (a): part of the ground-truth image no. 3 of the Kodak base, (b)-(c): parts of the reconstructed images with the method in [3] and the proposed method, respectively, from the mosaicked and noisy (σ = 20) version of (a). (d)-(e): same as (b)-(c) with a higher zoom factor, to emphasize the differences.

39.43 35.70 32.80 29.78Table 1 .

 1 Average CPSNR (in dB) over the 24 images of the Kodak test set for different demosaicking+denoising methods.

		σ	1	5	10	20
	Method	µ	0.5	0.45	0.4	0.35
	[3]		38.49 35.35 32.56 29.57
	[3]+mosaicking+[15] 39.42 35.66 32.73 29.68
	proposed		38.41 35.36 32.63 29.70
	prop.+mosaicking+[15]