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 [4]. 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 [5], 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 [6], 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 [7], to determine the edges in the various channels, a suitable Jacobian operator is applied. In [8], 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 [11] 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 [15] 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 [16], where the interpolation can be done also in the diagonal direction, while in [17] the weighted directional interpolation is used by means of a fuzzy membership assignment. A second order operator is employed as a correction term. [18]

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) [19], 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 [22]. 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 [20]. In particular, a constraint is imposed, to make the derivatives of all channels similar as soon as possible [20]. 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 [10]. 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.

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 [23]. 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.

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 [11]. 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⁽⁰⁾ 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 t₁, t₂, t₃, t₄ of the mean are computed by using the green and the blue channels. We define

(1)

In particular, we get

(2)

where k = 1, 2, 3, 4, ϕ is a suitable positive decreasing detection function and τ_k is defined by

and the pixels e_k, f_k, p_kand q_k are as in (1). When the differences between the green values on the pixels e_k and f_k (see the yellow arc in Fig. (2) for k = 1), and p_k (see the cyan arc in Fig. (2) for k = 1), and between the blue values on the pixels (i, j)and q_k (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 e_k 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 t₁ to the green value in the position e_k. Thus, when the value τ_k is small, the probability of having a discontinuity between the pixels (i, j) and e_k 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 t_4+k, k = 1, 2, 3, 4 of the mean are computed by using the green and the red channels. In particular, t_4+k = ϕ (τ_4+k), where

and e_k, f_k, p_k and q_k are as in (1). We argue analogously as in the computation of the weight t_k, 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⁽⁰⁾ 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

(3)

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

(4)

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 t_8+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, t_8+k is given by ϕ(τ_8+k), k = 1, 2, 3, 4, where ϕ is the detection function used in initializing the green channel, and

where e_k, f_k, p_kand q_k are as in (1). When the differences between the red values on the pixels e_k and f_k, e_k and p_k, and between the blue values on the pixels (i, j) and q_k, are sufficiently small, then we can suppose that there are no edges between the pixels (i, j) and e_k. So, in the calculus of the red value on the pixel (i, j), we have a large weight t_8+k, k = 1, 2, 3, 4, in correspondence with the red value in the position e_k.

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⁽⁰⁾ 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

(5)

As we know , we derive the value of from (5).

In the third case, similarly as above, we get

(6)

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 t_12+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 t_12+k = ϕ (τ_12+k), where

where e_k, f_k, p_kand q_k are as in (1).

A classical filter, often used to solve the demosaicing problem, is the Freeman filter also [10]. 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 [10]:

where

(7)

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.

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 [24], 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 [20] 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 ϕ₄. To evaluate whether the function ϕ₄ 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 ϕ₄.

From now on, we use the detection function ϕ₄ 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 [10] 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 [20], 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 [20], 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 [20] 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 [20], 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].