add LSIM metric to WebPPictureDistortion()

LSIM stands for "local similarity": before matching
a compressed pixel to the source, we search around in the source
and minimise the squared error. So, this is close to PSNR calculation,
but mitigates some of its limitations (pure translation and noise for instance).

There's a new -print_lsim option to cwebp too.

Change-Id: Ia38561034c7a90e71d2ea0f55bb1de527eda245b

README

@@ -158,8 +158,9 @@ options:
   -crop <x> <y> <w> <h> .. crop picture with the given rectangle
   -resize <w> <h> ........ resize picture (after any cropping)
   -map <int> ............. print map of extra info.
-  -print_ssim ............ prints averaged SSIM distortion.
   -print_psnr ............ prints averaged PSNR distortion.
+  -print_ssim ............ prints averaged SSIM distortion.
+  -print_lsim ............ prints local-similarity distortion.
   -d <file.pgm> .......... dump the compressed output (PGM file).
   -alpha_method <int> .... Transparency-compression method (0..1)
   -alpha_filter <string> . predictive filtering for alpha plane.

examples/cwebp.c

@@ -833,8 +833,9 @@ static void HelpLong(void) {
   printf("  -444 / -422 / -gray ..... Change colorspace\n");
 #endif
   printf("  -map <int> ............. print map of extra info.\n");
-  printf("  -print_ssim ............ prints averaged SSIM distortion.\n");
   printf("  -print_psnr ............ prints averaged PSNR distortion.\n");
+  printf("  -print_ssim ............ prints averaged SSIM distortion.\n");
+  printf("  -print_lsim ............ prints local-similarity distortion.\n");
   printf("  -d <file.pgm> .......... dump the compressed output (PGM file).\n");
   printf("  -alpha_method <int> .... Transparency-compression method (0..1)\n");
   printf("  -alpha_filter <string> . predictive filtering for alpha plane.\n");
@@ -898,7 +899,7 @@ int main(int argc, const char *argv[]) {
   int resize_w = 0, resize_h = 0;
   int show_progress = 0;
   WebPPicture picture;
-  int print_distortion = 0;        // 1=PSNR, 2=SSIM
+  int print_distortion = -1;        // -1=off, 0=PSNR, 1=SSIM, 2=LSIM
   WebPPicture original_picture;    // when PSNR or SSIM is requested
   WebPConfig config;
   WebPAuxStats stats;
@@ -928,12 +929,15 @@ int main(int argc, const char *argv[]) {
     } else if (!strcmp(argv[c], "-d") && c < argc - 1) {
       dump_file = argv[++c];
       config.show_compressed = 1;
-    } else if (!strcmp(argv[c], "-print_ssim")) {
-      config.show_compressed = 1;
-      print_distortion = 2;
     } else if (!strcmp(argv[c], "-print_psnr")) {
+      config.show_compressed = 1;
+      print_distortion = 0;
+    } else if (!strcmp(argv[c], "-print_ssim")) {
       config.show_compressed = 1;
       print_distortion = 1;
+    } else if (!strcmp(argv[c], "-print_lsim")) {
+      config.show_compressed = 1;
+      print_distortion = 2;
     } else if (!strcmp(argv[c], "-short")) {
       short_output++;
     } else if (!strcmp(argv[c], "-s") && c < argc - 2) {
@@ -1149,7 +1153,7 @@ int main(int argc, const char *argv[]) {
   if (picture.extra_info_type > 0) {
     AllocExtraInfo(&picture);
   }
-  if (print_distortion > 0) {  // Save original picture for later comparison
+  if (print_distortion >= 0) {  // Save original picture for later comparison
     WebPPictureCopy(&picture, &original_picture);
   }
   if (!WebPEncode(&config, &picture)) {
@@ -1179,12 +1183,13 @@ int main(int argc, const char *argv[]) {
       PrintExtraInfoLossy(&picture, short_output, in_file);
     }
   }
-  if (!quiet && !short_output && print_distortion > 0) {  // print distortion
+  if (!quiet && !short_output && print_distortion >= 0) {  // print distortion
+    static const char* distortion_names[] = { "PSNR", "SSIM", "LSIM" };
     float values[5];
     WebPPictureDistortion(&picture, &original_picture,
-                          (print_distortion == 1) ? 0 : 1, values);
+                          print_distortion, values);
     fprintf(stderr, "%s: Y:%.2f U:%.2f V:%.2f A:%.2f  Total:%.2f\n",
-            (print_distortion == 1) ? "PSNR" : "SSIM",
+            distortion_names[print_distortion],
             values[0], values[1], values[2], values[3], values[4]);
   }
   return_value = 0;

man/cwebp.1

@@ -178,6 +178,9 @@ Compute and report average PSNR (Peak-Signal-To-Noise ratio).
 .B \-print_ssim
 Compute and report average SSIM (structural similarity metric)
 .TP
+.B \-print_lsim
+Compute and report local similarity metric.
+.TP
 .B \-progress
 Report encoding progress in percent.
 .TP

src/enc/picture.c

@@ -906,54 +906,122 @@ void WebPCleanupTransparentArea(WebPPicture* pic) {
 #undef SIZE
 #undef SIZE2
 
+//------------------------------------------------------------------------------
+// local-min distortion
+//
+// For every pixel in the *reference* picture, we search for the local best
+// match in the compressed image. This is not a symmetrical measure.
+
+// search radius. Shouldn't be too large.
+#define RADIUS 2
+
+static double AccumulateLSIM(const uint8_t* src, int src_stride,
+                             const uint8_t* ref, int ref_stride,
+                             int w, int h) {
+  int x, y;
+  double total_sse = 0.;
+  for (y = 0; y < h; ++y) {
+    const int y0 = (y - RADIUS < 0) ? 0 : y - RADIUS;
+    const int y1 = (y + RADIUS + 1 >= h) ? h : y + RADIUS + 1;
+    for (x = 0; x < w; ++x) {
+      const int x0 = (x - RADIUS < 0) ? 0 : x - RADIUS;
+      const int x1 = (x + RADIUS + 1 >= w) ? w : x + RADIUS + 1;
+      double best_sse = 255. * 255.;
+      const double value = (double)ref[y * ref_stride + x];
+      int i, j;
+      for (j = y0; j < y1; ++j) {
+        const uint8_t* s = src + j * src_stride;
+        for (i = x0; i < x1; ++i) {
+          const double sse = (double)(s[i] - value) * (s[i] - value);
+          if (sse < best_sse) best_sse = sse;
+        }
+      }
+      total_sse += best_sse;
+    }
+  }
+  return total_sse;
+}
+#undef RADIUS
 
 //------------------------------------------------------------------------------
 // Distortion
 
 // Max value returned in case of exact similarity.
 static const double kMinDistortion_dB = 99.;
+static float GetPSNR(const double v) {
+  return (float)((v > 0.) ? -4.3429448 * log(v / (255 * 255.))
+                          : kMinDistortion_dB);
+}
 
-int WebPPictureDistortion(const WebPPicture* pic1, const WebPPicture* pic2,
+int WebPPictureDistortion(const WebPPicture* src, const WebPPicture* ref,
                           int type, float result[5]) {
-  int c;
   DistoStats stats[5];
   int has_alpha;
+  int uv_w, uv_h;
 
-  if (pic1 == NULL || pic2 == NULL ||
+  if (src == NULL || ref == NULL ||
-      pic1->width != pic2->width || pic1->height != pic2->height ||
+      src->width != ref->width || src->height != ref->height ||
-      pic1->y == NULL || pic2->y == NULL ||
+      src->y == NULL || ref->y == NULL ||
-      pic1->u == NULL || pic2->u == NULL ||
+      src->u == NULL || ref->u == NULL ||
-      pic1->v == NULL || pic2->v == NULL ||
+      src->v == NULL || ref->v == NULL ||
       result == NULL) {
     return 0;
   }
   // TODO(skal): provide distortion for ARGB too.
-  if (pic1->use_argb == 1 || pic1->use_argb != pic2->use_argb) {
+  if (src->use_argb == 1 || src->use_argb != ref->use_argb) {
     return 0;
   }
 
-  has_alpha = !!(pic1->colorspace & WEBP_CSP_ALPHA_BIT);
+  has_alpha = !!(src->colorspace & WEBP_CSP_ALPHA_BIT);
-  if (has_alpha != !!(pic2->colorspace & WEBP_CSP_ALPHA_BIT) ||
+  if (has_alpha != !!(ref->colorspace & WEBP_CSP_ALPHA_BIT) ||
-      (has_alpha && (pic1->a == NULL || pic2->a == NULL))) {
+      (has_alpha && (src->a == NULL || ref->a == NULL))) {
     return 0;
   }
 
   memset(stats, 0, sizeof(stats));
-  VP8SSIMAccumulatePlane(pic1->y, pic1->y_stride,
+
-                         pic2->y, pic2->y_stride,
+  uv_w = HALVE(src->width);
-                         pic1->width, pic1->height, &stats[0]);
+  uv_h = HALVE(src->height);
-  VP8SSIMAccumulatePlane(pic1->u, pic1->uv_stride,
+  if (type >= 2) {
-                         pic2->u, pic2->uv_stride,
+    float sse[4];
-                         (pic1->width + 1) >> 1, (pic1->height + 1) >> 1,
+    sse[0] = AccumulateLSIM(src->y, src->y_stride,
-                         &stats[1]);
+                            ref->y, ref->y_stride, src->width, src->height);
-  VP8SSIMAccumulatePlane(pic1->v, pic1->uv_stride,
+    sse[1] = AccumulateLSIM(src->u, src->uv_stride,
-                         pic2->v, pic2->uv_stride,
+                            ref->u, ref->uv_stride, uv_w, uv_h);
-                         (pic1->width + 1) >> 1, (pic1->height + 1) >> 1,
+    sse[2] = AccumulateLSIM(src->v, src->uv_stride,
-                         &stats[2]);
+                            ref->v, ref->uv_stride, uv_w, uv_h);
+    sse[3] = has_alpha ? AccumulateLSIM(src->a, src->a_stride,
+                                        ref->a, ref->a_stride,
+                                        src->width, src->height)
+                       : 0;
+    result[0] = GetPSNR(sse[0] / (src->width * src->height));
+    result[1] = GetPSNR(sse[1] / (uv_w * uv_h));
+    result[2] = GetPSNR(sse[2] / (uv_w * uv_h));
+    result[3] = GetPSNR(sse[3] / (src->width * src->height));
+    {
+      double total_sse = sse[0] + sse[1] + sse[2];
+      int total_pixels = src->width * src->height + 2 * uv_w * uv_h;
       if (has_alpha) {
-    VP8SSIMAccumulatePlane(pic1->a, pic1->a_stride,
+        total_pixels += src->width * src->height;
-                           pic2->a, pic2->a_stride,
+        total_sse += sse[3];
-                           pic1->width, pic1->height, &stats[3]);
+      }
+      result[4] = GetPSNR(total_sse / total_pixels);
+    }
+  } else {
+    int c;
+    VP8SSIMAccumulatePlane(src->y, src->y_stride,
+                           ref->y, ref->y_stride,
+                           src->width, src->height, &stats[0]);
+    VP8SSIMAccumulatePlane(src->u, src->uv_stride,
+                           ref->u, ref->uv_stride,
+                           uv_w, uv_h, &stats[1]);
+    VP8SSIMAccumulatePlane(src->v, src->uv_stride,
+                           ref->v, ref->uv_stride,
+                           uv_w, uv_h, &stats[2]);
+    if (has_alpha) {
+      VP8SSIMAccumulatePlane(src->a, src->a_stride,
+                             ref->a, ref->a_stride,
+                             src->width, src->height, &stats[3]);
     }
     for (c = 0; c <= 4; ++c) {
       if (type == 1) {
@@ -962,12 +1030,12 @@ int WebPPictureDistortion(const WebPPicture* pic1, const WebPPicture* pic2,
                                      : kMinDistortion_dB);
       } else {
         const double v = VP8SSIMGetSquaredError(&stats[c]);
-      result[c] = (float)((v > 0.) ? -4.3429448 * log(v / (255 * 255.))
+        result[c] = GetPSNR(v);
-                                   : kMinDistortion_dB);
       }
       // Accumulate forward
       if (c < 4) VP8SSIMAddStats(&stats[c], &stats[4]);
     }
+  }
   return 1;
 }
 

src/webp/encode.h

@@ -359,13 +359,13 @@ WEBP_EXTERN(void) WebPPictureFree(WebPPicture* picture);
 // Returns false in case of memory allocation error.
 WEBP_EXTERN(int) WebPPictureCopy(const WebPPicture* src, WebPPicture* dst);
 
-// Compute PSNR or SSIM distortion between two pictures.
+// Compute PSNR, SSIM or LSIM distortion metric between two pictures.
 // Result is in dB, stores in result[] in the Y/U/V/Alpha/All order.
-// Returns false in case of error (pic1 and pic2 don't have same dimension, ...)
+// Returns false in case of error (src and ref don't have same dimension, ...)
 // Warning: this function is rather CPU-intensive.
 WEBP_EXTERN(int) WebPPictureDistortion(
-    const WebPPicture* pic1, const WebPPicture* pic2,
+    const WebPPicture* src, const WebPPicture* ref,
-    int metric_type,           // 0 = PSNR, 1 = SSIM
+    int metric_type,           // 0 = PSNR, 1 = SSIM, 2 = LSIM
     float result[5]);
 
 // self-crops a picture to the rectangle defined by top/left/width/height.