BellMoonVassiliev

Back to Psych 221 Projects 2013

Background and Motivation

Most digital image sensors rely on a color filter array (CFA) in order to sense color, so demosaicing is an important and near-universal part of the imaging pipeline. Demosiacing is by nature error-prone, since it is an attempt to fill in missing data points with educated guesses. However, if there are multiple images of the same scene - whether captured as individual still frames or from a video - then color information which was lost in one image may be captured in another. By performing joint demosaicing on a set of images, we are able to more accurately reconstruct the color of the scene than with a single image.

A second research area within image processing involves combining a series of images of the same scene in order to achieve a higher signal-to-noise ratio (SNR), higher resolution, or both. Work on super-resolution has focused primarily on grayscale data, since RGB data adds several challenges. In particular, the image channels must remain aligned through any filtering operations, or else the colors may be distorted and leave unpleasant fringes. When these multi-image operations are performed after demosaicing, errors from the demosaicing process are treated as "signal" in the following steps. While denoising is often done on raw Bayer images, super-resolution cannot easily done before demosaicing, since the result. In this project, we treat demosaicing, super-resolution, and noise reduction as different aspects of the same problem, which is to use all of the information in the source images to produce a single result which most accurately represents the scene as it would have been captured by an ideal camera.

Even though modern image sensors have large pixel count (8MP and up), superresolution remains relevant. It enhances images obtained with digital zoom and adds another dimension of design flexibility: sensor resolution versus computation. Thanks to Moore's law, the energy price of computation is low and continues to fall, but the cost of moving data on and off a chip remains nearly constant constant. Using a sensor with fewer pixels (requiring less data movement) in combination with super resolution techniques may result in a more energy efficient system at the expense of small reduction of image quality. This argument is even stronger for computer vision systems where photographic quality is usually less important.

Method

Our approach consists of an image alignment step, which finds the relative horizontal and vertical offsets between the LR images, and a demosaicing step which combines the data from LR images into a single HR image. The quality of the reconstruction depends on the camera's sensor and optics, and on the relative alignment of the images.

The combination of the sensor and optics determines the optical point spread function, which in turn determines how difficult it is to obtain a sharp image result. Using multiple images does not reduce or remove the blur induced by the optics; it only provides a higher SNR image.

When multiple images are aligned and averaged, the result is as if the image were convolved with a PSF the size of a pixel. For hypothetical pixels with 100% fill factor, this PSF is a square of exactly the pixel size. Conversely, for pixels with 0% fill factor (i.e, sample points with no area), the PSF is a delta function - and no additional blur is induced. The obvious downside is that a 0% fill-factor pixel cannot capture any light.

In both of these cases, non-blind deconvolution can be applied to the image to sharpen it. However, deconvolution is a complex problem and we did not attempt it in this project. Any resolution improvement we obtain is by reconstructing aliased high frequencies, based on the fact that the sensor fill factor is less than 100%.

The relative positions of the LR images are also important to the reconstruction. In the ideal case, four images are spaced exactly one pixel apart, giving red and blue color samples at every point, and two green samples for every point. In the worst case, the images overlap perfectly, which allows us to average pixels together to reduce the SNR, but does not offer any additional color information. A burst of images from a real camera will generally fall somewhere between these two extremes. As long as the image registration returns the correct offsets, a larger number of images always has the potential for a better result, at the expense of additional computation time.

In this project, we assume ideal camera optics. This is the approach taken by most recent super-resolution and demosaicing works. TODO: cite

Image Alignment

Parabola fit to find subpixel offset

To align the images, we use phase correlation directly on the mosiaced data.

Some works low-pass filter the data in order to mask the effect of the Bayer pattern, but we did not find this to be necessary. If the image does not have any strong content, then high frequency induced by the alternating pixel colors will dominate the actual frequencies in the image, and always force a "zero" registration. We found that this was a problem for spatial cross-correlation, but not for phase correlation. Because phase correlation works in the frequency domain, it assumes that the image is circular (i.e, it is continuous from one edge back around to the opposite edge). Since the image is not circular in reality, it can be multiplied by a window function to reduce edge-induced artifacts. However, in our testing, windowing the images did not produce a more accurate registration.

The result of the correlation is a matrix of correlation values, which for a good registration form a single sharp peak. To obtain a subpixel-resolution offset, we fit a parabola to a 7x7 window around the peak, and then solve for the maximum of the parabola. We found that the (rather large) 7x7 window gave more accurate results than 5x5 or 3x3 windows.

We perform correlation between every unique pair of images, which gives us 1+2+...+(n-1) relative offsets. We fix the first image at (0, 0), and then find the image positions which best fit the measured offsets using a least-squares fit. Although our method only measures X and Y translations, it is straightforward to extend the process to handle rotation and scale between images by representing the images with polar coordinates.

Spacial x correlation

• More efficient for HW implementation.
• Works well for small image offsets.
• Size of correlation window is critical - fails if it's too small
• Can you several patches from different places in the image (4 corners + center) to handle images with repeating patterns.

Super resolution

After all the images are aligned we'll combine them to form a single higher resolution image. Depending on the actual offsets the result can be very different. For example if we have zero offset, we can't do anything. Another extreme would be 1 pixel offset which will give us the missing color channel(s) but not the extra resolution. Finally, if we have sub pixel offset, we can get higher resolution for all color channels, but the picture will require demosaicing. It should be noted that because of the Bayer pattern in the input images, any offset multiple to 2px in each direction is irrelevant. Another note is that the combined image will loose Bayer pattern.

Weights function

The first method for combining low resolution images treats each pixel as a point on a higher resolution grid and calculates the value on the grid using weighted average of the data. We used Gaussian distribution as a weight function. Here is how it works:

• Super resolution using Gaussian weight function
• 
  Start with the aligned LR images
• 
  A dd HR grid and weighting function at a point on the new grid
• 
  Weight function and data points in 3D

For every point of the new HR grid, take all close by data from LR images, multiply by the value of the weight function a the coordinate of the data, sum and divide by total weights. This guaranties that the coefficient for all data used in the calculation always sum to one. So that in case that there is only one data point in the vicinity, that value will be taken as is, without any division.

Area averaging

An alternative method with we've explored treats each pixel as a square and relies on the assumption that LR image registration is only accurate to the first digit after the decimal point. Every data point is replicated 100 times to cover 10x10 area on an intermediate grid. Then we combine all the images at every point by simple averaging followed by an interpolation to fill in the missing data. After that we subsample intermediate values using the final HR grid to get the final result.

• Super resolution using area averaging
• 
  Start with the aligned LR images
• 
  Replicate data samples to the 100x grid
• 
  Number of samples per point
• 
  Average and interpolated data at every point
• 
  Interpolated data in 3D
• 
  Result after subsampling using desired HR grid

• Treat data as overlapping squares
• Map to 10x grid in each direction
• Average
• Interpolate
• Subsample according to the desired grid

More efficient for hardware

Comparison of the interpolation methods

Both methods using averaged existing data to estimate values on the grid.

In our final code we used weighed function because it calculates the missing value in a single step and requires less memory. However it requires floating point operations to calculate the value of the weight function and floating point division to normalize the result.

In contrast, the Area averaging approach requires more intermediate memory and two steps (map to 100x grid with averaging and interpolation), but can use all integer math and no complicated division. For this reason, it will be better suited for hardware implementation.

It should should be noted that unlike matlab which allocates memory for the full 100x image array and thus runs slow, in hardware we would be working with a small sliding window of data and thus won't require large amounts of extra memory

Demosaicking

Principles from single-image demosaicing

Since single-image demosacing is a well-studied problem, we began by implementing a simple demosiacing algorithm.

Before directly working with multiple images for demosaicing, we first apply our demosaicing algorithms to a single image to figure out how effective it is within a single-image framework. Through this experiment, we expected to get some ideas of how to extend this algorithms to multiple-frame demosaicing.

Our algorithm is executed in two steps.

1. Interpolate Green channel by using gradient information of Red and Blue channels.

2. Interpolate Red and Blue by using inter-color correlation parameters, Red-to-Green and Blue-to-Green ratio.

Green Interpolation

Blue Interpolation

Since red and blue channels are down-sampled two times more than the green channel, it's reasonable to interpolate the green channel first then red and blue channel with the interpolated green channel. The simplest approach would be to estimate the unknown pixel values by simple linear interpolation. However, this approach will ignore important information about the correlation between color channels and will cause severe color artifacts. Another critical drawback of simple linear interpolation is that uncorrelated pixels across edges will be averaged and will result in blurred edges. Considering this fact, our algorithm uses second derivatives of the red and blue channel to give different weights for the horizontal and vertical interpolation of the green channel.

Weights are computed from the red or blue channel as follows.

Here is the equation for the second derivatives.

(The formula above is described for the green value estimation at blue points. The unknown values for the green channel at red points can be estimated by the same formular with replacing the blue channel gradients with red channel gradients.)

The intuition behind this algorithm is that second gradients in the blue or red channels imply possible existence of an edge so that we could use different weight for the horizontal and vertical direction interpolation in the green channel. The large value of second gradients implies existence of an edge while the small value implies either uniform color region or uniformly varying color region. Considering this fact, it's reasonable to giver larger weight to the direction with smaller gradient and smaller weight to the direction with larger gradient when interpolating unknown green pixels. This method will help prevent blurring across edges. If both gradients in the horizontal and vertical directions are zero, equal weights are used (0.5) to equally interpolate in both directions.

Now that the unknown pixels in the green channel have been interpolated, inter-color correlation could be taken into account for the unknown pixels in the blue and red channels. Here, we assume that there are correlation between the green and the blue (red) component. This assumption makes sense, because the green spectral power density overlaps with the blue and red spectral power density. After this point, the interpolation process will be explained based on the blue channel interpolation, but the same process is applied to the red channel interpolation.

The blue (red) channel interpolation involves the following two steps.

2-1 Filling in the unknown blue (red) values at the pixels on the green bayer pattern.

These pixels have adjacent known blue (red) values so that those values could be used to interpolate the pixel point. Also, the green value at the pixel point could be used to predict the blue (red) value at the point. The following equations take these two aspects into consideration.

The first term is an average of two neighboring blue values, and the second term uses the green value and the blue-to-green ratio to predict the blue value at the current pixel point. The blue-to-green ratio at the current point can be estimated from the blue-to-green ratio at the adjacent points. (If the adjacent blue values are located vertically, the first term of the equation should take two vertical values instead of horizontal ones)

The weight to the inter-color correlation term (the second term) varies across different images, and we got the following numbers through our own experiment on several test images.

2-2 Filling in the unknown blue (red) values at the pixels on the red (blue) bayer pattern.

At these locations, there are no adjacent blue values. While more intelligent algorithms could be applied to interpolate these pixels, we simply apply linear interpolation with the four adjacent blue values that have been interpolated in the previous phase.

Experimental Results

We compare the result of our demosaicing algorithm with the one that the decmosaic function in Matlab produces.

Test images

We use PSNR and S-CIELAB as the image quality metrics. The results shows that our algorithm usually gives better S-CIELAB number, and this motivates extending this algorithm to the multi-image demosaicing frame work.

File:Test images result.png

Matlab code for the single-image demosaicing

demosaic_single.m - Given a reference image, this function creates RGB bayer pattern out of it, apply our algorithm and Matlab demosaic separately to them, and give PSNR, S-CIELAB values.

test_demosaic_single.m - This script call demosaic_single function with five different test images to extract PSNR and S-CIELAB values.

getPSNR - compute PSNR given a reference image and demosaiced image.

getscielab - comput S-CIELAB value given a reference image and demosaiced image.

Generalization to multiple images

When multiple images with arbitrary translations are merged together, the resulting mesh of sample points no longer forms a Bayer pattern. If rotation and scale are also considered, the mesh is not even regular. Thus, we need to generalize the demosiacing methods to handle and arbitrary configuration of sample points.

Failure cases

Images with parallax

There is no good way to handle parallax or scene motion where parts of one frame are not visible in another and vice-versa. Not only is it impossible to find a global registration between the images, it is impossible to find a corresponding match for every pixel in an image. It would be possible to use a number of computer vision methods to segment the image and demosaic parts separately, but such algorithms typically do not have the precision necessary for our demosiacing to run well. A more effective way to handle parallax is to eliminate or reduce it it from the outset by capturing the scene at a very high frame rate, perhaps more than 1000 frames per second. At such speeds, the pixel-wise offsets between frames will be in the single digits, and so the global registration assumption is correct, or nearly correct.

Conclusions

Here is where you say what your results mean.

References

[1] S. Farsiu, D. M. Robinson, M. Elad, and P. Milanfar, "Dynamic demosaicing and color superresolution of video sequences", in Proceedings of SPIE, vol. 5562, 2004, pp. 169-178.

[2] S. Farsiu, M. Elad, and P. Milanfar, "Multiframe demosaicing and super-resolution of color images", Image Processing, IEEE Transactions on, vol. 15, no. 1, pp. 141-159, 2006.

[3] M. Guizar-Sicairos, S. T. Thurman, and J. R. Fienup, "Effcient subpixel image registration algorithms", Optics letters, vol. 33, no. 2, pp. 156-158, 2008.

[4] T. Q. Pham, L. J. Van Vliet, and K. Schutte, "Robust fusion of irregularly sampled data using adaptive normalized convolution", EURASIP Journal on Advances in Signal Processing, vol. 2006, 2006.

[5] H. Takeda, S. Farsiu, and P. Milanfar, "Kernel regression for image processing and reconstruction", Image Processing, IEEE Transactions on, vol. 16, no. 2, pp. 349-366, 2007.

[6] P. Vandewalle, K. Krichane, D. Alleysson, and S. Süsstrunk, "Joint demosaicing and super-resolution imaging from a set of unregistered aliased images", in Electronic Imaging 2007. International Society for Optics and Photonics, 2007, pp. 65 020A-65 020A.

Appendix I - Code and Data

All of our code is written in MATLAB and is available on GitHub: [1]

Test data files are also included in the git repository.

Appendix II - Work partition

TODO