1 /*
2 * Copyright 2011 The LibYuv Project Authors. All rights reserved.
3 *
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.
9 */
10
11 #include "libyuv/compare.h"
12
13 #include <float.h>
14 #include <math.h>
15 #ifdef _OPENMP
16 #include <omp.h>
17 #endif
18
19 #include "libyuv/basic_types.h"
20 #include "libyuv/compare_row.h"
21 #include "libyuv/cpu_id.h"
22 #include "libyuv/row.h"
23 #include "libyuv/video_common.h"
24
25 #ifdef __cplusplus
26 namespace libyuv {
27 extern "C" {
28 #endif
29
30 // hash seed of 5381 recommended.
31 LIBYUV_API
HashDjb2(const uint8_t * src,uint64_t count,uint32_t seed)32 uint32_t HashDjb2(const uint8_t* src, uint64_t count, uint32_t seed) {
33 const int kBlockSize = 1 << 15; // 32768;
34 int remainder;
35 uint32_t (*HashDjb2_SSE)(const uint8_t* src, int count, uint32_t seed) =
36 HashDjb2_C;
37 #if defined(HAS_HASHDJB2_SSE41)
38 if (TestCpuFlag(kCpuHasSSE41)) {
39 HashDjb2_SSE = HashDjb2_SSE41;
40 }
41 #endif
42 #if defined(HAS_HASHDJB2_AVX2)
43 if (TestCpuFlag(kCpuHasAVX2)) {
44 HashDjb2_SSE = HashDjb2_AVX2;
45 }
46 #endif
47
48 while (count >= (uint64_t)(kBlockSize)) {
49 seed = HashDjb2_SSE(src, kBlockSize, seed);
50 src += kBlockSize;
51 count -= kBlockSize;
52 }
53 remainder = (int)count & ~15;
54 if (remainder) {
55 seed = HashDjb2_SSE(src, remainder, seed);
56 src += remainder;
57 count -= remainder;
58 }
59 remainder = (int)count & 15;
60 if (remainder) {
61 seed = HashDjb2_C(src, remainder, seed);
62 }
63 return seed;
64 }
65
ARGBDetectRow_C(const uint8_t * argb,int width)66 static uint32_t ARGBDetectRow_C(const uint8_t* argb, int width) {
67 int x;
68 for (x = 0; x < width - 1; x += 2) {
69 if (argb[0] != 255) { // First byte is not Alpha of 255, so not ARGB.
70 return FOURCC_BGRA;
71 }
72 if (argb[3] != 255) { // 4th byte is not Alpha of 255, so not BGRA.
73 return FOURCC_ARGB;
74 }
75 if (argb[4] != 255) { // Second pixel first byte is not Alpha of 255.
76 return FOURCC_BGRA;
77 }
78 if (argb[7] != 255) { // Second pixel 4th byte is not Alpha of 255.
79 return FOURCC_ARGB;
80 }
81 argb += 8;
82 }
83 if (width & 1) {
84 if (argb[0] != 255) { // First byte is not Alpha of 255, so not ARGB.
85 return FOURCC_BGRA;
86 }
87 if (argb[3] != 255) { // 4th byte is not Alpha of 255, so not BGRA.
88 return FOURCC_ARGB;
89 }
90 }
91 return 0;
92 }
93
94 // Scan an opaque argb image and return fourcc based on alpha offset.
95 // Returns FOURCC_ARGB, FOURCC_BGRA, or 0 if unknown.
96 LIBYUV_API
ARGBDetect(const uint8_t * argb,int stride_argb,int width,int height)97 uint32_t ARGBDetect(const uint8_t* argb,
98 int stride_argb,
99 int width,
100 int height) {
101 uint32_t fourcc = 0;
102 int h;
103
104 // Coalesce rows.
105 if (stride_argb == width * 4) {
106 width *= height;
107 height = 1;
108 stride_argb = 0;
109 }
110 for (h = 0; h < height && fourcc == 0; ++h) {
111 fourcc = ARGBDetectRow_C(argb, width);
112 argb += stride_argb;
113 }
114 return fourcc;
115 }
116
117 // NEON version accumulates in 16 bit shorts which overflow at 65536 bytes.
118 // So actual maximum is 1 less loop, which is 64436 - 32 bytes.
119
120 LIBYUV_API
ComputeHammingDistance(const uint8_t * src_a,const uint8_t * src_b,int count)121 uint64_t ComputeHammingDistance(const uint8_t* src_a,
122 const uint8_t* src_b,
123 int count) {
124 const int kBlockSize = 1 << 15; // 32768;
125 const int kSimdSize = 64;
126 // SIMD for multiple of 64, and C for remainder
127 int remainder = count & (kBlockSize - 1) & ~(kSimdSize - 1);
128 uint64_t diff = 0;
129 int i;
130 uint32_t (*HammingDistance)(const uint8_t* src_a, const uint8_t* src_b,
131 int count) = HammingDistance_C;
132 #if defined(HAS_HAMMINGDISTANCE_NEON)
133 if (TestCpuFlag(kCpuHasNEON)) {
134 HammingDistance = HammingDistance_NEON;
135 }
136 #endif
137 #if defined(HAS_HAMMINGDISTANCE_SSSE3)
138 if (TestCpuFlag(kCpuHasSSSE3)) {
139 HammingDistance = HammingDistance_SSSE3;
140 }
141 #endif
142 #if defined(HAS_HAMMINGDISTANCE_SSE42)
143 if (TestCpuFlag(kCpuHasSSE42)) {
144 HammingDistance = HammingDistance_SSE42;
145 }
146 #endif
147 #if defined(HAS_HAMMINGDISTANCE_AVX2)
148 if (TestCpuFlag(kCpuHasAVX2)) {
149 HammingDistance = HammingDistance_AVX2;
150 }
151 #endif
152 #if defined(HAS_HAMMINGDISTANCE_MSA)
153 if (TestCpuFlag(kCpuHasMSA)) {
154 HammingDistance = HammingDistance_MSA;
155 }
156 #endif
157 #ifdef _OPENMP
158 #pragma omp parallel for reduction(+ : diff)
159 #endif
160 for (i = 0; i < (count - (kBlockSize - 1)); i += kBlockSize) {
161 diff += HammingDistance(src_a + i, src_b + i, kBlockSize);
162 }
163 src_a += count & ~(kBlockSize - 1);
164 src_b += count & ~(kBlockSize - 1);
165 if (remainder) {
166 diff += HammingDistance(src_a, src_b, remainder);
167 src_a += remainder;
168 src_b += remainder;
169 }
170 remainder = count & (kSimdSize - 1);
171 if (remainder) {
172 diff += HammingDistance_C(src_a, src_b, remainder);
173 }
174 return diff;
175 }
176
177 // TODO(fbarchard): Refactor into row function.
178 LIBYUV_API
ComputeSumSquareError(const uint8_t * src_a,const uint8_t * src_b,int count)179 uint64_t ComputeSumSquareError(const uint8_t* src_a,
180 const uint8_t* src_b,
181 int count) {
182 // SumSquareError returns values 0 to 65535 for each squared difference.
183 // Up to 65536 of those can be summed and remain within a uint32_t.
184 // After each block of 65536 pixels, accumulate into a uint64_t.
185 const int kBlockSize = 65536;
186 int remainder = count & (kBlockSize - 1) & ~31;
187 uint64_t sse = 0;
188 int i;
189 uint32_t (*SumSquareError)(const uint8_t* src_a, const uint8_t* src_b,
190 int count) = SumSquareError_C;
191 #if defined(HAS_SUMSQUAREERROR_NEON)
192 if (TestCpuFlag(kCpuHasNEON)) {
193 SumSquareError = SumSquareError_NEON;
194 }
195 #endif
196 #if defined(HAS_SUMSQUAREERROR_SSE2)
197 if (TestCpuFlag(kCpuHasSSE2)) {
198 // Note only used for multiples of 16 so count is not checked.
199 SumSquareError = SumSquareError_SSE2;
200 }
201 #endif
202 #if defined(HAS_SUMSQUAREERROR_AVX2)
203 if (TestCpuFlag(kCpuHasAVX2)) {
204 // Note only used for multiples of 32 so count is not checked.
205 SumSquareError = SumSquareError_AVX2;
206 }
207 #endif
208 #if defined(HAS_SUMSQUAREERROR_MSA)
209 if (TestCpuFlag(kCpuHasMSA)) {
210 SumSquareError = SumSquareError_MSA;
211 }
212 #endif
213 #ifdef _OPENMP
214 #pragma omp parallel for reduction(+ : sse)
215 #endif
216 for (i = 0; i < (count - (kBlockSize - 1)); i += kBlockSize) {
217 sse += SumSquareError(src_a + i, src_b + i, kBlockSize);
218 }
219 src_a += count & ~(kBlockSize - 1);
220 src_b += count & ~(kBlockSize - 1);
221 if (remainder) {
222 sse += SumSquareError(src_a, src_b, remainder);
223 src_a += remainder;
224 src_b += remainder;
225 }
226 remainder = count & 31;
227 if (remainder) {
228 sse += SumSquareError_C(src_a, src_b, remainder);
229 }
230 return sse;
231 }
232
233 LIBYUV_API
ComputeSumSquareErrorPlane(const uint8_t * src_a,int stride_a,const uint8_t * src_b,int stride_b,int width,int height)234 uint64_t ComputeSumSquareErrorPlane(const uint8_t* src_a,
235 int stride_a,
236 const uint8_t* src_b,
237 int stride_b,
238 int width,
239 int height) {
240 uint64_t sse = 0;
241 int h;
242 // Coalesce rows.
243 if (stride_a == width && stride_b == width) {
244 width *= height;
245 height = 1;
246 stride_a = stride_b = 0;
247 }
248 for (h = 0; h < height; ++h) {
249 sse += ComputeSumSquareError(src_a, src_b, width);
250 src_a += stride_a;
251 src_b += stride_b;
252 }
253 return sse;
254 }
255
256 LIBYUV_API
SumSquareErrorToPsnr(uint64_t sse,uint64_t count)257 double SumSquareErrorToPsnr(uint64_t sse, uint64_t count) {
258 double psnr;
259 if (sse > 0) {
260 double mse = (double)count / (double)sse;
261 psnr = 10.0 * log10(255.0 * 255.0 * mse);
262 } else {
263 psnr = kMaxPsnr; // Limit to prevent divide by 0
264 }
265
266 if (psnr > kMaxPsnr) {
267 psnr = kMaxPsnr;
268 }
269
270 return psnr;
271 }
272
273 LIBYUV_API
CalcFramePsnr(const uint8_t * src_a,int stride_a,const uint8_t * src_b,int stride_b,int width,int height)274 double CalcFramePsnr(const uint8_t* src_a,
275 int stride_a,
276 const uint8_t* src_b,
277 int stride_b,
278 int width,
279 int height) {
280 const uint64_t samples = (uint64_t)width * (uint64_t)height;
281 const uint64_t sse = ComputeSumSquareErrorPlane(src_a, stride_a, src_b,
282 stride_b, width, height);
283 return SumSquareErrorToPsnr(sse, samples);
284 }
285
286 LIBYUV_API
I420Psnr(const uint8_t * src_y_a,int stride_y_a,const uint8_t * src_u_a,int stride_u_a,const uint8_t * src_v_a,int stride_v_a,const uint8_t * src_y_b,int stride_y_b,const uint8_t * src_u_b,int stride_u_b,const uint8_t * src_v_b,int stride_v_b,int width,int height)287 double I420Psnr(const uint8_t* src_y_a,
288 int stride_y_a,
289 const uint8_t* src_u_a,
290 int stride_u_a,
291 const uint8_t* src_v_a,
292 int stride_v_a,
293 const uint8_t* src_y_b,
294 int stride_y_b,
295 const uint8_t* src_u_b,
296 int stride_u_b,
297 const uint8_t* src_v_b,
298 int stride_v_b,
299 int width,
300 int height) {
301 const uint64_t sse_y = ComputeSumSquareErrorPlane(
302 src_y_a, stride_y_a, src_y_b, stride_y_b, width, height);
303 const int width_uv = (width + 1) >> 1;
304 const int height_uv = (height + 1) >> 1;
305 const uint64_t sse_u = ComputeSumSquareErrorPlane(
306 src_u_a, stride_u_a, src_u_b, stride_u_b, width_uv, height_uv);
307 const uint64_t sse_v = ComputeSumSquareErrorPlane(
308 src_v_a, stride_v_a, src_v_b, stride_v_b, width_uv, height_uv);
309 const uint64_t samples = (uint64_t)width * (uint64_t)height +
310 2 * ((uint64_t)width_uv * (uint64_t)height_uv);
311 const uint64_t sse = sse_y + sse_u + sse_v;
312 return SumSquareErrorToPsnr(sse, samples);
313 }
314
315 static const int64_t cc1 = 26634; // (64^2*(.01*255)^2
316 static const int64_t cc2 = 239708; // (64^2*(.03*255)^2
317
Ssim8x8_C(const uint8_t * src_a,int stride_a,const uint8_t * src_b,int stride_b)318 static double Ssim8x8_C(const uint8_t* src_a,
319 int stride_a,
320 const uint8_t* src_b,
321 int stride_b) {
322 int64_t sum_a = 0;
323 int64_t sum_b = 0;
324 int64_t sum_sq_a = 0;
325 int64_t sum_sq_b = 0;
326 int64_t sum_axb = 0;
327
328 int i;
329 for (i = 0; i < 8; ++i) {
330 int j;
331 for (j = 0; j < 8; ++j) {
332 sum_a += src_a[j];
333 sum_b += src_b[j];
334 sum_sq_a += src_a[j] * src_a[j];
335 sum_sq_b += src_b[j] * src_b[j];
336 sum_axb += src_a[j] * src_b[j];
337 }
338
339 src_a += stride_a;
340 src_b += stride_b;
341 }
342
343 {
344 const int64_t count = 64;
345 // scale the constants by number of pixels
346 const int64_t c1 = (cc1 * count * count) >> 12;
347 const int64_t c2 = (cc2 * count * count) >> 12;
348
349 const int64_t sum_a_x_sum_b = sum_a * sum_b;
350
351 const int64_t ssim_n = (2 * sum_a_x_sum_b + c1) *
352 (2 * count * sum_axb - 2 * sum_a_x_sum_b + c2);
353
354 const int64_t sum_a_sq = sum_a * sum_a;
355 const int64_t sum_b_sq = sum_b * sum_b;
356
357 const int64_t ssim_d =
358 (sum_a_sq + sum_b_sq + c1) *
359 (count * sum_sq_a - sum_a_sq + count * sum_sq_b - sum_b_sq + c2);
360
361 if (ssim_d == 0.0) {
362 return DBL_MAX;
363 }
364 return ssim_n * 1.0 / ssim_d;
365 }
366 }
367
368 // We are using a 8x8 moving window with starting location of each 8x8 window
369 // on the 4x4 pixel grid. Such arrangement allows the windows to overlap
370 // block boundaries to penalize blocking artifacts.
371 LIBYUV_API
CalcFrameSsim(const uint8_t * src_a,int stride_a,const uint8_t * src_b,int stride_b,int width,int height)372 double CalcFrameSsim(const uint8_t* src_a,
373 int stride_a,
374 const uint8_t* src_b,
375 int stride_b,
376 int width,
377 int height) {
378 int samples = 0;
379 double ssim_total = 0;
380 double (*Ssim8x8)(const uint8_t* src_a, int stride_a, const uint8_t* src_b,
381 int stride_b) = Ssim8x8_C;
382
383 // sample point start with each 4x4 location
384 int i;
385 for (i = 0; i < height - 8; i += 4) {
386 int j;
387 for (j = 0; j < width - 8; j += 4) {
388 ssim_total += Ssim8x8(src_a + j, stride_a, src_b + j, stride_b);
389 samples++;
390 }
391
392 src_a += stride_a * 4;
393 src_b += stride_b * 4;
394 }
395
396 ssim_total /= samples;
397 return ssim_total;
398 }
399
400 LIBYUV_API
I420Ssim(const uint8_t * src_y_a,int stride_y_a,const uint8_t * src_u_a,int stride_u_a,const uint8_t * src_v_a,int stride_v_a,const uint8_t * src_y_b,int stride_y_b,const uint8_t * src_u_b,int stride_u_b,const uint8_t * src_v_b,int stride_v_b,int width,int height)401 double I420Ssim(const uint8_t* src_y_a,
402 int stride_y_a,
403 const uint8_t* src_u_a,
404 int stride_u_a,
405 const uint8_t* src_v_a,
406 int stride_v_a,
407 const uint8_t* src_y_b,
408 int stride_y_b,
409 const uint8_t* src_u_b,
410 int stride_u_b,
411 const uint8_t* src_v_b,
412 int stride_v_b,
413 int width,
414 int height) {
415 const double ssim_y =
416 CalcFrameSsim(src_y_a, stride_y_a, src_y_b, stride_y_b, width, height);
417 const int width_uv = (width + 1) >> 1;
418 const int height_uv = (height + 1) >> 1;
419 const double ssim_u = CalcFrameSsim(src_u_a, stride_u_a, src_u_b, stride_u_b,
420 width_uv, height_uv);
421 const double ssim_v = CalcFrameSsim(src_v_a, stride_v_a, src_v_b, stride_v_b,
422 width_uv, height_uv);
423 return ssim_y * 0.8 + 0.1 * (ssim_u + ssim_v);
424 }
425
426 #ifdef __cplusplus
427 } // extern "C"
428 } // namespace libyuv
429 #endif
430