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