2 * Copyright (c) 2010 The WebM project authors. All Rights Reserved.
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.
10 * This code was originally written by: Nathan E. Egge, at the Daala
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;
22 #define SSIM_C1 (255 * 255 * 0.01 * 0.01)
23 #define SSIM_C2 (255 * 255 * 0.03 * 0.03)
25 #define FS_MINI(_a, _b) ((_a) < (_b) ? (_a) : (_b))
26 #define FS_MAXI(_a, _b) ((_a) > (_b) ? (_a) : (_b))
42 static void fs_ctx_init(fs_ctx *_ctx, int _w, int _h, int _nlevels) {
50 data_size = _nlevels * sizeof(fs_level)
51 + 2 * (lw + 8) * 8 * sizeof(*_ctx->col_buf);
52 for (l = 0; l < _nlevels; l++) {
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;
65 data = (unsigned char *) malloc(data_size);
66 _ctx->level = (fs_level *) data;
67 _ctx->nlevels = _nlevels;
68 data += _nlevels * sizeof(*_ctx->level);
71 for (l = 0; l < _nlevels; l++) {
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;
84 _ctx->level[l].ssim = (double *) data;
85 data += im_size * sizeof(*_ctx->level[l].ssim);
89 _ctx->col_buf = (unsigned *) data;
92 static void fs_ctx_clear(fs_ctx *_ctx) {
96 static void fs_downsample_level(fs_ctx *_ctx, int _l) {
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++) {
119 j1offs = FS_MINI(2 * j + 1, h2) * w2;
120 for (i = 0; i < w; 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];
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) {
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++) {
150 j1 = FS_MINI(j0 + 1, _h);
151 for (i = 0; i < w; 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];
166 static void fs_apply_luminance(fs_ctx *_ctx, int _l) {
167 unsigned *col_sums_x;
168 unsigned *col_sums_y;
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];
196 ssim = _ctx->level[_l].ssim;
197 c1 = (double) (SSIM_C1 * 4096 * (1 << 4 * _l));
198 for (j = 0; j < h; j++) {
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];
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);
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];
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];
235 #define FS_COL_SET(_col, _joffs, _ioffs) \
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; \
247 #define FS_COL_ADD(_col, _joffs, _ioffs) \
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; \
259 #define FS_COL_SUB(_col, _joffs, _ioffs) \
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; \
271 #define FS_COL_COPY(_col1, _col2) \
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)]; \
279 #define FS_COL_HALVE(_col1, _col2) \
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; \
287 #define FS_COL_DOUBLE(_col1, _col2) \
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; \
295 static void fs_calc_structure(fs_ctx *_ctx, int _l) {
301 double col_sums_gx2[8];
302 double col_sums_gy2[8];
303 double col_sums_gxgy[8];
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;
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++) {
322 for (i = 0; i < w - 1; i++) {
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;
337 memset(gx_buf + (j & 7) * stride, 0, stride * sizeof(*gx_buf));
338 memset(gy_buf + (j & 7) * stride, 0, stride * sizeof(*gy_buf));
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);
349 for (k = 1; k < 8 - i; k++) {
351 FS_COL_ADD(i, -k - 1, 0);
355 for (i = 0; i < w; i++) {
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);
370 FS_COL_SET(0, -1, 1);
372 FS_COL_SUB(2, -3, 2);
375 FS_COL_SUB(3, -4, 3);
380 FS_COL_ADD(4, -4, 5);
383 FS_COL_ADD(5, -3, 6);
386 FS_COL_ADD(6, -2, 7);
388 FS_COL_SET(7, -1, 8);
396 #define FS_NLEVELS (4)
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.*/
402 static const double FS_WEIGHTS[FS_NLEVELS] = {0.2989654541015625,
403 0.3141326904296875, 0.2473602294921875, 0.1395416259765625};
405 static double fs_average(fs_ctx *_ctx, int _l) {
412 w = _ctx->level[_l].w;
413 h = _ctx->level[_l].h;
414 ssim = _ctx->level[_l].ssim;
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]);
422 static double calc_ssim(const unsigned char *_src, int _systride,
423 const unsigned char *_dst, int _dystride, int _w, int _h) {
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);
435 fs_calc_structure(&ctx, l);
436 fs_apply_luminance(&ctx, l);
437 ret *= fs_average(&ctx, l);
442 static double convert_ssim_db(double _ssim, double _weight) {
443 return 10 * (log10(_weight) - log10(_weight - _ssim));
446 double vp9_calc_fastssim(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
447 double *ssim_y, double *ssim_u, double *ssim_v) {
449 vp9_clear_system_state();
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);
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);
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));
464 return convert_ssim_db(ssimv, 1.0);