Ferréol Soulez

Éric Thiébaut

JOINT DECONVOLUTION AND DEMOSAICING

Keywords: demosaicing, deconvolution, image reconstruction

We present a new method to jointly perform deblurring and colordemosaicing of RGB images. Our method is derived following an inverse problem approach in a MAP framework. To avoid noise amplification and allow for interpolation of missing data, we make use of edge-preserving spatial regularization and spectral regularization. We demonstrate the improvements brought by our algorithm by processing both simulated and real RGB images obtained with a Bayer's color filter and with different types of blurring.

INTRODUCTION

Most digital color cameras use a single sensor with a color filter array (CFA) to acquire color images. The Bayer CFA [START_REF] Bayer | Color imaging array[END_REF] shown on Fig. 1 is the most common CFA. The color interpolation techniques to produce a color image are known as demosaicing (see [START_REF] Alleysson | 30 ans de démosaicage -30 years of demosaicing[END_REF][START_REF] Gunturk | Demosaicking: color filter array interpolation[END_REF] for a review). In the case of blurred images, the information of unmeasured pixels spread to measured neighbor pixels. The main idea in this paper is to use this phenomenon to elaborate a joint deconvolution and demosaicing algorithm. The problem of demosaicing blurred data was studied for the first by Na et al. [START_REF] Na | An image restoration system for a single-ccd color camcorder[END_REF] who design a electronic device for that purpose. It was based on Weiner filtering of decorelated channels (YCrCb). Taubman [START_REF] Taubman | Generalized wiener reconstruction of images from colour sensor data using a scale invariant prior[END_REF] and Trussel [START_REF] Trussell | A mmse estimate for demosaicking[END_REF] then proposed theoretical methods using same Wiener filtering approachs. More recently, Trimeche et al. [START_REF] Trimeche | Multichannel image deblurring of raw color components[END_REF] use deconvolution technics as preprocessing before demosaicing. Based of super-resolution framework Vega et al. [START_REF] Vega | A bayesian super-resolution approach to demosaicing of blurred images[END_REF] described a Bayesian framework for deconvolution of mosaicked data. We can also point out the work of Farsiu et al. [START_REF] Farsiu | A practical approach to superresolution[END_REF] about the super-resolution of multi-frame blurred and mosaicked data.

MODEL DESCRIPTION

The recording of a given scene using a Bayer camera can be broken down in several successive operations. First the scene get blurred by different physical phenomenons like optical aberrations (e.g. out of focus) and motion; then it is sampled by the Bayer's color filter array (CFA) and finally it is recorded by the sensor with additive noise. This acquisition scheme is sum up on Fig. 2. Fig. 2. Acquisition process with a Bayer camera.

Blurring

The blurred RGB image g(s) in a direction s is made of three color channels:

g(s) = {g R (s), g G (s), g B (s)} .

(1) In color channel c ∈ {R, G, B}, the blurred image g c (s) in a direction s is given by:

g c (s) = h c (s|s ′ ) x c (s ′ ) ds ′ , (2)
where x c (s ′ ) is the object brightness distribution in color channel c and direction s ′ , h c (s|s ′ ) is the point spread function (PSF). The PSF h c (s|s ′ ) is the observed relative brightness distribution in color channel c and direction s for a point source located in direction s ′ . If the PSF is isoplanatic (shift-invariant), then h c (s|s ′ ) = h c (ss ′ ) and Eq. ( 2) becomes a convolution product. For data sampled on N pixels, Eq. ( 2) takes a matrix form:

g c = H c • x c , (3)
where H c is the N × N response matrix, g c and x c are vectors of size N corresponding respectively to the blurred and incident brightness distributions. Taking the three color channels into account yields:

g = H • x (4)
where g and x are vectors of size 3 N obtained by packing the color components, e.g. g = ( g T R , g T G , g T B ) T ; likewise, H is a 3 N ×3 N sparse block diagonal matrix:

H =          H R 0 0 0 H G 0 0 0 H B          . (5)
In this study, we consider known spatial PSF's that are isoplanatic and identical in all color channels. This simplifies computations but has no other incidences: our algorithm can still be used if these assumptions have to be relaxed.

When a Bayer camera is used, the vector of available data y is a noisy and incomplete subset of g:

y = B • g + n = B • H • x + n , (6)
where B is an N × 3 N down-sampling matrix defined according to the color filter array and n is a N vector accounting for the noise (signal noise plus detector readout noise plus digitization noise and model approximations).

INVERSE PROBLEM APPROACH

To recover a color image from the measurements y, one has to solve a deblurring problem with incomplete data and additive noise.

To solve this difficult inverse problem which is both ill-posed and ill-conditioned, we use a penalized likelihood or maximum a posteriori [START_REF] Thiébaut | Strict a priori constraints for maximum likelihood blind deconvolution[END_REF] framework. The maximum a posteriori solution x MAP minizes the cost function ε(x):

x MAP = arg min x ε(x) (7)
where this penalty function ε(x) expands as:

ε(x) = Φ lkl (x; y) + µ Φ prior (x) , (8)
Minimizing the likelihood penalty Φ lkl (x; y) enforces the agreement between the model x and the data y, whereas minimizing the regularization penalty Φ prior (x) accounts for subjective a priori knowledge about the sought image. Hence minimizing the criterion in Eq. ( 8) yields a solution which is a trade-off between fitting the data and keeping the solution consistent with the priors. The balance of the trade-off is tuned by the hyper-parameter µ.

The Likelihood Penalty

For Gaussian noise, using the model defined in Eq. ( 6), the likelihood penalty writes:

Φ lkl (x; y) = [y -B • H • x] T • C -1 noise • [y -B • H • x] , (9)
where C noise = n • n T is the covariance matrix of the noise. Assuming uncorrelated noise, C noise is diagonal and Eq. ( 9) becomes:

Φ lkl (x; y) = k 1 σ 2 k (B • H • x) k -y k 2 , (10)
where σ 2 k is the noise variance for measured pixel k. For the results presented at the end of this paper, we assumed uncorrelated uniform noise with a statistic independent of the color of pixels, i.e. σ 2 k = σ 2 .

The Regularization Penalty

Restoration of the RGB image x involves interpolation of colors and deblurring. In this case, and at least to avoid noise amplification, the regularization term must enforce the solution to be somewhat smooth. Besides, since the spatial and spectral (color) dimensions of x are not homogeneous, we consider a separable regularization penalty:

µ Φ prior (x) = µ Φ spectral (x) + c∈{R,G,B} µ c Φ spatial (x c ) , (11)
where Φ spectral (x) is a spectral regularization and Φ spatial (x c ) is the spatial regularization for color channel c.

Spatial Regularization

As the noise mostly contaminates high frequencies, smoothness constraint yields the most effective regularization. We therefore choose:

Φ spatial (x c ) = k k ′ ∈V k ϕ x c,k -x c,k ′ ℓ k,k ′ , (12)
where ϕ(u) is some increasing function of |u|, x c,k is the value of k-th pixel in color channel c of the restored image, k ′ is a pixel in the neighborhood V k of pixel k (usually a V8 neighborhood) and ℓ k,k ′ is the distance between these pixels. Since the likelihood term is quadratic, using a quadratic regularization, i.e. ϕ(u) ∝ u 2 , would give an analytical expression for the solution of Eq. [START_REF] Trimeche | Multichannel image deblurring of raw color components[END_REF]. Quadratic regularization however has the drawback that it reintroduces some unwanted blur and ripples by over-smoothing sharp edges of the image. Non-quadratic norms can be used to preserve such sharp feature. Following this approach, we choose a so-called ℓ 2ℓ 1 norm [START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF] for ϕ:

ϕ(u; λ) = 2 λ 2 |u|/λ -log (1 + |u|/λ) , (13)
where λ > 0 is the threshold beyond which the difference between neighbor pixels is most certainly due to an edge and must not be too much smoothed. Indeed the norm in Eq. ( 13) is asymptotically quadratic (resp. linear) for small (resp. large) pixel differences compared to the threshold λ. In principle, one can choose a different threshold for each color channel and take:

Φ spatial (x c ; λ c ) = k k ′ ∈V k ϕ x c,k -x c,k ′ ℓ k,k ′ ; λ c . (14)

Spectral Regularization

Gunturk et al. [START_REF] Gunturk | Color plane interpolation using alternating projections[END_REF] have shown that there is strong correlation between color channels at high spatial frequencies. According to this, we design our spectral regularization penalty so as to minimize the variance between color channels at high frequencies:

Φ spectral (x) = P • d R,G 2
2 + P • d R,B 2
2 + P • d B,G 2 2 , (15)
where d c,c ′ = x cx c ′ is the pixel-wise difference between the monochromatic images in color channels c and c ′ and where P is a high-pass filter. For our tests, we took P = I -L where I is the identity and L is the flux-preserving low-pass filter obtained by convolving its argument by the 5 × 5 separable smoothing kernel:

K = 1 256                   1 4
                  . ( 16
)

ALGORITHM SUMMARY

In this Bayesian framework, restoring the de-blurred image is obtained by seeking for the RGB image x which minimizes the composite criterion:

ε(x; y, θ) = Φ lkl (x; y)+ c∈{R,G,B} µ c Φ spatial (x; λ c )+µ Φ spectral (x) , (17)
where θ = {µ, µ R , µ G , µ B , λ R , λ G , λ B } is the set of chosen hyperparameters. Given the observed image y, the restored image x only depends on θ, that is on four regularization levels (the µ's) and three edge-preserving thresholds (the λ's).

To derive proper hyper-parameter values, we assumed that the image dynamic is identical in each color channel. As a consequence, the spatial regularization weights and the edge-preserving thresholds are the same in every color channels: µ R = µ G = µ B = µ RGB , and λ R = λ G = λ B = λ RGB . Hence, only three hyper-parameters have to be chosen: θ = {µ, µ RGB , λ RGB }. Despite this simplification, choosing the optimal values of the hyper-parameters is cumbersome and difficult. Whether methods such as generalized cross-validation (GCV) [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF] or the L-curve [START_REF] Hansen | Analysis of discrete ill-posed problems by means of the l-curve[END_REF] are suitable for this task deserves an extensive study which is out of the scope of this paper. In the present work, we simply choose hyper-parameter values by trial and error and visual inspection of the resulting image.

Optimization Algorithm

In order to determine the optimal image x MAP in our inverse problem approach, one has to minimize a criterion with respect to a very large number of variables (all the pixel values for every color channels). To that end, we used the VMLM-B algorithm [START_REF] Thiébaut | Optimization issues in blind deconvolution algorithms[END_REF] which is a limited memory variant of the variable metric method with BFGS updates [START_REF] Nocedal | Theory of algorithms for unconstrained optimization[END_REF]. This algorithm, which can further account for bound constraints on the parameters, has proven effectiveness for image reconstruction and only requires the computation of the penalty function to be minimized and its gradient. The memory requirement is a few times the size of the problem.

RESULTS

Several experiments on both simulated and real data were carried out to assess the performance of the presented algorithm. To quantify the reconstruction quality, we use the Peak Signal to Noise Ratio (PSNR) which is classically used to measure improvements of digital image quality. This criterion corresponds to a mean squared error normalized by the maximum pixel value M:

E PSNR = -10 log 10         1 K M 2 c,k xc,k -x c,k 2         , (18)
where K is the number of pixels, x and x are respectively the restored and the true images. Here the signal is assigned to the original image and the noise corresponds to the reconstruction error. Thus, the larger is the PSNR, the better is the reconstruction. In the following, the results are presented in terms of PSNR improvement with respect to the PSNR between the blurred data before mosaicing. As convolution is done by means of Fourier Transforms, field of view aliasing corrupts a narrow border of the image which is thus not taken into account for this quality assessment.

Simulation

Simulations were performed on several images from the Kodak database (http://www.cipr.rpi.edu/resource/stills/kodak.html) and shown by Fig. 3. All these images are of size 512 × 728 and are digitized on 8 bits (256 quantification levels).

To assess the deblurring ability of our algorithm, we consider different types of blur: (i) Gaussian blur with full width at half maximum of 4 pixels, (ii) out of focus blur by a uniform disk PSF of radius 4 pixels, and (iii) vertical uniform motion blur of 5 pixels. The reconstruction using our method (column J in Tab. 1) is compared with a simple bilinear interpolation reconstruction (column B Two main tendencies can be retrieved from these results. First, as blur attenuates the high frequencies and therefore smooths the edges of the image, bilinear interpolation is quite close (-0.45 dB) to the blurred image before mosaicing. Second, except in one case, the presented method achieves better reconstruction than interpolation followed by deconvolution.

Experimental data

We carried out an experimental test of our algorithm on really blurred pictures. These pictures were taken using a Canon EOS 350D reflex digital camera which could record images in a mosaicked raw format. Several pictures of the same scene were taken with different level of defocus. We used a parametric myopic deconvolution to derive the approximate shape of the defocus PSF shown in Fig. 5(c). Deconvolution by the defocus PSF was then performed on a 1024 × 1024 pixels part of this scene. The result is shown on Fig. 5(b) and can be compared with bilinear reconstruction shown on Fig. 5(a) and with an in-focus and automatically processed (demosaicing, white balance, dynamics...) picture of the same scene on Fig. 5(d). Although some ringing artifacts are visible, the image was correctly restored (it is possible to read the writing on the reconstruction). This experimental result shows that the presented algorithm is robust to imperfection in the PSF knowledge.

CONCLUSION

In this paper we have presented a unified method to achieve combined deconvolution and demosaicing. On simulation, we show that joint deconvolution and demosaicing is more effective than a successive demosaicing then deconvolution. Experimental results validate the algorithm on real data's and show its robustness to partially known PSF.

Fig. 1 .

 1 Fig. 1. Bayer color filter array.

Fig. 3 .

 3 Fig. 3. Data set used in our simulation. From left to right: lighthouse, sail, statue and window.in Tab. 1) and the deconvolution of this bilinear reconstruction (column BD in Tab. 1). This deconvolution of the bilinear interpolation is performed using the same algorithm as the presented one but without down-sampling by the CFA matrix B. Results are shown in table1. A detail of the images processed for the 9-th row of this table are shown on Fig.4.

Fig. 4 .

 4 Fig. 4.Details of the statue picture with the motion blur as presented in Tab. 1. Let us point your attention to the crown and the wrinkle details.

Fig. 5 .

 5 Fig. 5. Experimentally defocused image. Let us point your attention to the inscription details.

  table are shown on Fig. 4.

	Image	Blur	Blurred	B	BD	J
	Light-	Gaussian	24.61	-0.28 +3.99 +5.60
	house	Defocus	22.97	+0.01 +6.17 +6.90
		Motion	27.37	-1.42 +0.82 +3.58
		Gaussian	27.96	-0.36 +5.00 +5.93
	Sail	Defocus	26.50	-0.05 +6.30 +6.46
		Motion	30.66	-0.98 +2.03 +4.33
		Gaussian	28.36	-0.29 +4.11 +4.86
	Statue	Defocus	27.20	-0.06 +5.60 +5.69
		Motion	29.83	-0.52 +2.94 +5.57
		Gaussian	28.10	-0.46 +5.69 +6.55
	Window	Defocus	26.59	-0.17 +6.81 +6.23
		Motion	30.10	-0.77 +4.04 +5.65
	average		-0.45 +4.46 +5.61
	Table 1. PSNR improvement in decibel of the bilinear interpolation
	(B), its deconvolution (BD) and the algorithm presented in this paper
	(J).