1 // Copyright 2014 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 // WebPPicture tools for measuring distortion
11 //
12 // Author: Skal (pascal.massimino@gmail.com)
13 
14 #include "src/webp/encode.h"
15 
16 #if !(defined(WEBP_DISABLE_STATS) || defined(WEBP_REDUCE_SIZE))
17 
18 #include <math.h>
19 #include <stdlib.h>
20 
21 #include "src/dsp/dsp.h"
22 #include "src/enc/vp8i_enc.h"
23 #include "src/utils/utils.h"
24 
25 typedef double (*AccumulateFunc)(const uint8_t* src, int src_stride,
26                                  const uint8_t* ref, int ref_stride,
27                                  int w, int h);
28 
29 //------------------------------------------------------------------------------
30 // local-min distortion
31 //
32 // For every pixel in the *reference* picture, we search for the local best
33 // match in the compressed image. This is not a symmetrical measure.
34 
35 #define RADIUS 2  // search radius. Shouldn't be too large.
36 
37 static double AccumulateLSIM(const uint8_t* src, int src_stride,
38                              const uint8_t* ref, int ref_stride,
39                              int w, int h) {
40   int x, y;
41   double total_sse = 0.;
42   for (y = 0; y < h; ++y) {
43     const int y_0 = (y - RADIUS < 0) ? 0 : y - RADIUS;
44     const int y_1 = (y + RADIUS + 1 >= h) ? h : y + RADIUS + 1;
45     for (x = 0; x < w; ++x) {
46       const int x_0 = (x - RADIUS < 0) ? 0 : x - RADIUS;
47       const int x_1 = (x + RADIUS + 1 >= w) ? w : x + RADIUS + 1;
48       double best_sse = 255. * 255.;
49       const double value = (double)ref[y * ref_stride + x];
50       int i, j;
51       for (j = y_0; j < y_1; ++j) {
52         const uint8_t* const s = src + j * src_stride;
53         for (i = x_0; i < x_1; ++i) {
54           const double diff = s[i] - value;
55           const double sse = diff * diff;
56           if (sse < best_sse) best_sse = sse;
57         }
58       }
59       total_sse += best_sse;
60     }
61   }
62   return total_sse;
63 }
64 #undef RADIUS
65 
66 static double AccumulateSSE(const uint8_t* src, int src_stride,
67                             const uint8_t* ref, int ref_stride,
68                             int w, int h) {
69   int y;
70   double total_sse = 0.;
71   for (y = 0; y < h; ++y) {
72     total_sse += VP8AccumulateSSE(src, ref, w);
73     src += src_stride;
74     ref += ref_stride;
75   }
76   return total_sse;
77 }
78 
79 //------------------------------------------------------------------------------
80 
81 static double AccumulateSSIM(const uint8_t* src, int src_stride,
82                              const uint8_t* ref, int ref_stride,
83                              int w, int h) {
84   const int w0 = (w < VP8_SSIM_KERNEL) ? w : VP8_SSIM_KERNEL;
85   const int w1 = w - VP8_SSIM_KERNEL - 1;
86   const int h0 = (h < VP8_SSIM_KERNEL) ? h : VP8_SSIM_KERNEL;
87   const int h1 = h - VP8_SSIM_KERNEL - 1;
88   int x, y;
89   double sum = 0.;
90   for (y = 0; y < h0; ++y) {
91     for (x = 0; x < w; ++x) {
92       sum += VP8SSIMGetClipped(src, src_stride, ref, ref_stride, x, y, w, h);
93     }
94   }
95   for (; y < h1; ++y) {
96     for (x = 0; x < w0; ++x) {
97       sum += VP8SSIMGetClipped(src, src_stride, ref, ref_stride, x, y, w, h);
98     }
99     for (; x < w1; ++x) {
100       const int off1 = x - VP8_SSIM_KERNEL + (y - VP8_SSIM_KERNEL) * src_stride;
101       const int off2 = x - VP8_SSIM_KERNEL + (y - VP8_SSIM_KERNEL) * ref_stride;
102       sum += VP8SSIMGet(src + off1, src_stride, ref + off2, ref_stride);
103     }
104     for (; x < w; ++x) {
105       sum += VP8SSIMGetClipped(src, src_stride, ref, ref_stride, x, y, w, h);
106     }
107   }
108   for (; y < h; ++y) {
109     for (x = 0; x < w; ++x) {
110       sum += VP8SSIMGetClipped(src, src_stride, ref, ref_stride, x, y, w, h);
111     }
112   }
113   return sum;
114 }
115 
116 //------------------------------------------------------------------------------
117 // Distortion
118 
119 // Max value returned in case of exact similarity.
120 static const double kMinDistortion_dB = 99.;
121 
122 static double GetPSNR(double v, double size) {
123   return (v > 0. && size > 0.) ? -4.3429448 * log(v / (size * 255 * 255.))
124                                : kMinDistortion_dB;
125 }
126 
127 static double GetLogSSIM(double v, double size) {
128   v = (size > 0.) ? v / size : 1.;
129   return (v < 1.) ? -10.0 * log10(1. - v) : kMinDistortion_dB;
130 }
131 
132 int WebPPlaneDistortion(const uint8_t* src, size_t src_stride,
133                         const uint8_t* ref, size_t ref_stride,
134                         int width, int height, size_t x_step,
135                         int type, float* distortion, float* result) {
136   uint8_t* allocated = NULL;
137   const AccumulateFunc metric = (type == 0) ? AccumulateSSE :
138                                 (type == 1) ? AccumulateSSIM :
139                                               AccumulateLSIM;
140   if (src == NULL || ref == NULL ||
141       src_stride < x_step * width || ref_stride < x_step * width ||
142       result == NULL || distortion == NULL) {
143     return 0;
144   }
145 
146   VP8SSIMDspInit();
147   if (x_step != 1) {   // extract a packed plane if needed
148     int x, y;
149     uint8_t* tmp1;
150     uint8_t* tmp2;
151     allocated =
152         (uint8_t*)WebPSafeMalloc(2ULL * width * height, sizeof(*allocated));
153     if (allocated == NULL) return 0;
154     tmp1 = allocated;
155     tmp2 = tmp1 + (size_t)width * height;
156     for (y = 0; y < height; ++y) {
157       for (x = 0; x < width; ++x) {
158         tmp1[x + y * width] = src[x * x_step + y * src_stride];
159         tmp2[x + y * width] = ref[x * x_step + y * ref_stride];
160       }
161     }
162     src = tmp1;
163     ref = tmp2;
164   }
165   *distortion = (float)metric(src, width, ref, width, width, height);
166   WebPSafeFree(allocated);
167 
168   *result = (type == 1) ? (float)GetLogSSIM(*distortion, (double)width * height)
169                         : (float)GetPSNR(*distortion, (double)width * height);
170   return 1;
171 }
172 
173 #ifdef WORDS_BIGENDIAN
174 #define BLUE_OFFSET 3   // uint32_t 0x000000ff is 0x00,00,00,ff in memory
175 #else
176 #define BLUE_OFFSET 0   // uint32_t 0x000000ff is 0xff,00,00,00 in memory
177 #endif
178 
179 int WebPPictureDistortion(const WebPPicture* src, const WebPPicture* ref,
180                           int type, float results[5]) {
181   int w, h, c;
182   int ok = 0;
183   WebPPicture p0, p1;
184   double total_size = 0., total_distortion = 0.;
185   if (src == NULL || ref == NULL ||
186       src->width != ref->width || src->height != ref->height ||
187       results == NULL) {
188     return 0;
189   }
190 
191   VP8SSIMDspInit();
192   if (!WebPPictureInit(&p0) || !WebPPictureInit(&p1)) return 0;
193   w = src->width;
194   h = src->height;
195   if (!WebPPictureView(src, 0, 0, w, h, &p0)) goto Error;
196   if (!WebPPictureView(ref, 0, 0, w, h, &p1)) goto Error;
197 
198   // We always measure distortion in ARGB space.
199   if (p0.use_argb == 0 && !WebPPictureYUVAToARGB(&p0)) goto Error;
200   if (p1.use_argb == 0 && !WebPPictureYUVAToARGB(&p1)) goto Error;
201   for (c = 0; c < 4; ++c) {
202     float distortion;
203     const size_t stride0 = 4 * (size_t)p0.argb_stride;
204     const size_t stride1 = 4 * (size_t)p1.argb_stride;
205     // results are reported as BGRA
206     const int offset = c ^ BLUE_OFFSET;
207     if (!WebPPlaneDistortion((const uint8_t*)p0.argb + offset, stride0,
208                              (const uint8_t*)p1.argb + offset, stride1,
209                              w, h, 4, type, &distortion, results + c)) {
210       goto Error;
211     }
212     total_distortion += distortion;
213     total_size += w * h;
214   }
215 
216   results[4] = (type == 1) ? (float)GetLogSSIM(total_distortion, total_size)
217                            : (float)GetPSNR(total_distortion, total_size);
218   ok = 1;
219 
220  Error:
221   WebPPictureFree(&p0);
222   WebPPictureFree(&p1);
223   return ok;
224 }
225 
226 #undef BLUE_OFFSET
227 
228 #else  // defined(WEBP_DISABLE_STATS)
229 int WebPPlaneDistortion(const uint8_t* src, size_t src_stride,
230                         const uint8_t* ref, size_t ref_stride,
231                         int width, int height, size_t x_step,
232                         int type, float* distortion, float* result) {
233   (void)src;
234   (void)src_stride;
235   (void)ref;
236   (void)ref_stride;
237   (void)width;
238   (void)height;
239   (void)x_step;
240   (void)type;
241   if (distortion == NULL || result == NULL) return 0;
242   *distortion = 0.f;
243   *result = 0.f;
244   return 1;
245 }
246 
247 int WebPPictureDistortion(const WebPPicture* src, const WebPPicture* ref,
248                           int type, float results[5]) {
249   int i;
250   (void)src;
251   (void)ref;
252   (void)type;
253   if (results == NULL) return 0;
254   for (i = 0; i < 5; ++i) results[i] = 0.f;
255   return 1;
256 }
257 
258 #endif  // !defined(WEBP_DISABLE_STATS)
259