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