Gradient-based Feature Extraction From Raw Bayer Pattern Images

Shown in Fig. 1 is an ISP pipeline from Adobe DNG converter [17]. Although the specific algorithms and their orders may vary for different manufacturers, the basic steps in Fig. 1 are usually covered. The details of the functionality of each step is illustrated in Table I.

Demosaicing is a crucial step to convert a single-channel Bayer pattern image to a three-channel color image by interpolating the other two missing color components at each pixel. It has a decisive effect on the final image quality. In order to minimize the color artifacts, sophisticated demosaicing algorithms are always computation hungry.

The problem of demosaicing a Bayer pattern image has been intensively studied in the past decades and a lot of algorithms have been proposed [18, 19, 20, 21]. All these algorithms can be grouped into two categories. The first category considers only the spatial correlation of the pixels and interpolates the missing color components separately using the same color channel. Although these single-channel interpolation algorithms may achieve fairly good results in the low frequency (smooth) regions, they always fail in the high-frequency regions, especially in the areas with rich texture information or along the edges[2].

To improve the demosaicing performance, the other category of algorithms takes the nature of the images’ high spectral inter-channel correlation into account. Almost all these algorithms are based on either the color ratio constancy assumption [20] or the color difference constancy assumption [21]. According to the color image model in [20], which is a result of viewing Lambertian non-flat surface patches, the three color channels can be expressed as

where $\rho$ is the reflection coefficient, $\overrightarrow{N}(x,y)$ is the surface’s normal vector at location $(x,y)$, $\overrightarrow{l}$ is the incident light vector, $I(x,y)$ is the intensity at location $(x,y)$ and $k\subseteq\left\{R,G,B\right\}$ indicates one of the three channels. Note that a Lambertian surface is equally bright from all viewing directions and does not absorb any incident light[22].

At a given pixel location, the ratio of any two color components, denoted by $k$ and $k^{\prime}$, is given by

Suppose that objects are made up of one single material, i.e., the reflection coefficient $\rho$ for each channel is a constant, the ratio of $\rho_{k}\left(x,y\right)/\rho_{{k}^{\prime}}\left(x,y\right)$ reduces to a constant, such that (2) can be simplified as

Equation (3) is referred to as color ratio constancy. In the same manner, the color difference constancy assumption is given by

Note that the direction and amplitude of the incident light are assumed to be locally constant, such that the color component difference $C(x,y)$ is also a constant within a neighborhood of $(x,y)$[2].

The color ratio and difference constancy assumptions are widely used in various demosaicing algorithms [23]. In practical applications, the color difference constancy assumption always is preferred due to its superior peak signal to noise ratio (PSNR) performance. As will be shown, in this work, the color difference constancy can be utilized to directly extract the gradient information from the Bayer pattern images.

In the past decades, many different feature descriptors such as Harr-like features [24], LBP [25], SIFT [26] and histograms of oriented gradients (HOG) [27] have been proposed for object detection. In this work, we mainly focus on the central difference gradient-based feature descriptors, study their applicability on Bayer pattern images and analyze the corresponding performances. Without loss of generality, HOG and SIFT are taken as examples in the analysis and experiments. The results can be extended to other descriptors, such as SURF [28], Color-SIFT [29], Affine-SIFT [30] and F-HOG [31], as long as the central difference is used for gradient computation.

SIFT is a local feature descriptor which detects key points in images. The computation of SIFT can be divided into five steps[32] as

1. 1.

   Scale space construction. The scale space is approximated by the difference-of-Gaussian (DoG) pyramid, which is computed as

   $$\begin{split}D\left(x,y,\sigma\right)&\!=\left(G\left(x,y,l^{i}\sigma\right)\!-\!G\left(x,y,l^{i-1}\sigma\right)\right)\!\ast\!I(x,y)\\ &\!=\!L\left(x,y,l^{i}\sigma\right)\!-\!L\left(x,y,l^{i-1}\sigma\right).\end{split}$$

   (5)

   Here, $G\left(x,y,\sigma\right)=\frac{1}{2\pi\sigma^{2}}e^{-\frac{(x^{2}+y^{2})}{2\sigma^{2}}}$ is the Gaussian function, $l=2^{\frac{1}{s}}$ is a constant multiplicative factor whose value is determined by the number of scales $s$, $\ast$ denotes the convolution operator, $i$ indicates the $i$-th layer in DoG pyramid and $L\left(x,y,l^{i}\sigma\right)$ is the convolution of the original image with the Gaussion function at scale $l^{i}\sigma$.

2. 2.

   Extremum detection. To detect the local maxima and minima by comparing each pixel with its neighbors in a $3\times 3$ neighbourhood among the current scale, scale above and scale below.

3. 3.

   Key point localization. To perform a refinement of key point candidates identified in the previous step. The unstable key points such as points with low contrast or poorly localized along an edge are rejected.

4. 4.

   Orientation determination. To assign one or more orientations to each key point. A histogram is created for a region centered on the key point with radius of $3\sigma_{0}$, where $\sigma_{0}$ is 1.5 times that of the scale of the key point. The direction with the highest bar in the histogram is regarded as the dominant direction and directions with heights of larger than 80% of the highest bar is regarded as the auxiliary directions.

5. 5.

   Key point description. To construct a descriptor vector for each key point. A gradient histogram with 8 bins is created for each $16\times 16$ pixel region around the key point. The key point descriptor is constructed by concatenating the histograms of a set of $4\times 4$ regions around the key point.

HOG is a feature descriptor initially proposed for pedestrian detection [27]. It counts the number of occurrences of gradient orientation in a detection window. The key steps of HOG feature generation are similar to steps 4 and 5 in the SIFT descriptor. The main difference is that orientation histograms in HOG are usually computed on an $8\times 8$ cell and summarized as a global feature by a sliding window.

Image gradient measures the change of intensity in specific directions. Mathematically, for a two-dimensional function $f(x,y)$, the gradients can be computed by the derivatives with respect to $x$ and $y$. For a digital image where $x$ and $y$ are discrete values, the derivatives can be approximated by finite differences.

There are different ways to define the difference of a digital image, as long as the following three conditions are satisfied: (i) zero in constant intensity area; (ii) non-zero along the ramps and (iii) nonzero at the onset of an intensity step or ramp [33]. One of the most commonly used image gradient computation is the central difference based approach as

Here $I(x,y)$ is the intensity at location $(x,y)$, $G_{x}$ and $G_{y}$ represent the gradients in the horizontal and vertical directions, respectively. The computation of (6) and (7) can be implemented by the convolution of the templates in Fig. 2 with the images.

The fundamental idea of the proposed Bayer pattern image based gradient extraction is illustrated in Fig. 3. Instead of demosaicing the Bayer pattern images before difference computation as shown in Fig. 3, we propose to take advantage of the color difference constancy assumption directly for gradient extraction based on Bayer pattern images. Note that by convolving the filter templates in Fig. 2 directly with a Bayer pattern image, all the three conditions for a valid difference definition mentioned are satisfied. To illustrate this, let us consider the example in Fig. 4.

As we can see, the two input pixels for coefficients 1 and -1 in the convolution templates are from the same channel, i.e., differences are always computed on homogeneous pixels. As shown in Fig. 4, applying the convolution templates at locations $(1,2)$ generates

where $G^{B}$ and $G^{R}$ are the gradients of the blue and red channels, respectively.

In the demosaicing tasks, it is a common practice to interpolate the G channel first followed by the R/B channels. This is because there are twice as many G channel pixels as R/B channel pixels in Bayer pattern images. The color difference constancy assumption in (4) can then be used to estimate the missing pixels of the R and B channels.

Here, $k$ represents either R or B channel, $C^{k}(x,y)$ is the difference between the R/B channel and the G channel at pixel location $(x,y)$, which needs to be estimated in demosaicing tasks[34].

Consider two pixels within a small neighborhood at locations $(x,y)$ and $(x^{\prime},y^{\prime})$, according to (10), we have

Where $\delta^{k}(x,y,x^{\prime},y^{\prime})=C^{k}(x,y)-C^{k}(x^{\prime},y^{\prime})$. The value of $\delta^{k}(x,y,x^{\prime},y^{\prime})$ is crucial in our analysis and will be discussed in detail.

Generally, there are flat areas (e.g. background) and texture areas (e.g., corners and edges) in a natural image, these two situations will be discussed separately.

For the flat areas, the difference between two pixels is negligible such that

i.e., $\delta^{k}(x,y,x^{\prime},y^{\prime})\approx 0$. This means the intensity difference between channels is approximately constant across nearby pixel locations, i.e., the color changes are small in the neighborhood of flat areas.

Importantly, a constant $C^{k}$ also means that some non-smooth color transitions (i.e., texture) are included as well. For example, for the two synthetic images in the top row of Fig. 5, suppose $(x,y)$ is a point in the background and $(x^{\prime},y^{\prime})$ is another point in the foreground. These two images correspond to the situation of $C^{k}(x,y)=C^{k}(x^{\prime},y^{\prime})$. Fig. 5-5 are the difference images of $G-R$ and $G-B$, respectively, while Fig. 5-5 are the corresponding gradient maps. Note that all the gradient maps and difference images are displayed as inverse images (1 $-$ original gray value), where 1 means difference or gradient is zero and the the corresponding location is displayed as white (likewise for Fig. 6 and 12). Then,

1. 1.

   Fig. 5 top-left: $I^{R}=I^{G}=I^{B}$ for both background and foreground. This results in $C^{k}(x,y)=0$, which further leads to $\delta^{k}(x,y,x^{\prime},y^{\prime})=0$, as shown in Fig. 5-5, top-left images.

2. 2.

   Fig. 5 top-right: $I^{R}=I^{B}$ for both background and foreground. This results in a constant $C^{k}(x,y)$ for both background and foreground, which further leads to $\delta^{k}(x,y,x^{\prime},y^{\prime})=0$, as shown in Fig. 5-5, top-right images.

For the above two cases, although there are obvious edges in the original images, we still have $\delta^{k}(x,y,x^{\prime},y^{\prime})=0$.

For the more extreme texture areas, the analysis is more complex. To analyze the areas with complex textures, (11) can be further rewritten as

Note that image’s gradients are always computed among a small neighborhood. Considering the central difference-based horizontal gradient computation at pixel location $(x,y)$, we have

Here, $G-k$ represents the difference image of the G channel and the R/B channel. It can be observed from (LABEL:equ:delta_gradient) that $\delta^{k}(x+1,y,x-1,y)$ is exactly the gradient of the difference image at location $(x,y)$. It has been shown in [35] that the difference images are slowly-varying over a spatial domain, meaning that the gradient $G_{x}^{G-k}(x,y)$ in (LABEL:equ:delta_gradient) is negligible, i.e., $\delta^{k}(x+1,y,x-1,y)$ approximates to zero. This can be illustrated by the bottom-left image in Fig. 5. Suppose $(x,y)$ is a point on the background border such that $(x+1,y)$ and $(x-1,y)$ are two points in the foreground and background, respectively. In this case, the following relationship holds (5, bottom-left).

Here, $Sgn(\cdot)$ is the signum function. Let $k=B$, the results of (LABEL:equ:delta_gradient) is

The corresponding result is illustrated in Fig. 5 bottom-left, where the edge is negligible.

Note that there are also failure cases. For example, the bottom-right image in Fig. 5 satisfies

In this case, the result of (LABEL:equ:delta_gradient) will be $\delta^{k}(x+1,y,x-1,y)=1$, which is illustrated in Fig. 5 and 5, bottom-right. Details of failure cases will be discussion in Section IV-C1.

Fig. 6 illustrates two situations: image with dull colors (top) and image with high saturation colors (bottom). Fig. 6 and Fig. 6 are the difference images of G channel $-$ R channel and G channel $-$ B channel ($C^{k}(x,y)$ in Eq. 2), respectively, while Fig. 6 and Fig. 6 are the corresponding gradient magnitude maps of the difference images $(\delta^{k}(x+1,y,x-1,y)$ in (LABEL:equ:delta_gradient)) computed using the central difference operator. The corresponding gray-level histograms of these images are shown in Fig. 7-7. It can be found that for images with dull colors, the differences between G and R channels are distributed in a much smaller range than that of high saturation color images, while the differences between G and B channels are distributed in a larger range. For all these images, most of $\delta^{k}(x+1,y,x-1,y)$ are distributed in a small range, as shown in Fig. 7 and 7. The distribution in Fig. 7 bottom is wider than the other three plots, which are caused by texture areas such as hat and hair edges in Fig. 6. As we can see, apart from these small exceptions, Fig. 6 and Fig. 6 are almost all white. Therefore, for most cases, $\delta^{k}(x,y,x^{\prime},y^{\prime})$ is small and can be ignored if pixel locations $(x,y)$ and $(x^{\prime},y^{\prime})$ are within a small neighborhood.

As a result of the above discussion, (11) can be rewritten as

Combining the gradient definition of (6) and (7) with (18) and (19), we have

meaning that the gradients of natural images can be computed using any one of the three channels as long as the color difference constancy holds. Combining (20) with (8) and (9), we have

Therefore, even though two color components are missing at each pixel, the gradients of location $(1,2)$ can be computed directly from the Bayer pattern image using the blue and red channel. The gradients of any other pixel locations can be computed in the same manner.

Generally, the above conclusion can be extended to other symmetrical first-order differential operators (with alternating zero and nonzero coefficients) on any kind of Bayer pattern. Let us take the Sobel operators in Fig. 2 as an example. Applying the Sobel operators in Fig. 2 to the pixel location $(1,2)$ of the Bayer pattern image results

As for gradient computation using (6) and (7), differences are always computed on homogeneous pixels for Sobel-based differential operations in (23) and (24), i.e., pixel values are always subtracted from pixel values of the same channel. Moreover, according to the color difference constancy assumption, (23) and (24) can be rewritten as

where $\widehat{R}$, $\widehat{G}$ and $\widehat{B}$ represent the missing color components at the corresponding locations. Therefore, the Sobel-based gradients can also be extracted directly from the Bayer pattern images as long as the color difference constancy holds.

In terms of different Bayer patterns, they are merely different arrangements of the RGB pixels, while the alternating pattern of $R$, $G$ and $B$ at each row and column are preserved. For example, discarding the first column of the Bayer pattern in Fig. 4 generates the GRBG Bayer pattern. Therefore, different Bayer patterns do not have any impact on the applicability of the discussed differential operators to Bayer pattern images. Moreover, the discussed gradient extraction method can be directly extended to other special CFA patterns with alternating color filter arrangements, e.g., RYYB, RGB-IR, as long as the arrangement of CFA patterns matches the gradient operators such that gradient operations are performed on the same color channel, i.e., subtract or add operations are performed on the same color channel, and the coefficients of the subtract or add terms in the gradient operator are equal such that the gradients compute from R/B channel can be approximated to G channel.

To validate the proposed Bayer pattern image-based gradient extraction, the differential operators in Fig. 2 are applied to true color images Kodim17 and Kodim04 from the Kodak image dataset[36] and the corresponding resampled Bayer version. The generated gradient maps are shown in Fig. 6. For display purpose, images in Fig. 6 are shown as color images while the gradient magnitude maps in Fig. 6 and 6 are computed from the corresponding gray-scale images generated using Fig. 6. The Bayer pattern images in Fig. 6 are presented as three-channel images to illustrate its Bayer “mosaic” structure. For the clearness of presentation, all the gradient maps and difference images in Fig. 6 are displayed as inverse images. As illustrated in Fig. 6, the gradient maps generated from the Bayer pattern images look almost the same as that generated from the true color version. To compare these gradient maps, two GMS maps are presented in Fig. 6 and 6, and the corresponding distributions of GMS values are presented in Fig.7 and 7. As can be seen, the GMS maps are almost pure white, and the histograms are distributed in a small range, meaning that the compared gradient maps are very close to each other. Overall, gradients of Bayer images yield a good approximation of the image gradients, except for a few pixels around certain color edges.

In SIFT, the scale-space is approximated by a DoG pyramid. The construction of the DoG pyramid can be divided into two parts: Gaussian blurring at different scales and resizing of the blurred images. Due to the special alternating pixel arrange of Bayer pattern images, directly Gaussian blur the images will destroy the “mosaic structure”. This phenomenon is illustrated in Fig. 8. If the Bayer pattern image is directly Gaussian filtered, the resulting image (after demosaicing) looks like a “three-channel grayscale image” as shown in Fig. 8, meaning that the Bayer image is treated as a single channel image, ignoring the “mosaic structure” from the beginning and, effectively, loosing/destroying the color information. Thus, smoothing on Bayer pattern image directly will blend the channels, which is akin to a RGB-to-gray conversion. Moreover, loss of the color information makes some of the algorithms in the SIFT family such as “C-SIFT” and “RGB-SIFT” no longer applicable.

To address the above mentioned problem, the super-pixel approach as illustrated in Fig. 8 is used in this work. A super-pixel is a compound pixel consisting of a complete Bayer pattern. The Bayer pattern image can therefore be regarded as a “continuous” image filling with super-pixels. Operating on the super-pixel structure preserves the Bayer pattern of the original images. Fig. 8 shows the Gaussian blurred image (after demosaicing) based on the super-pixel structure. As can be seen, it is close to that generated by the full color approach. Moreover, the super-pixel structure can also be used for resizing when constructing the scale space. The detailed comparison results will be presented in Section IV-C.

There are five datasets used for differents experiments in this work. Among these five datasets, four are commonly used benchmarks in different image processing and computer vision tasks such as demosaicing, pedestrian detection. A brief description of these datasets is presented in Table II.