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 
29 /* TODO(yaowu): The block_variance calls the unoptimized versions of variance()
30  * and highbd_8_variance(). It should not.
31  */
encoder_variance(const uint8_t * a,int a_stride,const uint8_t * b,int b_stride,int w,int h,unsigned int * sse,int * sum)32 static void encoder_variance(const uint8_t *a, int a_stride, const uint8_t *b,
33                              int b_stride, int w, int h, unsigned int *sse,
34                              int *sum) {
35   int i, j;
36 
37   *sum = 0;
38   *sse = 0;
39 
40   for (i = 0; i < h; i++) {
41     for (j = 0; j < w; j++) {
42       const int diff = a[j] - b[j];
43       *sum += diff;
44       *sse += diff * diff;
45     }
46 
47     a += a_stride;
48     b += b_stride;
49   }
50 }
51 
encoder_highbd_variance64(const uint8_t * a8,int a_stride,const uint8_t * b8,int b_stride,int w,int h,uint64_t * sse,int64_t * sum)52 static void encoder_highbd_variance64(const uint8_t *a8, int a_stride,
53                                       const uint8_t *b8, int b_stride, int w,
54                                       int h, uint64_t *sse, int64_t *sum) {
55   const uint16_t *a = CONVERT_TO_SHORTPTR(a8);
56   const uint16_t *b = CONVERT_TO_SHORTPTR(b8);
57   int64_t tsum = 0;
58   uint64_t tsse = 0;
59   for (int i = 0; i < h; ++i) {
60     int32_t lsum = 0;
61     for (int j = 0; j < w; ++j) {
62       const int diff = a[j] - b[j];
63       lsum += diff;
64       tsse += (uint32_t)(diff * diff);
65     }
66     tsum += lsum;
67     a += a_stride;
68     b += b_stride;
69   }
70   *sum = tsum;
71   *sse = tsse;
72 }
73 
encoder_highbd_8_variance(const uint8_t * a8,int a_stride,const uint8_t * b8,int b_stride,int w,int h,unsigned int * sse,int * sum)74 static void encoder_highbd_8_variance(const uint8_t *a8, int a_stride,
75                                       const uint8_t *b8, int b_stride, int w,
76                                       int h, unsigned int *sse, int *sum) {
77   uint64_t sse_long = 0;
78   int64_t sum_long = 0;
79   encoder_highbd_variance64(a8, a_stride, b8, b_stride, w, h, &sse_long,
80                             &sum_long);
81   *sse = (unsigned int)sse_long;
82   *sum = (int)sum_long;
83 }
84 
get_sse(const uint8_t * a,int a_stride,const uint8_t * b,int b_stride,int width,int height)85 static int64_t get_sse(const uint8_t *a, int a_stride, const uint8_t *b,
86                        int b_stride, int width, int height) {
87   const int dw = width % 16;
88   const int dh = height % 16;
89   int64_t total_sse = 0;
90   unsigned int sse = 0;
91   int sum = 0;
92   int x, y;
93 
94   if (dw > 0) {
95     encoder_variance(&a[width - dw], a_stride, &b[width - dw], b_stride, dw,
96                      height, &sse, &sum);
97     total_sse += sse;
98   }
99 
100   if (dh > 0) {
101     encoder_variance(&a[(height - dh) * a_stride], a_stride,
102                      &b[(height - dh) * b_stride], b_stride, width - dw, dh,
103                      &sse, &sum);
104     total_sse += sse;
105   }
106 
107   for (y = 0; y < height / 16; ++y) {
108     const uint8_t *pa = a;
109     const uint8_t *pb = b;
110     for (x = 0; x < width / 16; ++x) {
111       aom_mse16x16(pa, a_stride, pb, b_stride, &sse);
112       total_sse += sse;
113 
114       pa += 16;
115       pb += 16;
116     }
117 
118     a += 16 * a_stride;
119     b += 16 * b_stride;
120   }
121 
122   return total_sse;
123 }
124 
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)125 static int64_t highbd_get_sse_shift(const uint8_t *a8, int a_stride,
126                                     const uint8_t *b8, int b_stride, int width,
127                                     int height, unsigned int input_shift) {
128   const uint16_t *a = CONVERT_TO_SHORTPTR(a8);
129   const uint16_t *b = CONVERT_TO_SHORTPTR(b8);
130   int64_t total_sse = 0;
131   int x, y;
132   for (y = 0; y < height; ++y) {
133     for (x = 0; x < width; ++x) {
134       int64_t diff;
135       diff = (a[x] >> input_shift) - (b[x] >> input_shift);
136       total_sse += diff * diff;
137     }
138     a += a_stride;
139     b += b_stride;
140   }
141   return total_sse;
142 }
143 
highbd_get_sse(const uint8_t * a,int a_stride,const uint8_t * b,int b_stride,int width,int height)144 static int64_t highbd_get_sse(const uint8_t *a, int a_stride, const uint8_t *b,
145                               int b_stride, int width, int height) {
146   int64_t total_sse = 0;
147   int x, y;
148   const int dw = width % 16;
149   const int dh = height % 16;
150   unsigned int sse = 0;
151   int sum = 0;
152   if (dw > 0) {
153     encoder_highbd_8_variance(&a[width - dw], a_stride, &b[width - dw],
154                               b_stride, dw, height, &sse, &sum);
155     total_sse += sse;
156   }
157   if (dh > 0) {
158     encoder_highbd_8_variance(&a[(height - dh) * a_stride], a_stride,
159                               &b[(height - dh) * b_stride], b_stride,
160                               width - dw, dh, &sse, &sum);
161     total_sse += sse;
162   }
163   for (y = 0; y < height / 16; ++y) {
164     const uint8_t *pa = a;
165     const uint8_t *pb = b;
166     for (x = 0; x < width / 16; ++x) {
167       aom_highbd_8_mse16x16(pa, a_stride, pb, b_stride, &sse);
168       total_sse += sse;
169       pa += 16;
170       pb += 16;
171     }
172     a += 16 * a_stride;
173     b += 16 * b_stride;
174   }
175   return total_sse;
176 }
177 
aom_get_y_sse_part(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,int hstart,int width,int vstart,int height)178 int64_t aom_get_y_sse_part(const YV12_BUFFER_CONFIG *a,
179                            const YV12_BUFFER_CONFIG *b, int hstart, int width,
180                            int vstart, int height) {
181   return get_sse(a->y_buffer + vstart * a->y_stride + hstart, a->y_stride,
182                  b->y_buffer + vstart * b->y_stride + hstart, b->y_stride,
183                  width, height);
184 }
185 
aom_get_y_sse(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b)186 int64_t aom_get_y_sse(const YV12_BUFFER_CONFIG *a,
187                       const YV12_BUFFER_CONFIG *b) {
188   assert(a->y_crop_width == b->y_crop_width);
189   assert(a->y_crop_height == b->y_crop_height);
190 
191   return get_sse(a->y_buffer, a->y_stride, b->y_buffer, b->y_stride,
192                  a->y_crop_width, a->y_crop_height);
193 }
194 
aom_get_u_sse_part(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,int hstart,int width,int vstart,int height)195 int64_t aom_get_u_sse_part(const YV12_BUFFER_CONFIG *a,
196                            const YV12_BUFFER_CONFIG *b, int hstart, int width,
197                            int vstart, int height) {
198   return get_sse(a->u_buffer + vstart * a->uv_stride + hstart, a->uv_stride,
199                  b->u_buffer + vstart * b->uv_stride + hstart, b->uv_stride,
200                  width, height);
201 }
202 
aom_get_u_sse(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b)203 int64_t aom_get_u_sse(const YV12_BUFFER_CONFIG *a,
204                       const YV12_BUFFER_CONFIG *b) {
205   assert(a->uv_crop_width == b->uv_crop_width);
206   assert(a->uv_crop_height == b->uv_crop_height);
207 
208   return get_sse(a->u_buffer, a->uv_stride, b->u_buffer, b->uv_stride,
209                  a->uv_crop_width, a->uv_crop_height);
210 }
211 
aom_get_v_sse_part(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,int hstart,int width,int vstart,int height)212 int64_t aom_get_v_sse_part(const YV12_BUFFER_CONFIG *a,
213                            const YV12_BUFFER_CONFIG *b, int hstart, int width,
214                            int vstart, int height) {
215   return get_sse(a->v_buffer + vstart * a->uv_stride + hstart, a->uv_stride,
216                  b->v_buffer + vstart * b->uv_stride + hstart, b->uv_stride,
217                  width, height);
218 }
219 
aom_get_v_sse(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b)220 int64_t aom_get_v_sse(const YV12_BUFFER_CONFIG *a,
221                       const YV12_BUFFER_CONFIG *b) {
222   assert(a->uv_crop_width == b->uv_crop_width);
223   assert(a->uv_crop_height == b->uv_crop_height);
224 
225   return get_sse(a->v_buffer, a->uv_stride, b->v_buffer, b->uv_stride,
226                  a->uv_crop_width, a->uv_crop_height);
227 }
228 
aom_highbd_get_y_sse_part(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,int hstart,int width,int vstart,int height)229 int64_t aom_highbd_get_y_sse_part(const YV12_BUFFER_CONFIG *a,
230                                   const YV12_BUFFER_CONFIG *b, int hstart,
231                                   int width, int vstart, int height) {
232   return highbd_get_sse(
233       a->y_buffer + vstart * a->y_stride + hstart, a->y_stride,
234       b->y_buffer + vstart * b->y_stride + hstart, b->y_stride, width, height);
235 }
236 
aom_highbd_get_y_sse(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b)237 int64_t aom_highbd_get_y_sse(const YV12_BUFFER_CONFIG *a,
238                              const YV12_BUFFER_CONFIG *b) {
239   assert(a->y_crop_width == b->y_crop_width);
240   assert(a->y_crop_height == b->y_crop_height);
241   assert((a->flags & YV12_FLAG_HIGHBITDEPTH) != 0);
242   assert((b->flags & YV12_FLAG_HIGHBITDEPTH) != 0);
243 
244   return highbd_get_sse(a->y_buffer, a->y_stride, b->y_buffer, b->y_stride,
245                         a->y_crop_width, a->y_crop_height);
246 }
247 
aom_highbd_get_u_sse_part(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,int hstart,int width,int vstart,int height)248 int64_t aom_highbd_get_u_sse_part(const YV12_BUFFER_CONFIG *a,
249                                   const YV12_BUFFER_CONFIG *b, int hstart,
250                                   int width, int vstart, int height) {
251   return highbd_get_sse(a->u_buffer + vstart * a->uv_stride + hstart,
252                         a->uv_stride,
253                         b->u_buffer + vstart * b->uv_stride + hstart,
254                         b->uv_stride, width, height);
255 }
256 
aom_highbd_get_u_sse(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b)257 int64_t aom_highbd_get_u_sse(const YV12_BUFFER_CONFIG *a,
258                              const YV12_BUFFER_CONFIG *b) {
259   assert(a->uv_crop_width == b->uv_crop_width);
260   assert(a->uv_crop_height == b->uv_crop_height);
261   assert((a->flags & YV12_FLAG_HIGHBITDEPTH) != 0);
262   assert((b->flags & YV12_FLAG_HIGHBITDEPTH) != 0);
263 
264   return highbd_get_sse(a->u_buffer, a->uv_stride, b->u_buffer, b->uv_stride,
265                         a->uv_crop_width, a->uv_crop_height);
266 }
267 
aom_highbd_get_v_sse_part(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,int hstart,int width,int vstart,int height)268 int64_t aom_highbd_get_v_sse_part(const YV12_BUFFER_CONFIG *a,
269                                   const YV12_BUFFER_CONFIG *b, int hstart,
270                                   int width, int vstart, int height) {
271   return highbd_get_sse(a->v_buffer + vstart * a->uv_stride + hstart,
272                         a->uv_stride,
273                         b->v_buffer + vstart * b->uv_stride + hstart,
274                         b->uv_stride, width, height);
275 }
276 
aom_highbd_get_v_sse(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b)277 int64_t aom_highbd_get_v_sse(const YV12_BUFFER_CONFIG *a,
278                              const YV12_BUFFER_CONFIG *b) {
279   assert(a->uv_crop_width == b->uv_crop_width);
280   assert(a->uv_crop_height == b->uv_crop_height);
281   assert((a->flags & YV12_FLAG_HIGHBITDEPTH) != 0);
282   assert((b->flags & YV12_FLAG_HIGHBITDEPTH) != 0);
283 
284   return highbd_get_sse(a->v_buffer, a->uv_stride, b->v_buffer, b->uv_stride,
285                         a->uv_crop_width, a->uv_crop_height);
286 }
287 
aom_get_sse_plane(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,int plane,int highbd)288 int64_t aom_get_sse_plane(const YV12_BUFFER_CONFIG *a,
289                           const YV12_BUFFER_CONFIG *b, int plane, int highbd) {
290   if (highbd) {
291     switch (plane) {
292       case 0: return aom_highbd_get_y_sse(a, b);
293       case 1: return aom_highbd_get_u_sse(a, b);
294       case 2: return aom_highbd_get_v_sse(a, b);
295       default: assert(plane >= 0 && plane <= 2); return 0;
296     }
297   }
298   switch (plane) {
299     case 0: return aom_get_y_sse(a, b);
300     case 1: return aom_get_u_sse(a, b);
301     case 2: return aom_get_v_sse(a, b);
302     default: assert(plane >= 0 && plane <= 2); return 0;
303   }
304 }
305 
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)306 void aom_calc_highbd_psnr(const YV12_BUFFER_CONFIG *a,
307                           const YV12_BUFFER_CONFIG *b, PSNR_STATS *psnr,
308                           uint32_t bit_depth, uint32_t in_bit_depth) {
309   const int widths[3] = { a->y_crop_width, a->uv_crop_width, a->uv_crop_width };
310   const int heights[3] = { a->y_crop_height, a->uv_crop_height,
311                            a->uv_crop_height };
312   const int a_strides[3] = { a->y_stride, a->uv_stride, a->uv_stride };
313   const int b_strides[3] = { b->y_stride, b->uv_stride, b->uv_stride };
314   int i;
315   uint64_t total_sse = 0;
316   uint32_t total_samples = 0;
317   const double peak = (double)((1 << in_bit_depth) - 1);
318   const unsigned int input_shift = bit_depth - in_bit_depth;
319 
320   for (i = 0; i < 3; ++i) {
321     const int w = widths[i];
322     const int h = heights[i];
323     const uint32_t samples = w * h;
324     uint64_t sse;
325     if (a->flags & YV12_FLAG_HIGHBITDEPTH) {
326       if (input_shift) {
327         sse = highbd_get_sse_shift(a->buffers[i], a_strides[i], b->buffers[i],
328                                    b_strides[i], w, h, input_shift);
329       } else {
330         sse = highbd_get_sse(a->buffers[i], a_strides[i], b->buffers[i],
331                              b_strides[i], w, h);
332       }
333     } else {
334       sse = get_sse(a->buffers[i], a_strides[i], b->buffers[i], b_strides[i], w,
335                     h);
336     }
337     psnr->sse[1 + i] = sse;
338     psnr->samples[1 + i] = samples;
339     psnr->psnr[1 + i] = aom_sse_to_psnr(samples, peak, (double)sse);
340 
341     total_sse += sse;
342     total_samples += samples;
343   }
344 
345   psnr->sse[0] = total_sse;
346   psnr->samples[0] = total_samples;
347   psnr->psnr[0] =
348       aom_sse_to_psnr((double)total_samples, peak, (double)total_sse);
349 }
350 
aom_calc_psnr(const YV12_BUFFER_CONFIG * a,const YV12_BUFFER_CONFIG * b,PSNR_STATS * psnr)351 void aom_calc_psnr(const YV12_BUFFER_CONFIG *a, const YV12_BUFFER_CONFIG *b,
352                    PSNR_STATS *psnr) {
353   static const double peak = 255.0;
354   const int widths[3] = { a->y_crop_width, a->uv_crop_width, a->uv_crop_width };
355   const int heights[3] = { a->y_crop_height, a->uv_crop_height,
356                            a->uv_crop_height };
357   const int a_strides[3] = { a->y_stride, a->uv_stride, a->uv_stride };
358   const int b_strides[3] = { b->y_stride, b->uv_stride, b->uv_stride };
359   int i;
360   uint64_t total_sse = 0;
361   uint32_t total_samples = 0;
362 
363   for (i = 0; i < 3; ++i) {
364     const int w = widths[i];
365     const int h = heights[i];
366     const uint32_t samples = w * h;
367     const uint64_t sse =
368         get_sse(a->buffers[i], a_strides[i], b->buffers[i], b_strides[i], w, h);
369     psnr->sse[1 + i] = sse;
370     psnr->samples[1 + i] = samples;
371     psnr->psnr[1 + i] = aom_sse_to_psnr(samples, peak, (double)sse);
372 
373     total_sse += sse;
374     total_samples += samples;
375   }
376 
377   psnr->sse[0] = total_sse;
378   psnr->samples[0] = total_samples;
379   psnr->psnr[0] =
380       aom_sse_to_psnr((double)total_samples, peak, (double)total_sse);
381 }
382