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 | |
11 #include <math.h> | |
12 #include "./vp9_rtcd.h" | |
13 #include "vpx_ports/mem.h" | |
14 #include "vp9/encoder/vp9_ssim.h" | |
15 | |
16 void vp9_ssim_parms_16x16_c(uint8_t *s, int sp, uint8_t *r, | |
17 int rp, unsigned long *sum_s, unsigned long *sum_r, | |
18 unsigned long *sum_sq_s, unsigned long *sum_sq_r, | |
19 unsigned long *sum_sxr) { | |
20 int i, j; | |
21 for (i = 0; i < 16; i++, s += sp, r += rp) { | |
22 for (j = 0; j < 16; j++) { | |
23 *sum_s += s[j]; | |
24 *sum_r += r[j]; | |
25 *sum_sq_s += s[j] * s[j]; | |
26 *sum_sq_r += r[j] * r[j]; | |
27 *sum_sxr += s[j] * r[j]; | |
28 } | |
29 } | |
30 } | |
31 void vp9_ssim_parms_8x8_c(uint8_t *s, int sp, uint8_t *r, int rp, | |
32 unsigned long *sum_s, unsigned long *sum_r, | |
33 unsigned long *sum_sq_s, unsigned long *sum_sq_r, | |
34 unsigned long *sum_sxr) { | |
35 int i, j; | |
36 for (i = 0; i < 8; i++, s += sp, r += rp) { | |
37 for (j = 0; j < 8; j++) { | |
38 *sum_s += s[j]; | |
39 *sum_r += r[j]; | |
40 *sum_sq_s += s[j] * s[j]; | |
41 *sum_sq_r += r[j] * r[j]; | |
42 *sum_sxr += s[j] * r[j]; | |
43 } | |
44 } | |
45 } | |
46 | |
47 #if CONFIG_VP9_HIGHBITDEPTH | |
48 void vp9_highbd_ssim_parms_8x8_c(uint16_t *s, int sp, uint16_t *r, int rp, | |
49 uint32_t *sum_s, uint32_t *sum_r, | |
50 uint32_t *sum_sq_s, uint32_t *sum_sq_r, | |
51 uint32_t *sum_sxr) { | |
52 int i, j; | |
53 for (i = 0; i < 8; i++, s += sp, r += rp) { | |
54 for (j = 0; j < 8; j++) { | |
55 *sum_s += s[j]; | |
56 *sum_r += r[j]; | |
57 *sum_sq_s += s[j] * s[j]; | |
58 *sum_sq_r += r[j] * r[j]; | |
59 *sum_sxr += s[j] * r[j]; | |
60 } | |
61 } | |
62 } | |
63 #endif // CONFIG_VP9_HIGHBITDEPTH | |
64 | |
65 static const int64_t cc1 = 26634; // (64^2*(.01*255)^2 | |
66 static const int64_t cc2 = 239708; // (64^2*(.03*255)^2 | |
67 | |
68 static double similarity(unsigned long sum_s, unsigned long sum_r, | |
69 unsigned long sum_sq_s, unsigned long sum_sq_r, | |
70 unsigned long sum_sxr, int count) { | |
71 int64_t ssim_n, ssim_d; | |
72 int64_t c1, c2; | |
73 | |
74 // scale the constants by number of pixels | |
75 c1 = (cc1 * count * count) >> 12; | |
76 c2 = (cc2 * count * count) >> 12; | |
77 | |
78 ssim_n = (2 * sum_s * sum_r + c1) * ((int64_t) 2 * count * sum_sxr - | |
79 (int64_t) 2 * sum_s * sum_r + c2); | |
80 | |
81 ssim_d = (sum_s * sum_s + sum_r * sum_r + c1) * | |
82 ((int64_t)count * sum_sq_s - (int64_t)sum_s * sum_s + | |
83 (int64_t)count * sum_sq_r - (int64_t) sum_r * sum_r + c2); | |
84 | |
85 return ssim_n * 1.0 / ssim_d; | |
86 } | |
87 | |
88 static double ssim_8x8(uint8_t *s, int sp, uint8_t *r, int rp) { | |
89 unsigned long sum_s = 0, sum_r = 0, sum_sq_s = 0, sum_sq_r = 0, sum_sxr = 0; | |
90 vp9_ssim_parms_8x8(s, sp, r, rp, &sum_s, &sum_r, &sum_sq_s, &sum_sq_r, | |
91 &sum_sxr); | |
92 return similarity(sum_s, sum_r, sum_sq_s, sum_sq_r, sum_sxr, 64); | |
93 } | |
94 | |
95 #if CONFIG_VP9_HIGHBITDEPTH | |
96 static double highbd_ssim_8x8(uint16_t *s, int sp, uint16_t *r, int rp, | |
97 unsigned int bd) { | |
98 uint32_t sum_s = 0, sum_r = 0, sum_sq_s = 0, sum_sq_r = 0, sum_sxr = 0; | |
99 const int oshift = bd - 8; | |
100 vp9_highbd_ssim_parms_8x8(s, sp, r, rp, &sum_s, &sum_r, &sum_sq_s, &sum_sq_r, | |
101 &sum_sxr); | |
102 return similarity(sum_s >> oshift, | |
103 sum_r >> oshift, | |
104 sum_sq_s >> (2 * oshift), | |
105 sum_sq_r >> (2 * oshift), | |
106 sum_sxr >> (2 * oshift), | |
107 64); | |
108 } | |
109 #endif // CONFIG_VP9_HIGHBITDEPTH | |
110 | |
111 // We are using a 8x8 moving window with starting location of each 8x8 window | |
112 // on the 4x4 pixel grid. Such arrangement allows the windows to overlap | |
113 // block boundaries to penalize blocking artifacts. | |
114 double vp9_ssim2(uint8_t *img1, uint8_t *img2, int stride_img1, | |
115 int stride_img2, int width, int height) { | |
116 int i, j; | |
117 int samples = 0; | |
118 double ssim_total = 0; | |
119 | |
120 // sample point start with each 4x4 location | |
121 for (i = 0; i <= height - 8; | |
122 i += 4, img1 += stride_img1 * 4, img2 += stride_img2 * 4) { | |
123 for (j = 0; j <= width - 8; j += 4) { | |
124 double v = ssim_8x8(img1 + j, stride_img1, img2 + j, stride_img2); | |
125 ssim_total += v; | |
126 samples++; | |
127 } | |
128 } | |
129 ssim_total /= samples; | |
130 return ssim_total; | |
131 } | |
132 | |
133 #if CONFIG_VP9_HIGHBITDEPTH | |
134 double vp9_highbd_ssim2(uint8_t *img1, uint8_t *img2, int stride_img1, | |
135 int stride_img2, int width, int height, | |
136 unsigned int bd) { | |
137 int i, j; | |
138 int samples = 0; | |
139 double ssim_total = 0; | |
140 | |
141 // sample point start with each 4x4 location | |
142 for (i = 0; i <= height - 8; | |
143 i += 4, img1 += stride_img1 * 4, img2 += stride_img2 * 4) { | |
144 for (j = 0; j <= width - 8; j += 4) { | |
145 double v = highbd_ssim_8x8(CONVERT_TO_SHORTPTR(img1 + j), stride_img1, | |
146 CONVERT_TO_SHORTPTR(img2 + j), stride_img2, | |
147 bd); | |
148 ssim_total += v; | |
149 samples++; | |
150 } | |
151 } | |
152 ssim_total /= samples; | |
153 return ssim_total; | |
154 } | |
155 #endif // CONFIG_VP9_HIGHBITDEPTH | |
156 | |
157 double vp9_calc_ssim(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest, | |
158 double *weight) { | |
159 double a, b, c; | |
160 double ssimv; | |
161 | |
162 a = vp9_ssim2(source->y_buffer, dest->y_buffer, | |
163 source->y_stride, dest->y_stride, | |
164 source->y_crop_width, source->y_crop_height); | |
165 | |
166 b = vp9_ssim2(source->u_buffer, dest->u_buffer, | |
167 source->uv_stride, dest->uv_stride, | |
168 source->uv_crop_width, source->uv_crop_height); | |
169 | |
170 c = vp9_ssim2(source->v_buffer, dest->v_buffer, | |
171 source->uv_stride, dest->uv_stride, | |
172 source->uv_crop_width, source->uv_crop_height); | |
173 | |
174 ssimv = a * .8 + .1 * (b + c); | |
175 | |
176 *weight = 1; | |
177 | |
178 return ssimv; | |
179 } | |
180 | |
181 double vp9_calc_ssimg(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest, | |
182 double *ssim_y, double *ssim_u, double *ssim_v) { | |
183 double ssim_all = 0; | |
184 double a, b, c; | |
185 | |
186 a = vp9_ssim2(source->y_buffer, dest->y_buffer, | |
187 source->y_stride, dest->y_stride, | |
188 source->y_crop_width, source->y_crop_height); | |
189 | |
190 b = vp9_ssim2(source->u_buffer, dest->u_buffer, | |
191 source->uv_stride, dest->uv_stride, | |
192 source->uv_crop_width, source->uv_crop_height); | |
193 | |
194 c = vp9_ssim2(source->v_buffer, dest->v_buffer, | |
195 source->uv_stride, dest->uv_stride, | |
196 source->uv_crop_width, source->uv_crop_height); | |
197 *ssim_y = a; | |
198 *ssim_u = b; | |
199 *ssim_v = c; | |
200 ssim_all = (a * 4 + b + c) / 6; | |
201 | |
202 return ssim_all; | |
203 } | |
204 | |
205 // traditional ssim as per: http://en.wikipedia.org/wiki/Structural_similarity | |
206 // | |
207 // Re working out the math -> | |
208 // | |
209 // ssim(x,y) = (2*mean(x)*mean(y) + c1)*(2*cov(x,y)+c2) / | |
210 // ((mean(x)^2+mean(y)^2+c1)*(var(x)+var(y)+c2)) | |
211 // | |
212 // mean(x) = sum(x) / n | |
213 // | |
214 // cov(x,y) = (n*sum(xi*yi)-sum(x)*sum(y))/(n*n) | |
215 // | |
216 // var(x) = (n*sum(xi*xi)-sum(xi)*sum(xi))/(n*n) | |
217 // | |
218 // ssim(x,y) = | |
219 // (2*sum(x)*sum(y)/(n*n) + c1)*(2*(n*sum(xi*yi)-sum(x)*sum(y))/(n*n)+c2) / | |
220 // (((sum(x)*sum(x)+sum(y)*sum(y))/(n*n) +c1) * | |
221 // ((n*sum(xi*xi) - sum(xi)*sum(xi))/(n*n)+ | |
222 // (n*sum(yi*yi) - sum(yi)*sum(yi))/(n*n)+c2))) | |
223 // | |
224 // factoring out n*n | |
225 // | |
226 // ssim(x,y) = | |
227 // (2*sum(x)*sum(y) + n*n*c1)*(2*(n*sum(xi*yi)-sum(x)*sum(y))+n*n*c2) / | |
228 // (((sum(x)*sum(x)+sum(y)*sum(y)) + n*n*c1) * | |
229 // (n*sum(xi*xi)-sum(xi)*sum(xi)+n*sum(yi*yi)-sum(yi)*sum(yi)+n*n*c2)) | |
230 // | |
231 // Replace c1 with n*n * c1 for the final step that leads to this code: | |
232 // The final step scales by 12 bits so we don't lose precision in the constants. | |
233 | |
234 double ssimv_similarity(Ssimv *sv, int64_t n) { | |
235 // Scale the constants by number of pixels. | |
236 const int64_t c1 = (cc1 * n * n) >> 12; | |
237 const int64_t c2 = (cc2 * n * n) >> 12; | |
238 | |
239 const double l = 1.0 * (2 * sv->sum_s * sv->sum_r + c1) / | |
240 (sv->sum_s * sv->sum_s + sv->sum_r * sv->sum_r + c1); | |
241 | |
242 // Since these variables are unsigned sums, convert to double so | |
243 // math is done in double arithmetic. | |
244 const double v = (2.0 * n * sv->sum_sxr - 2 * sv->sum_s * sv->sum_r + c2) | |
245 / (n * sv->sum_sq_s - sv->sum_s * sv->sum_s + n * sv->sum_sq_r | |
246 - sv->sum_r * sv->sum_r + c2); | |
247 | |
248 return l * v; | |
249 } | |
250 | |
251 // The first term of the ssim metric is a luminance factor. | |
252 // | |
253 // (2*mean(x)*mean(y) + c1)/ (mean(x)^2+mean(y)^2+c1) | |
254 // | |
255 // This luminance factor is super sensitive to the dark side of luminance | |
256 // values and completely insensitive on the white side. check out 2 sets | |
257 // (1,3) and (250,252) the term gives ( 2*1*3/(1+9) = .60 | |
258 // 2*250*252/ (250^2+252^2) => .99999997 | |
259 // | |
260 // As a result in this tweaked version of the calculation in which the | |
261 // luminance is taken as percentage off from peak possible. | |
262 // | |
263 // 255 * 255 - (sum_s - sum_r) / count * (sum_s - sum_r) / count | |
264 // | |
265 double ssimv_similarity2(Ssimv *sv, int64_t n) { | |
266 // Scale the constants by number of pixels. | |
267 const int64_t c1 = (cc1 * n * n) >> 12; | |
268 const int64_t c2 = (cc2 * n * n) >> 12; | |
269 | |
270 const double mean_diff = (1.0 * sv->sum_s - sv->sum_r) / n; | |
271 const double l = (255 * 255 - mean_diff * mean_diff + c1) / (255 * 255 + c1); | |
272 | |
273 // Since these variables are unsigned, sums convert to double so | |
274 // math is done in double arithmetic. | |
275 const double v = (2.0 * n * sv->sum_sxr - 2 * sv->sum_s * sv->sum_r + c2) | |
276 / (n * sv->sum_sq_s - sv->sum_s * sv->sum_s + | |
277 n * sv->sum_sq_r - sv->sum_r * sv->sum_r + c2); | |
278 | |
279 return l * v; | |
280 } | |
281 void ssimv_parms(uint8_t *img1, int img1_pitch, uint8_t *img2, int img2_pitch, | |
282 Ssimv *sv) { | |
283 vp9_ssim_parms_8x8(img1, img1_pitch, img2, img2_pitch, | |
284 &sv->sum_s, &sv->sum_r, &sv->sum_sq_s, &sv->sum_sq_r, | |
285 &sv->sum_sxr); | |
286 } | |
287 | |
288 double vp9_get_ssim_metrics(uint8_t *img1, int img1_pitch, | |
289 uint8_t *img2, int img2_pitch, | |
290 int width, int height, | |
291 Ssimv *sv2, Metrics *m, | |
292 int do_inconsistency) { | |
293 double dssim_total = 0; | |
294 double ssim_total = 0; | |
295 double ssim2_total = 0; | |
296 double inconsistency_total = 0; | |
297 int i, j; | |
298 int c = 0; | |
299 double norm; | |
300 double old_ssim_total = 0; | |
301 vp9_clear_system_state(); | |
302 // We can sample points as frequently as we like start with 1 per 4x4. | |
303 for (i = 0; i < height; i += 4, | |
304 img1 += img1_pitch * 4, img2 += img2_pitch * 4) { | |
305 for (j = 0; j < width; j += 4, ++c) { | |
306 Ssimv sv = {0}; | |
307 double ssim; | |
308 double ssim2; | |
309 double dssim; | |
310 uint32_t var_new; | |
311 uint32_t var_old; | |
312 uint32_t mean_new; | |
313 uint32_t mean_old; | |
314 double ssim_new; | |
315 double ssim_old; | |
316 | |
317 // Not sure there's a great way to handle the edge pixels | |
318 // in ssim when using a window. Seems biased against edge pixels | |
319 // however you handle this. This uses only samples that are | |
320 // fully in the frame. | |
321 if (j + 8 <= width && i + 8 <= height) { | |
322 ssimv_parms(img1 + j, img1_pitch, img2 + j, img2_pitch, &sv); | |
323 } | |
324 | |
325 ssim = ssimv_similarity(&sv, 64); | |
326 ssim2 = ssimv_similarity2(&sv, 64); | |
327 | |
328 sv.ssim = ssim2; | |
329 | |
330 // dssim is calculated to use as an actual error metric and | |
331 // is scaled up to the same range as sum square error. | |
332 // Since we are subsampling every 16th point maybe this should be | |
333 // *16 ? | |
334 dssim = 255 * 255 * (1 - ssim2) / 2; | |
335 | |
336 // Here I introduce a new error metric: consistency-weighted | |
337 // SSIM-inconsistency. This metric isolates frames where the | |
338 // SSIM 'suddenly' changes, e.g. if one frame in every 8 is much | |
339 // sharper or blurrier than the others. Higher values indicate a | |
340 // temporally inconsistent SSIM. There are two ideas at work: | |
341 // | |
342 // 1) 'SSIM-inconsistency': the total inconsistency value | |
343 // reflects how much SSIM values are changing between this | |
344 // source / reference frame pair and the previous pair. | |
345 // | |
346 // 2) 'consistency-weighted': weights de-emphasize areas in the | |
347 // frame where the scene content has changed. Changes in scene | |
348 // content are detected via changes in local variance and local | |
349 // mean. | |
350 // | |
351 // Thus the overall measure reflects how inconsistent the SSIM | |
352 // values are, over consistent regions of the frame. | |
353 // | |
354 // The metric has three terms: | |
355 // | |
356 // term 1 -> uses change in scene Variance to weight error score | |
357 // 2 * var(Fi)*var(Fi-1) / (var(Fi)^2+var(Fi-1)^2) | |
358 // larger changes from one frame to the next mean we care | |
359 // less about consistency. | |
360 // | |
361 // term 2 -> uses change in local scene luminance to weight error | |
362 // 2 * avg(Fi)*avg(Fi-1) / (avg(Fi)^2+avg(Fi-1)^2) | |
363 // larger changes from one frame to the next mean we care | |
364 // less about consistency. | |
365 // | |
366 // term3 -> measures inconsistency in ssim scores between frames | |
367 // 1 - ( 2 * ssim(Fi)*ssim(Fi-1)/(ssim(Fi)^2+sssim(Fi-1)^2). | |
368 // | |
369 // This term compares the ssim score for the same location in 2 | |
370 // subsequent frames. | |
371 var_new = sv.sum_sq_s - sv.sum_s * sv.sum_s / 64; | |
372 var_old = sv2[c].sum_sq_s - sv2[c].sum_s * sv2[c].sum_s / 64; | |
373 mean_new = sv.sum_s; | |
374 mean_old = sv2[c].sum_s; | |
375 ssim_new = sv.ssim; | |
376 ssim_old = sv2[c].ssim; | |
377 | |
378 if (do_inconsistency) { | |
379 // We do the metric once for every 4x4 block in the image. Since | |
380 // we are scaling the error to SSE for use in a psnr calculation | |
381 // 1.0 = 4x4x255x255 the worst error we can possibly have. | |
382 static const double kScaling = 4. * 4 * 255 * 255; | |
383 | |
384 // The constants have to be non 0 to avoid potential divide by 0 | |
385 // issues other than that they affect kind of a weighting between | |
386 // the terms. No testing of what the right terms should be has been | |
387 // done. | |
388 static const double c1 = 1, c2 = 1, c3 = 1; | |
389 | |
390 // This measures how much consistent variance is in two consecutive | |
391 // source frames. 1.0 means they have exactly the same variance. | |
392 const double variance_term = (2.0 * var_old * var_new + c1) / | |
393 (1.0 * var_old * var_old + 1.0 * var_new * var_new + c1); | |
394 | |
395 // This measures how consistent the local mean are between two | |
396 // consecutive frames. 1.0 means they have exactly the same mean. | |
397 const double mean_term = (2.0 * mean_old * mean_new + c2) / | |
398 (1.0 * mean_old * mean_old + 1.0 * mean_new * mean_new + c2); | |
399 | |
400 // This measures how consistent the ssims of two | |
401 // consecutive frames is. 1.0 means they are exactly the same. | |
402 double ssim_term = pow((2.0 * ssim_old * ssim_new + c3) / | |
403 (ssim_old * ssim_old + ssim_new * ssim_new + c3), | |
404 5); | |
405 | |
406 double this_inconsistency; | |
407 | |
408 // Floating point math sometimes makes this > 1 by a tiny bit. | |
409 // We want the metric to scale between 0 and 1.0 so we can convert | |
410 // it to an snr scaled value. | |
411 if (ssim_term > 1) | |
412 ssim_term = 1; | |
413 | |
414 // This converts the consistency metric to an inconsistency metric | |
415 // ( so we can scale it like psnr to something like sum square error. | |
416 // The reason for the variance and mean terms is the assumption that | |
417 // if there are big changes in the source we shouldn't penalize | |
418 // inconsistency in ssim scores a bit less as it will be less visible | |
419 // to the user. | |
420 this_inconsistency = (1 - ssim_term) * variance_term * mean_term; | |
421 | |
422 this_inconsistency *= kScaling; | |
423 inconsistency_total += this_inconsistency; | |
424 } | |
425 sv2[c] = sv; | |
426 ssim_total += ssim; | |
427 ssim2_total += ssim2; | |
428 dssim_total += dssim; | |
429 | |
430 old_ssim_total += ssim_old; | |
431 } | |
432 old_ssim_total += 0; | |
433 } | |
434 | |
435 norm = 1. / (width / 4) / (height / 4); | |
436 ssim_total *= norm; | |
437 ssim2_total *= norm; | |
438 m->ssim2 = ssim2_total; | |
439 m->ssim = ssim_total; | |
440 if (old_ssim_total == 0) | |
441 inconsistency_total = 0; | |
442 | |
443 m->ssimc = inconsistency_total; | |
444 | |
445 m->dssim = dssim_total; | |
446 return inconsistency_total; | |
447 } | |
448 | |
449 | |
450 #if CONFIG_VP9_HIGHBITDEPTH | |
451 double vp9_highbd_calc_ssim(YV12_BUFFER_CONFIG *source, | |
452 YV12_BUFFER_CONFIG *dest, | |
453 double *weight, unsigned int bd) { | |
454 double a, b, c; | |
455 double ssimv; | |
456 | |
457 a = vp9_highbd_ssim2(source->y_buffer, dest->y_buffer, | |
458 source->y_stride, dest->y_stride, | |
459 source->y_crop_width, source->y_crop_height, bd); | |
460 | |
461 b = vp9_highbd_ssim2(source->u_buffer, dest->u_buffer, | |
462 source->uv_stride, dest->uv_stride, | |
463 source->uv_crop_width, source->uv_crop_height, bd); | |
464 | |
465 c = vp9_highbd_ssim2(source->v_buffer, dest->v_buffer, | |
466 source->uv_stride, dest->uv_stride, | |
467 source->uv_crop_width, source->uv_crop_height, bd); | |
468 | |
469 ssimv = a * .8 + .1 * (b + c); | |
470 | |
471 *weight = 1; | |
472 | |
473 return ssimv; | |
474 } | |
475 | |
476 double vp9_highbd_calc_ssimg(YV12_BUFFER_CONFIG *source, | |
477 YV12_BUFFER_CONFIG *dest, double *ssim_y, | |
478 double *ssim_u, double *ssim_v, unsigned int bd) { | |
479 double ssim_all = 0; | |
480 double a, b, c; | |
481 | |
482 a = vp9_highbd_ssim2(source->y_buffer, dest->y_buffer, | |
483 source->y_stride, dest->y_stride, | |
484 source->y_crop_width, source->y_crop_height, bd); | |
485 | |
486 b = vp9_highbd_ssim2(source->u_buffer, dest->u_buffer, | |
487 source->uv_stride, dest->uv_stride, | |
488 source->uv_crop_width, source->uv_crop_height, bd); | |
489 | |
490 c = vp9_highbd_ssim2(source->v_buffer, dest->v_buffer, | |
491 source->uv_stride, dest->uv_stride, | |
492 source->uv_crop_width, source->uv_crop_height, bd); | |
493 *ssim_y = a; | |
494 *ssim_u = b; | |
495 *ssim_v = c; | |
496 ssim_all = (a * 4 + b + c) / 6; | |
497 | |
498 return ssim_all; | |
499 } | |
500 #endif // CONFIG_VP9_HIGHBITDEPTH | |
OLD | NEW |