summaryrefslogtreecommitdiffstats
path: root/media/libyuv/source/compare.cc
blob: e3846bdfdd11521325cc4953cdfe29738b91a4ed (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
/*
 *  Copyright 2011 The LibYuv Project Authors. All rights reserved.
 *
 *  Use of this source code is governed by a BSD-style license
 *  that can be found in the LICENSE file in the root of the source
 *  tree. An additional intellectual property rights grant can be found
 *  in the file PATENTS. All contributing project authors may
 *  be found in the AUTHORS file in the root of the source tree.
 */

#include "libyuv/compare.h"

#include <float.h>
#include <math.h>
#ifdef _OPENMP
#include <omp.h>
#endif

#include "libyuv/basic_types.h"
#include "libyuv/compare_row.h"
#include "libyuv/cpu_id.h"
#include "libyuv/row.h"
#include "libyuv/video_common.h"

#ifdef __cplusplus
namespace libyuv {
extern "C" {
#endif

// hash seed of 5381 recommended.
LIBYUV_API
uint32 HashDjb2(const uint8* src, uint64 count, uint32 seed) {
  const int kBlockSize = 1 << 15;  // 32768;
  int remainder;
  uint32 (*HashDjb2_SSE)(const uint8* src, int count, uint32 seed) =
      HashDjb2_C;
#if defined(HAS_HASHDJB2_SSE41)
  if (TestCpuFlag(kCpuHasSSE41)) {
    HashDjb2_SSE = HashDjb2_SSE41;
  }
#endif
#if defined(HAS_HASHDJB2_AVX2)
  if (TestCpuFlag(kCpuHasAVX2)) {
    HashDjb2_SSE = HashDjb2_AVX2;
  }
#endif

  while (count >= (uint64)(kBlockSize)) {
    seed = HashDjb2_SSE(src, kBlockSize, seed);
    src += kBlockSize;
    count -= kBlockSize;
  }
  remainder = (int)(count) & ~15;
  if (remainder) {
    seed = HashDjb2_SSE(src, remainder, seed);
    src += remainder;
    count -= remainder;
  }
  remainder = (int)(count) & 15;
  if (remainder) {
    seed = HashDjb2_C(src, remainder, seed);
  }
  return seed;
}

static uint32 ARGBDetectRow_C(const uint8* argb, int width) {
  int x;
  for (x = 0; x < width - 1; x += 2) {
    if (argb[0] != 255) {  // First byte is not Alpha of 255, so not ARGB.
      return FOURCC_BGRA;
    }
    if (argb[3] != 255) {  // 4th byte is not Alpha of 255, so not BGRA.
      return FOURCC_ARGB;
    }
    if (argb[4] != 255) {  // Second pixel first byte is not Alpha of 255.
      return FOURCC_BGRA;
    }
    if (argb[7] != 255) {  // Second pixel 4th byte is not Alpha of 255.
      return FOURCC_ARGB;
    }
    argb += 8;
  }
  if (width & 1) {
    if (argb[0] != 255) {  // First byte is not Alpha of 255, so not ARGB.
      return FOURCC_BGRA;
    }
    if (argb[3] != 255) {  // 4th byte is not Alpha of 255, so not BGRA.
      return FOURCC_ARGB;
    }
  }
  return 0;
}

// Scan an opaque argb image and return fourcc based on alpha offset.
// Returns FOURCC_ARGB, FOURCC_BGRA, or 0 if unknown.
LIBYUV_API
uint32 ARGBDetect(const uint8* argb, int stride_argb, int width, int height) {
  uint32 fourcc = 0;
  int h;

  // Coalesce rows.
  if (stride_argb == width * 4) {
    width *= height;
    height = 1;
    stride_argb = 0;
  }
  for (h = 0; h < height && fourcc == 0; ++h) {
    fourcc = ARGBDetectRow_C(argb, width);
    argb += stride_argb;
  }
  return fourcc;
}

// TODO(fbarchard): Refactor into row function.
LIBYUV_API
uint64 ComputeSumSquareError(const uint8* src_a, const uint8* src_b,
                             int count) {
  // SumSquareError returns values 0 to 65535 for each squared difference.
  // Up to 65536 of those can be summed and remain within a uint32.
  // After each block of 65536 pixels, accumulate into a uint64.
  const int kBlockSize = 65536;
  int remainder = count & (kBlockSize - 1) & ~31;
  uint64 sse = 0;
  int i;
  uint32 (*SumSquareError)(const uint8* src_a, const uint8* src_b, int count) =
      SumSquareError_C;
#if defined(HAS_SUMSQUAREERROR_NEON)
  if (TestCpuFlag(kCpuHasNEON)) {
    SumSquareError = SumSquareError_NEON;
  }
#endif
#if defined(HAS_SUMSQUAREERROR_SSE2)
  if (TestCpuFlag(kCpuHasSSE2)) {
    // Note only used for multiples of 16 so count is not checked.
    SumSquareError = SumSquareError_SSE2;
  }
#endif
#if defined(HAS_SUMSQUAREERROR_AVX2)
  if (TestCpuFlag(kCpuHasAVX2)) {
    // Note only used for multiples of 32 so count is not checked.
    SumSquareError = SumSquareError_AVX2;
  }
#endif
#ifdef _OPENMP
#pragma omp parallel for reduction(+: sse)
#endif
  for (i = 0; i < (count - (kBlockSize - 1)); i += kBlockSize) {
    sse += SumSquareError(src_a + i, src_b + i, kBlockSize);
  }
  src_a += count & ~(kBlockSize - 1);
  src_b += count & ~(kBlockSize - 1);
  if (remainder) {
    sse += SumSquareError(src_a, src_b, remainder);
    src_a += remainder;
    src_b += remainder;
  }
  remainder = count & 31;
  if (remainder) {
    sse += SumSquareError_C(src_a, src_b, remainder);
  }
  return sse;
}

LIBYUV_API
uint64 ComputeSumSquareErrorPlane(const uint8* src_a, int stride_a,
                                  const uint8* src_b, int stride_b,
                                  int width, int height) {
  uint64 sse = 0;
  int h;
  // Coalesce rows.
  if (stride_a == width &&
      stride_b == width) {
    width *= height;
    height = 1;
    stride_a = stride_b = 0;
  }
  for (h = 0; h < height; ++h) {
    sse += ComputeSumSquareError(src_a, src_b, width);
    src_a += stride_a;
    src_b += stride_b;
  }
  return sse;
}

LIBYUV_API
double SumSquareErrorToPsnr(uint64 sse, uint64 count) {
  double psnr;
  if (sse > 0) {
    double mse = (double)(count) / (double)(sse);
    psnr = 10.0 * log10(255.0 * 255.0 * mse);
  } else {
    psnr = kMaxPsnr;      // Limit to prevent divide by 0
  }

  if (psnr > kMaxPsnr)
    psnr = kMaxPsnr;

  return psnr;
}

LIBYUV_API
double CalcFramePsnr(const uint8* src_a, int stride_a,
                     const uint8* src_b, int stride_b,
                     int width, int height) {
  const uint64 samples = width * height;
  const uint64 sse = ComputeSumSquareErrorPlane(src_a, stride_a,
                                                src_b, stride_b,
                                                width, height);
  return SumSquareErrorToPsnr(sse, samples);
}

LIBYUV_API
double I420Psnr(const uint8* src_y_a, int stride_y_a,
                const uint8* src_u_a, int stride_u_a,
                const uint8* src_v_a, int stride_v_a,
                const uint8* src_y_b, int stride_y_b,
                const uint8* src_u_b, int stride_u_b,
                const uint8* src_v_b, int stride_v_b,
                int width, int height) {
  const uint64 sse_y = ComputeSumSquareErrorPlane(src_y_a, stride_y_a,
                                                  src_y_b, stride_y_b,
                                                  width, height);
  const int width_uv = (width + 1) >> 1;
  const int height_uv = (height + 1) >> 1;
  const uint64 sse_u = ComputeSumSquareErrorPlane(src_u_a, stride_u_a,
                                                  src_u_b, stride_u_b,
                                                  width_uv, height_uv);
  const uint64 sse_v = ComputeSumSquareErrorPlane(src_v_a, stride_v_a,
                                                  src_v_b, stride_v_b,
                                                  width_uv, height_uv);
  const uint64 samples = width * height + 2 * (width_uv * height_uv);
  const uint64 sse = sse_y + sse_u + sse_v;
  return SumSquareErrorToPsnr(sse, samples);
}

static const int64 cc1 =  26634;  // (64^2*(.01*255)^2
static const int64 cc2 = 239708;  // (64^2*(.03*255)^2

static double Ssim8x8_C(const uint8* src_a, int stride_a,
                        const uint8* src_b, int stride_b) {
  int64 sum_a = 0;
  int64 sum_b = 0;
  int64 sum_sq_a = 0;
  int64 sum_sq_b = 0;
  int64 sum_axb = 0;

  int i;
  for (i = 0; i < 8; ++i) {
    int j;
    for (j = 0; j < 8; ++j) {
      sum_a += src_a[j];
      sum_b += src_b[j];
      sum_sq_a += src_a[j] * src_a[j];
      sum_sq_b += src_b[j] * src_b[j];
      sum_axb += src_a[j] * src_b[j];
    }

    src_a += stride_a;
    src_b += stride_b;
  }

  {
    const int64 count = 64;
    // scale the constants by number of pixels
    const int64 c1 = (cc1 * count * count) >> 12;
    const int64 c2 = (cc2 * count * count) >> 12;

    const int64 sum_a_x_sum_b = sum_a * sum_b;

    const int64 ssim_n = (2 * sum_a_x_sum_b + c1) *
                         (2 * count * sum_axb - 2 * sum_a_x_sum_b + c2);

    const int64 sum_a_sq = sum_a*sum_a;
    const int64 sum_b_sq = sum_b*sum_b;

    const int64 ssim_d = (sum_a_sq + sum_b_sq + c1) *
                         (count * sum_sq_a - sum_a_sq +
                          count * sum_sq_b - sum_b_sq + c2);

    if (ssim_d == 0.0) {
      return DBL_MAX;
    }
    return ssim_n * 1.0 / ssim_d;
  }
}

// We are using a 8x8 moving window with starting location of each 8x8 window
// on the 4x4 pixel grid. Such arrangement allows the windows to overlap
// block boundaries to penalize blocking artifacts.
LIBYUV_API
double CalcFrameSsim(const uint8* src_a, int stride_a,
                     const uint8* src_b, int stride_b,
                     int width, int height) {
  int samples = 0;
  double ssim_total = 0;
  double (*Ssim8x8)(const uint8* src_a, int stride_a,
                    const uint8* src_b, int stride_b) = Ssim8x8_C;

  // sample point start with each 4x4 location
  int i;
  for (i = 0; i < height - 8; i += 4) {
    int j;
    for (j = 0; j < width - 8; j += 4) {
      ssim_total += Ssim8x8(src_a + j, stride_a, src_b + j, stride_b);
      samples++;
    }

    src_a += stride_a * 4;
    src_b += stride_b * 4;
  }

  ssim_total /= samples;
  return ssim_total;
}

LIBYUV_API
double I420Ssim(const uint8* src_y_a, int stride_y_a,
                const uint8* src_u_a, int stride_u_a,
                const uint8* src_v_a, int stride_v_a,
                const uint8* src_y_b, int stride_y_b,
                const uint8* src_u_b, int stride_u_b,
                const uint8* src_v_b, int stride_v_b,
                int width, int height) {
  const double ssim_y = CalcFrameSsim(src_y_a, stride_y_a,
                                      src_y_b, stride_y_b, width, height);
  const int width_uv = (width + 1) >> 1;
  const int height_uv = (height + 1) >> 1;
  const double ssim_u = CalcFrameSsim(src_u_a, stride_u_a,
                                      src_u_b, stride_u_b,
                                      width_uv, height_uv);
  const double ssim_v = CalcFrameSsim(src_v_a, stride_v_a,
                                      src_v_b, stride_v_b,
                                      width_uv, height_uv);
  return ssim_y * 0.8 + 0.1 * (ssim_u + ssim_v);
}

#ifdef __cplusplus
}  // extern "C"
}  // namespace libyuv
#endif