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);