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 
od_bin_fdct8x8(tran_low_t * y,int ystride,const int16_t * x,int xstride)27 static 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.*/
36 static 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}};
53 static 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}};
70 static 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 
convert_score_db(double _score,double _weight)88 static double convert_score_db(double _score, double _weight) {
89   return 10 * (log10(255 * 255) - log10(_weight * _score));
90 }
91 
calc_psnrhvs(const unsigned char * _src,int _systride,const unsigned char * _dst,int _dystride,double _par,int _w,int _h,int _step,const float _csf[8][8])92 static 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 }
vpx_psnrhvs(const YV12_BUFFER_CONFIG * source,const YV12_BUFFER_CONFIG * dest,double * y_psnrhvs,double * u_psnrhvs,double * v_psnrhvs)206 double 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