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 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 88 static double convert_score_db(double _score, double _weight) { 89 return 10 * (log10(255 * 255) - log10(_weight * _score)); 90 } 91 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 } 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
