Digital Image Processing
Digital Image Processing ECE 6258
Popular in Course
verified elite notetaker
Popular in ELECTRICAL AND COMPUTER ENGINEERING
This 0 page Class Notes was uploaded by Cassidy Effertz on Monday November 2, 2015. The Class Notes belongs to ECE 6258 at Georgia Institute of Technology - Main Campus taught by Yucel Altunbasak in Fall. Since its upload, it has received 9 views. For similar materials see /class/233873/ece-6258-georgia-institute-of-technology-main-campus in ELECTRICAL AND COMPUTER ENGINEERING at Georgia Institute of Technology - Main Campus.
Reviews for Digital Image Processing
Report this Material
What is Karma?
Karma is the currency of StudySoup.
Date Created: 11/02/15
Demosaicking Color Filter Array Interpolation in SingleChip Digital Cameras B K Gunturk J Glotzbach Y Altunbasak R W Schafer and R M Mersereau Center for Signal and Image Processing Georgia Institute of Technology Atlanta7 Georgia7 30332 0250 bahadir glotz7 yucel7 rws7 rmm ecegatechedu Draft for the lEEE SPM Special Issue on Color Image Processing 1 Introduction Digital cameras have become popular and many people are choosing to take their pictures with digital cameras instead of lm cameras When a digital image is recorded7 the camera needs to perform a signi cant amount of processing to provide the user a viewable image This processing includes white balance adjustment7 gamma correction7 compression and more reference Color Image Processing Pipeline in Digital Still Cameras77 in this issue A very important part of this image processing chain is color lter array interpolation or demosaicking A color image requires at least three color samples at each pixel location Computer images often use red7 green7 and blue A camera would need three separate sensors to make these measurements To reduce size and cost7 many cameras use a single sensor array with a color lter array The color lter array allows only one part of the spectrum to pass to the sensor so that only one color is measured at each pixel This means that the camera must estimate the missing two color values at each pixel This process is known as demosaicking Several patterns exist for the lter array The most common array is the Bayer color lter array7 shown in Figure 1 The Bayer array measures the green image on a quincunx grid and the red and blue images on rectangular grids The green image is measured at a higher Figure 1 Bayer color lter array arrangement a b Figure 2 Bicubic interpolation used for color lter array interpolation results in numerous artifacts a Original image b Bicubic interpolation sampling rate because the peak sensitivity of the human visual system lies in the medium wavelengths7 corresponding to the green portion of the spectrum Although other patterns exist some using CMYG instead of RGB7 this article discusses the demosaicking problem with reference to this color lter array If the measured image is divided by measured color into three separate images7 this problem looks like a typical image interpolation problem Therefore7 one might try to apply standard image interpolation techniques Bi cubic interpolation is a common image interpolation tech nique that produces good interpolation results when applied to grayscale images However7 when bi cubic interpolation is used for this problem7 the resulting image shows many visible artifacts This is illustrated in Figure 2 This result motivates the need to nd a specialized algorithm for the demosaicking problem Bi cubic interpolation and other standard interpolation techniques treat the color image as three independent images However7 the three images in a color image are generally highly correlated Many algorithms have been published suggesting how to use this correlation This article surveys many of these algorithms and discusses the results in terms of objective and subjective measures 2 Demosaicking Methods We examine demosaicking methods in three groups The rst group consists of heuristic methods the second group uses restoration techniques and the last group is based on image formation modeling 21 Group I Heuristic Approaches Heuristic approaches do not try to solve a mathematically de ned problem7 but they are based on reasonable assumptions about color images Heuristic approaches are spatially adaptive7 and they may exploit correlation among the color channels We examine the heuristic approaches in six categories 211 EdgeDirected Interpolation Although non adaptive algorithms eg7 bilinear interpolation7 bi cubic interpolation can provide satisfactory results in smooth regions of an image7 they usually fail in textured regions and edges Edge directed interpolation is an adaptive approach7 where edge detection is performed for each pixel in question7 and interpolation is done along the edges rather than across them In the demosaicking problem7 edge directed interpolation is usually applied to the green channel7 which contains most of the spatial information andithereforeiis sampled more densely than the red and blue channels A simple way of performing edge detection is to compare the absolute difference among the neighboring pixels Referring to Figure 37 horizontal and vertical gradients at a missing green location can be calculated from the horizontally and vertically adjacent green pixels If the horizontal gradient is larger than the vertical gradient7 suggesting a possible edge in the horizontal direction7 interpolation is performed along the vertical direction If the vertical gradient is larger than the horizontal gradient7 interpolation is performed only in the horizontal direction When the horizontal and vertical gradients are equal7 the green value is obtained by averaging its four neighbors 1 Calculate horizontal gradient AH 1G2 7 G41 2 Calculate vertical gradient AV 1G1 7 G51 3 IfAH gt AV G3 7 G1 G52 Else ifAH lt AV G3 7 G2 G42 e G3 7 G1 G5 G2 G44 Figure 3 Edge directed interpolation in 1 is illustrated G1 G2 G4 and G5 are measured green values G3 is the estimated green value at pixel 3 1 Calculate horizontal gradient AH 1 R3 R72 7 R5 1 2 Calculate vertical gradient AV 1 R1 R92 7 R5 IfAH gt AV G5 7 G2 G82 Else ifAH lt AV G5 7 G4 G62 U e G5 7 G2 G8 G4 G 4 Figure 4 Edge directed interpolation in 2 is illustrated for estimating the green G value at pixel 5 The red R values are used to determine the edge direction When the missing green pixel is at a blue pixel the blue values are used to determine the edge direction It is also possible to compare the gradients against a predetermined threshold value The edge directed interpolation approach in 1 can be modi ed by using larger regions around the pixel in question with more complex predictors and by exploiting the texture similarity in different color channels In 2 the red and blue channels in the 5 gtlt 5 neighbor hood of the missing pixel are used instead of the green channel to determine the gradients In order to determine the horizontal and vertical gradients at a blue red sample second order derivatives of blue red values are computed in the corresponding direction This algorithm is illustrated in Figure 4 Once the missing samples of the green channel are found the red and blue channels are interpolated A typical approach for the redblue interpolation is constant hue based inter polation which is explained in the next section 212 ConstantHue Based Interpolation One commonly used assumption in demosaicking is that the hue color ratios within an object in an image is constant Although this is an oversimpli cation of image formation it is a decent assumption within small neighborhoods of an image This perfect inter channel correlation assumption is sometimes formulated such that the color differences within small neighborhoods are constant This constant color ratio or difference assumption prevents abrupt changes in color intensities and has been extensively used for the interpolation of the chrominance red and blue channels 3 4 5 2 6 7 The demosaicking algorithms that are based on this assumption are called constant hue based interpolation or smooth hue transition methods As a rst step these algorithms interpolate the luminance green channel which is done using bilinear or edge directed interpolation The chrominance red and blue channels are then estimated from the interpolated red hue red to green ratio and blue hue blue to green ratio To be more explicit the interpolated red hue and blue hue values are multiplied by the green value to determine the missing red and blue values at a particular pixel location The hues can be interpolated with any method bilinear bi cubic edge directed etc Instead of interpolating the color ratios it is also possible to interpolate the color differences or the logarithm of the color ratios The algorithm based on color difference interpolation is illustrated in Figure 5 The constant hue based interpolation is sometimes used with median ltering to reduce color artifacts For example in 6 the green channel is interpolated bilinearly and then the missing red and blue pixels are found from the median ltered red green difference and blue green difference Red a I Interpolated A Red Green Interpolate E E Figure 5 Constant difference based interpolation is illustrated 213 Weighted Sum In the edge directed interpolation the edge direction is found rst and then the missing sample is estimated by interpolating along the edge This is a hard decision process Instead the likelihood of an edge in a certain direction can be found and the interpolation can be done based on the edge likelihoods Such an algorithm was proposed by Kirnrnel in The algorithm de nes edge indicators in several directions as measures of edge likelihood in those directions and determines a missing pixel intensity as a weighted sum of its neighbors If the likelihood of an edge crossing in a particular direction is high the edge indicator returns a small value which results in less contribution from the neighboring pixel that direction The algorithm for one dirnensional signals is illustrated in Figure 6 The green channel is interpolated rst the red and blue channels are interpolated from the redgreen and bluegreen ratios The color channels are then updated iteratively to obey the color ratio rule The extension to two dirnensional images is straightforward A similar algorithm was proposed recently in 8 where edge indicators are determined in a 7 gtlt 7 window for the green and a 5 gtlt 5 window for the redblue channels In this case the edge indicator function is based on the L1 norrn absolute difference as opposed to the L2 norm of De ne 51F Sz391 Sz39 1 Object 1 Object 2 0 Interpolate the green at the missing locations e GO l ef GO l R R 9H ex1 60 0 Repeat for three times 0 Interpolate the red using the ratio rule 851 Rz 1 a Rz 1 Ga Gz391 l G G 9 ex1 Rz Gz 1 0 Correct the green to t the ratio rule er GO 1 GO 1 Rz 1 Rz 1 R 8 R ex1 G139Ri I 1 Figure 6 7 is illustrated for a one dimensional signal S is a generic symbol for red R and green dis is the gradient for channel S at location 239 and 615 is the corresponding edge indicator Calculate horizontal gradient AH G4 7 G6 RS 7 16 RS 7 R7 Calculate vertical gradient AV G2 7 G8 RS R1 RS 7 R9 IfAH gt AV G5 7 G2 G82 RS 7 R1 R5 7 R94 Else ifAH lt AV G5 7 G4 G62 RS 7 R3 R5 7 R74 Else G5G2G8 G4 G64R57R1R57R9R57R3R57R78 9 Figure 7 The Adams Hamilton method 9 is illustrated for estimating the green G value at pixel 5 The red R and green values are used to determine the edge direction and to estimate the missing value When the missing green pixel is at a blue pixel7 the blue and green values are used 214 SecondOrder Gradients As Correction Terms ln 97 Adams and Hamilton used edge directed interpolation for the green image This result was improved by using correction terms from the red and blue samples They computed the Laplacian for the red or blue samples along the interpolation row or column7 and used this to correct the simple averaging interpolation This correction term reduces aliasing and improves the mid frequency response for the green image Figure 7 illustrates this algorithm Assume that the horizontal direction has been chosen and consider the interpolation along a single row7 simplifying our analysis to one dimension The block diagram is shown in Figure 8 The output of this system is given by C3 Gs2H12 RsZH22 1 with the sampling equations given by G42 7 12 Gz12 672 2 1292 12 122 712 372 3 The aliasing terms have opposite signs because of the phase shift in the original sampling grid The lter kernel of H1z is 12 1 12 and the kernel of H2z is 714 0 12 0 7 14 The responses of these lters are shown in Figure 9 Because the response of H1 has a wide transition band7 components in the medium frequency range will be attenuated and the spectral copies of these components from the G7z term are passed into the output image causing aliasing distortions Figure 8 Block diagram for Adams Hamilton green interpolation Magnitude response Radian frequency w Figure 9 Filter responses for Adams Hamilton method solid H1z dotted The red lter H2 corrects this Because of the high correlation between channels it is assumed that the components that cause the distortion in the green image are also present in the red image The bandpass lter H2 isolates components in the transition band of H1 The expression for the aliasing terms 12 G72H12712 R72H22 shows that when the outputs of the lters are combined the aliasing distortions can be cancelled if H1z and H2z have the same magnitude at frequencies greater than 7r2 215 HomogeneityDirected Interpolation Instead of choosing the interpolation direction based on edge indicators it is possible to use different measures and to impose different constraints In 10 local homogeneity is used as an indicator to choose between horizontally and vertically interpolated intensities The homogeneity directed interpolation imposes the similarity ofthe luminance and chrominance values within small neighborhoods Referring to Figure 10 the RGB data is rst interpolated horizontally and vertically The green channel is interpolated using red and blue data as correction terms as in The red and blue channels are interpolated from the interpolated Horizontal RGB t0 Interpolation L a b Homo eneit Observed g y I Lab Fmal data based to RGB image selection Vertical RGB t0 Interpolation L a b Figure 10 Block diagram of the homogeneity directed interpolation in 10 red green difference and blue green difference as shown in Figure 5 The interpolated images are then converted to the CIELab space In the CIELab space7 horizontally or vertically interpolated intensities are chosen based on the local homogeneity To measure the local homogeneity7 the number of pixels that have similar luminance and chrominance values to the horizontal and vertical candidates are found and the candidate that is more similar to its neighborhood is chosen 216 Pattern Matching Several algorithms attempt to nd a pattern in the data or t the data to one of several templates A different interpolator is applied for each template This allows a different method to be used for edges and smooth regions In 3 Cok describes a pattern matching algorithm to be used on the green image Each missing green value is classi ed as a stripe7 edge7 or corner7 corresponding to the features ex pected to be found in natural images The human visual system is sensitive to these features so interpolating them correctly is important After classifying the pixel7 the appropriate interpolator is applied to estimate the missing value Wu et al uses the green image to detect patterns in the image in 11 The algorithm determines which directions have the least amount of change in the green image A direction with a large amount of change indicates that an edge is present in that region In this case7 the two directions with the smallest change are used for the interpolation This allows the 10 algorithm to detect edges in the horizontal vertical and diagonal directions In 12 Chang et al introduced a method using directional information and added the ability to use multiple directions This method uses eight possible horizontal vertical and diagonal interpolation directions A gradient is computed for each direction and then a threshold is computed based on these gradients to determine which directions are used The threshold is de ned as T klm k2M m where m is the minimum gradient in the set M is the maximum gradient of the the set and k1 and k2 are constants The authors suggest k1 15 and k2 05 For each direction included in the interpolation an average red green and blue value is computed For each of the missing colors at the current pixel the difference between the average of the missing color and the average of the color of the current pixel is calculated This color difference is added to the value of the current pixel to estimate the missing color value 22 Group II Reconstruction Approaches The second group of algorithms makes some assumptions about the inter channel correla tion or the prior image and solves a mathematical problem based on those assumptions In 13 Glotzbach et al assumes that the high frequency components of the red green and blue channels are identical and replaces the aliased components of the red and blue color planes with the high frequency components of the green plane which are less likely to be aliased In 14 Gunturk et al adds a data consistency constraint to the similar high frequency components assumption of 13 and develops an iterative reconstruction al gorithm In 15 Mukherjee et al proposes a Bayesian estimation algorithm that includes a spatial smoothness assumption as a regularization term 221 Alias Canceling Interpolation There are two observations that are important for the demosaicking problem The rst is that for natural images there is a high correlation among the red green and blue channels All three channels are very likely to have the same texture and edge locations Because of 11 the similar edge content7 we expect this inter channel correlation to be even higher when it is measured between the high frequency components of the channels The second observation is that digital cameras use a color lter array CFA in which the luminance green channel is sampled at a higher rate than the chrominance red and blue channels Therefore7 the green channel is less likely to be aliased7 and details are preserved better in the green channel than in the red and blue channels The high frequency components of the red and blue channels are affected the most in CFA sampling ln demosaicking7 it is the interpolation of the red and blue channels that is the limiting factor in performance Color artifacts7 which become severe in high frequency regions such as edges7 are caused primarily by aliasing in the red and blue channels Although this fact is acknowledged by the authors of most demosaicking algorithms7 inter channel correlation has not been used effectively to retrieve the aliased high frequency information in the red and blue channels In 137 the green image is used to add high frequency information and reduce aliasing in the red and blue images First7 the red and blue images are interpolated with a rectangular lowpass lter according to the rectangular sampling grid This lls in the missing values in the grid7 but allows aliasing distortions into the red and blue output images These output images are also missing the high frequency components needed to produce a sharp image However7 because the green image is sampled at a higher rate7 the high frequency information can be taken from the green image to improve an initial interpolation of the red and blue images A horizontal highpass lter and a vertical highpass lter are applied to the green image This provides the high frequency information that the low sampling rate of the red and blue images cannot preserve Aliasing occurs when high frequency components are shifted into the low frequency portion ofthe spectrum7 so ifthe outputs of the highpass lters are modulated into the low frequency regions7 an estimate of the aliasing in the red and blue images can be found This estimate is used to reduce the aliasing in the red and blue images7 as illustrated in Figure 11 This method relies on the assumption that the high frequency information in the red7 green7 and blue images is identical If this assumption does not hold7 the addition of the green information into the red and blue images can add unwanted distortions This method also makes the assumption that the input image is band limited 1 lowpass lm the sampled red rmage 2 lsolate the hgha equency compoan In the green rmage a Add the green hglqa equaqcy eomponamm to A Modulahe the green hghafxequency eompoe the red rmage nents to estemate alrasrng In the red rmage Frgure 11 ngh d ahasrng ln the red rmage L L A A L d N am tr r n ol the green When this assumptron falls the ahasrng artrlacts are enhanced rnstead ol reduced because the green rmage also contarns ahasrng Thrs system rs composed entrrely ol lrnear lters malnng rt emclent to lrnplernent 222 Projections onto Convex Sets POCS Approach The algorrthm presented ln the prevrous sectron proposes to decompose the green channel rnto rte lrequency components and then add the hlghrflequency components ol the green channel to the lowrpass ltered red and blue channels The ls based on the obsexwtlon that 13 the high frequency components of the red blue and green channels are similar and the fact that the green channel is less likely to be aliased There are a couple of problems with this approach First it is not possible to construct ideal low pass and high pass lters with nite spatial extent Secondly the high frequency components of the red green and blue channels may not be identical Therefore replacement of the high frequency components of the red and blue channels with those of the green channel may not work well In 14 Gunturk et al proposes an algorithm that ensures data consistency at the cost of higher computational complexity The algorithm de nes two constraint sets one irnposing sirnilar high frequency components in the color channels and the other ensuring data con sistency and reconstructs the color channels using a projections onto convex sets POCS technique A Constraint Sets Set theoretic reconstruction techniques produce solutions that are consistent with the infor mation arising from the observed data or prior knowledge about the solution Each piece of information is associated with a constraint set in the solution space and the intersection of these sets represents the space of acceptable solutions 16 In 14 two types of constraint sets are de ned for the dernosaicking problem The rst comes from the observed color sam ples The interpolated color channels must be consistent with the color sarnples captured by the digital camera We denote 0711 712 as this observed data which has red green and blue sarnples placed according to the CFA used 711712 are ordered pairs of integers denoting the pixel locations By de ning A3 A0 and AB as the set of pixel locations 711712 that have the samples of red green and blue channels respectively the observation77 constraint set 00 is written as follows 00 57117712 3 50117712 00117712 V 7117712 6 A57 S RvaB7 4 where S is a generic symbol for the interpolated color channels which can be R for the red channel G for the green channel and B for the blue channel The second constraint set is a result of the high inter channel correlation In 14 it was 14 shown that color channels have very similar detail high frequency subbands This infor mation would not be enough to de ne constraint sets if all channels lost the same amount of information in sampling However7 it was also pointed out that the red and blue channels lose more information details than the green channel when captured with a color lter array Therefore7 it is reasonable to de ne constraint sets on the red and blue channels that force their high frequency components to be similar to the high frequency components of the green channel In 147 Gunturk et al proposed to decompose the signals into their subbands using a lterbank and impose the high frequency similarity in the detail subbands Let h0 and h1 be the low pass and high pass analysis lters7 which form a perfect reconstruction lterbank with synthesis lters go and gl Then the four subbands of a two dimensional signal 011712 are W15 711712 h0n1 110712 Sn1 712 5 W25 711712 h1n1 110712 Sn1 712 6 W35 711712 h0n1 h1n2 Sn1 712 7 W45 711712 h1n1 111712 Sn1 712 8 where W15 is the approximation subband7 and W287 W35 W45 are the horizontal7 vertical7 and diagonal detail subbands7 respectively We can now de ne the detail constraint set Cd that forces the details high frequency components of the red and blue channels to be similar to the details of the green channel as follows 011712 711712 WkG 711712 3 Tn1n2 Cd 7 9 V 711712 6 A5 for k 234 and S RB where Tn1n2 is a positive threshold that quanti es the closeness77 of the detail subbands to each other If the color channels are highly correlated7 then this threshold should be small if the correlation is not high7 then the threshold should be larger Although Tn1n2 15 Observed Initial Projection onto Projection onto data gt interpolation T detail constraint gt observation set constraint set I I I I I Iterate Figure 12 Alternating projections algorithm is a function of image coordinates in general7 it is also possible to use a predetermined xed value for it One choice is to set Tn1n2 to zero for all 711712 which is a result of the high correlation assumption In 147 it was also discussed how to choose a non uniform threshold B Altemating Projections Algorithm The block diagram of the P008 algorithm is given in Figure 12 The rst constraint set that is used in the reconstruction is the observation77 constraint set given in 47 which states that the reconstructed signal must be consistent with the observed data Therefore7 projection onto the observation constraint set is simply insertion of the observed data into their corresponding locations in the current image This is illustrated in Figure 13 The other constraint set is the detail77 constraint set given in 97 which states that the high frequency portions of the redblue channels are similar to those of green channel The projection operation onto the detail constraint set is illustrated in Figure 14 The observation77 projection ensures that the interpolated channels are consistent with the observed data the detail77 projection reconstructs the high frequency information of the red and blue channels7 and imposes edge consistency between the channels By alternately applying these two projections onto the initial red and blue channel estimates7 these channels are enhanced Final I estimate Va 77 A V quot Il Current Updated estimate estimate Mask for the observed data Observed data Figure 13 Projection onto observation constraint set Current Updated est1mate estimate If XY gt T Z YT Else if XY lt T Z Y T Else Z X Green channel Figure 14 Projection onto detail constraint set 223 Bayesian Approach With the Bayesian estimation approach it is possible to incorporate prior knowledge about the solution such as spatial smoothness and constant color ratio and the noise statistics into the solution We start by de ning a GFA mask function fsn1n2 17 7117712 6 As fsn17 2 0 otherwise 10 where S is R G B for the red green and blue channels Denoting 11n1 ng as the additive sensor noise the observation process can be formulated as 07117712 Z fs711771257117712 Un17n2 11 SRGB In the maximum a posteriori probability MAP formulation the observed data Onln2 the full color channels Sn1n2 and the additive noise vn1n2 are all assumed to be random processes Denoting p SlO as the conditional probability density function PDF the MAP estimate 5 is given by S arggnax p 510 arggnaxp015p57 12 where Bayes7 rule and the fact that O is independent of S were used To nd the MAP estimate 5 the conditional PDF pOlS and the prior PDF 105 need to be modeled The conditional PDF pOlS is derived from the noise statistics which is usually assumed to be white Gaussian As for the prior PDF different models have been proposed In 15 and 17 Markov Random Field MRF models were used In MRF processing the conditional and prior PDFs can be modeled as Gibbs distributions The Gibbs distribution has an exponential form and it is characterized by an energy function and a temperature parameter A PDF with Gibbs distribution can be written as W 15W lt13 Z 7 where is the energy function T is the temperature parameter and Z is the normalization constant One feature of the MRF is that the total energy function U can be written as a 18 sum of local energy functions7 which allows for localized reconstruction 18 U Z ZakUkn1n2 14 711712 k where Ukn1n2 is kth type of energy function at location 711712 and 0 is the correspond ing weight of Ukn1n2 among the other energy function types In 157 three types of local energy functions are de ned at each pixel location The rst energy function is associated with the additive noise 2 011712 7 SOn1n2gt Un7n 17fn7n 7 15 1lt 1 2 lt slt 1 2 2020117712 lt gt where SOn1n2 is the initial estimate and 020117712 is the noise variance The second energy function imposes spatial smoothness7 and is de ned as follows 17 50117712 2 Z Z 1 f5lt7117712 50117712 5711 77117712 m12 m1 71 ms7p SRGB 16 The third energy function imposes constancy of cross color ratios F F Snin2 Sn1m1n2m1 Uamn2 1 7 fsn1n2 mg may 213 0mm 2 fTmgt 2Tm WW2 m1 T GB Once the local energy functions are de ned7 the solution minimizing the total energy can be found using a variety of techniques In 157 simulated annealing technique was used As an alternative7 17 proposed a prior based on the steerable wavelet decomposition With the steerable wavelet decomposition images can be represented as a sum of bandpass com ponents7 each of which can be decomposed into a set of oriented bands using steerable lters The idea is that instead of imposing an isotropic smoothness7 such a directional decompo sition allows edge oriented smoothness This means that across the edge averaging can be avoided Let Bkn1n2 be the kth scale in the Laplacian pyramid decomposition of the image 011712 The image 011712 can be constructed by summing of Bkn1n2 over all scales 011712 gBk0117712 19 Then the energy function Ukn1n2 at a particular scale k is 1011712 Z Z 9Ek9nmltn1n2gt Bkn17n2 7 Bkm17m227 18 1M mbmzeNWlM where is a monotonically increasing function Ek9nm is the directional energy at the kth scale and in the direction of 77117712 from 711712 and Nn1n2 is a neighborhood of 711712 The optimization algorithm tries to nd 011712 that minimizes the energy function Ukn1n2 at all scales The directional energy function EN is de ned so that it returns high values along directions with similar intensities 17 Therefore EM is small in directions where an edge is crossed Since the energy function Ukn1n2 neglects the difference between neighboring pixels the strength of the penalty is less in the cross edge directions A gradient descent procedure can be used to reconstruct 011712 17 See Figure 15 23 Group III Image Formation Modeling The last group of methods models the image formation process and formulates the de mosaicking problem based on this model To understand this approach we rst need to understand the imaging process This is usually modeled as a linear process between the light radiance arriving at the camera and the pixel intensities produced by the sensors Most digital cameras use charge coupled device COD sensors In a CCD camera there is a rect angular grid of electron collection sites laid over a silicon wafer to record the amount of light energy reaching each of them When photons strike these sensor sites electron hole pairs are generated and the electrons generated at each site are collected over a certain period of time The numbers of electrons are eventually converted to pixel values Each sensor type has a speci c spectral response L5 which is a function of the spectral wavelength A and a spatial response h5x y which results from optical blur and the spatial integration at each sensor site Typical spectral sensor sensitivities are illustrated in Figure 16 The imaging process can be formulated as Observed data l Chrominance Luminance Final estimate Figure 15 Block diagram of 17 21 Sxy LSh5x 7 uy 7 vruv Adu d1 dA N5zy 19 where Sxy is the pixel value at spatial location my ry A is the incident radiance and Nszy is the additive noise that is a result of thermalquantum effects and quanti zation There are couple of assumptions in this formulation the input output relation is assumed to be linear ii the spatial blur h5zy is assumed to be space invariant and independent of wavelength iii only the additive noise is considered These assumptions are reasonable for practical purposes Since we are dealing with digital data we need to have the discrete version of 19 011712 2 Z L5lh5n1 7 m1n2 7 m2rm1m2l N5n1n2 20 l mumz The color lters sample the signal Sn1 712 to produce a mosaic plane 0n1 712 which has red green and blue pixels Therefore the observation model is a linear system which can be written in the compact form OHrN 21 where I39 O and N are the stacked forms of rm1m2l Onln2 and N5n1n2 respec tively and H is the matrix that includes the combined effects of optical blur sensor blur spectral response and CFA sampling In 19 20 and 21 the minimum mean square error MMSE solution of 21 is given rMMSE E roT E ooT 1 o 22 where E is the expectation operation In 21 the point spread function is taken as an impulse function and I39 is represented as a weighted sum of spectral basis functions to reduce the dimensionality of the problem Later 22 extended 21 to include the PSP in the reconstruction In 20 adaptive reconstruction and ways to reduce computational complexity are discussed An interesting approach to solve 22 is provided in 19 The approach constructs a nite support lter that is optimal under a set of assumptions We now switch to frequency domain formulation to explain this approach as was done in 19 We rst de ne I39n as the stacked form of rn1 712 l where n is a short notation for 711712 With the application of principal component analysis on surface re ectance it has been found that the dimensionality of rn is typically less than 7 De ning f39w 0w and as the Discrete Space Fourier Transform DSFT of rn On and Nn respectively the frequency domain equivalent of 21 can be expressed as 0w Twfw Nw 23 where is the DSFT of a linear shift invariant matrix impulse response that includes the effects of the optical blur sensor spectral response and CFA sampling With this frequency domain observation formulation the corresponding MMSE solution in 22 can be written as fMMSEm Fr0wFOw 1Ow 24 where Pro and F0 are the power spectral densities of 10 and 0 respectively These two power spectral densities can be written in terms of the power spectral densities of the input signal and the additive noise 110W FrwTw7 25 F0w TwFrwTw FNW 26 where Fr is the power spectral density PSD of I39 and I N is the PSD of N The additive noise is often assumed as white noise process therefore the PSD I N has a diagonal form I N 021 On the other hand Fr is not obvious Taubman showed that one prior that is independent of scale at the which the image is formed has the following form 1 Frw7 7w9 7fweFo 27 w f where Frww9 is the prior PSD expressed in polar coordinates with w and we being the radial and angular frequencies The angular distribution fw9 is usually taken as constant and F0 is another constant arising from surface spectral re ectances In practice a nite 23 support lter whose DSFT is Fr0w1quot0w 1 in 24 is estimated To nd that lter Fro w and F0w 1 must be evaluated numerically Classical Wiener ltering results give the optimal solution F ZU l Z is a matrix whose row vectors are given by Zlklt in 28 U is a Toeplitz matrix with elements Uk de ned in 29 Zik Ern0nik P0weiw kciu k 6 7m 2s Uik E On0n i kl P0w67 tkdw k E 7Klz 29 Both P0w and Pro are unde ned at w 0 because Frw using the scale invariant prior is unde ned at w 0 Thus these integrals cannot be computed Taubman found that numerical integration of these integrals using a dense grid of frequencies excluding the point in 0 led to large numerical errors He suggests a solution to this problem numerically integrating the integrals in 28 and 29 on a coarse grid of frequencies to obtain estimates of the integrals These rough estimates are computed for each k 6 7 HP and assembled into matrices Z and U The solution is found by solving F Z SILPOF U sire 1 30 where H is a stacked matrix with 2 12 copies of H0 L is the dimension of F0 and s is a parameter chosen as large as possible without introducing numerical instability describes the DC response of the system similar to Tw excluding the sampling aliases This solution comes from solving the limit lin ZSUS 1 31 where the 8 modi ed matrices address the singularity at w 0 The details can be found in 19 Blue Green Red 400 500 600 700 Wavelength nm Figure 16 Typical spectral sensitivity of the sensors are shown The sensitivities may go beyond 700nm therefore7 digital cameras are equipped with infrared rejection lters The matrices Pro and F0 are different for each type of pixel R7 B7 GT and G5 Therefore7 a different lter kernel will be computed for each type of pixel Taubman combines this into a single problem by considering the overlying grid of superpixels7 GE BG7 and treating these as the fundamental sampling unit The lter kernel computed from this system generates an estimate of the superpixel instead of just one element of the group 3 Comparison In this section7 the algorithms are compared with objective measures mean square error and subjective image quality lmage results from each of the algorithms are provided For these experiments7 simulated sampling was used where full color images were sampled to match the CFA sampling process A better comparison would be to use real data from a camera system7 but we don7t have the capability of performing this comparison at this time Twenty four digital color images were used in the objective measure experiments These images are part of the Kodak color image database and include various scenes The images were sampled according to the Bayer CFA and reconstructed with each of the algorithms The mean square error was measured for each color plane on each of the output images to determine the difference between the original image and the reconstructed image The bar graph in Figure 17 shows the average MSE over the set of images7 along with error bars 25 showing the 25 75 range for the set of images The graph shows that the P008 method performs best on average in terms of MSE and the small range shown in the graph shows that it is also robust and performs well for all of the images Average MSE over 24 images a b C d e f g h 1 Figure 17 Average mean square error for different algorithms a Edge directed interpola tion in 1 b Constant hue based interpolation in 5 c Weighted sum in 7 d Second order gradients as correction terms in 9 e Bayesian approach in 15 f Homogeneity directed in 10 g Pattern matching Chang in 12 h Alias cancellation in 13 P008 in 14 The numbers can only provide part of the overall story An important evaluation is the subjective appearance of the output images For this7 several example images are presented The Sail image in Figure 18 shows how the algorithms respond to sharp transitions An example of Zipper effect artifacts can be seen in d and e Most of the algorithms do a decent job reconstructing this example Figure 19 shows the Lighthouse image This example includes a picket fence from a perspective that increases spatial frequency along the fence Aliasing is a prominent artifact in this image The homogeneity directed interpolation algorithm reconstructs this image best Very little aliasing is present in the output image The Boat image in Figure 20 contains lines at various angles across the image This is a good example to show how the algorithms respond to features at various orientations The P008 algorithm and the homogeneity directed interpolation algorithm show very few of the aliasing artifacts present in the other output images This shows that these algorithms are fairly robust to the orientation of various features According to the MSE measurements7 P008 is the best algorithm7 but the output images from the homogeneity directed method 26 have fewer artifacts This suggests the need to use subjective evaluations along with objective measures In 227 Longere et al provide a perceptual assessment of demosaicking algorithms They compared several algorithms with a subjective experiment The results of their rst experi ment showed that the subjects favored sharpness and the algorithms providing a sharp image were highly favored The experiment was repeated with the result images normalized for sharpness After this adjustment7 the results showed more variation and no one algorithm was highly favored 4 Conclusions and Future Directions The sensor size of digital cameras continues to decrease7 providing sensor arrays with larger numbers of pixels Today7 ve and six mega pixel cameras are common The increased sampling rate of these cameras reduces the probability of aliasing and other artifacts Also7 Foveon has invented an imaging sensor7 the X3 sensor7 that is able to capture red7 green7 and blue information at every pixel 7 eliminating the need for demosaicking in the digital camera pipeline However7 research into the demosaicking problem is still an important problem This research has provided an understanding of the image modelling process The estimation methods discussed in this article describe the image formation process7 describing how the natural scene is transformed into a digital image The correlation between the three color planes has also been explored This extends beyond three color planes into hyperspectral image processing Processing time is often an important measure for algorithms implemented in real time systems A photographer needs to be able to take pictures at a fast rate and the image processing can sometimes limit this Several cameras7 especially the more expensive digital SLR cameras7 provide access to the raw image data captured by the sensor With this data7 the images can be processed at a later time on a computer In this case7 processing time is 27 not critically important Therefore algorithms that perform well but are computationally complex can still be considered in off line processing applications References 1 R H Hibbard Apparatus and method for adaptively interpolating a full color image utilizing luminance gradients US Patent 5382976 January 1995 B C A Laroche and M A Prescott Apparatus and method for adaptively interpolating a full color image utilizing chrominance gradients US Patent 5373322 December 1994 E D R Cok Signal processing method and apparatus for producing interpolated chromi nance values in a sampled color image signal US Patent 4642678 February 1986 E J A Weldy Optimized design for a single sensor color electronic camera system SPIE vol 1071 pp 300307 1988 E J E Adams Interactions between color plane interpolation and other image processing functions in electronic photography Proc SPIE Cameras and Systems for Electronic Photography and Scienti c Imaging vol 2416 pp 1447151 February 1995 E W T Freeman Method and apparatus for reconstructing missing color samples US Patent 4774565 1988 7 R Kimmel Demosaicing image reconstruction from ccd samples Proc Trans Image Processing vol 8 pp 122171228 1999 8 W Lu and Y P Tan On new method and performance measures for color lter array to appear in IEEE Trans Image Processing 2003 E J E Adams and J F Hamilton Design of practical color lter array interpolation algorithms for digital cameras Proc SPIE Real Time Imaging II vol 3028 pp 1177 125 February 1997 llOl llll l H 3 Bl ll4l l H E ll6l ll7l ll8l ll9l 20l K Hirakawa and T W Parks7 Adaptive homogeneity directed demosaicing algorithm77 to appear in Proc IEEE Int Conf Image Processing7 2003 X Wu7 W K Choi7 and P Bao7 Color restoration from digital camera data by pattern matching77 SPIE7 vol 30187 pp 127177 1997 E Chang7 S Cheung7 and D Y Pan7 Color lter array recovery using a threshold based variable number of gradients77 SPIE7 vol 36507 pp 367437 1999 J W Glotzbach7 R W Schafer7 and K lllgner7 A method of color lter array interpo lation with alias cancellation properties77 in Proc IEEE Int Conf Image Processing7 2001 vol 1 pp 1417144 B K Gunturk7 Y Altunbasak7 and R M Mersereau7 Color plane interpolation using alternating projections77 IEEE Trans Image Processing7 vol 117 no 97 pp 997710137 September 2002 J Mukherjee7 R Parthasarathi7 and S Goyal7 Markov random eld processing for color demosaicing77 Pattern Recognition Letters7 vol 227 pp 33973517 2001 P L Combettes7 The foundations of set theoretic estimation7 Proc of the IEEE7 vol 817 no 27 pp 18272087 February 1993 Y Hel Or and D Keren7 lmage demosaicing method utilizing directional smoothing7 US Patent 6404918 July 2002 S Geman and D Geman7 Stochastic relaxation7 gibbs distributions and the bayesian distribution of images77 IEEE Trans Pattern Analysis and Machine Intelligence7 no 67 pp 72177417 1984 D Taubman7 Generalized Wiener reconstruction of images from colour sensor data using a scale invariant prior7 in Proc IEEE Int Conf Image Processing7 20007 vol 37 pp 8017804 H J Trussell and R E Hartwig7 Mathematics for demosaicking7 IEEE Trans Image Processing7 vol 37 no 117 pp 48574927 April 2002 29 21 D H Brainardl7 Bayesian method for reconstructing color images from trichromatic samples7 Prue IS 55 T47th Annual Meeting7 pp 37573807 1994 22 P Longere7 X Zhang7 P B Delahunt7 and D H Brainardl7 Perceptual assessment of demosaicing algorithm performance7 Proceedings of the IEEE7 vol 907 no 17 pp 12Z rli y27 January 2002 Figure 12 39 39 Weighted summl I 39 91g 15 h Homogeneityrduected m 10 1 Pattern matching Chang m 12 J Ahas cancellation m 13 k FOCS m 14 Figure 19 Result Images for example lighthouse Image a Ongmal Image 1 Bxlmear Interpolation c Edgerduected Interpolation m 1 d Constant huerbased Interpolation m 5 e Weighted sum m 7 Secondroxdex gradxems as eoneeuon terms In 9 g Bayesian approach In 15 h Homogeneityrduected m 10 1 Pattern matching Chang m 12 J Ahas caneeuauon m 13 k FOCS m 14