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