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