Zhang-Wu Directional LMMSE Image Demosaicking
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
dmzhangwu.c
Go to the documentation of this file.
1 
26 #include "basic.h"
27 #include "conv.h"
28 #include "dmzhangwu.h"
29 
30 
31 float DiagonalAverage(const float *Image, int Width, int Height, int x, int y);
32 float AxialAverage(const float *Image, int Width, int Height, int x, int y);
33 
34 
60 int ZhangWuDemosaic(float *Output, const float *Input,
61  int Width, int Height, int RedX, int RedY, int UseZhangCodeEst)
62 {
63  /* Window size for estimating LMMSE statistics */
64  const int M = 4;
65  /* Small value added in denominators to avoid divide-by-zero */
66  const float DivEpsilon = 0.1f/(255*255);
67  /* Interpolation filter used for horizontal and vertical interpolations */
68  static float InterpCoeff[5] = {-0.25f, 0.5f, 0.5f, 0.5f, -0.25f};
69  /* Approximately Gaussian smoothing filter used in the LMMSE denoising */
70  static float SmoothCoeff[9] = {0.03125f, 0.0703125f, 0.1171875f,
71  0.1796875f, 0.203125f, 0.1796875f, 0.1171875f, 0.0703125f, 0.03125f};
72  /* Use whole-sample symmeric boundary handling in convolutions */
73  boundaryext Boundary = GetBoundaryExt("symw");
74 
75  filter InterpFilter = MakeFilter(InterpCoeff, -2, 5);
76  filter SmoothFilter = MakeFilter(SmoothCoeff, -4, 9);
77  const int NumPixels = Width*Height;
78  const int Green = 1 - ((RedX + RedY) & 1);
79  float *OutputRed = Output;
80  float *OutputGreen = Output + NumPixels;
81  float *OutputBlue = Output + 2*NumPixels;
82  float *FilteredH = NULL, *FilteredV = NULL;
83  float *DiffH = NULL, *DiffV = NULL, *DiffGR, *DiffGB;
84  float mom1, h, H, mh, ph, Rh, v, V, mv, pv, Rv, Temp;
85  int x, y, i, m, m0, m1, Success = 0;
86 
87 
88  /* Allocate memory for workspace buffers */
89  if(!(FilteredH = (float *)Malloc(sizeof(float)*NumPixels))
90  || !(FilteredV = (float *)Malloc(sizeof(float)*NumPixels))
91  || !(DiffGR = DiffH = (float *)Malloc(sizeof(float)*NumPixels))
92  || !(DiffGB = DiffV = (float *)Malloc(sizeof(float)*NumPixels)))
93  goto Catch;
94 
95  /* Horizontal and vertical 1D interpolations */
96  for(y = 0; y < Height; y++)
97  Conv1D(FilteredH + Width*y, 1,
98  Input + Width*y, 1, InterpFilter, Boundary, Width);
99 
100  for(x = 0; x < Width; x++)
101  Conv1D(FilteredV + x, Width,
102  Input + x, Width, InterpFilter, Boundary, Height);
103 
104  /* Local noise estimation for LMMSE */
105  for(y = 0, i = 0; y < Height; y++)
106  for(x = 0; x < Width; x++, i++)
107  {
108  if(((x + y) & 1) == Green)
109  {
110  DiffH[i] = Input[i] - FilteredH[i];
111  DiffV[i] = Input[i] - FilteredV[i];
112  }
113  else
114  {
115  DiffH[i] = FilteredH[i] - Input[i];
116  DiffV[i] = FilteredV[i] - Input[i];
117  }
118  }
119 
120  /* Compute the smoothed signals for LMMSE */
121  for(y = 0; y < Height; y++)
122  Conv1D(FilteredH + Width*y, 1, DiffH + Width*y, 1,
123  SmoothFilter, Boundary, Width);
124 
125  for(x = 0; x < Width; x++)
126  Conv1D(FilteredV + x, Width, DiffV + x, Width,
127  SmoothFilter, Boundary, Height);
128 
129  /* LMMSE interpolation of the green channel */
130  for(y = 0, i = 0; y < Height; y++)
131  for(x = 0; x < Width; x++, i++)
132  {
133  if(((x + y) & 1) != Green) /* (x,y) is a red or blue location */
134  {
135  /* Adjust loop indices m = -M,...,M when necessary to
136  compensate for left and right boundaries. We effectively
137  do zero-padded boundary handling. */
138  m0 = (x >= M) ? -M : -x;
139  m1 = (x < Width - M) ? M : (Width - x - 1);
140 
141  /* The following computes
142  * ph = var FilteredH[i + m]
143  * m=-M,...,M
144  * Rh = mean (FilteredH[i + m] - DiffH[i + m])^2
145  * m=-M,...,M
146  * h = LMMSE estimate
147  * H = LMMSE estimate accuracy (estimated variance of h)
148  */
149 
150  for(m = m0, mom1 = ph = Rh = 0; m <= m1; m++)
151  {
152  Temp = FilteredH[i + m];
153  mom1 += Temp;
154  ph += Temp*Temp;
155  Temp -= DiffH[i + m];
156  Rh += Temp*Temp;
157  }
158 
159  if(!UseZhangCodeEst)
160  /* Compute mh = mean_m FilteredH[i + m] */
161  mh = mom1/(2*M + 1);
162  else
163  /* Compute mh as in Zhang's MATLAB code */
164  mh = FilteredH[i];
165 
166  ph = ph/(2*M) - mom1*mom1/(2*M*(2*M + 1));
167  Rh = Rh/(2*M + 1) + DivEpsilon;
168  h = mh + (ph/(ph + Rh))*(DiffH[i] - mh);
169  H = ph - (ph/(ph + Rh))*ph + DivEpsilon;
170 
171  /* Adjust loop indices for top and bottom boundaries. */
172  m0 = (y >= M) ? -M : -y;
173  m1 = (y < Height - M) ? M : (Height - y - 1);
174 
175  /* The following computes
176  * pv = var FilteredV[i + m]
177  * m=-M,...,M
178  * Rv = mean (FilteredV[i + m] - DiffV[i + m])^2
179  * m=-M,...,M
180  * v = LMMSE estimate
181  * V = LMMSE estimate accuracy (estimated variance of v)
182  */
183 
184  for(m = m0, mom1 = pv = Rv = 0; m <= m1; m++)
185  {
186  Temp = FilteredV[i + Width*m];
187  mom1 += Temp;
188  pv += Temp*Temp;
189  Temp -= DiffV[i + Width*m];
190  Rv += Temp*Temp;
191  }
192 
193  if(!UseZhangCodeEst)
194  /* Compute mv = mean_m FilteredV[i + m] */
195  mv = mom1/(2*M + 1);
196  else
197  /* Compute mv as in Zhang's MATLAB code */
198  mv = FilteredV[i];
199 
200  pv = pv/(2*M) - mom1*mom1/(2*M*(2*M + 1));
201  Rv = Rv/(2*M + 1) + DivEpsilon;
202  v = mv + (pv/(pv + Rv))*(DiffV[i] - mv);
203  V = pv - (pv/(pv + Rv))*pv + DivEpsilon;
204 
205  /* Fuse the directional estimates to obtain
206  the green component. */
207  OutputGreen[i] = Input[i] + (V*h + H*v) / (H + V);
208  }
209  else
210  OutputGreen[i] = Input[i];
211  }
212 
213  /* Compute the primary difference signals:
214  DiffGR = Green - Red (known at red locations)
215  DiffGB = Green - Blue (known at blue locations) */
216  for(y = 0, i = 0; y < Height; y++)
217  for(x = 0; x < Width; x++, i++)
218  if(((x + y) & 1) != Green)
219  (((y & 1) == RedY) ? DiffGR : DiffGB)[i]
220  = OutputGreen[i] - Input[i];
221 
222  /* Interpolate DiffGR at blue locations and DiffGB at red locations */
223  for(y = 0, i = 0; y < Height; y++)
224  for(x = 0; x < Width; x++, i++)
225  if(((x + y) & 1) != Green)
226  {
227  if((y & 1) == RedY)
228  DiffGB[i] = DiagonalAverage(DiffGB, Width, Height, x, y);
229  else
230  DiffGR[i] = DiagonalAverage(DiffGR, Width, Height, x, y);
231  }
232 
233  /* Interpolate DiffGR and DiffGB at green locations */
234  for(y = 0, i = 0; y < Height; y++)
235  for(x = 0; x < Width; x++, i++)
236  if(((x + y) & 1) == Green)
237  {
238  DiffGB[i] = AxialAverage(DiffGB, Width, Height, x, y);
239  DiffGR[i] = AxialAverage(DiffGR, Width, Height, x, y);
240  }
241 
242  /* Obtain the red and blue channel interpolations */
243  for(y = 0, i = 0; y < Height; y++)
244  for(x = 0; x < Width; x++, i++)
245  {
246  OutputRed[i] = OutputGreen[i] - DiffGR[i];
247  OutputBlue[i] = OutputGreen[i] - DiffGB[i];
248  }
249 
250  Success = 1;
251 Catch: /* This label is used for error handling. If something went wrong
252  above (which may be out of memory or a computation error), then
253  execution jumps to this point to clean up and exit. */
254  /* Free buffers */
255  Free(DiffV);
256  Free(DiffH);
257  Free(FilteredV);
258  Free(FilteredH);
259  return Success;
260 }
261 
262 
273 float DiagonalAverage(const float *Image, int Width, int Height, int x, int y)
274 {
275  Image += x + Width*y;
276 
277  if(y == 0)
278  {
279  if(x == 0)
280  return Image[1 + Width];
281  else if(x < Width - 1)
282  return (Image[-1 + Width] + Image[1 + Width]) / 2;
283  else
284  return Image[-1 + Width];
285  }
286  else if(y < Height - 1)
287  {
288  if(x == 0)
289  return (Image[1 - Width] + Image[1 + Width]) / 2;
290  else if(x < Width - 1)
291  return (Image[-1 - Width] + Image[1 - Width]
292  + Image[-1 + Width] + Image[1 + Width]) / 4;
293  else
294  return (Image[-1 - Width] + Image[-1 + Width]) / 2;
295  }
296  else
297  {
298  if(x == 0)
299  return Image[1 - Width];
300  else if(x < Width - 1)
301  return (Image[-1 - Width] + Image[1 - Width]) / 2;
302  else
303  return Image[-1 - Width];
304  }
305 }
306 
307 
318 float AxialAverage(const float *Image, int Width, int Height, int x, int y)
319 {
320  Image += x + Width*y;
321 
322  if(y == 0)
323  {
324  if(x == 0)
325  return (Image[1] + Image[Width]) / 2;
326  else if(x < Width - 1)
327  return (Image[-1] + Image[1] + 2*Image[Width]) / 4;
328  else
329  return (Image[-1] + Image[Width]) / 2;
330 
331  }
332  else if(y < Height - 1)
333  {
334  if(x == 0)
335  return (2*Image[1] + Image[-Width] + Image[Width]) / 4;
336  else if(x < Width - 1)
337  return (Image[-1] + Image[1] + Image[-Width] + Image[Width]) / 4;
338  else
339  return (2*Image[-1] + Image[-Width] + Image[Width]) / 4;
340  }
341  else
342  {
343  if(x == 0)
344  return (Image[1] + Image[-Width])/2;
345  else if(x < Width - 1)
346  return (Image[-1] + Image[1] + 2*Image[-Width]) / 4;
347  else
348  return (Image[-1] + Image[-Width]) / 2;
349  }
350 }
Generated on Fri Jul 5 2013 18:57:37 for Zhang-Wu Directional LMMSE Image Demosaicking by   doxygen 1.8.3.1