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