Zhang-Wu Directional LMMSE Image Demosaicking
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
imdiff.c
Go to the documentation of this file.
1 
33 #include <math.h>
34 #include <string.h>
35 #include <ctype.h>
36 
37 #include "imageio.h"
38 #include "conv.h"
39 
40 
42 #define DISPLAY_SCALING 255
43 
44 #define MSSIM_K1 0.01
45 #define MSSIM_K2 0.03
46 
47 #define MSSIM_C1 (MSSIM_K1*MSSIM_K1)
48 #define MSSIM_C2 (MSSIM_K2*MSSIM_K2)
49 
50 
54 
56 typedef struct
57 {
59  char *FileA;
61  char *FileB;
63  int JpegQuality;
69  int Pad;
70 
74  float D;
76 
77 
78 static int ParseParams(programparams *Param, int argc, char *argv[]);
79 void MakeDifferenceImage(float *A, const float *B,
80  int Width, int Height, int NumChannels, float D);
81 void BasicMetrics(float *Max, float *Mse, const float *A, const float *B,
82  int Width, int Height, int NumChannels, int Pad);
83 float ComputeMssim(const float *A, const float *B,
84  int Width, int Height, int NumChannels, int Pad);
85 
86 
89 {
90  puts("Image difference calculator, P. Getreuer 2010-2011\n");
91  puts("Usage: imdiff [options] <exact file> <distorted file>\n"
92  "Only " READIMAGE_FORMATS_SUPPORTED " images are supported.\n");
93  puts("Options:");
94  puts(" -m <metric> Metric to use for comparison, choices are");
95  puts(" max Maximum absolute difference, max_n |A_n - B_n|");
96  puts(" mse Mean squared error, 1/N sum |A_n - B_n|^2");
97  puts(" rmse Root mean squared error, (MSE)^1/2");
98  puts(" psnr Peak signal-to-noise ratio, -10 log10(MSE/255^2)");
99  puts(" mssim Mean structural similarity index\n");
100  puts(" -s Compute metric separately for each channel");
101  puts(" -p <pad> Remove a margin of <pad> pixels before comparison");
102  puts(" -D <number> D parameter for difference image\n");
103 #ifdef USE_LIBJPEG
104  puts(" -q <number> Quality for saving JPEG images (0 to 100)\n");
105 #endif
106  puts("Alternatively, a difference image is generated by the syntax\n"
107  " imdiff [-D <number>] <exact file> <distorted file> <output file>\n");
108  puts("The difference image is computed as\n"
109  " D_n = 255/D (A_n - B_n) + 255/2.\n"
110  "Values outside of the range [0,255] are saturated.\n");
111  puts("Example:\n"
112 #ifdef USE_LIBPNG
113  " imdiff -mpsnr frog-exact.png frog-4x.png");
114 #else
115  " imdiff -mpsnr frog-exact.bmp frog-4x.bmp");
116 #endif
117 }
118 
119 
120 int main(int argc, char *argv[])
121 {
122  struct
123  {
124  float *Data;
125  int Width;
126  int Height;
127  } A = {NULL, 0, 0}, B = {NULL, 0, 0};
128 
129  programparams Param;
130  float Max, MaxC[3], Mse, MseC[3], Mssim;
131  int Channel, Status = 1;
132 
133 
134  if(!ParseParams(&Param, argc, argv))
135  return 0;
136 
137  /* Read the exact image */
138  if(!(A.Data = (float *)ReadImage(&A.Width, &A.Height, Param.FileA,
140  goto Catch;
141 
142  /* Read the distorted image */
143  if(!(B.Data = (float *)ReadImage(&B.Width, &B.Height, Param.FileB,
145  goto Catch;
146 
147  if(A.Width != B.Width || A.Height != B.Height)
148  {
149  ErrorMessage("Image sizes don't match, %dx%d vs. %dx%d.\n",
150  A.Width, A.Height, B.Width, B.Height);
151  goto Catch;
152  }
153  else if(A.Width <= 2*Param.Pad || A.Height <= 2*Param.Pad)
154  {
155  ErrorMessage(
156  "Removal of %d-pixel padding removes entire %dx%d image.\n",
157  Param.Pad, A.Width, A.Height);
158  goto Catch;
159  }
160 
161  if(Param.DifferenceFile)
162  {
163  MakeDifferenceImage(A.Data, B.Data, A.Width, A.Height, 3, Param.D);
164 
165  if(!(WriteImage(A.Data, A.Width, A.Height, Param.DifferenceFile,
167  goto Catch;
168  }
169  else
170  {
171  Max = 0;
172  Mse = 0;
173 
174  for(Channel = 0; Channel < 3; Channel++)
175  {
176  BasicMetrics(&MaxC[Channel], &MseC[Channel],
177  A.Data + Channel*A.Width*A.Height,
178  B.Data + Channel*B.Width*B.Height,
179  A.Width, A.Height, 1, Param.Pad);
180 
181  if(MaxC[Channel] > Max)
182  Max = MaxC[Channel];
183 
184  Mse += MseC[Channel];
185  }
186 
187  Mse /= 3;
188 
189  switch(Param.Metric)
190  {
191  case DEFAULT_METRICS:
192  if(!Param.SeparateChannels)
193  {
194  printf("Maximum absolute difference: %g\n", DISPLAY_SCALING*Max);
195  printf("Peak signal-to-noise ratio: %.4f\n", -10*log10(Mse));
196  }
197  else
198  {
199  printf("Maximum absolute difference: %g %g %g\n",
200  DISPLAY_SCALING*MaxC[0], DISPLAY_SCALING*MaxC[1],
201  DISPLAY_SCALING*MaxC[2]);
202  printf("Peak signal-to-noise ratio: %.4f %.4f %.4f\n",
203  -10*log10(MseC[0]), -10*log10(MseC[1]), -10*log10(MseC[2]));
204  }
205 
206  if(A.Width <= 2*(5 + Param.Pad)
207  || A.Height <= 2*(5 + Param.Pad))
208  printf("Image size is too small to compute MSSIM.\n");
209  else
210  {
211 
212  Mssim = (Max == 0) ? 1 : ComputeMssim(A.Data, B.Data,
213  A.Width, A.Height, 3, Param.Pad);
214 
215  if(Mssim != -1)
216  printf("Mean structural similarity: %.4f\n", Mssim);
217  }
218  break;
219  case MAX_METRIC:
220  if(!Param.SeparateChannels)
221  printf("%g\n", DISPLAY_SCALING*Max);
222  else
223  printf("%g %g %g\n", DISPLAY_SCALING*MaxC[0],
224  DISPLAY_SCALING*MaxC[1], DISPLAY_SCALING*MaxC[2]);
225  break;
226  case MSE_METRIC:
227  if(!Param.SeparateChannels)
228  printf("%.4f\n", DISPLAY_SCALING*DISPLAY_SCALING*Mse);
229  else
230  printf("%.4f %.4f %.4f\n",
234  break;
235  case RMSE_METRIC:
236  if(!Param.SeparateChannels)
237  printf("%.4f\n", DISPLAY_SCALING*sqrt(Mse));
238  else
239  printf("%.4f %.4f %.4f\n",
240  DISPLAY_SCALING*sqrt(MseC[0]),
241  DISPLAY_SCALING*sqrt(MseC[1]),
242  DISPLAY_SCALING*sqrt(MseC[2]));
243  break;
244  case PSNR_METRIC:
245  if(!Param.SeparateChannels)
246  printf("%.4f\n", -10*log10(Mse));
247  else
248  printf("%.4f %.4f %.4f\n",
249  -10*log10(MseC[0]), -10*log10(MseC[1]), -10*log10(MseC[2]));
250  break;
251  case MSSIM_METRIC:
252  if(A.Width <= 2*(5 + Param.Pad)
253  || A.Height <= 2*(5 + Param.Pad))
254  printf("Image size is too small to compute MSSIM.\n");
255  else
256  {
257  Mssim = (Max == 0) ? 1 : ComputeMssim(A.Data, B.Data,
258  A.Width, A.Height, 3, Param.Pad);
259 
260  if(Mssim == -1)
261  ErrorMessage("Memory allocation failed.\n");
262  else
263  printf("%.4f\n", Mssim);
264  }
265  break;
266  }
267  }
268 
269  Status = 0; /* Finished successfully */
270 Catch:
271  Free(B.Data);
272  Free(A.Data);
273  return Status;
274 }
275 
276 
278 void MakeDifferenceImage(float *A, const float *B,
279  int Width, int Height, int NumChannels, float D)
280 {
281  const int NumEl = NumChannels*Width*Height;
282  int n;
283 
284  D /= 255;
285 
286  for(n = 0; n < NumEl; n++)
287  A[n] = (A[n] - B[n])/D + (float)0.5;
288 }
289 
290 
292 void BasicMetrics(float *Max, float *Mse, const float *A, const float *B,
293  int Width, int Height, int NumChannels, int Pad)
294 {
295  float Diff, CurMax = 0;
296  double AccumMse = 0;
297  int x, y, Channel, n;
298 
299 
300  for(Channel = 0; Channel < NumChannels; Channel++)
301  for(y = Pad; y < Height - Pad; y++)
302  for(x = Pad; x < Width - Pad; x++)
303  {
304  n = x + Width*(y + Height*Channel);
305  Diff = (float)fabs(A[n] - B[n]);
306 
307  if(CurMax < Diff)
308  CurMax = Diff;
309 
310  AccumMse += Diff*Diff;
311  }
312 
313  *Max = CurMax;
314  *Mse = (float)(AccumMse / (NumChannels*(Width - 2*Pad)*(Height - 2*Pad)));
315 }
316 
317 
319 float ComputeMssim(const float *A, const float *B,
320  int Width, int Height, int NumChannels, int Pad)
321 {
322  /* 11-tap Gaussian filter with standard deviation 1.5 */
323  const int R = 5;
324  filter Window = GaussianFilter(1.5, R);
325  /* Boundary does not matter, convolution is used only in the interior */
326  boundaryext Boundary = GetBoundaryExt("zpd");
327 
328  const int NumPixels = Width*Height;
329  const int NumEl = NumChannels*NumPixels;
330  float *Buffer = NULL, *MuA = NULL, *MuB = NULL,
331  *MuAA = NULL, *MuBB = NULL, *MuAB = NULL;
332  double MuASqr, MuBSqr, MuAMuB,
333  SigmaASqr, SigmaBSqr, SigmaAB, Mssim = -1;
334  int n, x, y, c;
335 
336 
337  if(IsNullFilter(Window)
338  || !(Buffer = (float *)Malloc(sizeof(float)*NumPixels))
339  || !(MuA = (float *)Malloc(sizeof(float)*NumEl))
340  || !(MuB = (float *)Malloc(sizeof(float)*NumEl))
341  || !(MuAA = (float *)Malloc(sizeof(float)*NumEl))
342  || !(MuBB = (float *)Malloc(sizeof(float)*NumEl))
343  || !(MuAB = (float *)Malloc(sizeof(float)*NumEl)))
344  goto Catch;
345 
346  SeparableConv2D(MuA, Buffer, A, Window, Window,
347  Boundary, Width, Height, NumChannels);
348  SeparableConv2D(MuB, Buffer, B, Window, Window,
349  Boundary, Width, Height, NumChannels);
350 
351  for(n = 0; n < NumEl; n++)
352  {
353  MuAA[n] = A[n]*A[n];
354  MuBB[n] = B[n]*B[n];
355  MuAB[n] = A[n]*B[n];
356  }
357 
358  SeparableConv2D(MuAA, Buffer, MuAA, Window, Window,
359  Boundary, Width, Height, NumChannels);
360  SeparableConv2D(MuBB, Buffer, MuBB, Window, Window,
361  Boundary, Width, Height, NumChannels);
362  SeparableConv2D(MuAB, Buffer, MuAB, Window, Window,
363  Boundary, Width, Height, NumChannels);
364  Mssim = 0;
365 
366  Pad += R;
367 
368  for(c = 0; c < NumChannels; c++)
369  for(y = Pad; y < Height - Pad; y++)
370  for(x = Pad; x < Width - Pad; x++)
371  {
372  n = x + Width*(y + Height*c);
373  MuASqr = MuA[n]*MuA[n];
374  MuBSqr = MuB[n]*MuB[n];
375  MuAMuB = MuA[n]*MuB[n];
376  SigmaASqr = MuAA[n] - MuASqr;
377  SigmaBSqr = MuBB[n] - MuBSqr;
378  SigmaAB = MuAB[n] - MuAMuB;
379 
380  Mssim += ((2*MuAMuB + MSSIM_C1)*(2*SigmaAB + MSSIM_C2))
381  / ((MuASqr + MuBSqr + MSSIM_C1)
382  *(SigmaASqr + SigmaBSqr + MSSIM_C2));
383  }
384 
385  Mssim /= NumChannels*(Width - 2*Pad)*(Height - 2*Pad);
386 
387 Catch:
388  FreeFilter(Window);
389  Free(MuAB);
390  Free(MuBB);
391  Free(MuAA);
392  Free(MuB);
393  Free(MuA);
394  Free(Buffer);
395  return (float)Mssim;
396 }
397 
398 
399 
400 static int ParseParams(programparams *Param, int argc, char *argv[])
401 {
402  char *OptionString;
403  char OptionChar;
404  int i;
405 
406 
407  if(argc < 2)
408  {
410  return 0;
411  }
412 
413  /* Set parameter defaults */
414  Param->FileA = NULL;
415  Param->FileB = NULL;
416  Param->Metric = DEFAULT_METRICS;
417  Param->SeparateChannels = 0;
418 
419  Param->Pad = 0;
420  Param->DifferenceFile = NULL;
421  Param->JpegQuality = 95;
422  Param->D = 20;
423 
424  for(i = 1; i < argc;)
425  {
426  if(argv[i] && argv[i][0] == '-')
427  {
428  if((OptionChar = argv[i][1]) == 0)
429  {
430  ErrorMessage("Invalid parameter format.\n");
431  return 0;
432  }
433 
434  if(argv[i][2])
435  OptionString = &argv[i][2];
436  else if(++i < argc)
437  OptionString = argv[i];
438  else
439  {
440  ErrorMessage("Invalid parameter format.\n");
441  return 0;
442  }
443 
444  switch(OptionChar)
445  {
446  case 'p':
447  Param->Pad = atoi(OptionString);
448 
449  if(Param->Pad < 0)
450  {
451  ErrorMessage("Pad must be nonnegative.\n");
452  return 0;
453  }
454  break;
455  case 's':
456  Param->SeparateChannels = 1;
457  i--;
458  break;
459  case 'D':
460  Param->D = (float)atof(OptionString);
461 
462  if(Param->D <= 0)
463  {
464  ErrorMessage("D must be positive.\n");
465  return 0;
466  }
467  break;
468  case 'm':
469  if(!strcmp(OptionString, "max"))
470  Param->Metric = MAX_METRIC;
471  else if(!strcmp(OptionString, "mse"))
472  Param->Metric = MSE_METRIC;
473  else if(!strcmp(OptionString, "rmse"))
474  Param->Metric = RMSE_METRIC;
475  else if(!strcmp(OptionString, "psnr"))
476  Param->Metric = PSNR_METRIC;
477  else if(!strcmp(OptionString, "mssim"))
478  Param->Metric = MSSIM_METRIC;
479  else
480  ErrorMessage("Unknown metric.\n");
481  break;
482 
483 #ifdef USE_LIBJPEG
484  case 'q':
485  Param->JpegQuality = atoi(OptionString);
486 
487  if(Param->JpegQuality <= 0 || Param->JpegQuality > 100)
488  {
489  ErrorMessage("JPEG quality must be between 0 and 100.\n");
490  return 0;
491  }
492  break;
493 #endif
494  case '-':
496  return 0;
497  default:
498  if(isprint(OptionChar))
499  ErrorMessage("Unknown option \"-%c\".\n", OptionChar);
500  else
501  ErrorMessage("Unknown option.\n");
502 
503  return 0;
504  }
505 
506  i++;
507  }
508  else
509  {
510  if(!Param->FileA)
511  Param->FileA = argv[i];
512  else if(!Param->FileB)
513  Param->FileB = argv[i];
514  else
515  Param->DifferenceFile = argv[i];
516 
517  i++;
518  }
519  }
520 
521  if(!Param->FileA || !Param->FileB)
522  {
524  return 0;
525  }
526 
527  return 1;
528 }