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