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