Department of Electrical Engineering, Stanford University
9:15 AM, March 12, 2003
A general image-processing pipeline in a digital camera can be mainly divided into two steps:
Sensor measurements à 1. Spatial Demosaicing è 2. Color Correction à Gamma correction
Depending on the method of each step, the final image quality can significantly vary.
However, image-processing algorithms in a digital camera is embedded in a chip, and we cannot modify them to test or improve the final image quality. Therefore, in this project, I manually applied each image-processing steps in order to understand, and potentially improve, the final image quality of a digital camera.
This report provides intuitive understandings on the design and trade-offs of spatial demosaicing and color correction (or color balancing) techniques, and their effect on the final image quality.
Note: PowerShot G1 from Canon, which has the raw sensor-data mode, was used. The CMYG sensor matrix was composed of approximately 2000 x 1500 pixels. More details are listed at http://www.dpreview.com/reviews/canong1/.
Figure 1 shows Color Filter Array (CFA) pattern for the digital camera used in this project. CFA covers the sensor matrix, and each pixel measures light with a color represented by the filter.
Figure 1. Color Filter Array (CFA) Pattern:
CMGY pattern used in this project,
and typical Bayer pattern for RGB.
The pattern is repeated covering the whole sensor matrix.
The recovery of full-color values at each pixel, i.e., the recovery of full-color image with the same spatial resolution of a pixel, requires a method called color interpolation or color demosaicing algorithms.
General purpose of a color correction matrix (or color balancing matrix), is to convert the measurements from sensor color-space to display color-space, which can be either sensor RGB to display RGB, or sensor CMYG to display RGB (see Fig. 2). The trade-off here will be between noise and color error in the transformed color-space . Here, our color correction matrix Cc will be a 3x4 matrix, which converts sensor CMYG values to display RGB values.
Figure 2. Typical CMYG and phosphor-RGB spectral waveform.
Equation 1. Color Correction or Color Balancing Equation.
( G, M, C, and Y are 2000 x 1500. R, G, and B are 2000 x 1500 )
In Eq. 1, if we think of m x n x l dimension of 3D matrix as vertical x horizontal1 x horizontal2, then the color correction can be thought of as performing linear combinations in vertical direction at each horizontal location, i.e., at each pixel. The spatial demosaicing can be thought of as linear operation (in fact, it is low-pass filtering; details will be discussed later) in the horizontal plane. Therefore, the order of color correction and spatial demosaicing operation can be interchanged.
For presentation purpose, I will first discuss “Step II – Color Correction,” whose results does not require sophisticated demosaicing algorithms to visualize, whereas “Step I – Spatial Demosaicing” requires decent color correction matrix to display the results.
Note: typical gamma value of 2.2 was used to do gamma correction for all results presented here.
To design color correction matrix Cc, we have to know the spectral waveform of the sensor (i.e., response-curve of the sensor CFA) and spectral waveform of display phosphors, similar to the ones shown in Fig. 2. Although, display phosphors may vary, the spectral waveform given in Fig. 2 can be a good approximation compared to the variations in the sensor response-curve among different digital cameras. Figure 3 shows the initial result compared with the direct TIFF output of the digital camera.
Figure 3. Macbeth Color-Checker Board:
(left) direct TIFF output from digital camera,
(right) raw-data color corrected with initial Cc;
Initial Cc. produced poor color contrast.
As shown in above figure, initial Cc, based on assuming simple CMYG response-curve of the sensor, could not be used. Therefore, it was necessary to measure of spectral response of CMYG CFA of the sensor. One common method of approximating the spectral response of the sensor, R, using the color-checker board is to solve below linear equation for R.
(Although we have to characterize the sensor noise and the light source, I focused on measuring the response since we can assume dominated noise as poisson and can use sun light to approximate the 6500K white illuminant D65.)
Equation 2. Linear equation for approximating the sensor response, R, using color patches. S represents spectrometer measurements of 24 color patches, and D represents the digital camera raw measurements of each CMYG sensors. The number of spectral samples was 101 in the visual frequency range.
We measured spectrums of Macbeth Color-Checker Board with 24 color patches on sunny day outdoor between 12 to 2 pm, in order to acquire approximately 6500K white illuminant D65, using our digital camera (Canon PowerShot G1) and the spectrometer that provides 101 spectral samples in the visual frequency range. For the digital camera, white balancing was set to daylight, and we acquired total of 21 pictures (7 pictures with 3 different exposure settings). 4 spectrometer measurements were acquired. Figure 4 shows the average of the measurements shown with the ideal responses from the simulation.
Figure 4. Spectral Responses of Macbeth Color-Checker Board:
(a) simulated spectral responses of ideal color-patches,
(b) actual measurements using spectrometer with the real color-checker,
(c) simulated CMYG sensor measurements with responses shown in Fig. 3,
(d) actual CMYG sensor measurements of the real color-checker board.
Solving Eq. 2 with pseudo-inverse gives a least norm solution. However, in the presence of the noise, we have to consider the condition number of k, which is the measure of the sensitivity of the solution to noise (see Eq. 3).
Equation 3. Condition Number of a linear equation y = Ax
and s(A) = singular value of A.
Large k(A) means the estimation x will be more sensitive to noise in the measurement y.
Note that the condition number is the measure of “geometric anisotropy” (thus k = 1 for an orthogonal matrix).
Our initial condition number k(S) is very large because the 24 color patches were illuminated with the same light, meaning that 24 vectors of 101-dimension were far from being orthogonal to each other. Even in simulations with the ideal S, the estimation of R was very sensitive to the noise (see Fig. 5). The k(S) was 2e3 for the simulation and 1e4 for the actual measurements.
Figure 5. Sensor Response Estimation.
(a) Simulated results with no noise, (b) with noise level of 0.003%, (c) with noise level of 0.03%, and (d) with noise level of 0.3%.
(d) actual results from the measurements;
the worn-down color patches may have produced the similarity in the response (e).
One way to reduce k(S) was to reduce the sample number of spectral wavelengths, i.e., reducing the column of S. Although this provided more robustness to noise in the simulation, no significant improvement was visible the actual estimation even with the spectral sample number of 3. The possible reason can come from the fact that noise level is much higher than the simulation, and that S also contained noise, which is not the case in the simulation.
For the comparison of spatial demosaicing algorithms, better color correction matrix Cc was still needed. The next method I tried was manually increasing the saturation of the colors, i.e., increasing the color contrast, in the intermediate color space, the standard XYZ color space .
CMGY color space of the sensor à XYZ color space à RGB color space of the display
After estimating the white points from one of the color patches, I amplified the distance of each color in the chromaticity plane (Fig. 6) with respect to the white point. (Note that doing this in regular XYZ space may fail because white point vector will have the largest norm, which means that amplifying the distance with respect to white points result in overall shift into one of the corners in the chromaticity plane.) The result is shown in Fig. 7.
Figure 6. Chromaticity Diagram.
Note that this is similar to projection of XYZ space onto X+Y+Z = 1 plane.
Figure 7. The result from Color Correction Matrix Estimation
using XYZ space method, and Least Squares (LS) on the difference method. The color contrast increased but at the same time, noise was amplified.
Even though XYZ space method improved our initial Cc, the poor result suggested for another approach.
Above two methods, solving linear equations from the spectral measurements of color patches and XYZ space method, can be implemented independent of existing image-processing pipelines. However, since we have already working image-processing pipeline that provided TIFF output, we can estimate the error term in Cc.
Equation 4. Linear equation for estimating the error term, DCc,e, which can be solved by simple Least Squares (LS); see Eq. 1.
The result by using the difference from the TIFF output was shown in Fig. 7.
Note that the result was using only one image of the color-checker from the digital camera.
Step I – Spatial Demosaicing
The previous class project by Chen  focused on implementing various demosaicing algorithms with RGB Bayer pattern . Typical RGB demosaicing algorithms can be categorized into two groups: non-adaptive and adaptive. Non-adaptive methods include Nearest Neighbor Replication (NNR), Bilinear Interpolation (BL), Cubic Convolution, Smooth Hue Transition, etc., and adaptive methods include Edge Sensing, Variable Number Gradients, Pattern Recognition, etc. However, the analysis of the spatial demosaicing algorithm has been based on image-domain perspective rather than Fourier domain perspective.
In general, any operation or processing performed can be analyzed in Fourier space.
Any spatial demosaicing or interpolation algorithms can be though of as
low-pass filtering the aliased spectrums, due to undersampling, in the Fourier domain.
Figure 8. 2D Spatial-Frequency Spectrum in Fourier-space:
(left) original image in image-space,
(middle) magnitude image of 2D-FT on fully-sampled image,
(right) magnitude image of 2D-FT on under-sampled image (i.e., Cyan pixels)
with typical low-pass filter (with cut-off frequency of ¼ Nyquist frequency) shown as rectangular box. Horizontal and vertical lines crossing the center are possibly due to line-pattern noise in the sensor.
Most non-adaptive demosaicing methods for RGB Bayer pattern are equivalent to using low-pass filters that has effectively low cut-off frequencies. However, NNR is equivalent to using rect kernel, which has sinc shaped low-pass filter that provides poor stop-band characteristics; moreover, uni-direction NNR has a low cut-off frequency in the direction of replication but no low-pass filtering is performed in the perpendicular direction, which results in severe high spatial frequency aliasing artifacts (see the results in ).
Adaptive method, which uses local spatial information to decide different demosaicing algorithms, is equivalent to using different-sized convolution kernels of low-pass filters at different spatial locations, namely, adaptive filtering.
The advantage of Fourier perspective is that we can directly control the trade-off between noise and spatial-resolutions. Also, various digital-filter-design techniques, which have already been extensively researched, can give us better understandings and more flexibility in designing demosaicing algorithms. The trade-off in computational cost can be regarded as the trade-off between FIR filter length and ringing (or error) in the pass-band of a low-pass filter.
Recently, Fourier domain approach was presented . They use normal low-pass filtering for the luminance channel but apply high-pass filter to use aliased spectral-islands instead of the central spectral-island for chrominance channel (see Fig. 2). However, the central spectral-island contains the information in the aliased spectral-islands because they are the replication of the central islands due to under-sampling.
Instead of converting the RGB-Bayer-pattern demosaicing methods to CMYG pattern in Fig. 1, which has not been not researched extensively, Fourier domain approach was used to solve CMYG demosaicing problem, i.e., solving the spatial demosaicing problem by designing a 2-dimensional low-pass FIR filter with well-known digital filter techniques . I applied the same low-pass filter for all four color (CMGY) channels.
The results from digital filter approach, and from Nearest Neighbor Replication and Bilinear Interpolation methods similar to the one with RGB Bayer pattern are shown in Fig. 9.
Figure 9. Spatial Demosaicing Results I.
Figure 10. Spatial Demosaicing Results II.
FIR filter with Kaiser Bessel window (b = 8) was using cut-off frequency of 0.25 (left figure) or 0.35 (middle figure) Nyquist frequency; note that the maximum cut-off frequency possible is 0.5 Nyquist frequency where Nyquist frequency represents 1/pixelwidth.
Nearest Neighbor Replication method with CMYG pattern cannot be done in uni-direction (see Fig. 1). Therefore, the result is better than the one shown in previous work .
The direct TIFF output looks sharper than FIR filter results (see left figures in Fig. 9 and 10). The possible reason is the sharpening operation in the embedded image-processing pipeline in the digital camera, which amplifies higher frequency components after the demosaicing (i.e., low-pass filtering) process. Note that the possible anti-aliasing filter in front of CFA may have suppressed higher frequency components before the sensor measurements because it is hard for physical anti-aliasing filter to have flat response in the low-frequency pass-band.
Optimal spatial demosaicing and color correction methods are essential for designing image-processing pipelines in digital cameras. As we have seen, depending on the color correction matrix and spatial demosaicing algorithms, the final image quality can significantly vary.
To achieve good color correction matrix, the spectral response information of the sensor is crucial. In the future work, using larger number of color patches or larger number of color checker measurements to reduce the condition number of the linear equation, and/or using brute-force adjustments of 12 = 3x4 elements in the color correction matrix can be done. Moreover, the accurate characterization of the sensor noise may be necessary.
To improve spatial demosaicing algorithms, the digital sharpening filter can be implemented. Moreover, it the computational time is important, one can directly use narrow Kaiser window or narrow Gaussian window as a convolution kernel instead of sinc kernel. The Kaiser or Gaussian window is widely used in various 1D- and 2D-filtering applications due to its controllability of both domains and the property of fast decay.
 Chen,T., “A study of spatial color interpolation algorithms for digital cameras,” Psych221/EE362 project, winter 1999.
 Bayer, B.E., “Color imaging array,” U.S. Patent 3,971,065
 Oppenheim, A.V., and Schafer, R.W., “Discrete-time signal processing,” Prentice Hall.
 Alleysson, D., et al., “Color demosaicing by estimating luminance and opponent chromatic signals in the Fourier domain,” IS&T/SID Tenth Color Imaging Conference, 2002.
 Barnhoefer, U., et al., “Color estimation error trade-offs,” ulrichb AT stanford DOT edu.
 Wandell, B.A., “Foundations of vision,” Sinauer.
MATLAB source codes
do_all.m : main project file.
crw600.exe : dos program for converint CRW raw-files to PPM files.
getraw4.m : reading PPM files.
LSdiff.m : forming linear equation variables for Least Squares on the difference method for color correction.
do_resp.m : extracting spectrometer data and solving Least Norm problem.
segment_patch.m : for region of interest segmentation.
MATLAB data files
whiteX0Z.mat : white patch data
xyz_phos_tutorial.mat : XYZ color-matching functions
new_color_corr.mat : the result of XYZ-space method for color correction.
Please contact me at below email address for other source codes and data-files.
THANK YOU FOR VISITING THIS WEBSITE !!!
Jong Buhm Park (piquant AT mrsrl DOT stanford DOT edu)
Research Assistant / PhD student
Magnetic Resonance Systems Research Laboratory
Department of Electrical Engineering