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