1// Copyright 2017 Google Inc. All Rights Reserved. 2// 3// Use of this source code is governed by a BSD-style license 4// that can be found in the COPYING file in the root of the source 5// tree. An additional intellectual property rights grant can be found 6// in the file PATENTS. All contributing project authors may 7// be found in the AUTHORS file in the root of the source tree. 8// ----------------------------------------------------------------------------- 9// 10// distortion calculation 11// 12// Author: Skal (pascal.massimino@gmail.com) 13 14#include <assert.h> 15#include <stdlib.h> // for abs() 16 17#include "src/dsp/dsp.h" 18 19#if !defined(WEBP_REDUCE_SIZE) 20 21//------------------------------------------------------------------------------ 22// SSIM / PSNR 23 24// hat-shaped filter. Sum of coefficients is equal to 16. 25static const uint32_t kWeight[2 * VP8_SSIM_KERNEL + 1] = { 26 1, 2, 3, 4, 3, 2, 1 27}; 28static const uint32_t kWeightSum = 16 * 16; // sum{kWeight}^2 29 30static WEBP_INLINE double SSIMCalculation( 31 const VP8DistoStats* const stats, uint32_t N /*num samples*/) { 32 const uint32_t w2 = N * N; 33 const uint32_t C1 = 20 * w2; 34 const uint32_t C2 = 60 * w2; 35 const uint32_t C3 = 8 * 8 * w2; // 'dark' limit ~= 6 36 const uint64_t xmxm = (uint64_t)stats->xm * stats->xm; 37 const uint64_t ymym = (uint64_t)stats->ym * stats->ym; 38 if (xmxm + ymym >= C3) { 39 const int64_t xmym = (int64_t)stats->xm * stats->ym; 40 const int64_t sxy = (int64_t)stats->xym * N - xmym; // can be negative 41 const uint64_t sxx = (uint64_t)stats->xxm * N - xmxm; 42 const uint64_t syy = (uint64_t)stats->yym * N - ymym; 43 // we descale by 8 to prevent overflow during the fnum/fden multiply. 44 const uint64_t num_S = (2 * (uint64_t)(sxy < 0 ? 0 : sxy) + C2) >> 8; 45 const uint64_t den_S = (sxx + syy + C2) >> 8; 46 const uint64_t fnum = (2 * xmym + C1) * num_S; 47 const uint64_t fden = (xmxm + ymym + C1) * den_S; 48 const double r = (double)fnum / fden; 49 assert(r >= 0. && r <= 1.0); 50 return r; 51 } 52 return 1.; // area is too dark to contribute meaningfully 53} 54 55double VP8SSIMFromStats(const VP8DistoStats* const stats) { 56 return SSIMCalculation(stats, kWeightSum); 57} 58 59double VP8SSIMFromStatsClipped(const VP8DistoStats* const stats) { 60 return SSIMCalculation(stats, stats->w); 61} 62 63static double SSIMGetClipped_C(const uint8_t* src1, int stride1, 64 const uint8_t* src2, int stride2, 65 int xo, int yo, int W, int H) { 66 VP8DistoStats stats = { 0, 0, 0, 0, 0, 0 }; 67 const int ymin = (yo - VP8_SSIM_KERNEL < 0) ? 0 : yo - VP8_SSIM_KERNEL; 68 const int ymax = (yo + VP8_SSIM_KERNEL > H - 1) ? H - 1 69 : yo + VP8_SSIM_KERNEL; 70 const int xmin = (xo - VP8_SSIM_KERNEL < 0) ? 0 : xo - VP8_SSIM_KERNEL; 71 const int xmax = (xo + VP8_SSIM_KERNEL > W - 1) ? W - 1 72 : xo + VP8_SSIM_KERNEL; 73 int x, y; 74 src1 += ymin * stride1; 75 src2 += ymin * stride2; 76 for (y = ymin; y <= ymax; ++y, src1 += stride1, src2 += stride2) { 77 for (x = xmin; x <= xmax; ++x) { 78 const uint32_t w = kWeight[VP8_SSIM_KERNEL + x - xo] 79 * kWeight[VP8_SSIM_KERNEL + y - yo]; 80 const uint32_t s1 = src1[x]; 81 const uint32_t s2 = src2[x]; 82 stats.w += w; 83 stats.xm += w * s1; 84 stats.ym += w * s2; 85 stats.xxm += w * s1 * s1; 86 stats.xym += w * s1 * s2; 87 stats.yym += w * s2 * s2; 88 } 89 } 90 return VP8SSIMFromStatsClipped(&stats); 91} 92 93static double SSIMGet_C(const uint8_t* src1, int stride1, 94 const uint8_t* src2, int stride2) { 95 VP8DistoStats stats = { 0, 0, 0, 0, 0, 0 }; 96 int x, y; 97 for (y = 0; y <= 2 * VP8_SSIM_KERNEL; ++y, src1 += stride1, src2 += stride2) { 98 for (x = 0; x <= 2 * VP8_SSIM_KERNEL; ++x) { 99 const uint32_t w = kWeight[x] * kWeight[y]; 100 const uint32_t s1 = src1[x]; 101 const uint32_t s2 = src2[x]; 102 stats.xm += w * s1; 103 stats.ym += w * s2; 104 stats.xxm += w * s1 * s1; 105 stats.xym += w * s1 * s2; 106 stats.yym += w * s2 * s2; 107 } 108 } 109 return VP8SSIMFromStats(&stats); 110} 111 112#endif // !defined(WEBP_REDUCE_SIZE) 113 114//------------------------------------------------------------------------------ 115 116#if !defined(WEBP_DISABLE_STATS) 117static uint32_t AccumulateSSE_C(const uint8_t* src1, 118 const uint8_t* src2, int len) { 119 int i; 120 uint32_t sse2 = 0; 121 assert(len <= 65535); // to ensure that accumulation fits within uint32_t 122 for (i = 0; i < len; ++i) { 123 const int32_t diff = src1[i] - src2[i]; 124 sse2 += diff * diff; 125 } 126 return sse2; 127} 128#endif 129 130//------------------------------------------------------------------------------ 131 132#if !defined(WEBP_REDUCE_SIZE) 133VP8SSIMGetFunc VP8SSIMGet; 134VP8SSIMGetClippedFunc VP8SSIMGetClipped; 135#endif 136#if !defined(WEBP_DISABLE_STATS) 137VP8AccumulateSSEFunc VP8AccumulateSSE; 138#endif 139 140extern void VP8SSIMDspInitSSE2(void); 141 142static volatile VP8CPUInfo ssim_last_cpuinfo_used = 143 (VP8CPUInfo)&ssim_last_cpuinfo_used; 144 145WEBP_TSAN_IGNORE_FUNCTION void VP8SSIMDspInit(void) { 146 if (ssim_last_cpuinfo_used == VP8GetCPUInfo) return; 147 148#if !defined(WEBP_REDUCE_SIZE) 149 VP8SSIMGetClipped = SSIMGetClipped_C; 150 VP8SSIMGet = SSIMGet_C; 151#endif 152 153#if !defined(WEBP_DISABLE_STATS) 154 VP8AccumulateSSE = AccumulateSSE_C; 155#endif 156 157 if (VP8GetCPUInfo != NULL) { 158#if defined(WEBP_USE_SSE2) 159 if (VP8GetCPUInfo(kSSE2)) { 160 VP8SSIMDspInitSSE2(); 161 } 162#endif 163 } 164 165 ssim_last_cpuinfo_used = VP8GetCPUInfo; 166} 167