blob: 0c21273714563cbb08ee08bdfe541f7907963b05 [file] [log] [blame]
Yaowu Xuc27fc142016-08-22 16:08:15 -07001/*
Krishna Rapaka7319db52021-09-28 20:35:29 -07002 * Copyright (c) 2021, Alliance for Open Media. All rights reserved
Yaowu Xuc27fc142016-08-22 16:08:15 -07003 *
Vibhoothi41c6dd72021-10-12 18:48:26 +00004 * This source code is subject to the terms of the BSD 3-Clause Clear License
5 * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear
6 * License was not distributed with this source code in the LICENSE file, you
7 * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. If the
8 * Alliance for Open Media Patent License 1.0 was not distributed with this
9 * source code in the PATENTS file, you can obtain it at
10 * aomedia.org/license/patent-license/.
Yaowu Xuc27fc142016-08-22 16:08:15 -070011 *
12 * This code was originally written by: Nathan E. Egge, at the Daala
13 * project.
14 */
15#include <assert.h>
16#include <math.h>
17#include <stdlib.h>
18#include <string.h>
Tom Finegan60e653d2018-05-22 11:34:58 -070019
20#include "config/aom_config.h"
Tom Finegan44702c82018-05-22 13:00:39 -070021#include "config/aom_dsp_rtcd.h"
Tom Finegan60e653d2018-05-22 11:34:58 -070022
Yaowu Xuc27fc142016-08-22 16:08:15 -070023#include "aom_dsp/ssim.h"
24#include "aom_ports/system_state.h"
25
26typedef struct fs_level fs_level;
27typedef struct fs_ctx fs_ctx;
28
29#define SSIM_C1 (255 * 255 * 0.01 * 0.01)
30#define SSIM_C2 (255 * 255 * 0.03 * 0.03)
Yaowu Xuc27fc142016-08-22 16:08:15 -070031#define SSIM_C1_10 (1023 * 1023 * 0.01 * 0.01)
32#define SSIM_C1_12 (4095 * 4095 * 0.01 * 0.01)
33#define SSIM_C2_10 (1023 * 1023 * 0.03 * 0.03)
34#define SSIM_C2_12 (4095 * 4095 * 0.03 * 0.03)
Yaowu Xud3e7c682017-12-21 14:08:25 -080035
Yaowu Xuc27fc142016-08-22 16:08:15 -070036#define FS_MINI(_a, _b) ((_a) < (_b) ? (_a) : (_b))
37#define FS_MAXI(_a, _b) ((_a) > (_b) ? (_a) : (_b))
38
39struct fs_level {
40 uint32_t *im1;
41 uint32_t *im2;
42 double *ssim;
43 int w;
44 int h;
45};
46
47struct fs_ctx {
48 fs_level *level;
49 int nlevels;
50 unsigned *col_buf;
51};
52
53static void fs_ctx_init(fs_ctx *_ctx, int _w, int _h, int _nlevels) {
54 unsigned char *data;
55 size_t data_size;
56 int lw;
57 int lh;
58 int l;
59 lw = (_w + 1) >> 1;
60 lh = (_h + 1) >> 1;
61 data_size =
62 _nlevels * sizeof(fs_level) + 2 * (lw + 8) * 8 * sizeof(*_ctx->col_buf);
63 for (l = 0; l < _nlevels; l++) {
64 size_t im_size;
65 size_t level_size;
66 im_size = lw * (size_t)lh;
67 level_size = 2 * im_size * sizeof(*_ctx->level[l].im1);
68 level_size += sizeof(*_ctx->level[l].ssim) - 1;
69 level_size /= sizeof(*_ctx->level[l].ssim);
70 level_size += im_size;
71 level_size *= sizeof(*_ctx->level[l].ssim);
72 data_size += level_size;
73 lw = (lw + 1) >> 1;
74 lh = (lh + 1) >> 1;
75 }
76 data = (unsigned char *)malloc(data_size);
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}
102
103static void fs_ctx_clear(fs_ctx *_ctx) { free(_ctx->level); }
104
105static void fs_downsample_level(fs_ctx *_ctx, int _l) {
106 const uint32_t *src1;
107 const uint32_t *src2;
108 uint32_t *dst1;
109 uint32_t *dst2;
110 int w2;
111 int h2;
112 int w;
113 int h;
114 int i;
115 int j;
116 w = _ctx->level[_l].w;
117 h = _ctx->level[_l].h;
118 dst1 = _ctx->level[_l].im1;
119 dst2 = _ctx->level[_l].im2;
120 w2 = _ctx->level[_l - 1].w;
121 h2 = _ctx->level[_l - 1].h;
122 src1 = _ctx->level[_l - 1].im1;
123 src2 = _ctx->level[_l - 1].im2;
124 for (j = 0; j < h; j++) {
125 int j0offs;
126 int j1offs;
127 j0offs = 2 * j * w2;
128 j1offs = FS_MINI(2 * j + 1, h2) * w2;
129 for (i = 0; i < w; i++) {
130 int i0;
131 int i1;
132 i0 = 2 * i;
133 i1 = FS_MINI(i0 + 1, w2);
134 dst1[j * w + i] = src1[j0offs + i0] + src1[j0offs + i1] +
135 src1[j1offs + i0] + src1[j1offs + i1];
136 dst2[j * w + i] = src2[j0offs + i0] + src2[j0offs + i1] +
137 src2[j1offs + i0] + src2[j1offs + i1];
138 }
139 }
140}
141
James Almer756717c2022-10-19 17:32:48 +0000142static void fs_downsample_level0(fs_ctx *_ctx, const uint16_t *src1s,
143 int _s1ystride, const uint16_t *src2s,
James Almer857e93f2022-05-25 16:44:43 +0000144 int _s2ystride, int _w, int _h,
145 uint32_t shift) {
Yaowu Xuc27fc142016-08-22 16:08:15 -0700146 uint32_t *dst1;
147 uint32_t *dst2;
148 int w;
149 int h;
150 int i;
151 int j;
152 w = _ctx->level[0].w;
153 h = _ctx->level[0].h;
154 dst1 = _ctx->level[0].im1;
155 dst2 = _ctx->level[0].im2;
156 for (j = 0; j < h; j++) {
157 int j0;
158 int j1;
159 j0 = 2 * j;
160 j1 = FS_MINI(j0 + 1, _h);
161 for (i = 0; i < w; i++) {
162 int i0;
163 int i1;
164 i0 = 2 * i;
165 i1 = FS_MINI(i0 + 1, _w);
James Almer857e93f2022-05-25 16:44:43 +0000166 dst1[j * w + i] = (src1s[j0 * _s1ystride + i0] >> shift) +
167 (src1s[j0 * _s1ystride + i1] >> shift) +
168 (src1s[j1 * _s1ystride + i0] >> shift) +
169 (src1s[j1 * _s1ystride + i1] >> shift);
170 dst2[j * w + i] = (src2s[j0 * _s2ystride + i0] >> shift) +
171 (src2s[j0 * _s2ystride + i1] >> shift) +
172 (src2s[j1 * _s2ystride + i0] >> shift) +
173 (src2s[j1 * _s2ystride + i1] >> shift);
Yaowu Xuc27fc142016-08-22 16:08:15 -0700174 }
175 }
176}
177
178static void fs_apply_luminance(fs_ctx *_ctx, int _l, int bit_depth) {
179 unsigned *col_sums_x;
180 unsigned *col_sums_y;
181 uint32_t *im1;
182 uint32_t *im2;
183 double *ssim;
184 double c1;
185 int w;
186 int h;
187 int j0offs;
188 int j1offs;
189 int i;
190 int j;
191 double ssim_c1 = SSIM_C1;
Yaowu Xud3e7c682017-12-21 14:08:25 -0800192
Yaowu Xuc27fc142016-08-22 16:08:15 -0700193 if (bit_depth == 10) ssim_c1 = SSIM_C1_10;
194 if (bit_depth == 12) ssim_c1 = SSIM_C1_12;
Yaowu Xud3e7c682017-12-21 14:08:25 -0800195
Yaowu Xuc27fc142016-08-22 16:08:15 -0700196 w = _ctx->level[_l].w;
197 h = _ctx->level[_l].h;
198 col_sums_x = _ctx->col_buf;
199 col_sums_y = col_sums_x + w;
200 im1 = _ctx->level[_l].im1;
201 im2 = _ctx->level[_l].im2;
202 for (i = 0; i < w; i++) col_sums_x[i] = 5 * im1[i];
203 for (i = 0; i < w; i++) col_sums_y[i] = 5 * im2[i];
204 for (j = 1; j < 4; j++) {
205 j1offs = FS_MINI(j, h - 1) * w;
206 for (i = 0; i < w; i++) col_sums_x[i] += im1[j1offs + i];
207 for (i = 0; i < w; i++) col_sums_y[i] += im2[j1offs + i];
208 }
209 ssim = _ctx->level[_l].ssim;
210 c1 = (double)(ssim_c1 * 4096 * (1 << 4 * _l));
211 for (j = 0; j < h; j++) {
212 unsigned mux;
213 unsigned muy;
214 int i0;
215 int i1;
216 mux = 5 * col_sums_x[0];
217 muy = 5 * col_sums_y[0];
218 for (i = 1; i < 4; i++) {
219 i1 = FS_MINI(i, w - 1);
220 mux += col_sums_x[i1];
221 muy += col_sums_y[i1];
222 }
223 for (i = 0; i < w; i++) {
224 ssim[j * w + i] *= (2 * mux * (double)muy + c1) /
225 (mux * (double)mux + muy * (double)muy + c1);
226 if (i + 1 < w) {
227 i0 = FS_MAXI(0, i - 4);
228 i1 = FS_MINI(i + 4, w - 1);
229 mux += col_sums_x[i1] - col_sums_x[i0];
230 muy += col_sums_x[i1] - col_sums_x[i0];
231 }
232 }
233 if (j + 1 < h) {
234 j0offs = FS_MAXI(0, j - 4) * w;
235 for (i = 0; i < w; i++) col_sums_x[i] -= im1[j0offs + i];
236 for (i = 0; i < w; i++) col_sums_y[i] -= im2[j0offs + i];
237 j1offs = FS_MINI(j + 4, h - 1) * w;
238 for (i = 0; i < w; i++) col_sums_x[i] += im1[j1offs + i];
239 for (i = 0; i < w; i++) col_sums_y[i] += im2[j1offs + i];
240 }
241 }
242}
243
244#define FS_COL_SET(_col, _joffs, _ioffs) \
245 do { \
246 unsigned gx; \
247 unsigned gy; \
248 gx = gx_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
249 gy = gy_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
250 col_sums_gx2[(_col)] = gx * (double)gx; \
251 col_sums_gy2[(_col)] = gy * (double)gy; \
252 col_sums_gxgy[(_col)] = gx * (double)gy; \
253 } while (0)
254
255#define FS_COL_ADD(_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_SUB(_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_COPY(_col1, _col2) \
278 do { \
279 col_sums_gx2[(_col1)] = col_sums_gx2[(_col2)]; \
280 col_sums_gy2[(_col1)] = col_sums_gy2[(_col2)]; \
281 col_sums_gxgy[(_col1)] = col_sums_gxgy[(_col2)]; \
282 } while (0)
283
284#define FS_COL_HALVE(_col1, _col2) \
285 do { \
286 col_sums_gx2[(_col1)] = col_sums_gx2[(_col2)] * 0.5; \
287 col_sums_gy2[(_col1)] = col_sums_gy2[(_col2)] * 0.5; \
288 col_sums_gxgy[(_col1)] = col_sums_gxgy[(_col2)] * 0.5; \
289 } while (0)
290
291#define FS_COL_DOUBLE(_col1, _col2) \
292 do { \
293 col_sums_gx2[(_col1)] = col_sums_gx2[(_col2)] * 2; \
294 col_sums_gy2[(_col1)] = col_sums_gy2[(_col2)] * 2; \
295 col_sums_gxgy[(_col1)] = col_sums_gxgy[(_col2)] * 2; \
296 } while (0)
297
298static void fs_calc_structure(fs_ctx *_ctx, int _l, int bit_depth) {
299 uint32_t *im1;
300 uint32_t *im2;
301 unsigned *gx_buf;
302 unsigned *gy_buf;
303 double *ssim;
304 double col_sums_gx2[8];
305 double col_sums_gy2[8];
306 double col_sums_gxgy[8];
307 double c2;
308 int stride;
309 int w;
310 int h;
311 int i;
312 int j;
313 double ssim_c2 = SSIM_C2;
Yaowu Xuc27fc142016-08-22 16:08:15 -0700314 if (bit_depth == 10) ssim_c2 = SSIM_C2_10;
315 if (bit_depth == 12) ssim_c2 = SSIM_C2_12;
Yaowu Xuc27fc142016-08-22 16:08:15 -0700316
317 w = _ctx->level[_l].w;
318 h = _ctx->level[_l].h;
319 im1 = _ctx->level[_l].im1;
320 im2 = _ctx->level[_l].im2;
321 ssim = _ctx->level[_l].ssim;
322 gx_buf = _ctx->col_buf;
323 stride = w + 8;
324 gy_buf = gx_buf + 8 * stride;
325 memset(gx_buf, 0, 2 * 8 * stride * sizeof(*gx_buf));
326 c2 = ssim_c2 * (1 << 4 * _l) * 16 * 104;
327 for (j = 0; j < h + 4; j++) {
328 if (j < h - 1) {
329 for (i = 0; i < w - 1; i++) {
330 unsigned g1;
331 unsigned g2;
332 unsigned gx;
333 unsigned gy;
334 g1 = abs((int)im1[(j + 1) * w + i + 1] - (int)im1[j * w + i]);
335 g2 = abs((int)im1[(j + 1) * w + i] - (int)im1[j * w + i + 1]);
336 gx = 4 * FS_MAXI(g1, g2) + FS_MINI(g1, g2);
337 g1 = abs((int)im2[(j + 1) * w + i + 1] - (int)im2[j * w + i]);
338 g2 = abs((int)im2[(j + 1) * w + i] - (int)im2[j * w + i + 1]);
339 gy = 4 * FS_MAXI(g1, g2) + FS_MINI(g1, g2);
340 gx_buf[(j & 7) * stride + i + 4] = gx;
341 gy_buf[(j & 7) * stride + i + 4] = gy;
342 }
343 } else {
344 memset(gx_buf + (j & 7) * stride, 0, stride * sizeof(*gx_buf));
345 memset(gy_buf + (j & 7) * stride, 0, stride * sizeof(*gy_buf));
346 }
347 if (j >= 4) {
348 int k;
349 col_sums_gx2[3] = col_sums_gx2[2] = col_sums_gx2[1] = col_sums_gx2[0] = 0;
350 col_sums_gy2[3] = col_sums_gy2[2] = col_sums_gy2[1] = col_sums_gy2[0] = 0;
351 col_sums_gxgy[3] = col_sums_gxgy[2] = col_sums_gxgy[1] =
352 col_sums_gxgy[0] = 0;
353 for (i = 4; i < 8; i++) {
354 FS_COL_SET(i, -1, 0);
355 FS_COL_ADD(i, 0, 0);
356 for (k = 1; k < 8 - i; k++) {
357 FS_COL_DOUBLE(i, i);
358 FS_COL_ADD(i, -k - 1, 0);
359 FS_COL_ADD(i, k, 0);
360 }
361 }
362 for (i = 0; i < w; i++) {
363 double mugx2;
364 double mugy2;
365 double mugxgy;
366 mugx2 = col_sums_gx2[0];
367 for (k = 1; k < 8; k++) mugx2 += col_sums_gx2[k];
368 mugy2 = col_sums_gy2[0];
369 for (k = 1; k < 8; k++) mugy2 += col_sums_gy2[k];
370 mugxgy = col_sums_gxgy[0];
371 for (k = 1; k < 8; k++) mugxgy += col_sums_gxgy[k];
372 ssim[(j - 4) * w + i] = (2 * mugxgy + c2) / (mugx2 + mugy2 + c2);
373 if (i + 1 < w) {
374 FS_COL_SET(0, -1, 1);
375 FS_COL_ADD(0, 0, 1);
376 FS_COL_SUB(2, -3, 2);
377 FS_COL_SUB(2, 2, 2);
378 FS_COL_HALVE(1, 2);
379 FS_COL_SUB(3, -4, 3);
380 FS_COL_SUB(3, 3, 3);
381 FS_COL_HALVE(2, 3);
382 FS_COL_COPY(3, 4);
383 FS_COL_DOUBLE(4, 5);
384 FS_COL_ADD(4, -4, 5);
385 FS_COL_ADD(4, 3, 5);
386 FS_COL_DOUBLE(5, 6);
387 FS_COL_ADD(5, -3, 6);
388 FS_COL_ADD(5, 2, 6);
389 FS_COL_DOUBLE(6, 7);
390 FS_COL_ADD(6, -2, 7);
391 FS_COL_ADD(6, 1, 7);
392 FS_COL_SET(7, -1, 8);
393 FS_COL_ADD(7, 0, 8);
394 }
395 }
396 }
397 }
398}
399
400#define FS_NLEVELS (4)
401
402/*These weights were derived from the default weights found in Wang's original
403 Matlab implementation: {0.0448, 0.2856, 0.2363, 0.1333}.
404 We drop the finest scale and renormalize the rest to sum to 1.*/
405
406static const double FS_WEIGHTS[FS_NLEVELS] = {
407 0.2989654541015625, 0.3141326904296875, 0.2473602294921875, 0.1395416259765625
408};
409
410static double fs_average(fs_ctx *_ctx, int _l) {
411 double *ssim;
412 double ret;
413 int w;
414 int h;
415 int i;
416 int j;
417 w = _ctx->level[_l].w;
418 h = _ctx->level[_l].h;
419 ssim = _ctx->level[_l].ssim;
420 ret = 0;
421 for (j = 0; j < h; j++)
422 for (i = 0; i < w; i++) ret += ssim[j * w + i];
423 return pow(ret / (w * h), FS_WEIGHTS[_l]);
424}
425
426static double convert_ssim_db(double _ssim, double _weight) {
427 assert(_weight >= _ssim);
428 if ((_weight - _ssim) < 1e-10) return MAX_SSIM_DB;
429 return 10 * (log10(_weight) - log10(_weight - _ssim));
430}
431
James Almer756717c2022-10-19 17:32:48 +0000432static double calc_ssim(const uint16_t *_src, int _systride,
433 const uint16_t *_dst, int _dystride, int _w, int _h,
434 uint32_t _bd, uint32_t _shift) {
Yaowu Xuc27fc142016-08-22 16:08:15 -0700435 fs_ctx ctx;
436 double ret;
437 int l;
438 ret = 1;
439 fs_ctx_init(&ctx, _w, _h, FS_NLEVELS);
James Almer857e93f2022-05-25 16:44:43 +0000440 fs_downsample_level0(&ctx, _src, _systride, _dst, _dystride, _w, _h, _shift);
Yaowu Xuc27fc142016-08-22 16:08:15 -0700441 for (l = 0; l < FS_NLEVELS - 1; l++) {
442 fs_calc_structure(&ctx, l, _bd);
443 ret *= fs_average(&ctx, l);
444 fs_downsample_level(&ctx, l + 1);
445 }
446 fs_calc_structure(&ctx, l, _bd);
447 fs_apply_luminance(&ctx, l, _bd);
448 ret *= fs_average(&ctx, l);
449 fs_ctx_clear(&ctx);
450 return ret;
451}
452
Yaowu Xuf883b422016-08-30 14:01:10 -0700453double aom_calc_fastssim(const YV12_BUFFER_CONFIG *source,
Yaowu Xuc27fc142016-08-22 16:08:15 -0700454 const YV12_BUFFER_CONFIG *dest, double *ssim_y,
455 double *ssim_u, double *ssim_v, uint32_t bd,
456 uint32_t in_bd) {
457 double ssimv;
458 uint32_t bd_shift = 0;
Yaowu Xuf883b422016-08-30 14:01:10 -0700459 aom_clear_system_state();
Yaowu Xuc27fc142016-08-22 16:08:15 -0700460 assert(bd >= in_bd);
Yaowu Xu0e1df782018-03-12 12:18:51 -0700461 assert(source->flags == dest->flags);
Yaowu Xuc27fc142016-08-22 16:08:15 -0700462 bd_shift = bd - in_bd;
463
464 *ssim_y = calc_ssim(source->y_buffer, source->y_stride, dest->y_buffer,
465 dest->y_stride, source->y_crop_width,
James Almer857e93f2022-05-25 16:44:43 +0000466 source->y_crop_height, in_bd, bd_shift);
Yaowu Xuc27fc142016-08-22 16:08:15 -0700467 *ssim_u = calc_ssim(source->u_buffer, source->uv_stride, dest->u_buffer,
468 dest->uv_stride, source->uv_crop_width,
James Almer857e93f2022-05-25 16:44:43 +0000469 source->uv_crop_height, in_bd, bd_shift);
Yaowu Xuc27fc142016-08-22 16:08:15 -0700470 *ssim_v = calc_ssim(source->v_buffer, source->uv_stride, dest->v_buffer,
471 dest->uv_stride, source->uv_crop_width,
James Almer857e93f2022-05-25 16:44:43 +0000472 source->uv_crop_height, in_bd, bd_shift);
Yaowu Xuc27fc142016-08-22 16:08:15 -0700473 ssimv = (*ssim_y) * .8 + .1 * ((*ssim_u) + (*ssim_v));
474 return convert_ssim_db(ssimv, 1.0);
475}