granicus.if.org Git - libvpx/blob - vp9/encoder/vp9_psnrhvs.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: Gregory Maxwell, at the Daala
  11  *  project.
  12  */
  13 #include <stdio.h>
  14 #include <stdlib.h>
  15 #include <math.h>
  16
  17 #include "./vpx_config.h"
  18 #include "./vp9_rtcd.h"
  19 #include "./vpx_dsp_rtcd.h"
  20 #include "vp9/encoder/vp9_ssim.h"
  21
  22 #if !defined(M_PI)
  23 # define M_PI (3.141592653589793238462643)
  24 #endif
  25 #include <string.h>
  26
  27 void od_bin_fdct8x8(tran_low_t *y, int ystride, const int16_t *x, int xstride) {
  28   (void) xstride;
  29   vpx_fdct8x8(x, y, ystride);
  30 }
  31
  32 /* Normalized inverse quantization matrix for 8x8 DCT at the point of
  33  * transparency. This is not the JPEG based matrix from the paper,
  34  this one gives a slightly higher MOS agreement.*/
  35 float csf_y[8][8] = {{1.6193873005, 2.2901594831, 2.08509755623, 1.48366094411,
  36     1.00227514334, 0.678296995242, 0.466224900598, 0.3265091542}, {2.2901594831,
  37     1.94321815382, 2.04793073064, 1.68731108984, 1.2305666963, 0.868920337363,
  38     0.61280991668, 0.436405793551}, {2.08509755623, 2.04793073064,
  39     1.34329019223, 1.09205635862, 0.875748795257, 0.670882927016,
  40     0.501731932449, 0.372504254596}, {1.48366094411, 1.68731108984,
  41     1.09205635862, 0.772819797575, 0.605636379554, 0.48309405692,
  42     0.380429446972, 0.295774038565}, {1.00227514334, 1.2305666963,
  43     0.875748795257, 0.605636379554, 0.448996256676, 0.352889268808,
  44     0.283006984131, 0.226951348204}, {0.678296995242, 0.868920337363,
  45     0.670882927016, 0.48309405692, 0.352889268808, 0.27032073436,
  46     0.215017739696, 0.17408067321}, {0.466224900598, 0.61280991668,
  47     0.501731932449, 0.380429446972, 0.283006984131, 0.215017739696,
  48     0.168869545842, 0.136153931001}, {0.3265091542, 0.436405793551,
  49     0.372504254596, 0.295774038565, 0.226951348204, 0.17408067321,
  50     0.136153931001, 0.109083846276}};
  51 float csf_cb420[8][8] = {
  52     {1.91113096927, 2.46074210438, 1.18284184739, 1.14982565193, 1.05017074788,
  53         0.898018824055, 0.74725392039, 0.615105596242}, {2.46074210438,
  54         1.58529308355, 1.21363250036, 1.38190029285, 1.33100189972,
  55         1.17428548929, 0.996404342439, 0.830890433625}, {1.18284184739,
  56         1.21363250036, 0.978712413627, 1.02624506078, 1.03145147362,
  57         0.960060382087, 0.849823426169, 0.731221236837}, {1.14982565193,
  58         1.38190029285, 1.02624506078, 0.861317501629, 0.801821139099,
  59         0.751437590932, 0.685398513368, 0.608694761374}, {1.05017074788,
  60         1.33100189972, 1.03145147362, 0.801821139099, 0.676555426187,
  61         0.605503172737, 0.55002013668, 0.495804539034}, {0.898018824055,
  62         1.17428548929, 0.960060382087, 0.751437590932, 0.605503172737,
  63         0.514674450957, 0.454353482512, 0.407050308965}, {0.74725392039,
  64         0.996404342439, 0.849823426169, 0.685398513368, 0.55002013668,
  65         0.454353482512, 0.389234902883, 0.342353999733}, {0.615105596242,
  66         0.830890433625, 0.731221236837, 0.608694761374, 0.495804539034,
  67         0.407050308965, 0.342353999733, 0.295530605237}};
  68 float csf_cr420[8][8] = {
  69     {2.03871978502, 2.62502345193, 1.26180942886, 1.11019789803, 1.01397751469,
  70         0.867069376285, 0.721500455585, 0.593906509971}, {2.62502345193,
  71         1.69112867013, 1.17180569821, 1.3342742857, 1.28513006198,
  72         1.13381474809, 0.962064122248, 0.802254508198}, {1.26180942886,
  73         1.17180569821, 0.944981930573, 0.990876405848, 0.995903384143,
  74         0.926972725286, 0.820534991409, 0.706020324706}, {1.11019789803,
  75         1.3342742857, 0.990876405848, 0.831632933426, 0.77418706195,
  76         0.725539939514, 0.661776842059, 0.587716619023}, {1.01397751469,
  77         1.28513006198, 0.995903384143, 0.77418706195, 0.653238524286,
  78         0.584635025748, 0.531064164893, 0.478717061273}, {0.867069376285,
  79         1.13381474809, 0.926972725286, 0.725539939514, 0.584635025748,
  80         0.496936637883, 0.438694579826, 0.393021669543}, {0.721500455585,
  81         0.962064122248, 0.820534991409, 0.661776842059, 0.531064164893,
  82         0.438694579826, 0.375820256136, 0.330555063063}, {0.593906509971,
  83         0.802254508198, 0.706020324706, 0.587716619023, 0.478717061273,
  84         0.393021669543, 0.330555063063, 0.285345396658}};
  85
  86 static double convert_score_db(double _score, double _weight) {
  87   return 10 * (log10(255 * 255) - log10(_weight * _score));
  88 }
  89
  90 static double calc_psnrhvs(const unsigned char *_src, int _systride,
  91                            const unsigned char *_dst, int _dystride,
  92                            double _par, int _w, int _h, int _step,
  93                            float _csf[8][8]) {
  94   float ret;
  95   int16_t dct_s[8 * 8], dct_d[8 * 8];
  96   tran_low_t dct_s_coef[8 * 8], dct_d_coef[8 * 8];
  97   float mask[8][8];
  98   int pixels;
  99   int x;
 100   int y;
 101   (void) _par;
 102   ret = pixels = 0;
 103   /*In the PSNR-HVS-M paper[1] the authors describe the construction of
 104    their masking table as "we have used the quantization table for the
 105    color component Y of JPEG [6] that has been also obtained on the
 106    basis of CSF. Note that the values in quantization table JPEG have
 107    been normalized and then squared." Their CSF matrix (from PSNR-HVS)
 108    was also constructed from the JPEG matrices. I can not find any obvious
 109    scheme of normalizing to produce their table, but if I multiply their
 110    CSF by 0.38857 and square the result I get their masking table.
 111    I have no idea where this constant comes from, but deviating from it
 112    too greatly hurts MOS agreement.
 113
 114    [1] Nikolay Ponomarenko, Flavia Silvestri, Karen Egiazarian, Marco Carli,
 115    Jaakko Astola, Vladimir Lukin, "On between-coefficient contrast masking
 116    of DCT basis functions", CD-ROM Proceedings of the Third
 117    International Workshop on Video Processing and Quality Metrics for Consumer
 118    Electronics VPQM-07, Scottsdale, Arizona, USA, 25-26 January, 2007, 4 p.*/
 119   for (x = 0; x < 8; x++)
 120     for (y = 0; y < 8; y++)
 121       mask[x][y] = (_csf[x][y] * 0.3885746225901003)
 122           * (_csf[x][y] * 0.3885746225901003);
 123   for (y = 0; y < _h - 7; y += _step) {
 124     for (x = 0; x < _w - 7; x += _step) {
 125       int i;
 126       int j;
 127       float s_means[4];
 128       float d_means[4];
 129       float s_vars[4];
 130       float d_vars[4];
 131       float s_gmean = 0;
 132       float d_gmean = 0;
 133       float s_gvar = 0;
 134       float d_gvar = 0;
 135       float s_mask = 0;
 136       float d_mask = 0;
 137       for (i = 0; i < 4; i++)
 138         s_means[i] = d_means[i] = s_vars[i] = d_vars[i] = 0;
 139       for (i = 0; i < 8; i++) {
 140         for (j = 0; j < 8; j++) {
 141           int sub = ((i & 12) >> 2) + ((j & 12) >> 1);
 142           dct_s[i * 8 + j] = _src[(y + i) * _systride + (j + x)];
 143           dct_d[i * 8 + j] = _dst[(y + i) * _dystride + (j + x)];
 144           s_gmean += dct_s[i * 8 + j];
 145           d_gmean += dct_d[i * 8 + j];
 146           s_means[sub] += dct_s[i * 8 + j];
 147           d_means[sub] += dct_d[i * 8 + j];
 148         }
 149       }
 150       s_gmean /= 64.f;
 151       d_gmean /= 64.f;
 152       for (i = 0; i < 4; i++)
 153         s_means[i] /= 16.f;
 154       for (i = 0; i < 4; i++)
 155         d_means[i] /= 16.f;
 156       for (i = 0; i < 8; i++) {
 157         for (j = 0; j < 8; j++) {
 158           int sub = ((i & 12) >> 2) + ((j & 12) >> 1);
 159           s_gvar += (dct_s[i * 8 + j] - s_gmean) * (dct_s[i * 8 + j] - s_gmean);
 160           d_gvar += (dct_d[i * 8 + j] - d_gmean) * (dct_d[i * 8 + j] - d_gmean);
 161           s_vars[sub] += (dct_s[i * 8 + j] - s_means[sub])
 162               * (dct_s[i * 8 + j] - s_means[sub]);
 163           d_vars[sub] += (dct_d[i * 8 + j] - d_means[sub])
 164               * (dct_d[i * 8 + j] - d_means[sub]);
 165         }
 166       }
 167       s_gvar *= 1 / 63.f * 64;
 168       d_gvar *= 1 / 63.f * 64;
 169       for (i = 0; i < 4; i++)
 170         s_vars[i] *= 1 / 15.f * 16;
 171       for (i = 0; i < 4; i++)
 172         d_vars[i] *= 1 / 15.f * 16;
 173       if (s_gvar > 0)
 174         s_gvar = (s_vars[0] + s_vars[1] + s_vars[2] + s_vars[3]) / s_gvar;
 175       if (d_gvar > 0)
 176         d_gvar = (d_vars[0] + d_vars[1] + d_vars[2] + d_vars[3]) / d_gvar;
 177       od_bin_fdct8x8(dct_s_coef, 8, dct_s, 8);
 178       od_bin_fdct8x8(dct_d_coef, 8, dct_d, 8);
 179       for (i = 0; i < 8; i++)
 180         for (j = (i == 0); j < 8; j++)
 181           s_mask += dct_s_coef[i * 8 + j] * dct_s_coef[i * 8 + j] * mask[i][j];
 182       for (i = 0; i < 8; i++)
 183         for (j = (i == 0); j < 8; j++)
 184           d_mask += dct_d_coef[i * 8 + j] * dct_d_coef[i * 8 + j] * mask[i][j];
 185       s_mask = sqrt(s_mask * s_gvar) / 32.f;
 186       d_mask = sqrt(d_mask * d_gvar) / 32.f;
 187       if (d_mask > s_mask)
 188         s_mask = d_mask;
 189       for (i = 0; i < 8; i++) {
 190         for (j = 0; j < 8; j++) {
 191           float err;
 192           err = fabs(dct_s_coef[i * 8 + j] - dct_d_coef[i * 8 + j]);
 193           if (i != 0 || j != 0)
 194             err = err < s_mask / mask[i][j] ? 0 : err - s_mask / mask[i][j];
 195           ret += (err * _csf[i][j]) * (err * _csf[i][j]);
 196           pixels++;
 197         }
 198       }
 199     }
 200   }
 201   ret /= pixels;
 202   return ret;
 203 }
 204 double vp9_psnrhvs(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
 205                    double *y_psnrhvs, double *u_psnrhvs, double *v_psnrhvs) {
 206   double psnrhvs;
 207   double par = 1.0;
 208   int step = 7;
 209   vp9_clear_system_state();
 210   *y_psnrhvs = calc_psnrhvs(source->y_buffer, source->y_stride, dest->y_buffer,
 211                             dest->y_stride, par, source->y_crop_width,
 212                             source->y_crop_height, step, csf_y);
 213
 214   *u_psnrhvs = calc_psnrhvs(source->u_buffer, source->uv_stride, dest->u_buffer,
 215                             dest->uv_stride, par, source->uv_crop_width,
 216                             source->uv_crop_height, step, csf_cb420);
 217
 218   *v_psnrhvs = calc_psnrhvs(source->v_buffer, source->uv_stride, dest->v_buffer,
 219                             dest->uv_stride, par, source->uv_crop_width,
 220                             source->uv_crop_height, step, csf_cr420);
 221   psnrhvs = (*y_psnrhvs) * .8 + .1 * ((*u_psnrhvs) + (*v_psnrhvs));
 222
 223   return convert_score_db(psnrhvs, 1.0);
 224 }