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