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