A Fast Algorithm for the Demosaicing Problem Concerning the Bayer Pattern
Antonio Boccuto1, Ivan Gerace1, *, Valentina Giorgetti1, 2, Matteo Rinaldi1
Identifiers and Pagination:Year: 2019
First Page: 1
Last Page: 14
Publisher Id: TOSIGPJ-6-1
Article History:Received Date: 29/6/2018
Revision Received Date: 6/12/2018
Acceptance Date: 18/12/2018
Electronic publication date: 29/1/2019
Collection year: 2019
open-access license: This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 International Public License (CC-BY 4.0), a copy of which is available at: https://creativecommons.org/licenses/by/4.0/legalcode. This license permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
In this paper, we deal with the demosaicing problem when the Bayer pattern is used. We propose a fast heuristic algorithm, consisting of three parts.
In the first one, we initialize the green channel by means of an edge-directed and weighted average technique. In the second part, the red and blue channels are updated, thanks to an equality constraint on the second derivatives. The third part consists of a constant-hue-based interpolation.
We show experimentally how the proposed algorithm gives in mean better reconstructions than more computationally expensive algorithms.
The demosaicing problem is related to the acquisition of RGB color images by means of CCD digital cameras. In the RGB model, each pixel of a digital color image is associated to a triple of numbers, which indicate the light intensity of the red, green and blue channel, respectively. However, most cameras use a single sensor, associated with a color filter that allows only the measure at each pixel of the reflectance of the scene at one of the three colors, according to a given scheme or pattern, called Color Filter Array (CFA). For this reason, at each pixel, the other two missing colors should be estimated. Different CFA’s are proposed for the acquisition [1-3]. The most common is the Bayer pattern . In this scheme, the numbers of pixels in which the green color is sampled are double with respect to those associated with the red and blue channels, because of the higher sensibility of the human eye to the green wavelengths. If we decompose the acquired image into three channels, we obtain three downsampled grayscale images, so that demosaicing could be interpreted as interpolating grayscale images from sparse data. In most cameras, demosaicing is a part of the processing required to obtain a visible images. The camera’s built-in-firmware is substantially based on fast local interpolation algorithms.
The heuristic approaches, which do not try to solve an optimization problem defined in mathematical terms, are widely used in the literature. These methods, in general, are very fast. Our proposed technique is of heuristic kind. In general, the heuristic techniques consist of filtering operations, which are formulated by means of suitable observations on color images. The nonadaptive algorithms, among which bilinear and bicubic interpolation, yield satisfactory results in smooth regions of an image, but they can fail in textured or edge areas. Edge-directed interpolation is an adaptive approach, where, by ana-lyzing the area around each pixel, we choose the possible interpolation direction. In practice, the interpolation direction is chosen to avoid interpolating across the edges. In , for each pixel the horizontal and vertical gradients are compared with a constant threshold. If the gradient in one direction is greater than the threshold, then interpolation is not performed along this direction. Some other direct interpolation methods use larger neighborhoods by examining different color channels. In , to determine the edges of the green channels, the red and blue channels are employed. On the other hand, to determine the edges of the red and blue channels, some discrete derivation operators of the second order are used, while in , to determine the edges in the various channels, a suitable Jacobian operator is applied. In , local homogeneity is used as an indicator to choose horizontally or vertically interpolated intensities. Thanks to homogeneity-directed interpolation, the luminance and chrominance values have to be similar in a suitable neighborhood. In demosaicing it is often assumed that the differences or the ratios of the intensity values in different channels are locally constant [5, 6, 9-14]. In  the probability of having an edge in a certain direction is determined and used to find the weights relative to the weighted average employed as an interpolation operator. In this algorithm, the color channels are updated iteratively according to the constant color ratio condition. In  a similar algorithm is proposed, where -size neighborhoods are employed to find the edges of the green channel, and -size neighborhoods are used to determine the edges of the red and blue channels. An analogous algorithm is defined in , where the interpolation can be done also in the diagonal direction, while in  the weighted directional interpolation is used by means of a fuzzy membership assignment. A second order operator is employed as a correction term. 
To have more accurate results, several techniques, which use iterative methods, are proposed. However, they have a higher computational cost with respect to the heuristic techniques. One of well-known techniques is the algorithm of Alternate Projections (AP) , which uses the strong correlation between the high frequences of the three colored components, by projecting alternately the estimated image in a constraint of observation and in a constraint which imposes similarity between the red and green edges and between the blue and green edges, until a fixed point is found. Another widely used technique is regularization [20, 21]. The algorithm is based on interpolation in a residual domain . The residuals are the differences between the observed and estimated pixel values which minimize a Laplacian energy.
The algorithm here presented consists of three steps. The first two ones are initialization steps, while the third one is an iterative steps. In the first one, the missing valued in the green component are determined, in particular a weighted average-type technique is used. The weights are determined in an edge-directed approach, in which we consider also the possible edges in the red and blue components. In the second step, we determine the missing values in the red and blue components. In this case we use two alternative techniques, according to the position of the involved pixel in the Bayer pattern. In the first technique, the missing value is determined by imposing that the second derivative of the intensity value of the red/blue channel is equal to the second derivative of the intensity values of the green channel. This is done according to the proposed approaches in the AP algorithm and the regularization algorithm given in . In particular, a constraint is imposed, to make the derivatives of all channels similar as soon as possible . At the third step, all values of the three channels are recursively updated, by means of a constant-hue-based technique. In particular, we assume the constant color difference. The technique we propose at this step is similar to that used by W. T. Freeman . Indeed, even here a median filter is employed, in order to correct small spurious imperfections. We repeat iteratively the third step. However, to avoid increasing excessively the computational cost, we experimentally estimate that only four iterations are necessary to obtain an accurate demosaicing. We call our technique as Local Edge Preserving (LEP) algorithm.
The paper is structured as follows. In Section 2 we give a mathematical formulation of the demosaicing problem. In Section 3 we describe the initialization of the proposed algorithm, which consists of the two first steps aforementioned. In Section 4 we give the third iterative step of our algorithm, highlighting the differences with the Freeman filter. In Section 5 our experimental results are presented. This section consists of two parts. In the first one, we determine the best detection function which can be used in order to evaluate the edges. In the second one, we compare our algorithm with some other techniques recently proposed in the literature and we show how the LEP method gives in mean more accurate reconstructions than the other considered algorithms.
2. THE DEMOSAICING PROBLEM
An RGB (Red-Green-Blue) color image with height n and width m is a vector of the type
where are the red, green and blue channels according to the lexicographic order, respectively. We consider the problem of acquisition of data from a digital camera, and call it "mosaicing problem". Given an ideal image , the acquired or "mosaiced image" is defined by
where and is a linear operator defined by setting
is the null matrix, and are diagonal matrices whose principal entries, if we use the Bayer pattern Fig. 1, are given by
where the symbol indicates that i-j is even. The corresponding demosaicing problem is the associated inverse problem, that is to determine the ideal color image x, knowing the mosaiced image and the linear operator M. An inverse problem is said to be well-posed (in the sense of Hadamard) if and only if the solution exists, is unique and stable with respect to data variation. A not well-posed problem is said to be ill-posed . Note that the demosaicing problem is ill-posed, since the matrix M in (2) is singular, as it is readily seen, and so there are infinitely many solutions.
|Fig. (1). Bayer pattern.|
3. THE INITIALIZATION OF THE PROPOSED ALGORITHM
In the initialization phase we proceed as follows: first we initialize the green channel, since in this channel we have more data than in the other ones, and successively we update the other two.
3.1. The Initialization of the Green Channel
We refer to a clique as a pair of adjacent pixels. Every missing value of the green channel is initialized by a weighted mean of the known green values in its neighborhood. The weights of the considered mean take into account possible discontinuities in a set of adjacent cliques. We consider cliques both in the blue and in the red channel, since it is well-known that there is a correlation between the discontinuities in the various channels related to edges, such as object borders and textures . Here we distinguish three cases: the first one is when we have the value of the green light intensity on a pixel; the second one is when we the blue value of the involved pixel is known, that is when i and j are both odd; the third one is when the red value on the considered pixel is known, namely when i and j are both even.
The first approximation g(0) of the green ideal image g is given by
Note that, in the first case, we keep the value we already have. In the second case, we do a weighted mean of the intensity values taken on the adjacent pixels where the green value is known. The weights t1, t2, t3, t4 of the mean are computed by using the green and the blue channels. We define
In particular, we get
where k = 1, 2, 3, 4, ϕ is a suitable positive decreasing detection function and τk is defined by
and the pixels ek, fk, pkand qk are as in (1). When the differences between the green values on the pixels ek and fk (see the yellow arc in Fig. (2) for k = 1), and pk (see the cyan arc in Fig. (2) for k = 1), and between the blue values on the pixels (i, j)and qk (see the brown arc in Fig. (2) for k = 1) are small enough, then we can assume that there are no discontinuities between the pixels ek and (i, j) (see the red line in Fig. (2) for k = 1). So, in the calculus of the green value on the pixel (i, j), we give a large weight t1 to the green value in the position ek. Thus, when the value τk is small, the probability of having a discontinuity between the pixels (i, j) and ek in the green channel is large, and vice versa. The computation of τk is illustrated in Fig. (2). For k = 1, 2, 3 the computation of τk can be described by an appropriately rotated similar figure.
|Fig. (2). Computation of for τk for k = 1.|
Even in the third case, we compute the weighted mean of the intensity values taken on the adjacent pixels where the green value is known. The weights t4+k, k = 1, 2, 3, 4 of the mean are computed by using the green and the red channels. In particular, t4+k = ϕ (τ4+k), where
and ek, fk, pk and qk are as in (1). We argue analogously as in the computation of the weight tk, k = 1, 2, 3, 4, where the role of the blue channel is played by the red component.
3.2. The Initialization of the Red Values
Here we distinguish four cases: the first one is when we already know the red value of a pixel; the second one is when we know the red values in the two adjacent pixels in the same column, that is i is odd and j is even (Fig. 3a); the third one is when we know the red values in the two adjacent pixels in the same row, namely i is even and j is odd (Fig. 3b); the fourth one is when we know the red values of the pixels adjacent in the corners of the involved pixel, that is i and j are both odd (Fig. 3c). In the second and in the third case we equalize the second derivatives of the red and the green channels previously computed. In the last case we use the computed values of the red channel to determine the weights of a suitable mean.
So, we define the initial estimate r(0) of the red ideal image r by
Note that, in the first case, we keep the value which we already have. In the second case, we pose that the finite difference of the second order in the vertical direction of the red channel coincides with that of the green channel, which we have already initialized, namely
Since we know , we can deduce the value of from (3).
|Fig. (3). Different cases in the initialization of the red channel.|
In the third case, we impose that the finite difference of the second order in the horizontal direction of the red channel coincides with that of the green channel, just already initialized, that is
By proceeding analogously as above, we obtain the value of from (4).
In the fourth case, we do a weighted mean of the intensity values taken on the adjacent pixels where the red value has just been computed. The weights t8+k, k = 1, 2, 3, 4, of the mean are calculated by using the just initialized red channel and the observed blue channel. In particular, t8+k is given by ϕ(τ8+k), k = 1, 2, 3, 4, where ϕ is the detection function used in initializing the green channel, and
where ek, fk, pkand qk are as in (1). When the differences between the red values on the pixels ek and fk, ek and pk, and between the blue values on the pixels (i, j) and qk, are sufficiently small, then we can suppose that there are no edges between the pixels (i, j) and ek. So, in the calculus of the red value on the pixel (i, j), we have a large weight t8+k, k = 1, 2, 3, 4, in correspondence with the red value in the position ek.
3.3. The Initialization of the Blue Values
Also in this setting, we distinguish four cases: the first one is given when we know the blue value of a pixel; the second one is when we know the blue values in the two adjacent pixels in the same column, that is i is even and j is odd (Fig. 4a); the third one is when we know the blue values in the two adjacent pixels in the same row, namely i is odd and j is even (Fig. 4b); the fourth one is when we know the blue values of the pixels adjacent in the corners of the involved pixel, that is i and j are both even (Fig. 4c). In the second and third cases we equalize the second derivatives of the blue and the green channels previously calculated. In the last case we use the computed values of the blue channel to determine the weights of a suitable mean.
|Fig. (4). Different cases in the initialization of the blue channel.|
Thus, we define the estimate b(0) of the blue ideal image b by
Note that, in the first case, we keep the value which we already have.
In the second case, analogously as before, we impose
As we know , we derive the value of from (5).
In the third case, similarly as above, we get
By arguing as in the previous section, we deduce the value of from (6).
In the fourth case, we do a weighted mean of the intensity values of the adjacent pixels where the blue value has just been computed. The weights t12+k, k = 1, 2, 3, 4, of the mean are calculated by using the observed red channel and the just initialized blue channel.
Analogously as before, we obtain t12+k = ϕ (τ12+k), where
where ek, fk, pkand qk are as in (1).
4. THE ITERATIVE PHASE OF THE PROPOSED ALGORITHM
A classical filter, often used to solve the demosaicing problem, is the Freeman filter also . The phase described in this section is a suitable modification of this filter. The Freeman filter performes the initialization phase by means of the bilinear filter, which works as follows. When the value of a certain color of a pixel is not available, such a value is computed by the arithmetic mean of the values of that color, which are assumed in the neighborhood of this pixel, that is the bilinear estimation is given as
Moreover, the following values are defined, by means of the median of the color differences of the channels red-green and blue-green :
with . The median turns out to be very useful to correctly preserve the edges which are in the images. Indeed, the median filter is often used to restore images corrupted by salt-and-pepper noise, namely by the noise present only in a few pixels not adjacent each other.
In the Freeman filter it is assumed that the color differences are constant in a suitable subarea. Thus, the Freeman estimation is defined as follows:
In this paper we modify the Freeman filter as follows.
From the initial estimation we define the following variables:
where . So, we define the estimates for s = 1, 2,... as follows:
We pose our final estimate as
We saw experimentally that a good approximation is given by . We call the technique associated to this estimate as Local Edge Preserving (LEP) algorithm.
5. EXPERIMENTAL RESULTS AND DISCUSSION
In this section we present the experimental results obtained from the implementation of the proposed algorithm, which was tested for the Bayer CFA on the set of Kodak sample images , of size, shown in Fig. (5). This dataset represents the typical benchmark images used in the literature to compare the various demosaicing algorithms. These high quality images have been acquired as raw images, in order to minimize the compression. We have implemented our algorithm in the C language on an Ubuntu operating system by means of a computer having a processor of speed of 3.40 GHz.
To define a specific LEP method, we assume that the radius t of the neighborhood of the median filter in the equation (7) is equal to, and we experimentally choose the detection function in (2). In particular, the tested functions are
Observe that the detection functions ϕj, j = 1,...,6 are decreasing and continuous. Moreover, we get ϕj(0) = 2 and .
In Table 1 there are the errors of the LEP algorithm in terms of Mean Square Error MSE, also  in reconstructing the images of the Kodak set as the detection function varies. The values in bold are related to the best reconstruction of a specific image. In the last line there are the means of the MSE obtained in the reconstruction of the Kodak sample images, as the detection function varies. Note that the best result can be obtained by different detection functions, but, if one takes the means, then the minimal error corresponds to the detection function ϕ4. To evaluate whether the function ϕ4 is actually the best detection function, we proceed as follows. For each sample image we give five points to the detection function which allows to obtain an estimate with the minimal error; four points to the detection function which obtain the second best minimal error; three points in correspondence with the third minimal error, and so on.
|Fig. (5). Kodak image set.|
In Table 2 there are the results obtained by the all detection functions on the single images, and in the last line there is the global score. Observe that, even in this case, the highest score is obtained by the detection function ϕ4.
From now on, we use the detection function ϕ4 in the LEP algorithm. In Fig. (6a) the reconstruction of Image 02 is shown. If we evaluate the results visually, it is very difficult to perceive the errors present during the reconstruction. Thus, in Fig. (6b) we present the image of errors, which is given by
|Fig. (6). LEP reconstruction of Figure 02.|
where is the estimate obtained by the LEP algorithm, x is the ideal image and e is the column vector belonging to , whose entries are equal to one. Again, it is difficult to note visually the errors of the algorithm. So, in Fig. (6c) we show the image of the enlarged errors, that is
where it is possible to see in detail the errors of the algorithm.
Since most algorithms existing in the literature do not allow to see easily the errors related to the reconstructions, because they seem to be perfect, then, to compare our algorithm with some of those proposed in the literature, we use the table of the errors in the reconstruction of the Kodak dataset. In Tables 3 and 4, we compare the LEP method with the original Freeman filter  and with some other recently published algorithms [8, 12, 19, 22, 26-33]. Although the proposed algorithm gives the best reconstruction of two images, the total mean of the errors obtained with the LEP algorithm is the smallest of the selected methods.
In the literature there exist several other algorithms, for instance that proposed , which is one of the best performed algorithms, obtaining a MSE mean equal to 6.11. However, in order to reach this goal, the needed mean computation time is equal to 27 minutes and 4 seconds, while the mean computation time for the LEP algorithm is equal to 0.16 seconds. The aim of the LEP algorithm is to obtain good results in a very short period of time. This method can be used as an initialization algorithm for the technique proposed , obtaining meaningful reductions of its computational cost.
In Fig. (7a) the reconstruction of Image 08 by LEP is presented. Its MSE, between the original Image 08 is about 19.99, obtained in a computational time of 0.16 seconds. The relative enlarged error image is presented in Fig. (7b). In Fig. (7c) the reconstruction of Image 08 by the algorithm proposed  is illustrated. Its MSE between the original Image 08 is about 12.33, obtained in a computational time of 27 minutes and 54.74 seconds. The relative enlarged error image is given in Fig. (7d). From the enlarged error images it is possible to notice how the algorithm proposed specially refines the reconstruction of the buildings on the left part of the image , however it does not allow an immediate processing of the image.
|Fig. (7). Reconstruction of Image 08.|
We investigated the demosaicing problem and proposed a heuristic technique, in order to obtain a very fast algorithm. In particular, we proposed an algorithm consisting of three steps. In the first one, the green channel was updated by means of an edge-directed and weighted average technique. In the second one, the red and blue channels were updated, by using also the constraint of equality of the second derivatives in the various channels. In the third step, we proposed an iterative algorithm, assuming the constant color difference. this choice allowed to obtain accurate reconstructions in very short calculation times. Moreover, similarly as in the Freeman technique, in this phase we employed a median filter. We fixed a maximum number of iterative steps as four, in order to obtain low computational costs. We called our algorithm Local Edge Preserving (LEP).
The choice to impose an edge-preserving reconstruction, to impose a correlation between the channels on the derivatives of the second order, and to impose that the difference between the channels is constant, was the result of an extensive experimentation. Such an experimentation showed how this choice allowed to obtain accurate reconstructions in very short calculation times. Thus the proposed technique turns out to be very competitive when compared with some other methods known in the literature also [8, 12, 19, 22, 26-33].
CONSENT FOR PUBLICATION
CONFLICT OF INTEREST
The authors declare no conflict of interest, financial or otherwise.
The first author was partially supported by the G. N. A. M. P. A. (the Italian National Group of Mathematical Analysis, Probability and Applications). The second and the third author were partially supported by the G. N. C. S. (the Italian National Group for the Scientific Computations).
This research was partially supported by the projects “Ricerca di Base 2017” (Metodi di Teoria dell'Approssimazione e di Analisi Reale per problemi di approssimazione ed applicazioni) and “Ricerca di Base 2018” (Metodi di Teoria dell'Approssimazione, Analisi Reale, Analisi Nonlineare e loro applicazioni).
|||Bai C, Lin Z, Yu J, Chen YW, Chen YW. Penrose demosaicking. IEEE Trans Image Process 2015; 24(5): 1672-84.
|||Hamilton JF, Compton JT. Processing color and panchromatic pixels. Patent US 0024879 2007.|
|||Hirakawa K, Wolfe P. Spatio-spectral color filter array design for enhanced image fidelity. Image Processing (ICIP) 2007; 81-4.
|||Bayer BE. Color imaging array. Patent US 3971065 1976.|
|||Hibbard RH. Apparatus and method for adaptively interpolating a full color image utilizing luminance gradients. Patent US 5382976 1995.|
|||Laroche CA, Prescott MA. Apparatus and method for adaptively interpolating a full color image utilizing chrominance gradients. Patent US 537332 1994.|
|||Kakarala R, Baharav Z. Adaptive demosaicking with the principal vector method. IEEE Trans Consum Electron 2002; 48: 932-7.
|||Buades A, Coll B, Morel JM, Sbert C. Self-similarity driven color demosaicking. IEEE Trans Image Process 2009; 18(6): 1192-202.
|||Adams JE. Interactions between color plane interpolation and other image processing functions in electronic photography. Proc SPIE 1995; 2416: 144-51.
|||Freeman W T. Median filter for reconstructing missing color samples. Patent US 4774565 1988.|
|||Kimmel R. Demosaicing: image reconstruction from color CCD samples. IEEE Trans Image Process 1999; 8(9): 1221-8.
|||Pei SC, Tam IK. Effective color interpolation in CCD color filter arrays using signal correlation. IEEE Trans Circ Syst Video Tech 2003; 13: 503-13.
|||Tao B, Tastl I, Cooper T, Blasgen M, Edwards E. Demosaicing using human visual properties and wavelet interpolation filtering. Proceedings of Color Imaging Conference: Color Science, Systems, Applications 1999; 252-6.|
|||Weldy JA. Optimized design for a single-sensor color electronic camera system. Proc SPIE 1988; 1071: 300-7.
|||Lu W, Tan YP. Color filter array demosaicking: New method and performance measures. IEEE Trans Image Process 2003; 12(10): 1194-210.
|||Wu X, Choi WK, Bao P. Color restoration from digital camera data by pattern matching. Proc SPIE 1997; 3018: 12-7.
|||Tsai PS, Acharya T, Ray AK. Adaptive fuzzy color interpolation. J Electron Imaging 2002; 11: 293-305.
|||Adams JE, Hamilton JF. Design of practical color filter array interpolation algorithms for digital cameras. Proc SPIE 1997; 3028: 117-25.
|||Gunturk BK, Altunbasak Y, Mersereau RM. Color plane interpolation using alternating projections. IEEE Trans Image Process 2002; 11(9): 997-1013.
|||Gerace I, Martinelli F, Tonazzini A. Demosaicing of noisy color images through edge-preserving regularization. Proceeding of 2014 International Workshop on Computational Intelligence for Multimedia Understanding (IWCIM) 2014; 1-5.
|||Menon D, Calvagno G. Regularization approaches to demosaicking. IEEE Trans Image Process 2009; 18(10): 2209-20.
|||Kiku D, Monno Y, Tanaka M, Okutomi M. Beyond color difference: Residual interpolation for color image demosaicking. IEEE Trans Image Process 2016; 25(3): 1288-300.
|||Hadamard J. Lectures on cauchy’s Problem in linear partial differential equations Yale 1923.|
|||Kodak lossless true color image suiter. Available on: http:// r0k. us/ graphics/ kodak/|
|||Hosam O. Side-informed image watermarking scheme based on dither modulation in the frequency domain. Open Signal Process J 2013; 5: 1-6.
|||Ferradans S, Bertalmío M, Caselles V. Geometry-based demosaicking. IEEE Trans Image Process 2009; 18(3): 665-70.
|||Hamilton JF, Adams JE. Adaptive color plane interpolation in single sensor color electronic camera. Patent US 5629734 1997.|
|||Hirakawa K, Parks TW. Adaptive homogeneity-directed demosaicing algorithm. IEEE Trans Image Process 2005; 14(3): 360-9.
|||Li X. Demosaicing by successive approximation. IEEE Trans Image Process 2005; 14(3): 370-9.
|||Lukac R, Plataniotis KN. Data-adaptive filters for demosaicing: A framework. IEEE Trans Consum Electron 2005; 51: 560-70.
|||Muresan DD, Parks TW. Demosaicing using optimal recovery. IEEE Trans Image Process 2005; 14(2): 267-78.
|||Sakamoto T, Nakanishi C, Hase T. Software pixel interpolation for digital style camera suitable for a 32-bit MCU. IEEE Trans Consum Electron 1998; 44: 1342-52.
|||Wu X, Zhang N. Primary-consistent soft-decision color demosaicking for digital cameras (patent pending). IEEE Trans Image Process 2004; 13(9): 1263-74.
|||Moghadam AA, Aghagolzadeh M, Kumar M, Radha H. Compressive framework for demosaicing of natural images. IEEE Trans Image Process 2013; 22(6): 2356-71.