granicus.if.org Git - libvpx/blob - tools/tiny_ssim.c

   1 /*
   2  *  Copyright (c) 2016 The WebM project authors. All Rights Reserved.
   3  *
   4  *  Use of this source code is governed by a BSD-style license
   5  *  that can be found in the LICENSE file in the root of the source
   6  *  tree. An additional intellectual property rights grant can be found
   7  *  in the file PATENTS.  All contributing project authors may
   8  *  be found in the AUTHORS file in the root of the source tree.
   9  */
  10
  11 #include <errno.h>
  12 #include <math.h>
  13 #include <stdio.h>
  14 #include <stdlib.h>
  15 #include <string.h>
  16 #include "vpx/vpx_codec.h"
  17 #include "vpx/vpx_integer.h"
  18 #include "./y4minput.h"
  19
  20 void vp8_ssim_parms_8x8_c(unsigned char *s, int sp, unsigned char *r, int rp,
  21                           uint32_t *sum_s, uint32_t *sum_r, uint32_t *sum_sq_s,
  22                           uint32_t *sum_sq_r, uint32_t *sum_sxr) {
  23   int i, j;
  24   for (i = 0; i < 8; i++, s += sp, r += rp) {
  25     for (j = 0; j < 8; j++) {
  26       *sum_s += s[j];
  27       *sum_r += r[j];
  28       *sum_sq_s += s[j] * s[j];
  29       *sum_sq_r += r[j] * r[j];
  30       *sum_sxr += s[j] * r[j];
  31     }
  32   }
  33 }
  34
  35 static const int64_t cc1 = 26634;   // (64^2*(.01*255)^2
  36 static const int64_t cc2 = 239708;  // (64^2*(.03*255)^2
  37
  38 static double similarity(uint32_t sum_s, uint32_t sum_r, uint32_t sum_sq_s,
  39                          uint32_t sum_sq_r, uint32_t sum_sxr, int count) {
  40   int64_t ssim_n, ssim_d;
  41   int64_t c1, c2;
  42
  43   // scale the constants by number of pixels
  44   c1 = (cc1 * count * count) >> 12;
  45   c2 = (cc2 * count * count) >> 12;
  46
  47   ssim_n = (2 * sum_s * sum_r + c1) *
  48            ((int64_t)2 * count * sum_sxr - (int64_t)2 * sum_s * sum_r + c2);
  49
  50   ssim_d = (sum_s * sum_s + sum_r * sum_r + c1) *
  51            ((int64_t)count * sum_sq_s - (int64_t)sum_s * sum_s +
  52             (int64_t)count * sum_sq_r - (int64_t)sum_r * sum_r + c2);
  53
  54   return ssim_n * 1.0 / ssim_d;
  55 }
  56
  57 static double ssim_8x8(unsigned char *s, int sp, unsigned char *r, int rp) {
  58   uint32_t sum_s = 0, sum_r = 0, sum_sq_s = 0, sum_sq_r = 0, sum_sxr = 0;
  59   vp8_ssim_parms_8x8_c(s, sp, r, rp, &sum_s, &sum_r, &sum_sq_s, &sum_sq_r,
  60                        &sum_sxr);
  61   return similarity(sum_s, sum_r, sum_sq_s, sum_sq_r, sum_sxr, 64);
  62 }
  63
  64 // We are using a 8x8 moving window with starting location of each 8x8 window
  65 // on the 4x4 pixel grid. Such arrangement allows the windows to overlap
  66 // block boundaries to penalize blocking artifacts.
  67 double vp8_ssim2(unsigned char *img1, unsigned char *img2, int stride_img1,
  68                  int stride_img2, int width, int height) {
  69   int i, j;
  70   int samples = 0;
  71   double ssim_total = 0;
  72
  73   // sample point start with each 4x4 location
  74   for (i = 0; i <= height - 8;
  75        i += 4, img1 += stride_img1 * 4, img2 += stride_img2 * 4) {
  76     for (j = 0; j <= width - 8; j += 4) {
  77       double v = ssim_8x8(img1 + j, stride_img1, img2 + j, stride_img2);
  78       ssim_total += v;
  79       samples++;
  80     }
  81   }
  82   ssim_total /= samples;
  83   return ssim_total;
  84 }
  85
  86 static uint64_t calc_plane_error(uint8_t *orig, int orig_stride, uint8_t *recon,
  87                                  int recon_stride, unsigned int cols,
  88                                  unsigned int rows) {
  89   unsigned int row, col;
  90   uint64_t total_sse = 0;
  91   int diff;
  92
  93   for (row = 0; row < rows; row++) {
  94     for (col = 0; col < cols; col++) {
  95       diff = orig[col] - recon[col];
  96       total_sse += diff * diff;
  97     }
  98
  99     orig += orig_stride;
 100     recon += recon_stride;
 101   }
 102   return total_sse;
 103 }
 104
 105 #define MAX_PSNR 100
 106 double vp9_mse2psnr(double samples, double peak, double mse) {
 107   double psnr;
 108
 109   if (mse > 0.0)
 110     psnr = 10.0 * log10(peak * peak * samples / mse);
 111   else
 112     psnr = MAX_PSNR;  // Limit to prevent / 0
 113
 114   if (psnr > MAX_PSNR) psnr = MAX_PSNR;
 115
 116   return psnr;
 117 }
 118
 119 typedef enum { RAW_YUV, Y4M } input_file_type;
 120
 121 typedef struct input_file {
 122   FILE *file;
 123   input_file_type type;
 124   unsigned char *buf;
 125   y4m_input y4m;
 126   vpx_image_t img;
 127   int w;
 128   int h;
 129 } input_file_t;
 130
 131 // Open a file and determine if its y4m or raw.  If y4m get the header.
 132 int open_input_file(const char *file_name, input_file_t *input, int w, int h) {
 133   char y4m_buf[4];
 134   size_t r1;
 135   input->type = RAW_YUV;
 136   input->buf = NULL;
 137   input->file = strcmp(file_name, "-") ? fopen(file_name, "rb") : stdin;
 138   if (input->file == NULL) return -1;
 139   r1 = fread(y4m_buf, 1, 4, input->file);
 140   if (r1 == 4) {
 141     if (memcmp(y4m_buf, "YUV4", 4) == 0) input->type = Y4M;
 142     switch (input->type) {
 143       case Y4M:
 144         y4m_input_open(&input->y4m, input->file, y4m_buf, 4, 0);
 145         input->w = input->y4m.pic_w;
 146         input->h = input->y4m.pic_h;
 147         // Y4M alloc's its own buf. Init this to avoid problems if we never
 148         // read frames.
 149         memset(&input->img, 0, sizeof(input->img));
 150         break;
 151       case RAW_YUV:
 152         fseek(input->file, 0, SEEK_SET);
 153         input->w = w;
 154         input->h = h;
 155         input->buf = malloc(w * h * 3 / 2);
 156         break;
 157     }
 158   }
 159   return 0;
 160 }
 161
 162 void close_input_file(input_file_t *in) {
 163   if (in->file) fclose(in->file);
 164   if (in->type == Y4M) {
 165     vpx_img_free(&in->img);
 166   } else {
 167     free(in->buf);
 168   }
 169 }
 170
 171 size_t read_input_file(input_file_t *in, unsigned char **y, unsigned char **u,
 172                        unsigned char **v) {
 173   size_t r1 = 0;
 174   switch (in->type) {
 175     case Y4M:
 176       r1 = y4m_input_fetch_frame(&in->y4m, in->file, &in->img);
 177       *y = in->img.planes[0];
 178       *u = in->img.planes[1];
 179       *v = in->img.planes[2];
 180       break;
 181     case RAW_YUV:
 182       r1 = fread(in->buf, in->w * in->h * 3 / 2, 1, in->file);
 183       *y = in->buf;
 184       *u = in->buf + in->w * in->h;
 185       *v = in->buf + 5 * in->w * in->h / 4;
 186       break;
 187   }
 188
 189   return r1;
 190 }
 191
 192 int main(int argc, char *argv[]) {
 193   FILE *framestats = NULL;
 194   int w = 0, h = 0, tl_skip = 0, tl_skips_remaining = 0;
 195   double ssimavg = 0, ssimyavg = 0, ssimuavg = 0, ssimvavg = 0;
 196   double psnrglb = 0, psnryglb = 0, psnruglb = 0, psnrvglb = 0;
 197   double psnravg = 0, psnryavg = 0, psnruavg = 0, psnrvavg = 0;
 198   double *ssimy = NULL, *ssimu = NULL, *ssimv = NULL;
 199   uint64_t *psnry = NULL, *psnru = NULL, *psnrv = NULL;
 200   size_t i, n_frames = 0, allocated_frames = 0;
 201   int return_value = 0;
 202   input_file_t in[2];
 203
 204   if (argc < 2) {
 205     fprintf(stderr,
 206             "Usage: %s file1.{yuv|y4m} file2.{yuv|y4m}"
 207             "[WxH tl_skip={0,1,3}]\n",
 208             argv[0]);
 209     return_value = 1;
 210     goto clean_up;
 211   }
 212
 213   if (argc > 3) {
 214     sscanf(argv[3], "%dx%d", &w, &h);
 215   }
 216
 217   if (open_input_file(argv[1], &in[0], w, h) < 0) {
 218     fprintf(stderr, "File %s can't be opened or parsed!\n", argv[2]);
 219     goto clean_up;
 220   }
 221
 222   if (w == 0 && h == 0) {
 223     // If a y4m is the first file and w, h is not set grab from first file.
 224     w = in[0].w;
 225     h = in[0].h;
 226   }
 227
 228   if (open_input_file(argv[2], &in[1], w, h) < 0) {
 229     fprintf(stderr, "File %s can't be opened or parsed!\n", argv[2]);
 230     goto clean_up;
 231   }
 232
 233   if (in[0].w != in[1].w || in[0].h != in[1].h || in[0].w != w ||
 234       in[0].h != h || w == 0 || h == 0) {
 235     fprintf(stderr,
 236             "Failing: Image dimensions don't match or are unspecified!\n");
 237     return_value = 1;
 238     goto clean_up;
 239   }
 240
 241   // Number of frames to skip from file1.yuv for every frame used. Normal values
 242   // 0, 1 and 3 correspond to TL2, TL1 and TL0 respectively for a 3TL encoding
 243   // in mode 10. 7 would be reasonable for comparing TL0 of a 4-layer encoding.
 244   if (argc > 4) {
 245     sscanf(argv[4], "%d", &tl_skip);
 246     if (argc > 5) {
 247       framestats = fopen(argv[5], "w");
 248       if (!framestats) {
 249         fprintf(stderr, "Could not open \"%s\" for writing: %s\n", argv[5],
 250                 strerror(errno));
 251         return_value = 1;
 252         goto clean_up;
 253       }
 254     }
 255   }
 256
 257   if (w & 1 || h & 1) {
 258     fprintf(stderr, "Invalid size %dx%d\n", w, h);
 259     return_value = 1;
 260     goto clean_up;
 261   }
 262
 263   while (1) {
 264     size_t r1, r2;
 265     unsigned char *y[2], *u[2], *v[2];
 266
 267     r1 = read_input_file(&in[0], &y[0], &u[0], &v[0]);
 268
 269     if (r1) {
 270       // Reading parts of file1.yuv that were not used in temporal layer.
 271       if (tl_skips_remaining > 0) {
 272         --tl_skips_remaining;
 273         continue;
 274       }
 275       // Use frame, but skip |tl_skip| after it.
 276       tl_skips_remaining = tl_skip;
 277     }
 278
 279     r2 = read_input_file(&in[1], &y[1], &u[1], &v[1]);
 280
 281     if (r1 && r2 && r1 != r2) {
 282       fprintf(stderr, "Failed to read data: %s [%d/%d]\n", strerror(errno),
 283               (int)r1, (int)r2);
 284       return_value = 1;
 285       goto clean_up;
 286     } else if (r1 == 0 || r2 == 0) {
 287       break;
 288     }
 289 #define psnr_and_ssim(ssim, psnr, buf0, buf1, w, h) \
 290   ssim = vp8_ssim2(buf0, buf1, w, w, w, h);         \
 291   psnr = calc_plane_error(buf0, w, buf1, w, w, h);
 292
 293     if (n_frames == allocated_frames) {
 294       allocated_frames = allocated_frames == 0 ? 1024 : allocated_frames * 2;
 295       ssimy = realloc(ssimy, allocated_frames * sizeof(*ssimy));
 296       ssimu = realloc(ssimu, allocated_frames * sizeof(*ssimu));
 297       ssimv = realloc(ssimv, allocated_frames * sizeof(*ssimv));
 298       psnry = realloc(psnry, allocated_frames * sizeof(*psnry));
 299       psnru = realloc(psnru, allocated_frames * sizeof(*psnru));
 300       psnrv = realloc(psnrv, allocated_frames * sizeof(*psnrv));
 301     }
 302     psnr_and_ssim(ssimy[n_frames], psnry[n_frames], y[0], y[1], w, h);
 303     psnr_and_ssim(ssimu[n_frames], psnru[n_frames], u[0], u[1], w / 2, h / 2);
 304     psnr_and_ssim(ssimv[n_frames], psnrv[n_frames], v[0], v[1], w / 2, h / 2);
 305
 306     n_frames++;
 307   }
 308
 309   if (framestats) {
 310     fprintf(framestats,
 311             "ssim,ssim-y,ssim-u,ssim-v,psnr,psnr-y,psnr-u,psnr-v\n");
 312   }
 313
 314   for (i = 0; i < n_frames; ++i) {
 315     double frame_ssim;
 316     double frame_psnr, frame_psnry, frame_psnru, frame_psnrv;
 317
 318     frame_ssim = 0.8 * ssimy[i] + 0.1 * (ssimu[i] + ssimv[i]);
 319     ssimavg += frame_ssim;
 320     ssimyavg += ssimy[i];
 321     ssimuavg += ssimu[i];
 322     ssimvavg += ssimv[i];
 323
 324     frame_psnr = vp9_mse2psnr(w * h * 6 / 4, 255.0,
 325                               (double)psnry[i] + psnru[i] + psnrv[i]);
 326     frame_psnry = vp9_mse2psnr(w * h * 4 / 4, 255.0, (double)psnry[i]);
 327     frame_psnru = vp9_mse2psnr(w * h * 1 / 4, 255.0, (double)psnru[i]);
 328     frame_psnrv = vp9_mse2psnr(w * h * 1 / 4, 255.0, (double)psnrv[i]);
 329
 330     psnravg += frame_psnr;
 331     psnryavg += frame_psnry;
 332     psnruavg += frame_psnru;
 333     psnrvavg += frame_psnrv;
 334
 335     psnryglb += psnry[i];
 336     psnruglb += psnru[i];
 337     psnrvglb += psnrv[i];
 338
 339     if (framestats) {
 340       fprintf(framestats, "%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf\n", frame_ssim,
 341               ssimy[i], ssimu[i], ssimv[i], frame_psnr, frame_psnry,
 342               frame_psnru, frame_psnrv);
 343     }
 344   }
 345
 346   ssimavg /= n_frames;
 347   ssimyavg /= n_frames;
 348   ssimuavg /= n_frames;
 349   ssimvavg /= n_frames;
 350
 351   printf("VpxSSIM: %lf\n", 100 * pow(ssimavg, 8.0));
 352   printf("SSIM: %lf\n", ssimavg);
 353   printf("SSIM-Y: %lf\n", ssimyavg);
 354   printf("SSIM-U: %lf\n", ssimuavg);
 355   printf("SSIM-V: %lf\n", ssimvavg);
 356   puts("");
 357
 358   psnravg /= n_frames;
 359   psnryavg /= n_frames;
 360   psnruavg /= n_frames;
 361   psnrvavg /= n_frames;
 362
 363   printf("AvgPSNR: %lf\n", psnravg);
 364   printf("AvgPSNR-Y: %lf\n", psnryavg);
 365   printf("AvgPSNR-U: %lf\n", psnruavg);
 366   printf("AvgPSNR-V: %lf\n", psnrvavg);
 367   puts("");
 368
 369   psnrglb = psnryglb + psnruglb + psnrvglb;
 370   psnrglb = vp9_mse2psnr((double)n_frames * w * h * 6 / 4, 255.0, psnrglb);
 371   psnryglb = vp9_mse2psnr((double)n_frames * w * h * 4 / 4, 255.0, psnryglb);
 372   psnruglb = vp9_mse2psnr((double)n_frames * w * h * 1 / 4, 255.0, psnruglb);
 373   psnrvglb = vp9_mse2psnr((double)n_frames * w * h * 1 / 4, 255.0, psnrvglb);
 374
 375   printf("GlbPSNR: %lf\n", psnrglb);
 376   printf("GlbPSNR-Y: %lf\n", psnryglb);
 377   printf("GlbPSNR-U: %lf\n", psnruglb);
 378   printf("GlbPSNR-V: %lf\n", psnrvglb);
 379   puts("");
 380
 381   printf("Nframes: %d\n", (int)n_frames);
 382
 383 clean_up:
 384
 385   close_input_file(&in[0]);
 386   close_input_file(&in[1]);
 387
 388   if (framestats) fclose(framestats);
 389
 390   free(ssimy);
 391   free(ssimu);
 392   free(ssimv);
 393
 394   free(psnry);
 395   free(psnru);
 396   free(psnrv);
 397
 398   return return_value;
 399 }