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.
25 static const uint32_t kWeight[2 * VP8_SSIM_KERNEL + 1] = {
26 1, 2, 3, 4, 3, 2, 1
27 };
28 static const uint32_t kWeightSum = 16 * 16; // sum{kWeight}^2
29
SSIMCalculation(const VP8DistoStats * const stats,uint32_t N)30 static 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
VP8SSIMFromStats(const VP8DistoStats * const stats)55 double VP8SSIMFromStats(const VP8DistoStats* const stats) {
56 return SSIMCalculation(stats, kWeightSum);
57 }
58
VP8SSIMFromStatsClipped(const VP8DistoStats * const stats)59 double VP8SSIMFromStatsClipped(const VP8DistoStats* const stats) {
60 return SSIMCalculation(stats, stats->w);
61 }
62
SSIMGetClipped_C(const uint8_t * src1,int stride1,const uint8_t * src2,int stride2,int xo,int yo,int W,int H)63 static 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
SSIMGet_C(const uint8_t * src1,int stride1,const uint8_t * src2,int stride2)93 static 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)
AccumulateSSE_C(const uint8_t * src1,const uint8_t * src2,int len)117 static 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)
133 VP8SSIMGetFunc VP8SSIMGet;
134 VP8SSIMGetClippedFunc VP8SSIMGetClipped;
135 #endif
136 #if !defined(WEBP_DISABLE_STATS)
137 VP8AccumulateSSEFunc VP8AccumulateSSE;
138 #endif
139
140 extern void VP8SSIMDspInitSSE2(void);
141
WEBP_DSP_INIT_FUNC(VP8SSIMDspInit)142 WEBP_DSP_INIT_FUNC(VP8SSIMDspInit) {
143 #if !defined(WEBP_REDUCE_SIZE)
144 VP8SSIMGetClipped = SSIMGetClipped_C;
145 VP8SSIMGet = SSIMGet_C;
146 #endif
147
148 #if !defined(WEBP_DISABLE_STATS)
149 VP8AccumulateSSE = AccumulateSSE_C;
150 #endif
151
152 if (VP8GetCPUInfo != NULL) {
153 #if defined(WEBP_HAVE_SSE2)
154 if (VP8GetCPUInfo(kSSE2)) {
155 VP8SSIMDspInitSSE2();
156 }
157 #endif
158 }
159 }
160