vpx_dsp/fastssim.c

   1 /*
   2  *  Copyright (c) 2010 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  *  This code was originally written by: Nathan E. Egge, at the Daala
  11  *  project.
  12  */
  13 #include <math.h>
  14 #include <stdlib.h>
  15 #include <string.h>
  16 #include "./vpx_config.h"
  17 #include "./vpx_dsp_rtcd.h"
  18 #include "vpx_dsp/ssim.h"
  19 #include "vpx_ports/system_state.h"
  20 /* TODO(jbb): High bit depth version of this code needed */
  21 typedef struct fs_level fs_level;
  22 typedef struct fs_ctx fs_ctx;
  23
  24 #define SSIM_C1 (255 * 255 * 0.01 * 0.01)
  25 #define SSIM_C2 (255 * 255 * 0.03 * 0.03)
  26
  27 #define FS_MINI(_a, _b) ((_a) < (_b) ? (_a) : (_b))
  28 #define FS_MAXI(_a, _b) ((_a) > (_b) ? (_a) : (_b))
  29
  30 struct fs_level {
  31   uint16_t *im1;
  32   uint16_t *im2;
  33   double *ssim;
  34   int w;
  35   int h;
  36 };
  37
  38 struct fs_ctx {
  39   fs_level *level;
  40   int nlevels;
  41   unsigned *col_buf;
  42 };
  43
  44 static void fs_ctx_init(fs_ctx *_ctx, int _w, int _h, int _nlevels) {
  45   unsigned char *data;
  46   size_t data_size;
  47   int lw;
  48   int lh;
  49   int l;
  50   lw = (_w + 1) >> 1;
  51   lh = (_h + 1) >> 1;
  52   data_size = _nlevels * sizeof(fs_level)
  53       + 2 * (lw + 8) * 8 * sizeof(*_ctx->col_buf);
  54   for (l = 0; l < _nlevels; l++) {
  55     size_t im_size;
  56     size_t level_size;
  57     im_size = lw * (size_t) lh;
  58     level_size = 2 * im_size * sizeof(*_ctx->level[l].im1);
  59     level_size += sizeof(*_ctx->level[l].ssim) - 1;
  60     level_size /= sizeof(*_ctx->level[l].ssim);
  61     level_size += im_size;
  62     level_size *= sizeof(*_ctx->level[l].ssim);
  63     data_size += level_size;
  64     lw = (lw + 1) >> 1;
  65     lh = (lh + 1) >> 1;
  66   }
  67   data = (unsigned char *) malloc(data_size);
  68   _ctx->level = (fs_level *) data;
  69   _ctx->nlevels = _nlevels;
  70   data += _nlevels * sizeof(*_ctx->level);
  71   lw = (_w + 1) >> 1;
  72   lh = (_h + 1) >> 1;
  73   for (l = 0; l < _nlevels; l++) {
  74     size_t im_size;
  75     size_t level_size;
  76     _ctx->level[l].w = lw;
  77     _ctx->level[l].h = lh;
  78     im_size = lw * (size_t) lh;
  79     level_size = 2 * im_size * sizeof(*_ctx->level[l].im1);
  80     level_size += sizeof(*_ctx->level[l].ssim) - 1;
  81     level_size /= sizeof(*_ctx->level[l].ssim);
  82     level_size *= sizeof(*_ctx->level[l].ssim);
  83     _ctx->level[l].im1 = (uint16_t *) data;
  84     _ctx->level[l].im2 = _ctx->level[l].im1 + im_size;
  85     data += level_size;
  86     _ctx->level[l].ssim = (double *) data;
  87     data += im_size * sizeof(*_ctx->level[l].ssim);
  88     lw = (lw + 1) >> 1;
  89     lh = (lh + 1) >> 1;
  90   }
  91   _ctx->col_buf = (unsigned *) data;
  92 }
  93
  94 static void fs_ctx_clear(fs_ctx *_ctx) {
  95   free(_ctx->level);
  96 }
  97
  98 static void fs_downsample_level(fs_ctx *_ctx, int _l) {
  99   const uint16_t *src1;
 100   const uint16_t *src2;
 101   uint16_t *dst1;
 102   uint16_t *dst2;
 103   int w2;
 104   int h2;
 105   int w;
 106   int h;
 107   int i;
 108   int j;
 109   w = _ctx->level[_l].w;
 110   h = _ctx->level[_l].h;
 111   dst1 = _ctx->level[_l].im1;
 112   dst2 = _ctx->level[_l].im2;
 113   w2 = _ctx->level[_l - 1].w;
 114   h2 = _ctx->level[_l - 1].h;
 115   src1 = _ctx->level[_l - 1].im1;
 116   src2 = _ctx->level[_l - 1].im2;
 117   for (j = 0; j < h; j++) {
 118     int j0offs;
 119     int j1offs;
 120     j0offs = 2 * j * w2;
 121     j1offs = FS_MINI(2 * j + 1, h2) * w2;
 122     for (i = 0; i < w; i++) {
 123       int i0;
 124       int i1;
 125       i0 = 2 * i;
 126       i1 = FS_MINI(i0 + 1, w2);
 127       dst1[j * w + i] = src1[j0offs + i0] + src1[j0offs + i1]
 128           + src1[j1offs + i0] + src1[j1offs + i1];
 129       dst2[j * w + i] = src2[j0offs + i0] + src2[j0offs + i1]
 130           + src2[j1offs + i0] + src2[j1offs + i1];
 131     }
 132   }
 133 }
 134
 135 static void fs_downsample_level0(fs_ctx *_ctx, const unsigned char *_src1,
 136                                  int _s1ystride, const unsigned char *_src2,
 137                                  int _s2ystride, int _w, int _h) {
 138   uint16_t *dst1;
 139   uint16_t *dst2;
 140   int w;
 141   int h;
 142   int i;
 143   int j;
 144   w = _ctx->level[0].w;
 145   h = _ctx->level[0].h;
 146   dst1 = _ctx->level[0].im1;
 147   dst2 = _ctx->level[0].im2;
 148   for (j = 0; j < h; j++) {
 149     int j0;
 150     int j1;
 151     j0 = 2 * j;
 152     j1 = FS_MINI(j0 + 1, _h);
 153     for (i = 0; i < w; i++) {
 154       int i0;
 155       int i1;
 156       i0 = 2 * i;
 157       i1 = FS_MINI(i0 + 1, _w);
 158       dst1[j * w + i] = _src1[j0 * _s1ystride + i0]
 159           + _src1[j0 * _s1ystride + i1] + _src1[j1 * _s1ystride + i0]
 160           + _src1[j1 * _s1ystride + i1];
 161       dst2[j * w + i] = _src2[j0 * _s2ystride + i0]
 162           + _src2[j0 * _s2ystride + i1] + _src2[j1 * _s2ystride + i0]
 163           + _src2[j1 * _s2ystride + i1];
 164     }
 165   }
 166 }
 167
 168 static void fs_apply_luminance(fs_ctx *_ctx, int _l) {
 169   unsigned *col_sums_x;
 170   unsigned *col_sums_y;
 171   uint16_t *im1;
 172   uint16_t *im2;
 173   double *ssim;
 174   double c1;
 175   int w;
 176   int h;
 177   int j0offs;
 178   int j1offs;
 179   int i;
 180   int j;
 181   w = _ctx->level[_l].w;
 182   h = _ctx->level[_l].h;
 183   col_sums_x = _ctx->col_buf;
 184   col_sums_y = col_sums_x + w;
 185   im1 = _ctx->level[_l].im1;
 186   im2 = _ctx->level[_l].im2;
 187   for (i = 0; i < w; i++)
 188     col_sums_x[i] = 5 * im1[i];
 189   for (i = 0; i < w; i++)
 190     col_sums_y[i] = 5 * im2[i];
 191   for (j = 1; j < 4; j++) {
 192     j1offs = FS_MINI(j, h - 1) * w;
 193     for (i = 0; i < w; i++)
 194       col_sums_x[i] += im1[j1offs + i];
 195     for (i = 0; i < w; i++)
 196       col_sums_y[i] += im2[j1offs + i];
 197   }
 198   ssim = _ctx->level[_l].ssim;
 199   c1 = (double) (SSIM_C1 * 4096 * (1 << 4 * _l));
 200   for (j = 0; j < h; j++) {
 201     unsigned mux;
 202     unsigned muy;
 203     int i0;
 204     int i1;
 205     mux = 5 * col_sums_x[0];
 206     muy = 5 * col_sums_y[0];
 207     for (i = 1; i < 4; i++) {
 208       i1 = FS_MINI(i, w - 1);
 209       mux += col_sums_x[i1];
 210       muy += col_sums_y[i1];
 211     }
 212     for (i = 0; i < w; i++) {
 213       ssim[j * w + i] *= (2 * mux * (double) muy + c1)
 214           / (mux * (double) mux + muy * (double) muy + c1);
 215       if (i + 1 < w) {
 216         i0 = FS_MAXI(0, i - 4);
 217         i1 = FS_MINI(i + 4, w - 1);
 218         mux += col_sums_x[i1] - col_sums_x[i0];
 219         muy += col_sums_x[i1] - col_sums_x[i0];
 220       }
 221     }
 222     if (j + 1 < h) {
 223       j0offs = FS_MAXI(0, j - 4) * w;
 224       for (i = 0; i < w; i++)
 225         col_sums_x[i] -= im1[j0offs + i];
 226       for (i = 0; i < w; i++)
 227         col_sums_y[i] -= im2[j0offs + i];
 228       j1offs = FS_MINI(j + 4, h - 1) * w;
 229       for (i = 0; i < w; i++)
 230         col_sums_x[i] += im1[j1offs + i];
 231       for (i = 0; i < w; i++)
 232         col_sums_y[i] += im2[j1offs + i];
 233     }
 234   }
 235 }
 236
 237 #define FS_COL_SET(_col, _joffs, _ioffs) \
 238   do { \
 239     unsigned gx; \
 240     unsigned gy; \
 241     gx = gx_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
 242     gy = gy_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
 243     col_sums_gx2[(_col)] = gx * (double)gx; \
 244     col_sums_gy2[(_col)] = gy * (double)gy; \
 245     col_sums_gxgy[(_col)] = gx * (double)gy; \
 246   } \
 247   while (0)
 248
 249 #define FS_COL_ADD(_col, _joffs, _ioffs) \
 250   do { \
 251     unsigned gx; \
 252     unsigned gy; \
 253     gx = gx_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
 254     gy = gy_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
 255     col_sums_gx2[(_col)] += gx * (double)gx; \
 256     col_sums_gy2[(_col)] += gy * (double)gy; \
 257     col_sums_gxgy[(_col)] += gx * (double)gy; \
 258   } \
 259   while (0)
 260
 261 #define FS_COL_SUB(_col, _joffs, _ioffs) \
 262   do { \
 263     unsigned gx; \
 264     unsigned gy; \
 265     gx = gx_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
 266     gy = gy_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
 267     col_sums_gx2[(_col)] -= gx * (double)gx; \
 268     col_sums_gy2[(_col)] -= gy * (double)gy; \
 269     col_sums_gxgy[(_col)] -= gx * (double)gy; \
 270   } \
 271   while (0)
 272
 273 #define FS_COL_COPY(_col1, _col2) \
 274   do { \
 275     col_sums_gx2[(_col1)] = col_sums_gx2[(_col2)]; \
 276     col_sums_gy2[(_col1)] = col_sums_gy2[(_col2)]; \
 277     col_sums_gxgy[(_col1)] = col_sums_gxgy[(_col2)]; \
 278   } \
 279   while (0)
 280
 281 #define FS_COL_HALVE(_col1, _col2) \
 282   do { \
 283     col_sums_gx2[(_col1)] = col_sums_gx2[(_col2)] * 0.5; \
 284     col_sums_gy2[(_col1)] = col_sums_gy2[(_col2)] * 0.5; \
 285     col_sums_gxgy[(_col1)] = col_sums_gxgy[(_col2)] * 0.5; \
 286   } \
 287   while (0)
 288
 289 #define FS_COL_DOUBLE(_col1, _col2) \
 290   do { \
 291     col_sums_gx2[(_col1)] = col_sums_gx2[(_col2)] * 2; \
 292     col_sums_gy2[(_col1)] = col_sums_gy2[(_col2)] * 2; \
 293     col_sums_gxgy[(_col1)] = col_sums_gxgy[(_col2)] * 2; \
 294   } \
 295   while (0)
 296
 297 static void fs_calc_structure(fs_ctx *_ctx, int _l) {
 298   uint16_t *im1;
 299   uint16_t *im2;
 300   unsigned *gx_buf;
 301   unsigned *gy_buf;
 302   double *ssim;
 303   double col_sums_gx2[8];
 304   double col_sums_gy2[8];
 305   double col_sums_gxgy[8];
 306   double c2;
 307   int stride;
 308   int w;
 309   int h;
 310   int i;
 311   int j;
 312   w = _ctx->level[_l].w;
 313   h = _ctx->level[_l].h;
 314   im1 = _ctx->level[_l].im1;
 315   im2 = _ctx->level[_l].im2;
 316   ssim = _ctx->level[_l].ssim;
 317   gx_buf = _ctx->col_buf;
 318   stride = w + 8;
 319   gy_buf = gx_buf + 8 * stride;
 320   memset(gx_buf, 0, 2 * 8 * stride * sizeof(*gx_buf));
 321   c2 = SSIM_C2 * (1 << 4 * _l) * 16 * 104;
 322   for (j = 0; j < h + 4; j++) {
 323     if (j < h - 1) {
 324       for (i = 0; i < w - 1; i++) {
 325         unsigned g1;
 326         unsigned g2;
 327         unsigned gx;
 328         unsigned gy;
 329         g1 = abs(im1[(j + 1) * w + i + 1] - im1[j * w + i]);
 330         g2 = abs(im1[(j + 1) * w + i] - im1[j * w + i + 1]);
 331         gx = 4 * FS_MAXI(g1, g2) + FS_MINI(g1, g2);
 332         g1 = abs(im2[(j + 1) * w + i + 1] - im2[j * w + i]);
 333         g2 = abs(im2[(j + 1) * w + i] - im2[j * w + i + 1]);
 334         gy = 4 * FS_MAXI(g1, g2) + FS_MINI(g1, g2);
 335         gx_buf[(j & 7) * stride + i + 4] = gx;
 336         gy_buf[(j & 7) * stride + i + 4] = gy;
 337       }
 338     } else {
 339       memset(gx_buf + (j & 7) * stride, 0, stride * sizeof(*gx_buf));
 340       memset(gy_buf + (j & 7) * stride, 0, stride * sizeof(*gy_buf));
 341     }
 342     if (j >= 4) {
 343       int k;
 344       col_sums_gx2[3] = col_sums_gx2[2] = col_sums_gx2[1] = col_sums_gx2[0] = 0;
 345       col_sums_gy2[3] = col_sums_gy2[2] = col_sums_gy2[1] = col_sums_gy2[0] = 0;
 346       col_sums_gxgy[3] = col_sums_gxgy[2] = col_sums_gxgy[1] =
 347           col_sums_gxgy[0] = 0;
 348       for (i = 4; i < 8; i++) {
 349         FS_COL_SET(i, -1, 0);
 350         FS_COL_ADD(i, 0, 0);
 351         for (k = 1; k < 8 - i; k++) {
 352           FS_COL_DOUBLE(i, i);
 353           FS_COL_ADD(i, -k - 1, 0);
 354           FS_COL_ADD(i, k, 0);
 355         }
 356       }
 357       for (i = 0; i < w; i++) {
 358         double mugx2;
 359         double mugy2;
 360         double mugxgy;
 361         mugx2 = col_sums_gx2[0];
 362         for (k = 1; k < 8; k++)
 363           mugx2 += col_sums_gx2[k];
 364         mugy2 = col_sums_gy2[0];
 365         for (k = 1; k < 8; k++)
 366           mugy2 += col_sums_gy2[k];
 367         mugxgy = col_sums_gxgy[0];
 368         for (k = 1; k < 8; k++)
 369           mugxgy += col_sums_gxgy[k];
 370         ssim[(j - 4) * w + i] = (2 * mugxgy + c2) / (mugx2 + mugy2 + c2);
 371         if (i + 1 < w) {
 372           FS_COL_SET(0, -1, 1);
 373           FS_COL_ADD(0, 0, 1);
 374           FS_COL_SUB(2, -3, 2);
 375           FS_COL_SUB(2, 2, 2);
 376           FS_COL_HALVE(1, 2);
 377           FS_COL_SUB(3, -4, 3);
 378           FS_COL_SUB(3, 3, 3);
 379           FS_COL_HALVE(2, 3);
 380           FS_COL_COPY(3, 4);
 381           FS_COL_DOUBLE(4, 5);
 382           FS_COL_ADD(4, -4, 5);
 383           FS_COL_ADD(4, 3, 5);
 384           FS_COL_DOUBLE(5, 6);
 385           FS_COL_ADD(5, -3, 6);
 386           FS_COL_ADD(5, 2, 6);
 387           FS_COL_DOUBLE(6, 7);
 388           FS_COL_ADD(6, -2, 7);
 389           FS_COL_ADD(6, 1, 7);
 390           FS_COL_SET(7, -1, 8);
 391           FS_COL_ADD(7, 0, 8);
 392         }
 393       }
 394     }
 395   }
 396 }
 397
 398 #define FS_NLEVELS (4)
 399
 400 /*These weights were derived from the default weights found in Wang's original
 401  Matlab implementation: {0.0448, 0.2856, 0.2363, 0.1333}.
 402  We drop the finest scale and renormalize the rest to sum to 1.*/
 403
 404 static const double FS_WEIGHTS[FS_NLEVELS] = {0.2989654541015625,
 405     0.3141326904296875, 0.2473602294921875, 0.1395416259765625};
 406
 407 static double fs_average(fs_ctx *_ctx, int _l) {
 408   double *ssim;
 409   double ret;
 410   int w;
 411   int h;
 412   int i;
 413   int j;
 414   w = _ctx->level[_l].w;
 415   h = _ctx->level[_l].h;
 416   ssim = _ctx->level[_l].ssim;
 417   ret = 0;
 418   for (j = 0; j < h; j++)
 419     for (i = 0; i < w; i++)
 420       ret += ssim[j * w + i];
 421   return pow(ret / (w * h), FS_WEIGHTS[_l]);
 422 }
 423
 424 static double calc_ssim(const unsigned char *_src, int _systride,
 425                  const unsigned char *_dst, int _dystride, int _w, int _h) {
 426   fs_ctx ctx;
 427   double ret;
 428   int l;
 429   ret = 1;
 430   fs_ctx_init(&ctx, _w, _h, FS_NLEVELS);
 431   fs_downsample_level0(&ctx, _src, _systride, _dst, _dystride, _w, _h);
 432   for (l = 0; l < FS_NLEVELS - 1; l++) {
 433     fs_calc_structure(&ctx, l);
 434     ret *= fs_average(&ctx, l);
 435     fs_downsample_level(&ctx, l + 1);
 436   }
 437   fs_calc_structure(&ctx, l);
 438   fs_apply_luminance(&ctx, l);
 439   ret *= fs_average(&ctx, l);
 440   fs_ctx_clear(&ctx);
 441   return ret;
 442 }
 443
 444 static double convert_ssim_db(double _ssim, double _weight) {
 445   return 10 * (log10(_weight) - log10(_weight - _ssim));
 446 }
 447
 448 double vpx_calc_fastssim(const YV12_BUFFER_CONFIG *source,
 449                          const YV12_BUFFER_CONFIG *dest,
 450                          double *ssim_y, double *ssim_u, double *ssim_v) {
 451   double ssimv;
 452   vpx_clear_system_state();
 453
 454   *ssim_y = calc_ssim(source->y_buffer, source->y_stride, dest->y_buffer,
 455                       dest->y_stride, source->y_crop_width,
 456                       source->y_crop_height);
 457
 458   *ssim_u = calc_ssim(source->u_buffer, source->uv_stride, dest->u_buffer,
 459                       dest->uv_stride, source->uv_crop_width,
 460                       source->uv_crop_height);
 461
 462   *ssim_v = calc_ssim(source->v_buffer, source->uv_stride, dest->v_buffer,
 463                       dest->uv_stride, source->uv_crop_width,
 464                       source->uv_crop_height);
 465   ssimv = (*ssim_y) * .8 + .1 * ((*ssim_u) + (*ssim_v));
 466
 467   return convert_ssim_db(ssimv, 1.0);
 468 }