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.
12 #include "./vpx_dsp_rtcd.h"
13 #include "vpx_dsp/ssim.h"
14 #include "vpx_ports/mem.h"
15 #include "vpx_ports/system_state.h"
17 void vpx_ssim_parms_16x16_c(const uint8_t *s
, int sp
, const uint8_t *r
,
18 int rp
, uint32_t *sum_s
, uint32_t *sum_r
,
19 uint32_t *sum_sq_s
, uint32_t *sum_sq_r
,
22 for (i
= 0; i
< 16; i
++, s
+= sp
, r
+= rp
) {
23 for (j
= 0; j
< 16; j
++) {
26 *sum_sq_s
+= s
[j
] * s
[j
];
27 *sum_sq_r
+= r
[j
] * r
[j
];
28 *sum_sxr
+= s
[j
] * r
[j
];
32 void vpx_ssim_parms_8x8_c(const uint8_t *s
, int sp
, const uint8_t *r
, int rp
,
33 uint32_t *sum_s
, uint32_t *sum_r
,
34 uint32_t *sum_sq_s
, uint32_t *sum_sq_r
,
37 for (i
= 0; i
< 8; i
++, s
+= sp
, r
+= rp
) {
38 for (j
= 0; j
< 8; j
++) {
41 *sum_sq_s
+= s
[j
] * s
[j
];
42 *sum_sq_r
+= r
[j
] * r
[j
];
43 *sum_sxr
+= s
[j
] * r
[j
];
48 #if CONFIG_VP9_HIGHBITDEPTH
49 void vpx_highbd_ssim_parms_8x8_c(const uint16_t *s
, int sp
,
50 const uint16_t *r
, int rp
,
51 uint32_t *sum_s
, uint32_t *sum_r
,
52 uint32_t *sum_sq_s
, uint32_t *sum_sq_r
,
55 for (i
= 0; i
< 8; i
++, s
+= sp
, r
+= rp
) {
56 for (j
= 0; j
< 8; j
++) {
59 *sum_sq_s
+= s
[j
] * s
[j
];
60 *sum_sq_r
+= r
[j
] * r
[j
];
61 *sum_sxr
+= s
[j
] * r
[j
];
65 #endif // CONFIG_VP9_HIGHBITDEPTH
67 static const int64_t cc1
= 26634; // (64^2*(.01*255)^2
68 static const int64_t cc2
= 239708; // (64^2*(.03*255)^2
70 static double similarity(uint32_t sum_s
, uint32_t sum_r
,
71 uint32_t sum_sq_s
, uint32_t sum_sq_r
,
72 uint32_t sum_sxr
, int count
) {
73 int64_t ssim_n
, ssim_d
;
76 // scale the constants by number of pixels
77 c1
= (cc1
* count
* count
) >> 12;
78 c2
= (cc2
* count
* count
) >> 12;
80 ssim_n
= (2 * sum_s
* sum_r
+ c1
) * ((int64_t) 2 * count
* sum_sxr
-
81 (int64_t) 2 * sum_s
* sum_r
+ c2
);
83 ssim_d
= (sum_s
* sum_s
+ sum_r
* sum_r
+ c1
) *
84 ((int64_t)count
* sum_sq_s
- (int64_t)sum_s
* sum_s
+
85 (int64_t)count
* sum_sq_r
- (int64_t) sum_r
* sum_r
+ c2
);
87 return ssim_n
* 1.0 / ssim_d
;
90 static double ssim_8x8(const uint8_t *s
, int sp
, const uint8_t *r
, int rp
) {
91 uint32_t sum_s
= 0, sum_r
= 0, sum_sq_s
= 0, sum_sq_r
= 0, sum_sxr
= 0;
92 vpx_ssim_parms_8x8(s
, sp
, r
, rp
, &sum_s
, &sum_r
, &sum_sq_s
, &sum_sq_r
,
94 return similarity(sum_s
, sum_r
, sum_sq_s
, sum_sq_r
, sum_sxr
, 64);
97 #if CONFIG_VP9_HIGHBITDEPTH
98 static double highbd_ssim_8x8(const uint16_t *s
, int sp
, const uint16_t *r
,
99 int rp
, unsigned int bd
) {
100 uint32_t sum_s
= 0, sum_r
= 0, sum_sq_s
= 0, sum_sq_r
= 0, sum_sxr
= 0;
101 const int oshift
= bd
- 8;
102 vpx_highbd_ssim_parms_8x8(s
, sp
, r
, rp
, &sum_s
, &sum_r
, &sum_sq_s
, &sum_sq_r
,
104 return similarity(sum_s
>> oshift
,
106 sum_sq_s
>> (2 * oshift
),
107 sum_sq_r
>> (2 * oshift
),
108 sum_sxr
>> (2 * oshift
),
111 #endif // CONFIG_VP9_HIGHBITDEPTH
113 // We are using a 8x8 moving window with starting location of each 8x8 window
114 // on the 4x4 pixel grid. Such arrangement allows the windows to overlap
115 // block boundaries to penalize blocking artifacts.
116 static double vpx_ssim2(const uint8_t *img1
, const uint8_t *img2
,
117 int stride_img1
, int stride_img2
, int width
,
121 double ssim_total
= 0;
123 // sample point start with each 4x4 location
124 for (i
= 0; i
<= height
- 8;
125 i
+= 4, img1
+= stride_img1
* 4, img2
+= stride_img2
* 4) {
126 for (j
= 0; j
<= width
- 8; j
+= 4) {
127 double v
= ssim_8x8(img1
+ j
, stride_img1
, img2
+ j
, stride_img2
);
132 ssim_total
/= samples
;
136 #if CONFIG_VP9_HIGHBITDEPTH
137 static double vpx_highbd_ssim2(const uint8_t *img1
, const uint8_t *img2
,
138 int stride_img1
, int stride_img2
, int width
,
139 int height
, unsigned int bd
) {
142 double ssim_total
= 0;
144 // sample point start with each 4x4 location
145 for (i
= 0; i
<= height
- 8;
146 i
+= 4, img1
+= stride_img1
* 4, img2
+= stride_img2
* 4) {
147 for (j
= 0; j
<= width
- 8; j
+= 4) {
148 double v
= highbd_ssim_8x8(CONVERT_TO_SHORTPTR(img1
+ j
), stride_img1
,
149 CONVERT_TO_SHORTPTR(img2
+ j
), stride_img2
,
155 ssim_total
/= samples
;
158 #endif // CONFIG_VP9_HIGHBITDEPTH
160 double vpx_calc_ssim(const YV12_BUFFER_CONFIG
*source
,
161 const YV12_BUFFER_CONFIG
*dest
,
166 a
= vpx_ssim2(source
->y_buffer
, dest
->y_buffer
,
167 source
->y_stride
, dest
->y_stride
,
168 source
->y_crop_width
, source
->y_crop_height
);
170 b
= vpx_ssim2(source
->u_buffer
, dest
->u_buffer
,
171 source
->uv_stride
, dest
->uv_stride
,
172 source
->uv_crop_width
, source
->uv_crop_height
);
174 c
= vpx_ssim2(source
->v_buffer
, dest
->v_buffer
,
175 source
->uv_stride
, dest
->uv_stride
,
176 source
->uv_crop_width
, source
->uv_crop_height
);
178 ssimv
= a
* .8 + .1 * (b
+ c
);
185 double vpx_calc_ssimg(const YV12_BUFFER_CONFIG
*source
,
186 const YV12_BUFFER_CONFIG
*dest
,
187 double *ssim_y
, double *ssim_u
, double *ssim_v
) {
191 a
= vpx_ssim2(source
->y_buffer
, dest
->y_buffer
,
192 source
->y_stride
, dest
->y_stride
,
193 source
->y_crop_width
, source
->y_crop_height
);
195 b
= vpx_ssim2(source
->u_buffer
, dest
->u_buffer
,
196 source
->uv_stride
, dest
->uv_stride
,
197 source
->uv_crop_width
, source
->uv_crop_height
);
199 c
= vpx_ssim2(source
->v_buffer
, dest
->v_buffer
,
200 source
->uv_stride
, dest
->uv_stride
,
201 source
->uv_crop_width
, source
->uv_crop_height
);
205 ssim_all
= (a
* 4 + b
+ c
) / 6;
210 // traditional ssim as per: http://en.wikipedia.org/wiki/Structural_similarity
212 // Re working out the math ->
214 // ssim(x,y) = (2*mean(x)*mean(y) + c1)*(2*cov(x,y)+c2) /
215 // ((mean(x)^2+mean(y)^2+c1)*(var(x)+var(y)+c2))
217 // mean(x) = sum(x) / n
219 // cov(x,y) = (n*sum(xi*yi)-sum(x)*sum(y))/(n*n)
221 // var(x) = (n*sum(xi*xi)-sum(xi)*sum(xi))/(n*n)
224 // (2*sum(x)*sum(y)/(n*n) + c1)*(2*(n*sum(xi*yi)-sum(x)*sum(y))/(n*n)+c2) /
225 // (((sum(x)*sum(x)+sum(y)*sum(y))/(n*n) +c1) *
226 // ((n*sum(xi*xi) - sum(xi)*sum(xi))/(n*n)+
227 // (n*sum(yi*yi) - sum(yi)*sum(yi))/(n*n)+c2)))
232 // (2*sum(x)*sum(y) + n*n*c1)*(2*(n*sum(xi*yi)-sum(x)*sum(y))+n*n*c2) /
233 // (((sum(x)*sum(x)+sum(y)*sum(y)) + n*n*c1) *
234 // (n*sum(xi*xi)-sum(xi)*sum(xi)+n*sum(yi*yi)-sum(yi)*sum(yi)+n*n*c2))
236 // Replace c1 with n*n * c1 for the final step that leads to this code:
237 // The final step scales by 12 bits so we don't lose precision in the constants.
239 static double ssimv_similarity(const Ssimv
*sv
, int64_t n
) {
240 // Scale the constants by number of pixels.
241 const int64_t c1
= (cc1
* n
* n
) >> 12;
242 const int64_t c2
= (cc2
* n
* n
) >> 12;
244 const double l
= 1.0 * (2 * sv
->sum_s
* sv
->sum_r
+ c1
) /
245 (sv
->sum_s
* sv
->sum_s
+ sv
->sum_r
* sv
->sum_r
+ c1
);
247 // Since these variables are unsigned sums, convert to double so
248 // math is done in double arithmetic.
249 const double v
= (2.0 * n
* sv
->sum_sxr
- 2 * sv
->sum_s
* sv
->sum_r
+ c2
)
250 / (n
* sv
->sum_sq_s
- sv
->sum_s
* sv
->sum_s
+ n
* sv
->sum_sq_r
251 - sv
->sum_r
* sv
->sum_r
+ c2
);
256 // The first term of the ssim metric is a luminance factor.
258 // (2*mean(x)*mean(y) + c1)/ (mean(x)^2+mean(y)^2+c1)
260 // This luminance factor is super sensitive to the dark side of luminance
261 // values and completely insensitive on the white side. check out 2 sets
262 // (1,3) and (250,252) the term gives ( 2*1*3/(1+9) = .60
263 // 2*250*252/ (250^2+252^2) => .99999997
265 // As a result in this tweaked version of the calculation in which the
266 // luminance is taken as percentage off from peak possible.
268 // 255 * 255 - (sum_s - sum_r) / count * (sum_s - sum_r) / count
270 static double ssimv_similarity2(const Ssimv
*sv
, int64_t n
) {
271 // Scale the constants by number of pixels.
272 const int64_t c1
= (cc1
* n
* n
) >> 12;
273 const int64_t c2
= (cc2
* n
* n
) >> 12;
275 const double mean_diff
= (1.0 * sv
->sum_s
- sv
->sum_r
) / n
;
276 const double l
= (255 * 255 - mean_diff
* mean_diff
+ c1
) / (255 * 255 + c1
);
278 // Since these variables are unsigned, sums convert to double so
279 // math is done in double arithmetic.
280 const double v
= (2.0 * n
* sv
->sum_sxr
- 2 * sv
->sum_s
* sv
->sum_r
+ c2
)
281 / (n
* sv
->sum_sq_s
- sv
->sum_s
* sv
->sum_s
+
282 n
* sv
->sum_sq_r
- sv
->sum_r
* sv
->sum_r
+ c2
);
286 static void ssimv_parms(uint8_t *img1
, int img1_pitch
, uint8_t *img2
,
287 int img2_pitch
, Ssimv
*sv
) {
288 vpx_ssim_parms_8x8(img1
, img1_pitch
, img2
, img2_pitch
,
289 &sv
->sum_s
, &sv
->sum_r
, &sv
->sum_sq_s
, &sv
->sum_sq_r
,
293 double vpx_get_ssim_metrics(uint8_t *img1
, int img1_pitch
,
294 uint8_t *img2
, int img2_pitch
,
295 int width
, int height
,
296 Ssimv
*sv2
, Metrics
*m
,
297 int do_inconsistency
) {
298 double dssim_total
= 0;
299 double ssim_total
= 0;
300 double ssim2_total
= 0;
301 double inconsistency_total
= 0;
305 double old_ssim_total
= 0;
306 vpx_clear_system_state();
307 // We can sample points as frequently as we like start with 1 per 4x4.
308 for (i
= 0; i
< height
; i
+= 4,
309 img1
+= img1_pitch
* 4, img2
+= img2_pitch
* 4) {
310 for (j
= 0; j
< width
; j
+= 4, ++c
) {
322 // Not sure there's a great way to handle the edge pixels
323 // in ssim when using a window. Seems biased against edge pixels
324 // however you handle this. This uses only samples that are
325 // fully in the frame.
326 if (j
+ 8 <= width
&& i
+ 8 <= height
) {
327 ssimv_parms(img1
+ j
, img1_pitch
, img2
+ j
, img2_pitch
, &sv
);
330 ssim
= ssimv_similarity(&sv
, 64);
331 ssim2
= ssimv_similarity2(&sv
, 64);
335 // dssim is calculated to use as an actual error metric and
336 // is scaled up to the same range as sum square error.
337 // Since we are subsampling every 16th point maybe this should be
339 dssim
= 255 * 255 * (1 - ssim2
) / 2;
341 // Here I introduce a new error metric: consistency-weighted
342 // SSIM-inconsistency. This metric isolates frames where the
343 // SSIM 'suddenly' changes, e.g. if one frame in every 8 is much
344 // sharper or blurrier than the others. Higher values indicate a
345 // temporally inconsistent SSIM. There are two ideas at work:
347 // 1) 'SSIM-inconsistency': the total inconsistency value
348 // reflects how much SSIM values are changing between this
349 // source / reference frame pair and the previous pair.
351 // 2) 'consistency-weighted': weights de-emphasize areas in the
352 // frame where the scene content has changed. Changes in scene
353 // content are detected via changes in local variance and local
356 // Thus the overall measure reflects how inconsistent the SSIM
357 // values are, over consistent regions of the frame.
359 // The metric has three terms:
361 // term 1 -> uses change in scene Variance to weight error score
362 // 2 * var(Fi)*var(Fi-1) / (var(Fi)^2+var(Fi-1)^2)
363 // larger changes from one frame to the next mean we care
364 // less about consistency.
366 // term 2 -> uses change in local scene luminance to weight error
367 // 2 * avg(Fi)*avg(Fi-1) / (avg(Fi)^2+avg(Fi-1)^2)
368 // larger changes from one frame to the next mean we care
369 // less about consistency.
371 // term3 -> measures inconsistency in ssim scores between frames
372 // 1 - ( 2 * ssim(Fi)*ssim(Fi-1)/(ssim(Fi)^2+sssim(Fi-1)^2).
374 // This term compares the ssim score for the same location in 2
375 // subsequent frames.
376 var_new
= sv
.sum_sq_s
- sv
.sum_s
* sv
.sum_s
/ 64;
377 var_old
= sv2
[c
].sum_sq_s
- sv2
[c
].sum_s
* sv2
[c
].sum_s
/ 64;
379 mean_old
= sv2
[c
].sum_s
;
381 ssim_old
= sv2
[c
].ssim
;
383 if (do_inconsistency
) {
384 // We do the metric once for every 4x4 block in the image. Since
385 // we are scaling the error to SSE for use in a psnr calculation
386 // 1.0 = 4x4x255x255 the worst error we can possibly have.
387 static const double kScaling
= 4. * 4 * 255 * 255;
389 // The constants have to be non 0 to avoid potential divide by 0
390 // issues other than that they affect kind of a weighting between
391 // the terms. No testing of what the right terms should be has been
393 static const double c1
= 1, c2
= 1, c3
= 1;
395 // This measures how much consistent variance is in two consecutive
396 // source frames. 1.0 means they have exactly the same variance.
397 const double variance_term
= (2.0 * var_old
* var_new
+ c1
) /
398 (1.0 * var_old
* var_old
+ 1.0 * var_new
* var_new
+ c1
);
400 // This measures how consistent the local mean are between two
401 // consecutive frames. 1.0 means they have exactly the same mean.
402 const double mean_term
= (2.0 * mean_old
* mean_new
+ c2
) /
403 (1.0 * mean_old
* mean_old
+ 1.0 * mean_new
* mean_new
+ c2
);
405 // This measures how consistent the ssims of two
406 // consecutive frames is. 1.0 means they are exactly the same.
407 double ssim_term
= pow((2.0 * ssim_old
* ssim_new
+ c3
) /
408 (ssim_old
* ssim_old
+ ssim_new
* ssim_new
+ c3
),
411 double this_inconsistency
;
413 // Floating point math sometimes makes this > 1 by a tiny bit.
414 // We want the metric to scale between 0 and 1.0 so we can convert
415 // it to an snr scaled value.
419 // This converts the consistency metric to an inconsistency metric
420 // ( so we can scale it like psnr to something like sum square error.
421 // The reason for the variance and mean terms is the assumption that
422 // if there are big changes in the source we shouldn't penalize
423 // inconsistency in ssim scores a bit less as it will be less visible
425 this_inconsistency
= (1 - ssim_term
) * variance_term
* mean_term
;
427 this_inconsistency
*= kScaling
;
428 inconsistency_total
+= this_inconsistency
;
432 ssim2_total
+= ssim2
;
433 dssim_total
+= dssim
;
435 old_ssim_total
+= ssim_old
;
440 norm
= 1. / (width
/ 4) / (height
/ 4);
443 m
->ssim2
= ssim2_total
;
444 m
->ssim
= ssim_total
;
445 if (old_ssim_total
== 0)
446 inconsistency_total
= 0;
448 m
->ssimc
= inconsistency_total
;
450 m
->dssim
= dssim_total
;
451 return inconsistency_total
;
455 #if CONFIG_VP9_HIGHBITDEPTH
456 double vpx_highbd_calc_ssim(const YV12_BUFFER_CONFIG
*source
,
457 const YV12_BUFFER_CONFIG
*dest
,
458 double *weight
, unsigned int bd
) {
462 a
= vpx_highbd_ssim2(source
->y_buffer
, dest
->y_buffer
,
463 source
->y_stride
, dest
->y_stride
,
464 source
->y_crop_width
, source
->y_crop_height
, bd
);
466 b
= vpx_highbd_ssim2(source
->u_buffer
, dest
->u_buffer
,
467 source
->uv_stride
, dest
->uv_stride
,
468 source
->uv_crop_width
, source
->uv_crop_height
, bd
);
470 c
= vpx_highbd_ssim2(source
->v_buffer
, dest
->v_buffer
,
471 source
->uv_stride
, dest
->uv_stride
,
472 source
->uv_crop_width
, source
->uv_crop_height
, bd
);
474 ssimv
= a
* .8 + .1 * (b
+ c
);
481 double vpx_highbd_calc_ssimg(const YV12_BUFFER_CONFIG
*source
,
482 const YV12_BUFFER_CONFIG
*dest
, double *ssim_y
,
483 double *ssim_u
, double *ssim_v
, unsigned int bd
) {
487 a
= vpx_highbd_ssim2(source
->y_buffer
, dest
->y_buffer
,
488 source
->y_stride
, dest
->y_stride
,
489 source
->y_crop_width
, source
->y_crop_height
, bd
);
491 b
= vpx_highbd_ssim2(source
->u_buffer
, dest
->u_buffer
,
492 source
->uv_stride
, dest
->uv_stride
,
493 source
->uv_crop_width
, source
->uv_crop_height
, bd
);
495 c
= vpx_highbd_ssim2(source
->v_buffer
, dest
->v_buffer
,
496 source
->uv_stride
, dest
->uv_stride
,
497 source
->uv_crop_width
, source
->uv_crop_height
, bd
);
501 ssim_all
= (a
* 4 + b
+ c
) / 6;
505 #endif // CONFIG_VP9_HIGHBITDEPTH