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 <math.h>
15
16#include "./vp8enci.h"
17
18//------------------------------------------------------------------------------
19// local-min distortion
20//
21// For every pixel in the *reference* picture, we search for the local best
22// match in the compressed image. This is not a symmetrical measure.
23
24#define RADIUS 2  // search radius. Shouldn't be too large.
25
26static float AccumulateLSIM(const uint8_t* src, int src_stride,
27                            const uint8_t* ref, int ref_stride,
28                            int w, int h) {
29  int x, y;
30  double total_sse = 0.;
31  for (y = 0; y < h; ++y) {
32    const int y_0 = (y - RADIUS < 0) ? 0 : y - RADIUS;
33    const int y_1 = (y + RADIUS + 1 >= h) ? h : y + RADIUS + 1;
34    for (x = 0; x < w; ++x) {
35      const int x_0 = (x - RADIUS < 0) ? 0 : x - RADIUS;
36      const int x_1 = (x + RADIUS + 1 >= w) ? w : x + RADIUS + 1;
37      double best_sse = 255. * 255.;
38      const double value = (double)ref[y * ref_stride + x];
39      int i, j;
40      for (j = y_0; j < y_1; ++j) {
41        const uint8_t* s = src + j * src_stride;
42        for (i = x_0; i < x_1; ++i) {
43          const double sse = (double)(s[i] - value) * (s[i] - value);
44          if (sse < best_sse) best_sse = sse;
45        }
46      }
47      total_sse += best_sse;
48    }
49  }
50  return (float)total_sse;
51}
52#undef RADIUS
53
54//------------------------------------------------------------------------------
55// Distortion
56
57// Max value returned in case of exact similarity.
58static const double kMinDistortion_dB = 99.;
59static float GetPSNR(const double v) {
60  return (float)((v > 0.) ? -4.3429448 * log(v / (255 * 255.))
61                          : kMinDistortion_dB);
62}
63
64int WebPPictureDistortion(const WebPPicture* src, const WebPPicture* ref,
65                          int type, float result[5]) {
66  DistoStats stats[5];
67  int has_alpha;
68  int uv_w, uv_h;
69
70  if (src == NULL || ref == NULL ||
71      src->width != ref->width || src->height != ref->height ||
72      src->y == NULL || ref->y == NULL ||
73      src->u == NULL || ref->u == NULL ||
74      src->v == NULL || ref->v == NULL ||
75      result == NULL) {
76    return 0;
77  }
78  // TODO(skal): provide distortion for ARGB too.
79  if (src->use_argb == 1 || src->use_argb != ref->use_argb) {
80    return 0;
81  }
82
83  has_alpha = !!(src->colorspace & WEBP_CSP_ALPHA_BIT);
84  if (has_alpha != !!(ref->colorspace & WEBP_CSP_ALPHA_BIT) ||
85      (has_alpha && (src->a == NULL || ref->a == NULL))) {
86    return 0;
87  }
88
89  memset(stats, 0, sizeof(stats));
90
91  uv_w = (src->width + 1) >> 1;
92  uv_h = (src->height + 1) >> 1;
93  if (type >= 2) {
94    float sse[4];
95    sse[0] = AccumulateLSIM(src->y, src->y_stride,
96                            ref->y, ref->y_stride, src->width, src->height);
97    sse[1] = AccumulateLSIM(src->u, src->uv_stride,
98                            ref->u, ref->uv_stride, uv_w, uv_h);
99    sse[2] = AccumulateLSIM(src->v, src->uv_stride,
100                            ref->v, ref->uv_stride, uv_w, uv_h);
101    sse[3] = has_alpha ? AccumulateLSIM(src->a, src->a_stride,
102                                        ref->a, ref->a_stride,
103                                        src->width, src->height)
104                       : 0.f;
105    result[0] = GetPSNR(sse[0] / (src->width * src->height));
106    result[1] = GetPSNR(sse[1] / (uv_w * uv_h));
107    result[2] = GetPSNR(sse[2] / (uv_w * uv_h));
108    result[3] = GetPSNR(sse[3] / (src->width * src->height));
109    {
110      double total_sse = sse[0] + sse[1] + sse[2];
111      int total_pixels = src->width * src->height + 2 * uv_w * uv_h;
112      if (has_alpha) {
113        total_pixels += src->width * src->height;
114        total_sse += sse[3];
115      }
116      result[4] = GetPSNR(total_sse / total_pixels);
117    }
118  } else {
119    int c;
120    VP8SSIMAccumulatePlane(src->y, src->y_stride,
121                           ref->y, ref->y_stride,
122                           src->width, src->height, &stats[0]);
123    VP8SSIMAccumulatePlane(src->u, src->uv_stride,
124                           ref->u, ref->uv_stride,
125                           uv_w, uv_h, &stats[1]);
126    VP8SSIMAccumulatePlane(src->v, src->uv_stride,
127                           ref->v, ref->uv_stride,
128                           uv_w, uv_h, &stats[2]);
129    if (has_alpha) {
130      VP8SSIMAccumulatePlane(src->a, src->a_stride,
131                             ref->a, ref->a_stride,
132                             src->width, src->height, &stats[3]);
133    }
134    for (c = 0; c <= 4; ++c) {
135      if (type == 1) {
136        const double v = VP8SSIMGet(&stats[c]);
137        result[c] = (float)((v < 1.) ? -10.0 * log10(1. - v)
138                                     : kMinDistortion_dB);
139      } else {
140        const double v = VP8SSIMGetSquaredError(&stats[c]);
141        result[c] = GetPSNR(v);
142      }
143      // Accumulate forward
144      if (c < 4) VP8SSIMAddStats(&stats[c], &stats[4]);
145    }
146  }
147  return 1;
148}
149
150//------------------------------------------------------------------------------
151