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 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 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 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 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 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 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 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 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 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 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 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 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