OLD | NEW |
| (Empty) |
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 <math.h> | |
14 #include <string.h> | |
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; | |
21 | |
22 #define SSIM_C1 (255 * 255 * 0.01 * 0.01) | |
23 #define SSIM_C2 (255 * 255 * 0.03 * 0.03) | |
24 | |
25 #define FS_MINI(_a, _b) ((_a) < (_b) ? (_a) : (_b)) | |
26 #define FS_MAXI(_a, _b) ((_a) > (_b) ? (_a) : (_b)) | |
27 | |
28 struct fs_level { | |
29 uint16_t *im1; | |
30 uint16_t *im2; | |
31 double *ssim; | |
32 int w; | |
33 int h; | |
34 }; | |
35 | |
36 struct fs_ctx { | |
37 fs_level *level; | |
38 int nlevels; | |
39 unsigned *col_buf; | |
40 }; | |
41 | |
42 static void fs_ctx_init(fs_ctx *_ctx, int _w, int _h, int _nlevels) { | |
43 unsigned char *data; | |
44 size_t data_size; | |
45 int lw; | |
46 int lh; | |
47 int l; | |
48 lw = (_w + 1) >> 1; | |
49 lh = (_h + 1) >> 1; | |
50 data_size = _nlevels * sizeof(fs_level) | |
51 + 2 * (lw + 8) * 8 * sizeof(*_ctx->col_buf); | |
52 for (l = 0; l < _nlevels; l++) { | |
53 size_t im_size; | |
54 size_t level_size; | |
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; | |
62 lw = (lw + 1) >> 1; | |
63 lh = (lh + 1) >> 1; | |
64 } | |
65 data = (unsigned char *) malloc(data_size); | |
66 _ctx->level = (fs_level *) data; | |
67 _ctx->nlevels = _nlevels; | |
68 data += _nlevels * sizeof(*_ctx->level); | |
69 lw = (_w + 1) >> 1; | |
70 lh = (_h + 1) >> 1; | |
71 for (l = 0; l < _nlevels; l++) { | |
72 size_t im_size; | |
73 size_t level_size; | |
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; | |
83 data += level_size; | |
84 _ctx->level[l].ssim = (double *) data; | |
85 data += im_size * sizeof(*_ctx->level[l].ssim); | |
86 lw = (lw + 1) >> 1; | |
87 lh = (lh + 1) >> 1; | |
88 } | |
89 _ctx->col_buf = (unsigned *) data; | |
90 } | |
91 | |
92 static void fs_ctx_clear(fs_ctx *_ctx) { | |
93 free(_ctx->level); | |
94 } | |
95 | |
96 static void fs_downsample_level(fs_ctx *_ctx, int _l) { | |
97 const uint16_t *src1; | |
98 const uint16_t *src2; | |
99 uint16_t *dst1; | |
100 uint16_t *dst2; | |
101 int w2; | |
102 int h2; | |
103 int w; | |
104 int h; | |
105 int i; | |
106 int j; | |
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++) { | |
116 int j0offs; | |
117 int j1offs; | |
118 j0offs = 2 * j * w2; | |
119 j1offs = FS_MINI(2 * j + 1, h2) * w2; | |
120 for (i = 0; i < w; i++) { | |
121 int i0; | |
122 int i1; | |
123 i0 = 2 * 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]; | |
129 } | |
130 } | |
131 } | |
132 | |
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) { | |
136 uint16_t *dst1; | |
137 uint16_t *dst2; | |
138 int w; | |
139 int h; | |
140 int i; | |
141 int j; | |
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++) { | |
147 int j0; | |
148 int j1; | |
149 j0 = 2 * j; | |
150 j1 = FS_MINI(j0 + 1, _h); | |
151 for (i = 0; i < w; i++) { | |
152 int i0; | |
153 int i1; | |
154 i0 = 2 * 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]; | |
162 } | |
163 } | |
164 } | |
165 | |
166 static void fs_apply_luminance(fs_ctx *_ctx, int _l) { | |
167 unsigned *col_sums_x; | |
168 unsigned *col_sums_y; | |
169 uint16_t *im1; | |
170 uint16_t *im2; | |
171 double *ssim; | |
172 double c1; | |
173 int w; | |
174 int h; | |
175 int j0offs; | |
176 int j1offs; | |
177 int i; | |
178 int j; | |
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]; | |
195 } | |
196 ssim = _ctx->level[_l].ssim; | |
197 c1 = (double) (SSIM_C1 * 4096 * (1 << 4 * _l)); | |
198 for (j = 0; j < h; j++) { | |
199 unsigned mux; | |
200 unsigned muy; | |
201 int i0; | |
202 int i1; | |
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]; | |
209 } | |
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); | |
213 if (i + 1 < w) { | |
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]; | |
218 } | |
219 } | |
220 if (j + 1 < h) { | |
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]; | |
231 } | |
232 } | |
233 } | |
234 | |
235 #define FS_COL_SET(_col, _joffs, _ioffs) \ | |
236 do { \ | |
237 unsigned gx; \ | |
238 unsigned gy; \ | |
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; \ | |
244 } \ | |
245 while (0) | |
246 | |
247 #define FS_COL_ADD(_col, _joffs, _ioffs) \ | |
248 do { \ | |
249 unsigned gx; \ | |
250 unsigned gy; \ | |
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; \ | |
256 } \ | |
257 while (0) | |
258 | |
259 #define FS_COL_SUB(_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 } \ | |
269 while (0) | |
270 | |
271 #define FS_COL_COPY(_col1, _col2) \ | |
272 do { \ | |
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)]; \ | |
276 } \ | |
277 while (0) | |
278 | |
279 #define FS_COL_HALVE(_col1, _col2) \ | |
280 do { \ | |
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; \ | |
284 } \ | |
285 while (0) | |
286 | |
287 #define FS_COL_DOUBLE(_col1, _col2) \ | |
288 do { \ | |
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; \ | |
292 } \ | |
293 while (0) | |
294 | |
295 static void fs_calc_structure(fs_ctx *_ctx, int _l) { | |
296 uint16_t *im1; | |
297 uint16_t *im2; | |
298 unsigned *gx_buf; | |
299 unsigned *gy_buf; | |
300 double *ssim; | |
301 double col_sums_gx2[8]; | |
302 double col_sums_gy2[8]; | |
303 double col_sums_gxgy[8]; | |
304 double c2; | |
305 int stride; | |
306 int w; | |
307 int h; | |
308 int i; | |
309 int j; | |
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; | |
316 stride = w + 8; | |
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++) { | |
321 if (j < h - 1) { | |
322 for (i = 0; i < w - 1; i++) { | |
323 unsigned g1; | |
324 unsigned g2; | |
325 unsigned gx; | |
326 unsigned gy; | |
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; | |
335 } | |
336 } else { | |
337 memset(gx_buf + (j & 7) * stride, 0, stride * sizeof(*gx_buf)); | |
338 memset(gy_buf + (j & 7) * stride, 0, stride * sizeof(*gy_buf)); | |
339 } | |
340 if (j >= 4) { | |
341 int k; | |
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); | |
348 FS_COL_ADD(i, 0, 0); | |
349 for (k = 1; k < 8 - i; k++) { | |
350 FS_COL_DOUBLE(i, i); | |
351 FS_COL_ADD(i, -k - 1, 0); | |
352 FS_COL_ADD(i, k, 0); | |
353 } | |
354 } | |
355 for (i = 0; i < w; i++) { | |
356 double mugx2; | |
357 double mugy2; | |
358 double mugxgy; | |
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); | |
369 if (i + 1 < w) { | |
370 FS_COL_SET(0, -1, 1); | |
371 FS_COL_ADD(0, 0, 1); | |
372 FS_COL_SUB(2, -3, 2); | |
373 FS_COL_SUB(2, 2, 2); | |
374 FS_COL_HALVE(1, 2); | |
375 FS_COL_SUB(3, -4, 3); | |
376 FS_COL_SUB(3, 3, 3); | |
377 FS_COL_HALVE(2, 3); | |
378 FS_COL_COPY(3, 4); | |
379 FS_COL_DOUBLE(4, 5); | |
380 FS_COL_ADD(4, -4, 5); | |
381 FS_COL_ADD(4, 3, 5); | |
382 FS_COL_DOUBLE(5, 6); | |
383 FS_COL_ADD(5, -3, 6); | |
384 FS_COL_ADD(5, 2, 6); | |
385 FS_COL_DOUBLE(6, 7); | |
386 FS_COL_ADD(6, -2, 7); | |
387 FS_COL_ADD(6, 1, 7); | |
388 FS_COL_SET(7, -1, 8); | |
389 FS_COL_ADD(7, 0, 8); | |
390 } | |
391 } | |
392 } | |
393 } | |
394 } | |
395 | |
396 #define FS_NLEVELS (4) | |
397 | |
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.*/ | |
401 | |
402 static const double FS_WEIGHTS[FS_NLEVELS] = {0.2989654541015625, | |
403 0.3141326904296875, 0.2473602294921875, 0.1395416259765625}; | |
404 | |
405 static double fs_average(fs_ctx *_ctx, int _l) { | |
406 double *ssim; | |
407 double ret; | |
408 int w; | |
409 int h; | |
410 int i; | |
411 int j; | |
412 w = _ctx->level[_l].w; | |
413 h = _ctx->level[_l].h; | |
414 ssim = _ctx->level[_l].ssim; | |
415 ret = 0; | |
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]); | |
420 } | |
421 | |
422 static double calc_ssim(const unsigned char *_src, int _systride, | |
423 const unsigned char *_dst, int _dystride, int _w, int _h) { | |
424 fs_ctx ctx; | |
425 double ret; | |
426 int l; | |
427 ret = 1; | |
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); | |
434 } | |
435 fs_calc_structure(&ctx, l); | |
436 fs_apply_luminance(&ctx, l); | |
437 ret *= fs_average(&ctx, l); | |
438 fs_ctx_clear(&ctx); | |
439 return ret; | |
440 } | |
441 | |
442 static double convert_ssim_db(double _ssim, double _weight) { | |
443 return 10 * (log10(_weight) - log10(_weight - _ssim)); | |
444 } | |
445 | |
446 double vp9_calc_fastssim(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest, | |
447 double *ssim_y, double *ssim_u, double *ssim_v) { | |
448 double ssimv; | |
449 vp9_clear_system_state(); | |
450 | |
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); | |
454 | |
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); | |
458 | |
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)); | |
463 | |
464 return convert_ssim_db(ssimv, 1.0); | |
465 } | |
OLD | NEW |