1 // Copyright 2016 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 // Simple tool to load two webp/png/jpg/tiff files and compute PSNR/SSIM.
11 // This is mostly a wrapper around WebPPictureDistortion().
12 //
13 /*
14 gcc -o get_disto get_disto.c -O3 -I../ -L../examples -L../imageio \
15 -lexample_util -limageio_util -limagedec -lwebp -L/opt/local/lib \
16 -lpng -lz -ljpeg -ltiff -lm -lpthread
17 */
18 //
19 // Author: Skal (pascal.massimino@gmail.com)
20
21 #include <assert.h>
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <string.h>
25
26 #include "webp/encode.h"
27 #include "imageio/image_dec.h"
28 #include "imageio/imageio_util.h"
29
ReadPicture(const char * const filename,WebPPicture * const pic,int keep_alpha)30 static size_t ReadPicture(const char* const filename, WebPPicture* const pic,
31 int keep_alpha) {
32 const uint8_t* data = NULL;
33 size_t data_size = 0;
34 WebPImageReader reader = NULL;
35 int ok = ImgIoUtilReadFile(filename, &data, &data_size);
36 if (!ok) goto End;
37
38 pic->use_argb = 1; // force ARGB
39
40 #ifdef HAVE_WINCODEC_H
41 // Try to decode the file using WIC falling back to the other readers for
42 // e.g., WebP.
43 ok = ReadPictureWithWIC(filename, pic, keep_alpha, NULL);
44 if (ok) goto End;
45 #endif
46 reader = WebPGuessImageReader(data, data_size);
47 ok = reader(data, data_size, pic, keep_alpha, NULL);
48
49 End:
50 if (!ok) {
51 fprintf(stderr, "Error! Could not process file %s\n", filename);
52 }
53 free((void*)data);
54 return ok ? data_size : 0;
55 }
56
RescalePlane(uint8_t * plane,int width,int height,int x_stride,int y_stride,int max)57 static void RescalePlane(uint8_t* plane, int width, int height,
58 int x_stride, int y_stride, int max) {
59 const uint32_t factor = (max > 0) ? (255u << 16) / max : 0;
60 int x, y;
61 for (y = 0; y < height; ++y) {
62 uint8_t* const ptr = plane + y * y_stride;
63 for (x = 0; x < width * x_stride; x += x_stride) {
64 const uint32_t diff = (ptr[x] * factor + (1 << 15)) >> 16;
65 ptr[x] = diff;
66 }
67 }
68 }
69
70 // Return the max absolute difference.
DiffScaleChannel(uint8_t * src1,int stride1,const uint8_t * src2,int stride2,int x_stride,int w,int h,int do_scaling)71 static int DiffScaleChannel(uint8_t* src1, int stride1,
72 const uint8_t* src2, int stride2,
73 int x_stride, int w, int h, int do_scaling) {
74 int x, y;
75 int max = 0;
76 for (y = 0; y < h; ++y) {
77 uint8_t* const ptr1 = src1 + y * stride1;
78 const uint8_t* const ptr2 = src2 + y * stride2;
79 for (x = 0; x < w * x_stride; x += x_stride) {
80 const int diff = abs(ptr1[x] - ptr2[x]);
81 if (diff > max) max = diff;
82 ptr1[x] = diff;
83 }
84 }
85
86 if (do_scaling) RescalePlane(src1, w, h, x_stride, stride1, max);
87 return max;
88 }
89
90 //------------------------------------------------------------------------------
91 // SSIM calculation. We re-implement these functions here, out of dsp/, to avoid
92 // breaking the library's hidden visibility. This code duplication avoids the
93 // bigger annoyance of having to open up internal details of libdsp...
94
95 #define SSIM_KERNEL 3 // total size of the kernel: 2 * SSIM_KERNEL + 1
96
97 // struct for accumulating statistical moments
98 typedef struct {
99 uint32_t w; // sum(w_i) : sum of weights
100 uint32_t xm, ym; // sum(w_i * x_i), sum(w_i * y_i)
101 uint32_t xxm, xym, yym; // sum(w_i * x_i * x_i), etc.
102 } DistoStats;
103
104 // hat-shaped filter. Sum of coefficients is equal to 16.
105 static const uint32_t kWeight[2 * SSIM_KERNEL + 1] = { 1, 2, 3, 4, 3, 2, 1 };
106
SSIMCalculation(const DistoStats * const stats)107 static WEBP_INLINE double SSIMCalculation(const DistoStats* const stats) {
108 const uint32_t N = stats->w;
109 const uint32_t w2 = N * N;
110 const uint32_t C1 = 20 * w2;
111 const uint32_t C2 = 60 * w2;
112 const uint32_t C3 = 8 * 8 * w2; // 'dark' limit ~= 6
113 const uint64_t xmxm = (uint64_t)stats->xm * stats->xm;
114 const uint64_t ymym = (uint64_t)stats->ym * stats->ym;
115 if (xmxm + ymym >= C3) {
116 const int64_t xmym = (int64_t)stats->xm * stats->ym;
117 const int64_t sxy = (int64_t)stats->xym * N - xmym; // can be negative
118 const uint64_t sxx = (uint64_t)stats->xxm * N - xmxm;
119 const uint64_t syy = (uint64_t)stats->yym * N - ymym;
120 // we descale by 8 to prevent overflow during the fnum/fden multiply.
121 const uint64_t num_S = (2 * (uint64_t)(sxy < 0 ? 0 : sxy) + C2) >> 8;
122 const uint64_t den_S = (sxx + syy + C2) >> 8;
123 const uint64_t fnum = (2 * xmym + C1) * num_S;
124 const uint64_t fden = (xmxm + ymym + C1) * den_S;
125 const double r = (double)fnum / fden;
126 assert(r >= 0. && r <= 1.0);
127 return r;
128 }
129 return 1.; // area is too dark to contribute meaningfully
130 }
131
SSIMGetClipped(const uint8_t * src1,int stride1,const uint8_t * src2,int stride2,int xo,int yo,int W,int H)132 static double SSIMGetClipped(const uint8_t* src1, int stride1,
133 const uint8_t* src2, int stride2,
134 int xo, int yo, int W, int H) {
135 DistoStats stats = { 0, 0, 0, 0, 0, 0 };
136 const int ymin = (yo - SSIM_KERNEL < 0) ? 0 : yo - SSIM_KERNEL;
137 const int ymax = (yo + SSIM_KERNEL > H - 1) ? H - 1 : yo + SSIM_KERNEL;
138 const int xmin = (xo - SSIM_KERNEL < 0) ? 0 : xo - SSIM_KERNEL;
139 const int xmax = (xo + SSIM_KERNEL > W - 1) ? W - 1 : xo + SSIM_KERNEL;
140 int x, y;
141 src1 += ymin * stride1;
142 src2 += ymin * stride2;
143 for (y = ymin; y <= ymax; ++y, src1 += stride1, src2 += stride2) {
144 for (x = xmin; x <= xmax; ++x) {
145 const uint32_t w = kWeight[SSIM_KERNEL + x - xo]
146 * kWeight[SSIM_KERNEL + y - yo];
147 const uint32_t s1 = src1[x];
148 const uint32_t s2 = src2[x];
149 stats.w += w;
150 stats.xm += w * s1;
151 stats.ym += w * s2;
152 stats.xxm += w * s1 * s1;
153 stats.xym += w * s1 * s2;
154 stats.yym += w * s2 * s2;
155 }
156 }
157 return SSIMCalculation(&stats);
158 }
159
160 // Compute SSIM-score map. Return -1 in case of error, max diff otherwise.
SSIMScaleChannel(uint8_t * src1,int stride1,const uint8_t * src2,int stride2,int x_stride,int w,int h,int do_scaling)161 static int SSIMScaleChannel(uint8_t* src1, int stride1,
162 const uint8_t* src2, int stride2,
163 int x_stride, int w, int h, int do_scaling) {
164 int x, y;
165 int max = 0;
166 uint8_t* const plane1 = (uint8_t*)malloc(2 * w * h * sizeof(*plane1));
167 uint8_t* const plane2 = plane1 + w * h;
168 if (plane1 == NULL) return -1;
169
170 // extract plane
171 for (y = 0; y < h; ++y) {
172 for (x = 0; x < w; ++x) {
173 plane1[x + y * w] = src1[x * x_stride + y * stride1];
174 plane2[x + y * w] = src2[x * x_stride + y * stride2];
175 }
176 }
177 for (y = 0; y < h; ++y) {
178 for (x = 0; x < w; ++x) {
179 const double ssim = SSIMGetClipped(plane1, w, plane2, w, x, y, w, h);
180 int diff = (int)(255 * (1. - ssim));
181 if (diff < 0) {
182 diff = 0;
183 } else if (diff > max) {
184 max = diff;
185 }
186 src1[x * x_stride + y * stride1] = (diff > 255) ? 255u : (uint8_t)diff;
187 }
188 }
189 free(plane1);
190
191 if (do_scaling) RescalePlane(src1, w, h, x_stride, stride1, max);
192 return max;
193 }
194
195 // Convert an argb picture to luminance.
ConvertToGray(WebPPicture * const pic)196 static void ConvertToGray(WebPPicture* const pic) {
197 int x, y;
198 assert(pic != NULL);
199 assert(pic->use_argb);
200 for (y = 0; y < pic->height; ++y) {
201 uint32_t* const row = &pic->argb[y * pic->argb_stride];
202 for (x = 0; x < pic->width; ++x) {
203 const uint32_t argb = row[x];
204 const uint32_t r = (argb >> 16) & 0xff;
205 const uint32_t g = (argb >> 8) & 0xff;
206 const uint32_t b = (argb >> 0) & 0xff;
207 // We use BT.709 for converting to luminance.
208 const uint32_t Y = (uint32_t)(0.2126 * r + 0.7152 * g + 0.0722 * b + .5);
209 row[x] = (argb & 0xff000000u) | (Y * 0x010101u);
210 }
211 }
212 }
213
Help(void)214 static void Help(void) {
215 fprintf(stderr,
216 "Usage: get_disto [-ssim][-psnr][-alpha] compressed.webp orig.webp\n"
217 " -ssim ..... print SSIM distortion\n"
218 " -psnr ..... print PSNR distortion (default)\n"
219 " -alpha .... preserve alpha plane\n"
220 " -h ........ this message\n"
221 " -o <file> . save the diff map as a WebP lossless file\n"
222 " -scale .... scale the difference map to fit [0..255] range\n"
223 " -gray ..... use grayscale for difference map (-scale)\n"
224 " Also handles PNG, JPG and TIFF files, in addition to WebP.\n");
225 }
226
main(int argc,const char * argv[])227 int main(int argc, const char *argv[]) {
228 WebPPicture pic1, pic2;
229 size_t size1 = 0, size2 = 0;
230 int ret = 1;
231 float disto[5];
232 int type = 0;
233 int c;
234 int help = 0;
235 int keep_alpha = 0;
236 int scale = 0;
237 int use_gray = 0;
238 const char* name1 = NULL;
239 const char* name2 = NULL;
240 const char* output = NULL;
241
242 if (!WebPPictureInit(&pic1) || !WebPPictureInit(&pic2)) {
243 fprintf(stderr, "Can't init pictures\n");
244 return 1;
245 }
246
247 for (c = 1; c < argc; ++c) {
248 if (!strcmp(argv[c], "-ssim")) {
249 type = 1;
250 } else if (!strcmp(argv[c], "-psnr")) {
251 type = 0;
252 } else if (!strcmp(argv[c], "-alpha")) {
253 keep_alpha = 1;
254 } else if (!strcmp(argv[c], "-scale")) {
255 scale = 1;
256 } else if (!strcmp(argv[c], "-gray")) {
257 use_gray = 1;
258 } else if (!strcmp(argv[c], "-h")) {
259 help = 1;
260 ret = 0;
261 } else if (!strcmp(argv[c], "-o")) {
262 if (++c == argc) {
263 fprintf(stderr, "missing file name after %s option.\n", argv[c - 1]);
264 goto End;
265 }
266 output = argv[c];
267 } else if (name1 == NULL) {
268 name1 = argv[c];
269 } else {
270 name2 = argv[c];
271 }
272 }
273 if (help || name1 == NULL || name2 == NULL) {
274 if (!help) {
275 fprintf(stderr, "Error: missing arguments.\n");
276 }
277 Help();
278 goto End;
279 }
280 size1 = ReadPicture(name1, &pic1, 1);
281 size2 = ReadPicture(name2, &pic2, 1);
282 if (size1 == 0 || size2 == 0) goto End;
283
284 if (!keep_alpha) {
285 WebPBlendAlpha(&pic1, 0x00000000);
286 WebPBlendAlpha(&pic2, 0x00000000);
287 }
288
289 if (!WebPPictureDistortion(&pic1, &pic2, type, disto)) {
290 fprintf(stderr, "Error while computing the distortion.\n");
291 goto End;
292 }
293 printf("%u %.2f %.2f %.2f %.2f %.2f [ %.2f bpp ]\n",
294 (unsigned int)size1,
295 disto[4], disto[0], disto[1], disto[2], disto[3],
296 8.f * size1 / pic1.width / pic1.height);
297
298 if (output != NULL) {
299 uint8_t* data = NULL;
300 size_t data_size = 0;
301 if (pic1.use_argb != pic2.use_argb) {
302 fprintf(stderr, "Pictures are not in the same argb format. "
303 "Can't save the difference map.\n");
304 goto End;
305 }
306 if (pic1.use_argb) {
307 int n;
308 fprintf(stderr, "max differences per channel: ");
309 for (n = 0; n < 3; ++n) { // skip the alpha channel
310 const int range = (type == 1) ?
311 SSIMScaleChannel((uint8_t*)pic1.argb + n, pic1.argb_stride * 4,
312 (const uint8_t*)pic2.argb + n, pic2.argb_stride * 4,
313 4, pic1.width, pic1.height, scale) :
314 DiffScaleChannel((uint8_t*)pic1.argb + n, pic1.argb_stride * 4,
315 (const uint8_t*)pic2.argb + n, pic2.argb_stride * 4,
316 4, pic1.width, pic1.height, scale);
317 if (range < 0) fprintf(stderr, "\nError computing diff map\n");
318 fprintf(stderr, "[%d]", range);
319 }
320 fprintf(stderr, "\n");
321 if (use_gray) ConvertToGray(&pic1);
322 } else {
323 fprintf(stderr, "Can only compute the difference map in ARGB format.\n");
324 goto End;
325 }
326 #if !defined(WEBP_REDUCE_CSP)
327 data_size = WebPEncodeLosslessBGRA((const uint8_t*)pic1.argb,
328 pic1.width, pic1.height,
329 pic1.argb_stride * 4,
330 &data);
331 if (data_size == 0) {
332 fprintf(stderr, "Error during lossless encoding.\n");
333 goto End;
334 }
335 ret = ImgIoUtilWriteFile(output, data, data_size) ? 0 : 1;
336 WebPFree(data);
337 if (ret) goto End;
338 #else
339 (void)data;
340 (void)data_size;
341 fprintf(stderr, "Cannot save the difference map. Please recompile "
342 "without the WEBP_REDUCE_CSP flag.\n");
343 #endif // WEBP_REDUCE_CSP
344 }
345 ret = 0;
346
347 End:
348 WebPPictureFree(&pic1);
349 WebPPictureFree(&pic2);
350 return ret;
351 }
352