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
fs_ctx_init(fs_ctx * _ctx,int _w,int _h,int _nlevels)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
fs_ctx_clear(fs_ctx * _ctx)104 static void fs_ctx_clear(fs_ctx *_ctx) { free(_ctx->level); }
105
fs_downsample_level(fs_ctx * _ctx,int _l)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
fs_downsample_level0(fs_ctx * _ctx,const uint8_t * _src1,int _s1ystride,const uint8_t * _src2,int _s2ystride,int _w,int _h,uint32_t shift,int buf_is_hbd)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
fs_apply_luminance(fs_ctx * _ctx,int _l,int bit_depth)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
fs_calc_structure(fs_ctx * _ctx,int _l,int bit_depth)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
fs_average(fs_ctx * _ctx,int _l)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
convert_ssim_db(double _ssim,double _weight)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
calc_ssim(const uint8_t * _src,int _systride,const uint8_t * _dst,int _dystride,int _w,int _h,uint32_t _bd,uint32_t _shift,int buf_is_hbd)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
aom_calc_fastssim(const YV12_BUFFER_CONFIG * source,const YV12_BUFFER_CONFIG * dest,double * ssim_y,double * ssim_u,double * ssim_v,uint32_t bd,uint32_t in_bd)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 }
489