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