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: Nathan E. Egge, at the Daala 11 * project. 12 */ 13#include <assert.h> 14#include <math.h> 15#include <stdlib.h> 16#include <string.h> 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 22typedef struct fs_level fs_level; 23typedef struct fs_ctx fs_ctx; 24 25#define SSIM_C1 (255 * 255 * 0.01 * 0.01) 26#define SSIM_C2 (255 * 255 * 0.03 * 0.03) 27#if CONFIG_VP9_HIGHBITDEPTH 28#define SSIM_C1_10 (1023 * 1023 * 0.01 * 0.01) 29#define SSIM_C1_12 (4095 * 4095 * 0.01 * 0.01) 30#define SSIM_C2_10 (1023 * 1023 * 0.03 * 0.03) 31#define SSIM_C2_12 (4095 * 4095 * 0.03 * 0.03) 32#endif 33#define FS_MINI(_a, _b) ((_a) < (_b) ? (_a) : (_b)) 34#define FS_MAXI(_a, _b) ((_a) > (_b) ? (_a) : (_b)) 35 36struct fs_level { 37 uint32_t *im1; 38 uint32_t *im2; 39 double *ssim; 40 int w; 41 int h; 42}; 43 44struct fs_ctx { 45 fs_level *level; 46 int nlevels; 47 unsigned *col_buf; 48}; 49 50static void fs_ctx_init(fs_ctx *_ctx, int _w, int _h, int _nlevels) { 51 unsigned char *data; 52 size_t data_size; 53 int lw; 54 int lh; 55 int l; 56 lw = (_w + 1) >> 1; 57 lh = (_h + 1) >> 1; 58 data_size = 59 _nlevels * sizeof(fs_level) + 2 * (lw + 8) * 8 * sizeof(*_ctx->col_buf); 60 for (l = 0; l < _nlevels; l++) { 61 size_t im_size; 62 size_t level_size; 63 im_size = lw * (size_t)lh; 64 level_size = 2 * im_size * sizeof(*_ctx->level[l].im1); 65 level_size += sizeof(*_ctx->level[l].ssim) - 1; 66 level_size /= sizeof(*_ctx->level[l].ssim); 67 level_size += im_size; 68 level_size *= sizeof(*_ctx->level[l].ssim); 69 data_size += level_size; 70 lw = (lw + 1) >> 1; 71 lh = (lh + 1) >> 1; 72 } 73 data = (unsigned char *)malloc(data_size); 74 _ctx->level = (fs_level *)data; 75 _ctx->nlevels = _nlevels; 76 data += _nlevels * sizeof(*_ctx->level); 77 lw = (_w + 1) >> 1; 78 lh = (_h + 1) >> 1; 79 for (l = 0; l < _nlevels; l++) { 80 size_t im_size; 81 size_t level_size; 82 _ctx->level[l].w = lw; 83 _ctx->level[l].h = lh; 84 im_size = lw * (size_t)lh; 85 level_size = 2 * im_size * sizeof(*_ctx->level[l].im1); 86 level_size += sizeof(*_ctx->level[l].ssim) - 1; 87 level_size /= sizeof(*_ctx->level[l].ssim); 88 level_size *= sizeof(*_ctx->level[l].ssim); 89 _ctx->level[l].im1 = (uint32_t *)data; 90 _ctx->level[l].im2 = _ctx->level[l].im1 + im_size; 91 data += level_size; 92 _ctx->level[l].ssim = (double *)data; 93 data += im_size * sizeof(*_ctx->level[l].ssim); 94 lw = (lw + 1) >> 1; 95 lh = (lh + 1) >> 1; 96 } 97 _ctx->col_buf = (unsigned *)data; 98} 99 100static void fs_ctx_clear(fs_ctx *_ctx) { free(_ctx->level); } 101 102static void fs_downsample_level(fs_ctx *_ctx, int _l) { 103 const uint32_t *src1; 104 const uint32_t *src2; 105 uint32_t *dst1; 106 uint32_t *dst2; 107 int w2; 108 int h2; 109 int w; 110 int h; 111 int i; 112 int j; 113 w = _ctx->level[_l].w; 114 h = _ctx->level[_l].h; 115 dst1 = _ctx->level[_l].im1; 116 dst2 = _ctx->level[_l].im2; 117 w2 = _ctx->level[_l - 1].w; 118 h2 = _ctx->level[_l - 1].h; 119 src1 = _ctx->level[_l - 1].im1; 120 src2 = _ctx->level[_l - 1].im2; 121 for (j = 0; j < h; j++) { 122 int j0offs; 123 int j1offs; 124 j0offs = 2 * j * w2; 125 j1offs = FS_MINI(2 * j + 1, h2) * w2; 126 for (i = 0; i < w; i++) { 127 int i0; 128 int i1; 129 i0 = 2 * i; 130 i1 = FS_MINI(i0 + 1, w2); 131 dst1[j * w + i] = src1[j0offs + i0] + src1[j0offs + i1] + 132 src1[j1offs + i0] + src1[j1offs + i1]; 133 dst2[j * w + i] = src2[j0offs + i0] + src2[j0offs + i1] + 134 src2[j1offs + i0] + src2[j1offs + i1]; 135 } 136 } 137} 138 139static void fs_downsample_level0(fs_ctx *_ctx, const uint8_t *_src1, 140 int _s1ystride, const uint8_t *_src2, 141 int _s2ystride, int _w, int _h, uint32_t bd, 142 uint32_t shift) { 143 uint32_t *dst1; 144 uint32_t *dst2; 145 int w; 146 int h; 147 int i; 148 int j; 149 w = _ctx->level[0].w; 150 h = _ctx->level[0].h; 151 dst1 = _ctx->level[0].im1; 152 dst2 = _ctx->level[0].im2; 153 for (j = 0; j < h; j++) { 154 int j0; 155 int j1; 156 j0 = 2 * j; 157 j1 = FS_MINI(j0 + 1, _h); 158 for (i = 0; i < w; i++) { 159 int i0; 160 int i1; 161 i0 = 2 * i; 162 i1 = FS_MINI(i0 + 1, _w); 163 if (bd == 8 && shift == 0) { 164 dst1[j * w + i] = 165 _src1[j0 * _s1ystride + i0] + _src1[j0 * _s1ystride + i1] + 166 _src1[j1 * _s1ystride + i0] + _src1[j1 * _s1ystride + i1]; 167 dst2[j * w + i] = 168 _src2[j0 * _s2ystride + i0] + _src2[j0 * _s2ystride + i1] + 169 _src2[j1 * _s2ystride + i0] + _src2[j1 * _s2ystride + i1]; 170 } else { 171 uint16_t *src1s = CONVERT_TO_SHORTPTR(_src1); 172 uint16_t *src2s = CONVERT_TO_SHORTPTR(_src2); 173 dst1[j * w + i] = (src1s[j0 * _s1ystride + i0] >> shift) + 174 (src1s[j0 * _s1ystride + i1] >> shift) + 175 (src1s[j1 * _s1ystride + i0] >> shift) + 176 (src1s[j1 * _s1ystride + i1] >> shift); 177 dst2[j * w + i] = (src2s[j0 * _s2ystride + i0] >> shift) + 178 (src2s[j0 * _s2ystride + i1] >> shift) + 179 (src2s[j1 * _s2ystride + i0] >> shift) + 180 (src2s[j1 * _s2ystride + i1] >> shift); 181 } 182 } 183 } 184} 185 186static void fs_apply_luminance(fs_ctx *_ctx, int _l, int bit_depth) { 187 unsigned *col_sums_x; 188 unsigned *col_sums_y; 189 uint32_t *im1; 190 uint32_t *im2; 191 double *ssim; 192 double c1; 193 int w; 194 int h; 195 int j0offs; 196 int j1offs; 197 int i; 198 int j; 199 double ssim_c1 = SSIM_C1; 200#if CONFIG_VP9_HIGHBITDEPTH 201 if (bit_depth == 10) ssim_c1 = SSIM_C1_10; 202 if (bit_depth == 12) ssim_c1 = SSIM_C1_12; 203#else 204 assert(bit_depth == 8); 205 (void)bit_depth; 206#endif 207 w = _ctx->level[_l].w; 208 h = _ctx->level[_l].h; 209 col_sums_x = _ctx->col_buf; 210 col_sums_y = col_sums_x + w; 211 im1 = _ctx->level[_l].im1; 212 im2 = _ctx->level[_l].im2; 213 for (i = 0; i < w; i++) col_sums_x[i] = 5 * im1[i]; 214 for (i = 0; i < w; i++) col_sums_y[i] = 5 * im2[i]; 215 for (j = 1; j < 4; j++) { 216 j1offs = FS_MINI(j, h - 1) * w; 217 for (i = 0; i < w; i++) col_sums_x[i] += im1[j1offs + i]; 218 for (i = 0; i < w; i++) col_sums_y[i] += im2[j1offs + i]; 219 } 220 ssim = _ctx->level[_l].ssim; 221 c1 = (double)(ssim_c1 * 4096 * (1 << 4 * _l)); 222 for (j = 0; j < h; j++) { 223 unsigned mux; 224 unsigned muy; 225 int i0; 226 int i1; 227 mux = 5 * col_sums_x[0]; 228 muy = 5 * col_sums_y[0]; 229 for (i = 1; i < 4; i++) { 230 i1 = FS_MINI(i, w - 1); 231 mux += col_sums_x[i1]; 232 muy += col_sums_y[i1]; 233 } 234 for (i = 0; i < w; i++) { 235 ssim[j * w + i] *= (2 * mux * (double)muy + c1) / 236 (mux * (double)mux + muy * (double)muy + c1); 237 if (i + 1 < w) { 238 i0 = FS_MAXI(0, i - 4); 239 i1 = FS_MINI(i + 4, w - 1); 240 mux += col_sums_x[i1] - col_sums_x[i0]; 241 muy += col_sums_x[i1] - col_sums_x[i0]; 242 } 243 } 244 if (j + 1 < h) { 245 j0offs = FS_MAXI(0, j - 4) * w; 246 for (i = 0; i < w; i++) col_sums_x[i] -= im1[j0offs + i]; 247 for (i = 0; i < w; i++) col_sums_y[i] -= im2[j0offs + i]; 248 j1offs = FS_MINI(j + 4, h - 1) * w; 249 for (i = 0; i < w; i++) col_sums_x[i] += im1[j1offs + i]; 250 for (i = 0; i < w; i++) col_sums_y[i] += im2[j1offs + i]; 251 } 252 } 253} 254 255#define FS_COL_SET(_col, _joffs, _ioffs) \ 256 do { \ 257 unsigned gx; \ 258 unsigned gy; \ 259 gx = gx_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \ 260 gy = gy_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \ 261 col_sums_gx2[(_col)] = gx * (double)gx; \ 262 col_sums_gy2[(_col)] = gy * (double)gy; \ 263 col_sums_gxgy[(_col)] = gx * (double)gy; \ 264 } while (0) 265 266#define FS_COL_ADD(_col, _joffs, _ioffs) \ 267 do { \ 268 unsigned gx; \ 269 unsigned gy; \ 270 gx = gx_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \ 271 gy = gy_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \ 272 col_sums_gx2[(_col)] += gx * (double)gx; \ 273 col_sums_gy2[(_col)] += gy * (double)gy; \ 274 col_sums_gxgy[(_col)] += gx * (double)gy; \ 275 } while (0) 276 277#define FS_COL_SUB(_col, _joffs, _ioffs) \ 278 do { \ 279 unsigned gx; \ 280 unsigned gy; \ 281 gx = gx_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \ 282 gy = gy_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \ 283 col_sums_gx2[(_col)] -= gx * (double)gx; \ 284 col_sums_gy2[(_col)] -= gy * (double)gy; \ 285 col_sums_gxgy[(_col)] -= gx * (double)gy; \ 286 } while (0) 287 288#define FS_COL_COPY(_col1, _col2) \ 289 do { \ 290 col_sums_gx2[(_col1)] = col_sums_gx2[(_col2)]; \ 291 col_sums_gy2[(_col1)] = col_sums_gy2[(_col2)]; \ 292 col_sums_gxgy[(_col1)] = col_sums_gxgy[(_col2)]; \ 293 } while (0) 294 295#define FS_COL_HALVE(_col1, _col2) \ 296 do { \ 297 col_sums_gx2[(_col1)] = col_sums_gx2[(_col2)] * 0.5; \ 298 col_sums_gy2[(_col1)] = col_sums_gy2[(_col2)] * 0.5; \ 299 col_sums_gxgy[(_col1)] = col_sums_gxgy[(_col2)] * 0.5; \ 300 } while (0) 301 302#define FS_COL_DOUBLE(_col1, _col2) \ 303 do { \ 304 col_sums_gx2[(_col1)] = col_sums_gx2[(_col2)] * 2; \ 305 col_sums_gy2[(_col1)] = col_sums_gy2[(_col2)] * 2; \ 306 col_sums_gxgy[(_col1)] = col_sums_gxgy[(_col2)] * 2; \ 307 } while (0) 308 309static void fs_calc_structure(fs_ctx *_ctx, int _l, int bit_depth) { 310 uint32_t *im1; 311 uint32_t *im2; 312 unsigned *gx_buf; 313 unsigned *gy_buf; 314 double *ssim; 315 double col_sums_gx2[8]; 316 double col_sums_gy2[8]; 317 double col_sums_gxgy[8]; 318 double c2; 319 int stride; 320 int w; 321 int h; 322 int i; 323 int j; 324 double ssim_c2 = SSIM_C2; 325#if CONFIG_VP9_HIGHBITDEPTH 326 if (bit_depth == 10) ssim_c2 = SSIM_C2_10; 327 if (bit_depth == 12) ssim_c2 = SSIM_C2_12; 328#else 329 assert(bit_depth == 8); 330 (void)bit_depth; 331#endif 332 333 w = _ctx->level[_l].w; 334 h = _ctx->level[_l].h; 335 im1 = _ctx->level[_l].im1; 336 im2 = _ctx->level[_l].im2; 337 ssim = _ctx->level[_l].ssim; 338 gx_buf = _ctx->col_buf; 339 stride = w + 8; 340 gy_buf = gx_buf + 8 * stride; 341 memset(gx_buf, 0, 2 * 8 * stride * sizeof(*gx_buf)); 342 c2 = ssim_c2 * (1 << 4 * _l) * 16 * 104; 343 for (j = 0; j < h + 4; j++) { 344 if (j < h - 1) { 345 for (i = 0; i < w - 1; i++) { 346 unsigned g1; 347 unsigned g2; 348 unsigned gx; 349 unsigned gy; 350 g1 = abs((int)im1[(j + 1) * w + i + 1] - (int)im1[j * w + i]); 351 g2 = abs((int)im1[(j + 1) * w + i] - (int)im1[j * w + i + 1]); 352 gx = 4 * FS_MAXI(g1, g2) + FS_MINI(g1, g2); 353 g1 = abs((int)im2[(j + 1) * w + i + 1] - (int)im2[j * w + i]); 354 g2 = abs((int)im2[(j + 1) * w + i] - (int)im2[j * w + i + 1]); 355 gy = 4 * FS_MAXI(g1, g2) + FS_MINI(g1, g2); 356 gx_buf[(j & 7) * stride + i + 4] = gx; 357 gy_buf[(j & 7) * stride + i + 4] = gy; 358 } 359 } else { 360 memset(gx_buf + (j & 7) * stride, 0, stride * sizeof(*gx_buf)); 361 memset(gy_buf + (j & 7) * stride, 0, stride * sizeof(*gy_buf)); 362 } 363 if (j >= 4) { 364 int k; 365 col_sums_gx2[3] = col_sums_gx2[2] = col_sums_gx2[1] = col_sums_gx2[0] = 0; 366 col_sums_gy2[3] = col_sums_gy2[2] = col_sums_gy2[1] = col_sums_gy2[0] = 0; 367 col_sums_gxgy[3] = col_sums_gxgy[2] = col_sums_gxgy[1] = 368 col_sums_gxgy[0] = 0; 369 for (i = 4; i < 8; i++) { 370 FS_COL_SET(i, -1, 0); 371 FS_COL_ADD(i, 0, 0); 372 for (k = 1; k < 8 - i; k++) { 373 FS_COL_DOUBLE(i, i); 374 FS_COL_ADD(i, -k - 1, 0); 375 FS_COL_ADD(i, k, 0); 376 } 377 } 378 for (i = 0; i < w; i++) { 379 double mugx2; 380 double mugy2; 381 double mugxgy; 382 mugx2 = col_sums_gx2[0]; 383 for (k = 1; k < 8; k++) mugx2 += col_sums_gx2[k]; 384 mugy2 = col_sums_gy2[0]; 385 for (k = 1; k < 8; k++) mugy2 += col_sums_gy2[k]; 386 mugxgy = col_sums_gxgy[0]; 387 for (k = 1; k < 8; k++) mugxgy += col_sums_gxgy[k]; 388 ssim[(j - 4) * w + i] = (2 * mugxgy + c2) / (mugx2 + mugy2 + c2); 389 if (i + 1 < w) { 390 FS_COL_SET(0, -1, 1); 391 FS_COL_ADD(0, 0, 1); 392 FS_COL_SUB(2, -3, 2); 393 FS_COL_SUB(2, 2, 2); 394 FS_COL_HALVE(1, 2); 395 FS_COL_SUB(3, -4, 3); 396 FS_COL_SUB(3, 3, 3); 397 FS_COL_HALVE(2, 3); 398 FS_COL_COPY(3, 4); 399 FS_COL_DOUBLE(4, 5); 400 FS_COL_ADD(4, -4, 5); 401 FS_COL_ADD(4, 3, 5); 402 FS_COL_DOUBLE(5, 6); 403 FS_COL_ADD(5, -3, 6); 404 FS_COL_ADD(5, 2, 6); 405 FS_COL_DOUBLE(6, 7); 406 FS_COL_ADD(6, -2, 7); 407 FS_COL_ADD(6, 1, 7); 408 FS_COL_SET(7, -1, 8); 409 FS_COL_ADD(7, 0, 8); 410 } 411 } 412 } 413 } 414} 415 416#define FS_NLEVELS (4) 417 418/*These weights were derived from the default weights found in Wang's original 419 Matlab implementation: {0.0448, 0.2856, 0.2363, 0.1333}. 420 We drop the finest scale and renormalize the rest to sum to 1.*/ 421 422static const double FS_WEIGHTS[FS_NLEVELS] = { 423 0.2989654541015625, 0.3141326904296875, 0.2473602294921875, 0.1395416259765625 424}; 425 426static double fs_average(fs_ctx *_ctx, int _l) { 427 double *ssim; 428 double ret; 429 int w; 430 int h; 431 int i; 432 int j; 433 w = _ctx->level[_l].w; 434 h = _ctx->level[_l].h; 435 ssim = _ctx->level[_l].ssim; 436 ret = 0; 437 for (j = 0; j < h; j++) 438 for (i = 0; i < w; i++) ret += ssim[j * w + i]; 439 return pow(ret / (w * h), FS_WEIGHTS[_l]); 440} 441 442static double convert_ssim_db(double _ssim, double _weight) { 443 assert(_weight >= _ssim); 444 if ((_weight - _ssim) < 1e-10) return MAX_SSIM_DB; 445 return 10 * (log10(_weight) - log10(_weight - _ssim)); 446} 447 448static double calc_ssim(const uint8_t *_src, int _systride, const uint8_t *_dst, 449 int _dystride, int _w, int _h, uint32_t _bd, 450 uint32_t _shift) { 451 fs_ctx ctx; 452 double ret; 453 int l; 454 ret = 1; 455 fs_ctx_init(&ctx, _w, _h, FS_NLEVELS); 456 fs_downsample_level0(&ctx, _src, _systride, _dst, _dystride, _w, _h, _bd, 457 _shift); 458 for (l = 0; l < FS_NLEVELS - 1; l++) { 459 fs_calc_structure(&ctx, l, _bd); 460 ret *= fs_average(&ctx, l); 461 fs_downsample_level(&ctx, l + 1); 462 } 463 fs_calc_structure(&ctx, l, _bd); 464 fs_apply_luminance(&ctx, l, _bd); 465 ret *= fs_average(&ctx, l); 466 fs_ctx_clear(&ctx); 467 return ret; 468} 469 470double vpx_calc_fastssim(const YV12_BUFFER_CONFIG *source, 471 const YV12_BUFFER_CONFIG *dest, double *ssim_y, 472 double *ssim_u, double *ssim_v, uint32_t bd, 473 uint32_t in_bd) { 474 double ssimv; 475 uint32_t bd_shift = 0; 476 vpx_clear_system_state(); 477 assert(bd >= in_bd); 478 bd_shift = bd - in_bd; 479 480 *ssim_y = calc_ssim(source->y_buffer, source->y_stride, dest->y_buffer, 481 dest->y_stride, source->y_crop_width, 482 source->y_crop_height, in_bd, bd_shift); 483 *ssim_u = calc_ssim(source->u_buffer, source->uv_stride, dest->u_buffer, 484 dest->uv_stride, source->uv_crop_width, 485 source->uv_crop_height, in_bd, bd_shift); 486 *ssim_v = calc_ssim(source->v_buffer, source->uv_stride, dest->v_buffer, 487 dest->uv_stride, source->uv_crop_width, 488 source->uv_crop_height, in_bd, bd_shift); 489 490 ssimv = (*ssim_y) * .8 + .1 * ((*ssim_u) + (*ssim_v)); 491 return convert_ssim_db(ssimv, 1.0); 492} 493