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 | |
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 | |
55 | double VP8SSIMFromStats(const VP8DistoStats* const stats) { |
56 | return SSIMCalculation(stats, kWeightSum); |
57 | } |
58 | |
59 | double VP8SSIMFromStatsClipped(const VP8DistoStats* const stats) { |
60 | return SSIMCalculation(stats, stats->w); |
61 | } |
62 | |
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 | |
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) |
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 | |
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_USE_SSE2) |
154 | if (VP8GetCPUInfo(kSSE2)) { |
155 | VP8SSIMDspInitSSE2(); |
156 | } |
157 | #endif |
158 | } |
159 | } |
160 | |