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