1 /*
2 * Copyright (c) 2016, Alliance for Open Media. All rights reserved.
3 *
4 * This source code is subject to the terms of the BSD 2 Clause License and
5 * the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License
6 * was not distributed with this source code in the LICENSE file, you can
7 * obtain it at www.aomedia.org/license/software. If the Alliance for Open
8 * Media Patent License 1.0 was not distributed with this source code in the
9 * PATENTS file, you can obtain it at www.aomedia.org/license/patent.
10 */
11
12 #include <assert.h>
13 #include <math.h>
14
15 #include "config/aom_config.h"
16 #include "config/aom_dsp_rtcd.h"
17
18 #include "aom_dsp/psnr.h"
19 #include "aom_scale/yv12config.h"
20
21 #if CONFIG_INTERNAL_STATS
22 #define STATIC
23 #else
24 #define STATIC static
25 #endif // CONFIG_INTERNAL_STATS
26
aom_sse_to_psnr(double samples,double peak,double sse)27 STATIC double aom_sse_to_psnr(double samples, double peak, double sse) {
28 if (sse > 0.0) {
29 const double psnr = 10.0 * log10(samples * peak * peak / sse);
30 return psnr > MAX_PSNR ? MAX_PSNR : psnr;
31 } else {
32 return MAX_PSNR;
33 }
34 }
35
36 #undef STATIC
37
encoder_sse(const uint8_t * a,int a_stride,const uint8_t * b,int b_stride,int w,int h)38 static int64_t encoder_sse(const uint8_t *a, int a_stride, const uint8_t *b,
39 int b_stride, int w, int h) {
40 int i, j;
41 int64_t sse = 0;
42
43 for (i = 0; i < h; i++) {
44 for (j = 0; j < w; j++) {
45 const int diff = a[j] - b[j];
46 sse += diff * diff;
47 }
48
49 a += a_stride;
50 b += b_stride;
51 }
52 return sse;
53 }
54
55 #if CONFIG_AV1_HIGHBITDEPTH
encoder_highbd_sse(const uint8_t * a8,int a_stride,const uint8_t * b8,int b_stride,int w,int h)56 static int64_t encoder_highbd_sse(const uint8_t *a8, int a_stride,
57 const uint8_t *b8, int b_stride, int w,
58 int h) {
59 const uint16_t *a = CONVERT_TO_SHORTPTR(a8);
60 const uint16_t *b = CONVERT_TO_SHORTPTR(b8);
61 int64_t sse = 0;
62 for (int i = 0; i < h; ++i) {
63 for (int j = 0; j < w; ++j) {
64 const int diff = a[j] - b[j];
65 sse += diff * diff;
66 }
67 a += a_stride;
68 b += b_stride;
69 }
70 return sse;
71 }
72
73 #endif // CONFIG_AV1_HIGHBITDEPTH
74
get_sse(const uint8_t * a,int a_stride,const uint8_t * b,int b_stride,int width,int height)75 static int64_t get_sse(const uint8_t *a, int a_stride, const uint8_t *b,
76 int b_stride, int width, int height) {
77 const int dw = width % 16;
78 const int dh = height % 16;
79 int64_t total_sse = 0;
80 int x, y;
81
82 if (dw > 0) {
83 total_sse += encoder_sse(&a[width - dw], a_stride, &b[width - dw], b_stride,
84 dw, height);
85 }
86
87 if (dh > 0) {
88 total_sse +=
89 encoder_sse(&a[(height - dh) * a_stride], a_stride,
90 &b[(height - dh) * b_stride], b_stride, width - dw, dh);
91 }
92
93 for (y = 0; y < height / 16; ++y) {
94 const uint8_t *pa = a;
95 const uint8_t *pb = b;
96 for (x = 0; x < width / 16; ++x) {
97 total_sse += aom_sse(pa, a_stride, pb, b_stride, 16, 16);
98
99 pa += 16;
100 pb += 16;
101 }
102
103 a += 16 * a_stride;
104 b += 16 * b_stride;
105 }
106
107 return total_sse;
108 }
109
110 #if CONFIG_AV1_HIGHBITDEPTH
highbd_get_sse_shift(const uint8_t * a8,int a_stride,const uint8_t * b8,int b_stride,int width,int height,unsigned int input_shift)111 static int64_t highbd_get_sse_shift(const uint8_t *a8, int a_stride,
112 const uint8_t *b8, int b_stride, int width,
113 int height, unsigned int input_shift) {
114 const uint16_t *a = CONVERT_TO_SHORTPTR(a8);
115 const uint16_t *b = CONVERT_TO_SHORTPTR(b8);
116 int64_t total_sse = 0;
117 int x, y;
118 for (y = 0; y < height; ++y) {
119 for (x = 0; x < width; ++x) {
120 int64_t diff;
121 diff = (a[x] >> input_shift) - (b[x] >> input_shift);
122 total_sse += diff * diff;
123 }
124 a += a_stride;
125 b += b_stride;
126 }
127 return total_sse;
128 }
129
highbd_get_sse(const uint8_t * a,int a_stride,const uint8_t * b,int b_stride,int width,int height)130 static int64_t highbd_get_sse(const uint8_t *a, int a_stride, const uint8_t *b,
131 int b_stride, int width, int height) {
132 int64_t total_sse = 0;
133 int x, y;
134 const int dw = width % 16;
135 const int dh = height % 16;
136
137 if (dw > 0) {
138 total_sse += encoder_highbd_sse(&a[width - dw], a_stride, &b[width - dw],
139 b_stride, dw, height);
140 }
141 if (dh > 0) {
142 total_sse += encoder_highbd_sse(&a[(height - dh) * a_stride], a_stride,
143 &b[(height - dh) * b_stride], b_stride,
144 width - dw, dh);
145 }
146
147 for (y = 0; y < height / 16; ++y) {
148 const uint8_t *pa = a;
149 const uint8_t *pb = b;
150 for (x = 0; x < width / 16; ++x) {
151 total_sse += aom_highbd_sse(pa, a_stride, pb, b_stride, 16, 16);
152 pa += 16;
153 pb += 16;
154 }
155 a += 16 * a_stride;
156 b += 16 * b_stride;
157 }
158 return total_sse;
159 }
160 #endif // CONFIG_AV1_HIGHBITDEPTH
161
aom_get_y_var(const YV12_BUFFER_CONFIG * a,int hstart,int width,int vstart,int height)162 uint64_t aom_get_y_var(const YV12_BUFFER_CONFIG *a, int hstart, int width,
163 int vstart, int height) {
164 return aom_var_2d_u8(a->y_buffer + vstart * a->y_stride + hstart, a->y_stride,
165 width, height) /
166 (width * height);
167 }
168
aom_get_u_var(const YV12_BUFFER_CONFIG * a,int hstart,int width,int vstart,int height)169 uint64_t aom_get_u_var(const YV12_BUFFER_CONFIG *a, int hstart, int width,
170 int vstart, int height) {
171 return aom_var_2d_u8(a->u_buffer + vstart * a->uv_stride + hstart,
172 a->uv_stride, width, height) /
173 (width * height);
174 }
175
aom_get_v_var(const YV12_BUFFER_CONFIG * a,int hstart,int width,int vstart,int height)176 uint64_t aom_get_v_var(const YV12_BUFFER_CONFIG *a, int hstart, int width,
177 int vstart, int height) {
178 return aom_var_2d_u8(a->v_buffer + vstart * a->uv_stride + hstart,
179 a->uv_stride, width, height) /
180 (width * height);
181 }
182
aom_get_y_sse_part(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,int hstart,int width,int vstart,int height)183 int64_t aom_get_y_sse_part(const YV12_BUFFER_CONFIG *a,
184 const YV12_BUFFER_CONFIG *b, int hstart, int width,
185 int vstart, int height) {
186 return get_sse(a->y_buffer + vstart * a->y_stride + hstart, a->y_stride,
187 b->y_buffer + vstart * b->y_stride + hstart, b->y_stride,
188 width, height);
189 }
190
aom_get_y_sse(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b)191 int64_t aom_get_y_sse(const YV12_BUFFER_CONFIG *a,
192 const YV12_BUFFER_CONFIG *b) {
193 assert(a->y_crop_width == b->y_crop_width);
194 assert(a->y_crop_height == b->y_crop_height);
195
196 return get_sse(a->y_buffer, a->y_stride, b->y_buffer, b->y_stride,
197 a->y_crop_width, a->y_crop_height);
198 }
199
aom_get_u_sse_part(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,int hstart,int width,int vstart,int height)200 int64_t aom_get_u_sse_part(const YV12_BUFFER_CONFIG *a,
201 const YV12_BUFFER_CONFIG *b, int hstart, int width,
202 int vstart, int height) {
203 return get_sse(a->u_buffer + vstart * a->uv_stride + hstart, a->uv_stride,
204 b->u_buffer + vstart * b->uv_stride + hstart, b->uv_stride,
205 width, height);
206 }
207
aom_get_u_sse(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b)208 int64_t aom_get_u_sse(const YV12_BUFFER_CONFIG *a,
209 const YV12_BUFFER_CONFIG *b) {
210 assert(a->uv_crop_width == b->uv_crop_width);
211 assert(a->uv_crop_height == b->uv_crop_height);
212
213 return get_sse(a->u_buffer, a->uv_stride, b->u_buffer, b->uv_stride,
214 a->uv_crop_width, a->uv_crop_height);
215 }
216
aom_get_v_sse_part(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,int hstart,int width,int vstart,int height)217 int64_t aom_get_v_sse_part(const YV12_BUFFER_CONFIG *a,
218 const YV12_BUFFER_CONFIG *b, int hstart, int width,
219 int vstart, int height) {
220 return get_sse(a->v_buffer + vstart * a->uv_stride + hstart, a->uv_stride,
221 b->v_buffer + vstart * b->uv_stride + hstart, b->uv_stride,
222 width, height);
223 }
224
aom_get_v_sse(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b)225 int64_t aom_get_v_sse(const YV12_BUFFER_CONFIG *a,
226 const YV12_BUFFER_CONFIG *b) {
227 assert(a->uv_crop_width == b->uv_crop_width);
228 assert(a->uv_crop_height == b->uv_crop_height);
229
230 return get_sse(a->v_buffer, a->uv_stride, b->v_buffer, b->uv_stride,
231 a->uv_crop_width, a->uv_crop_height);
232 }
233
234 #if CONFIG_AV1_HIGHBITDEPTH
aom_highbd_get_y_var(const YV12_BUFFER_CONFIG * a,int hstart,int width,int vstart,int height)235 uint64_t aom_highbd_get_y_var(const YV12_BUFFER_CONFIG *a, int hstart,
236 int width, int vstart, int height) {
237 return aom_var_2d_u16(a->y_buffer + vstart * a->y_stride + hstart,
238 a->y_stride, width, height) /
239 (width * height);
240 }
241
aom_highbd_get_u_var(const YV12_BUFFER_CONFIG * a,int hstart,int width,int vstart,int height)242 uint64_t aom_highbd_get_u_var(const YV12_BUFFER_CONFIG *a, int hstart,
243 int width, int vstart, int height) {
244 return aom_var_2d_u16(a->u_buffer + vstart * a->uv_stride + hstart,
245 a->uv_stride, width, height) /
246 (width * height);
247 }
248
aom_highbd_get_v_var(const YV12_BUFFER_CONFIG * a,int hstart,int width,int vstart,int height)249 uint64_t aom_highbd_get_v_var(const YV12_BUFFER_CONFIG *a, int hstart,
250 int width, int vstart, int height) {
251 return aom_var_2d_u16(a->v_buffer + vstart * a->uv_stride + hstart,
252 a->uv_stride, width, height) /
253 (width * height);
254 }
255
aom_highbd_get_y_sse_part(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,int hstart,int width,int vstart,int height)256 int64_t aom_highbd_get_y_sse_part(const YV12_BUFFER_CONFIG *a,
257 const YV12_BUFFER_CONFIG *b, int hstart,
258 int width, int vstart, int height) {
259 return highbd_get_sse(
260 a->y_buffer + vstart * a->y_stride + hstart, a->y_stride,
261 b->y_buffer + vstart * b->y_stride + hstart, b->y_stride, width, height);
262 }
263
aom_highbd_get_y_sse(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b)264 int64_t aom_highbd_get_y_sse(const YV12_BUFFER_CONFIG *a,
265 const YV12_BUFFER_CONFIG *b) {
266 assert(a->y_crop_width == b->y_crop_width);
267 assert(a->y_crop_height == b->y_crop_height);
268 assert((a->flags & YV12_FLAG_HIGHBITDEPTH) != 0);
269 assert((b->flags & YV12_FLAG_HIGHBITDEPTH) != 0);
270
271 return highbd_get_sse(a->y_buffer, a->y_stride, b->y_buffer, b->y_stride,
272 a->y_crop_width, a->y_crop_height);
273 }
274
aom_highbd_get_u_sse_part(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,int hstart,int width,int vstart,int height)275 int64_t aom_highbd_get_u_sse_part(const YV12_BUFFER_CONFIG *a,
276 const YV12_BUFFER_CONFIG *b, int hstart,
277 int width, int vstart, int height) {
278 return highbd_get_sse(a->u_buffer + vstart * a->uv_stride + hstart,
279 a->uv_stride,
280 b->u_buffer + vstart * b->uv_stride + hstart,
281 b->uv_stride, width, height);
282 }
283
aom_highbd_get_u_sse(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b)284 int64_t aom_highbd_get_u_sse(const YV12_BUFFER_CONFIG *a,
285 const YV12_BUFFER_CONFIG *b) {
286 assert(a->uv_crop_width == b->uv_crop_width);
287 assert(a->uv_crop_height == b->uv_crop_height);
288 assert((a->flags & YV12_FLAG_HIGHBITDEPTH) != 0);
289 assert((b->flags & YV12_FLAG_HIGHBITDEPTH) != 0);
290
291 return highbd_get_sse(a->u_buffer, a->uv_stride, b->u_buffer, b->uv_stride,
292 a->uv_crop_width, a->uv_crop_height);
293 }
294
aom_highbd_get_v_sse_part(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,int hstart,int width,int vstart,int height)295 int64_t aom_highbd_get_v_sse_part(const YV12_BUFFER_CONFIG *a,
296 const YV12_BUFFER_CONFIG *b, int hstart,
297 int width, int vstart, int height) {
298 return highbd_get_sse(a->v_buffer + vstart * a->uv_stride + hstart,
299 a->uv_stride,
300 b->v_buffer + vstart * b->uv_stride + hstart,
301 b->uv_stride, width, height);
302 }
303
aom_highbd_get_v_sse(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b)304 int64_t aom_highbd_get_v_sse(const YV12_BUFFER_CONFIG *a,
305 const YV12_BUFFER_CONFIG *b) {
306 assert(a->uv_crop_width == b->uv_crop_width);
307 assert(a->uv_crop_height == b->uv_crop_height);
308 assert((a->flags & YV12_FLAG_HIGHBITDEPTH) != 0);
309 assert((b->flags & YV12_FLAG_HIGHBITDEPTH) != 0);
310
311 return highbd_get_sse(a->v_buffer, a->uv_stride, b->v_buffer, b->uv_stride,
312 a->uv_crop_width, a->uv_crop_height);
313 }
314 #endif // CONFIG_AV1_HIGHBITDEPTH
315
aom_get_sse_plane(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,int plane,int highbd)316 int64_t aom_get_sse_plane(const YV12_BUFFER_CONFIG *a,
317 const YV12_BUFFER_CONFIG *b, int plane, int highbd) {
318 #if CONFIG_AV1_HIGHBITDEPTH
319 if (highbd) {
320 switch (plane) {
321 case 0: return aom_highbd_get_y_sse(a, b);
322 case 1: return aom_highbd_get_u_sse(a, b);
323 case 2: return aom_highbd_get_v_sse(a, b);
324 default: assert(plane >= 0 && plane <= 2); return 0;
325 }
326 } else {
327 switch (plane) {
328 case 0: return aom_get_y_sse(a, b);
329 case 1: return aom_get_u_sse(a, b);
330 case 2: return aom_get_v_sse(a, b);
331 default: assert(plane >= 0 && plane <= 2); return 0;
332 }
333 }
334 #else
335 (void)highbd;
336 switch (plane) {
337 case 0: return aom_get_y_sse(a, b);
338 case 1: return aom_get_u_sse(a, b);
339 case 2: return aom_get_v_sse(a, b);
340 default: assert(plane >= 0 && plane <= 2); return 0;
341 }
342 #endif
343 }
344
345 #if CONFIG_AV1_HIGHBITDEPTH
aom_calc_highbd_psnr(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,PSNR_STATS * psnr,uint32_t bit_depth,uint32_t in_bit_depth)346 void aom_calc_highbd_psnr(const YV12_BUFFER_CONFIG *a,
347 const YV12_BUFFER_CONFIG *b, PSNR_STATS *psnr,
348 uint32_t bit_depth, uint32_t in_bit_depth) {
349 assert(a->y_crop_width == b->y_crop_width);
350 assert(a->y_crop_height == b->y_crop_height);
351 assert(a->uv_crop_width == b->uv_crop_width);
352 assert(a->uv_crop_height == b->uv_crop_height);
353 const int widths[3] = { a->y_crop_width, a->uv_crop_width, a->uv_crop_width };
354 const int heights[3] = { a->y_crop_height, a->uv_crop_height,
355 a->uv_crop_height };
356 const int a_strides[3] = { a->y_stride, a->uv_stride, a->uv_stride };
357 const int b_strides[3] = { b->y_stride, b->uv_stride, b->uv_stride };
358 int i;
359 uint64_t total_sse = 0;
360 uint32_t total_samples = 0;
361 #if CONFIG_LIBVMAF_PSNR_PEAK
362 double peak = (double)(255 << (in_bit_depth - 8));
363 #else
364 double peak = (double)((1 << in_bit_depth) - 1);
365 #endif // CONFIG_LIBVMAF_PSNR_PEAK
366 const unsigned int input_shift = bit_depth - in_bit_depth;
367
368 for (i = 0; i < 3; ++i) {
369 const int w = widths[i];
370 const int h = heights[i];
371 const uint32_t samples = w * h;
372 uint64_t sse;
373 if (a->flags & YV12_FLAG_HIGHBITDEPTH) {
374 if (input_shift) {
375 sse = highbd_get_sse_shift(a->buffers[i], a_strides[i], b->buffers[i],
376 b_strides[i], w, h, input_shift);
377 } else {
378 sse = highbd_get_sse(a->buffers[i], a_strides[i], b->buffers[i],
379 b_strides[i], w, h);
380 }
381 } else {
382 sse = get_sse(a->buffers[i], a_strides[i], b->buffers[i], b_strides[i], w,
383 h);
384 }
385 psnr->sse[1 + i] = sse;
386 psnr->samples[1 + i] = samples;
387 psnr->psnr[1 + i] = aom_sse_to_psnr(samples, peak, (double)sse);
388
389 total_sse += sse;
390 total_samples += samples;
391 }
392
393 psnr->sse[0] = total_sse;
394 psnr->samples[0] = total_samples;
395 psnr->psnr[0] =
396 aom_sse_to_psnr((double)total_samples, peak, (double)total_sse);
397
398 // Compute PSNR based on stream bit depth
399 if ((a->flags & YV12_FLAG_HIGHBITDEPTH) && (in_bit_depth < bit_depth)) {
400 #if CONFIG_LIBVMAF_PSNR_PEAK
401 peak = (double)(255 << (bit_depth - 8));
402 #else
403 peak = (double)((1 << bit_depth) - 1);
404 #endif // CONFIG_LIBVMAF_PSNR_PEAK
405 total_sse = 0;
406 total_samples = 0;
407 for (i = 0; i < 3; ++i) {
408 const int w = widths[i];
409 const int h = heights[i];
410 const uint32_t samples = w * h;
411 uint64_t sse;
412 sse = highbd_get_sse(a->buffers[i], a_strides[i], b->buffers[i],
413 b_strides[i], w, h);
414 psnr->sse_hbd[1 + i] = sse;
415 psnr->samples_hbd[1 + i] = samples;
416 psnr->psnr_hbd[1 + i] = aom_sse_to_psnr(samples, peak, (double)sse);
417 total_sse += sse;
418 total_samples += samples;
419 }
420
421 psnr->sse_hbd[0] = total_sse;
422 psnr->samples_hbd[0] = total_samples;
423 psnr->psnr_hbd[0] =
424 aom_sse_to_psnr((double)total_samples, peak, (double)total_sse);
425 }
426 }
427 #endif
428
aom_calc_psnr(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,PSNR_STATS * psnr)429 void aom_calc_psnr(const YV12_BUFFER_CONFIG *a, const YV12_BUFFER_CONFIG *b,
430 PSNR_STATS *psnr) {
431 assert(a->y_crop_width == b->y_crop_width);
432 assert(a->y_crop_height == b->y_crop_height);
433 assert(a->uv_crop_width == b->uv_crop_width);
434 assert(a->uv_crop_height == b->uv_crop_height);
435 static const double peak = 255.0;
436 const int widths[3] = { a->y_crop_width, a->uv_crop_width, a->uv_crop_width };
437 const int heights[3] = { a->y_crop_height, a->uv_crop_height,
438 a->uv_crop_height };
439 const int a_strides[3] = { a->y_stride, a->uv_stride, a->uv_stride };
440 const int b_strides[3] = { b->y_stride, b->uv_stride, b->uv_stride };
441 int i;
442 uint64_t total_sse = 0;
443 uint32_t total_samples = 0;
444
445 for (i = 0; i < 3; ++i) {
446 const int w = widths[i];
447 const int h = heights[i];
448 const uint32_t samples = w * h;
449 const uint64_t sse =
450 get_sse(a->buffers[i], a_strides[i], b->buffers[i], b_strides[i], w, h);
451 psnr->sse[1 + i] = sse;
452 psnr->samples[1 + i] = samples;
453 psnr->psnr[1 + i] = aom_sse_to_psnr(samples, peak, (double)sse);
454
455 total_sse += sse;
456 total_samples += samples;
457 }
458
459 psnr->sse[0] = total_sse;
460 psnr->samples[0] = total_samples;
461 psnr->psnr[0] =
462 aom_sse_to_psnr((double)total_samples, peak, (double)total_sse);
463 }
464