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 "./vp8i_enc.h"
18 #include "../utils/utils.h"
19
20 typedef double (*AccumulateFunc)(const uint8_t* src, int src_stride,
21 const uint8_t* ref, int ref_stride,
22 int w, int h);
23
24 //------------------------------------------------------------------------------
25 // local-min distortion
26 //
27 // For every pixel in the *reference* picture, we search for the local best
28 // match in the compressed image. This is not a symmetrical measure.
29
30 #define RADIUS 2 // search radius. Shouldn't be too large.
31
AccumulateLSIM(const uint8_t * src,int src_stride,const uint8_t * ref,int ref_stride,int w,int h)32 static double AccumulateLSIM(const uint8_t* src, int src_stride,
33 const uint8_t* ref, int ref_stride,
34 int w, int h) {
35 int x, y;
36 double total_sse = 0.;
37 for (y = 0; y < h; ++y) {
38 const int y_0 = (y - RADIUS < 0) ? 0 : y - RADIUS;
39 const int y_1 = (y + RADIUS + 1 >= h) ? h : y + RADIUS + 1;
40 for (x = 0; x < w; ++x) {
41 const int x_0 = (x - RADIUS < 0) ? 0 : x - RADIUS;
42 const int x_1 = (x + RADIUS + 1 >= w) ? w : x + RADIUS + 1;
43 double best_sse = 255. * 255.;
44 const double value = (double)ref[y * ref_stride + x];
45 int i, j;
46 for (j = y_0; j < y_1; ++j) {
47 const uint8_t* const s = src + j * src_stride;
48 for (i = x_0; i < x_1; ++i) {
49 const double diff = s[i] - value;
50 const double sse = diff * diff;
51 if (sse < best_sse) best_sse = sse;
52 }
53 }
54 total_sse += best_sse;
55 }
56 }
57 return total_sse;
58 }
59 #undef RADIUS
60
AccumulateSSE(const uint8_t * src,int src_stride,const uint8_t * ref,int ref_stride,int w,int h)61 static double AccumulateSSE(const uint8_t* src, int src_stride,
62 const uint8_t* ref, int ref_stride,
63 int w, int h) {
64 int y;
65 double total_sse = 0.;
66 for (y = 0; y < h; ++y) {
67 total_sse += VP8AccumulateSSE(src, ref, w);
68 src += src_stride;
69 ref += ref_stride;
70 }
71 return total_sse;
72 }
73
74 //------------------------------------------------------------------------------
75
AccumulateSSIM(const uint8_t * src,int src_stride,const uint8_t * ref,int ref_stride,int w,int h)76 static double AccumulateSSIM(const uint8_t* src, int src_stride,
77 const uint8_t* ref, int ref_stride,
78 int w, int h) {
79 const int w0 = (w < VP8_SSIM_KERNEL) ? w : VP8_SSIM_KERNEL;
80 const int w1 = w - VP8_SSIM_KERNEL - 1;
81 const int h0 = (h < VP8_SSIM_KERNEL) ? h : VP8_SSIM_KERNEL;
82 const int h1 = h - VP8_SSIM_KERNEL - 1;
83 int x, y;
84 double sum = 0.;
85 for (y = 0; y < h0; ++y) {
86 for (x = 0; x < w; ++x) {
87 sum += VP8SSIMGetClipped(src, src_stride, ref, ref_stride, x, y, w, h);
88 }
89 }
90 for (; y < h1; ++y) {
91 for (x = 0; x < w0; ++x) {
92 sum += VP8SSIMGetClipped(src, src_stride, ref, ref_stride, x, y, w, h);
93 }
94 for (; x < w1; ++x) {
95 const int off1 = x - VP8_SSIM_KERNEL + (y - VP8_SSIM_KERNEL) * src_stride;
96 const int off2 = x - VP8_SSIM_KERNEL + (y - VP8_SSIM_KERNEL) * ref_stride;
97 sum += VP8SSIMGet(src + off1, src_stride, ref + off2, ref_stride);
98 }
99 for (; x < w; ++x) {
100 sum += VP8SSIMGetClipped(src, src_stride, ref, ref_stride, x, y, w, h);
101 }
102 }
103 for (; y < h; ++y) {
104 for (x = 0; x < w; ++x) {
105 sum += VP8SSIMGetClipped(src, src_stride, ref, ref_stride, x, y, w, h);
106 }
107 }
108 return sum;
109 }
110
111 //------------------------------------------------------------------------------
112 // Distortion
113
114 // Max value returned in case of exact similarity.
115 static const double kMinDistortion_dB = 99.;
116
GetPSNR(double v,double size)117 static double GetPSNR(double v, double size) {
118 return (v > 0. && size > 0.) ? -4.3429448 * log(v / (size * 255 * 255.))
119 : kMinDistortion_dB;
120 }
121
GetLogSSIM(double v,double size)122 static double GetLogSSIM(double v, double size) {
123 v = (size > 0.) ? v / size : 1.;
124 return (v < 1.) ? -10.0 * log10(1. - v) : kMinDistortion_dB;
125 }
126
WebPPlaneDistortion(const uint8_t * src,size_t src_stride,const uint8_t * ref,size_t ref_stride,int width,int height,size_t x_step,int type,float * distortion,float * result)127 int WebPPlaneDistortion(const uint8_t* src, size_t src_stride,
128 const uint8_t* ref, size_t ref_stride,
129 int width, int height, size_t x_step,
130 int type, float* distortion, float* result) {
131 uint8_t* allocated = NULL;
132 const AccumulateFunc metric = (type == 0) ? AccumulateSSE :
133 (type == 1) ? AccumulateSSIM :
134 AccumulateLSIM;
135 if (src == NULL || ref == NULL ||
136 src_stride < x_step * width || ref_stride < x_step * width ||
137 result == NULL || distortion == NULL) {
138 return 0;
139 }
140
141 VP8SSIMDspInit();
142 if (x_step != 1) { // extract a packed plane if needed
143 int x, y;
144 uint8_t* tmp1;
145 uint8_t* tmp2;
146 allocated =
147 (uint8_t*)WebPSafeMalloc(2ULL * width * height, sizeof(*allocated));
148 if (allocated == NULL) return 0;
149 tmp1 = allocated;
150 tmp2 = tmp1 + (size_t)width * height;
151 for (y = 0; y < height; ++y) {
152 for (x = 0; x < width; ++x) {
153 tmp1[x + y * width] = src[x * x_step + y * src_stride];
154 tmp2[x + y * width] = ref[x * x_step + y * ref_stride];
155 }
156 }
157 src = tmp1;
158 ref = tmp2;
159 }
160 *distortion = (float)metric(src, width, ref, width, width, height);
161 WebPSafeFree(allocated);
162
163 *result = (type == 1) ? (float)GetLogSSIM(*distortion, (double)width * height)
164 : (float)GetPSNR(*distortion, (double)width * height);
165 return 1;
166 }
167
WebPPictureDistortion(const WebPPicture * src,const WebPPicture * ref,int type,float results[5])168 int WebPPictureDistortion(const WebPPicture* src, const WebPPicture* ref,
169 int type, float results[5]) {
170 int w, h, c;
171 int ok = 0;
172 WebPPicture p0, p1;
173 double total_size = 0., total_distortion = 0.;
174 if (src == NULL || ref == NULL ||
175 src->width != ref->width || src->height != ref->height ||
176 results == NULL) {
177 return 0;
178 }
179
180 VP8SSIMDspInit();
181 if (!WebPPictureInit(&p0) || !WebPPictureInit(&p1)) return 0;
182 w = src->width;
183 h = src->height;
184 if (!WebPPictureView(src, 0, 0, w, h, &p0)) goto Error;
185 if (!WebPPictureView(ref, 0, 0, w, h, &p1)) goto Error;
186
187 // We always measure distortion in ARGB space.
188 if (p0.use_argb == 0 && !WebPPictureYUVAToARGB(&p0)) goto Error;
189 if (p1.use_argb == 0 && !WebPPictureYUVAToARGB(&p1)) goto Error;
190 for (c = 0; c < 4; ++c) {
191 float distortion;
192 const size_t stride0 = 4 * (size_t)p0.argb_stride;
193 const size_t stride1 = 4 * (size_t)p1.argb_stride;
194 if (!WebPPlaneDistortion((const uint8_t*)p0.argb + c, stride0,
195 (const uint8_t*)p1.argb + c, stride1,
196 w, h, 4, type, &distortion, results + c)) {
197 goto Error;
198 }
199 total_distortion += distortion;
200 total_size += w * h;
201 }
202
203 results[4] = (type == 1) ? (float)GetLogSSIM(total_distortion, total_size)
204 : (float)GetPSNR(total_distortion, total_size);
205 ok = 1;
206
207 Error:
208 WebPPictureFree(&p0);
209 WebPPictureFree(&p1);
210 return ok;
211 }
212
213 //------------------------------------------------------------------------------
214