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