blob: 08f8e92a03c238af310acb02b4c398a8ac86d97f
1 | /* |
2 | * Copyright (c) 2003-2013 Loren Merritt |
3 | * |
4 | * This program is free software; you can redistribute it and/or modify |
5 | * it under the terms of the GNU General Public License as published by |
6 | * the Free Software Foundation; either version 2 of the License, or |
7 | * (at your option) any later version. |
8 | * |
9 | * This program is distributed in the hope that it will be useful, |
10 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
11 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
12 | * GNU General Public License for more details. |
13 | * |
14 | * You should have received a copy of the GNU General Public License |
15 | * along with this program; if not, write to the Free Software |
16 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110 USA |
17 | */ |
18 | /* |
19 | * tiny_ssim.c |
20 | * Computes the Structural Similarity Metric between two rawYV12 video files. |
21 | * original algorithm: |
22 | * Z. Wang, A. C. Bovik, H. R. Sheikh and E. P. Simoncelli, |
23 | * "Image quality assessment: From error visibility to structural similarity," |
24 | * IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600-612, Apr. 2004. |
25 | * |
26 | * To improve speed, this implementation uses the standard approximation of |
27 | * overlapped 8x8 block sums, rather than the original gaussian weights. |
28 | */ |
29 | |
30 | #include "config.h" |
31 | #include <inttypes.h> |
32 | #include <limits.h> |
33 | #include <math.h> |
34 | #include <stdio.h> |
35 | #include <stdlib.h> |
36 | |
37 | #define FFSWAP(type,a,b) do{type SWAP_tmp= b; b= a; a= SWAP_tmp;}while(0) |
38 | #define FFMIN(a,b) ((a) > (b) ? (b) : (a)) |
39 | |
40 | #define BIT_DEPTH 8 |
41 | #define PIXEL_MAX ((1 << BIT_DEPTH)-1) |
42 | typedef uint8_t pixel; |
43 | |
44 | /**************************************************************************** |
45 | * structural similarity metric |
46 | ****************************************************************************/ |
47 | static void ssim_4x4x2_core( const pixel *pix1, intptr_t stride1, |
48 | const pixel *pix2, intptr_t stride2, |
49 | int sums[2][4] ) |
50 | { |
51 | int x,y,z; |
52 | |
53 | for( z = 0; z < 2; z++ ) |
54 | { |
55 | uint32_t s1 = 0, s2 = 0, ss = 0, s12 = 0; |
56 | for( y = 0; y < 4; y++ ) |
57 | for( x = 0; x < 4; x++ ) |
58 | { |
59 | int a = pix1[x+y*stride1]; |
60 | int b = pix2[x+y*stride2]; |
61 | s1 += a; |
62 | s2 += b; |
63 | ss += a*a; |
64 | ss += b*b; |
65 | s12 += a*b; |
66 | } |
67 | sums[z][0] = s1; |
68 | sums[z][1] = s2; |
69 | sums[z][2] = ss; |
70 | sums[z][3] = s12; |
71 | pix1 += 4; |
72 | pix2 += 4; |
73 | } |
74 | } |
75 | |
76 | static float ssim_end1( int s1, int s2, int ss, int s12 ) |
77 | { |
78 | /* Maximum value for 10-bit is: ss*64 = (2^10-1)^2*16*4*64 = 4286582784, which will overflow in some cases. |
79 | * s1*s1, s2*s2, and s1*s2 also obtain this value for edge cases: ((2^10-1)*16*4)^2 = 4286582784. |
80 | * Maximum value for 9-bit is: ss*64 = (2^9-1)^2*16*4*64 = 1069551616, which will not overflow. */ |
81 | #if BIT_DEPTH > 9 |
82 | typedef float type; |
83 | static const float ssim_c1 = .01*.01*PIXEL_MAX*PIXEL_MAX*64; |
84 | static const float ssim_c2 = .03*.03*PIXEL_MAX*PIXEL_MAX*64*63; |
85 | #else |
86 | typedef int type; |
87 | static const int ssim_c1 = (int)(.01*.01*PIXEL_MAX*PIXEL_MAX*64 + .5); |
88 | static const int ssim_c2 = (int)(.03*.03*PIXEL_MAX*PIXEL_MAX*64*63 + .5); |
89 | #endif |
90 | type fs1 = s1; |
91 | type fs2 = s2; |
92 | type fss = ss; |
93 | type fs12 = s12; |
94 | type vars = fss*64 - fs1*fs1 - fs2*fs2; |
95 | type covar = fs12*64 - fs1*fs2; |
96 | return (float)(2*fs1*fs2 + ssim_c1) * (float)(2*covar + ssim_c2) |
97 | / ((float)(fs1*fs1 + fs2*fs2 + ssim_c1) * (float)(vars + ssim_c2)); |
98 | } |
99 | |
100 | static float ssim_end4( int sum0[5][4], int sum1[5][4], int width ) |
101 | { |
102 | float ssim = 0.0; |
103 | int i; |
104 | |
105 | for( i = 0; i < width; i++ ) |
106 | ssim += ssim_end1( sum0[i][0] + sum0[i+1][0] + sum1[i][0] + sum1[i+1][0], |
107 | sum0[i][1] + sum0[i+1][1] + sum1[i][1] + sum1[i+1][1], |
108 | sum0[i][2] + sum0[i+1][2] + sum1[i][2] + sum1[i+1][2], |
109 | sum0[i][3] + sum0[i+1][3] + sum1[i][3] + sum1[i+1][3] ); |
110 | return ssim; |
111 | } |
112 | |
113 | float ssim_plane( |
114 | pixel *pix1, intptr_t stride1, |
115 | pixel *pix2, intptr_t stride2, |
116 | int width, int height, void *buf, int *cnt ) |
117 | { |
118 | int z = 0; |
119 | int x, y; |
120 | float ssim = 0.0; |
121 | int (*sum0)[4] = buf; |
122 | int (*sum1)[4] = sum0 + (width >> 2) + 3; |
123 | width >>= 2; |
124 | height >>= 2; |
125 | for( y = 1; y < height; y++ ) |
126 | { |
127 | for( ; z <= y; z++ ) |
128 | { |
129 | FFSWAP( void*, sum0, sum1 ); |
130 | for( x = 0; x < width; x+=2 ) |
131 | ssim_4x4x2_core( &pix1[4*(x+z*stride1)], stride1, &pix2[4*(x+z*stride2)], stride2, &sum0[x] ); |
132 | } |
133 | for( x = 0; x < width-1; x += 4 ) |
134 | ssim += ssim_end4( sum0+x, sum1+x, FFMIN(4,width-x-1) ); |
135 | } |
136 | // *cnt = (height-1) * (width-1); |
137 | return ssim / ((height-1) * (width-1)); |
138 | } |
139 | |
140 | |
141 | uint64_t ssd_plane( const uint8_t *pix1, const uint8_t *pix2, int size ) |
142 | { |
143 | uint64_t ssd = 0; |
144 | int i; |
145 | for( i=0; i<size; i++ ) |
146 | { |
147 | int d = pix1[i] - pix2[i]; |
148 | ssd += d*d; |
149 | } |
150 | return ssd; |
151 | } |
152 | |
153 | static double ssd_to_psnr( uint64_t ssd, uint64_t denom ) |
154 | { |
155 | return -10*log((double)ssd/(denom*255*255))/log(10); |
156 | } |
157 | |
158 | static double ssim_db( double ssim, double weight ) |
159 | { |
160 | return 10*(log(weight)/log(10)-log(weight-ssim)/log(10)); |
161 | } |
162 | |
163 | static void print_results(uint64_t ssd[3], double ssim[3], int frames, int w, int h) |
164 | { |
165 | printf( "PSNR Y:%.3f U:%.3f V:%.3f All:%.3f | ", |
166 | ssd_to_psnr( ssd[0], (uint64_t)frames*w*h ), |
167 | ssd_to_psnr( ssd[1], (uint64_t)frames*w*h/4 ), |
168 | ssd_to_psnr( ssd[2], (uint64_t)frames*w*h/4 ), |
169 | ssd_to_psnr( ssd[0] + ssd[1] + ssd[2], (uint64_t)frames*w*h*3/2 ) ); |
170 | printf( "SSIM Y:%.5f U:%.5f V:%.5f All:%.5f (%.5f)", |
171 | ssim[0] / frames, |
172 | ssim[1] / frames, |
173 | ssim[2] / frames, |
174 | (ssim[0]*4 + ssim[1] + ssim[2]) / (frames*6), |
175 | ssim_db(ssim[0] * 4 + ssim[1] + ssim[2], frames*6)); |
176 | } |
177 | |
178 | int main(int argc, char* argv[]) |
179 | { |
180 | FILE *f[2]; |
181 | uint8_t *buf[2], *plane[2][3]; |
182 | int *temp; |
183 | uint64_t ssd[3] = {0,0,0}; |
184 | double ssim[3] = {0,0,0}; |
185 | int frame_size, w, h; |
186 | int frames, seek; |
187 | int i; |
188 | |
189 | if( argc<4 || 2 != sscanf(argv[3], "%dx%d", &w, &h) ) |
190 | { |
191 | printf("tiny_ssim <file1.yuv> <file2.yuv> <width>x<height> [<seek>]\n"); |
192 | return -1; |
193 | } |
194 | |
195 | f[0] = fopen(argv[1], "rb"); |
196 | f[1] = fopen(argv[2], "rb"); |
197 | sscanf(argv[3], "%dx%d", &w, &h); |
198 | |
199 | if (w<=0 || h<=0 || w*(int64_t)h >= INT_MAX/3 || 2LL*w+12 >= INT_MAX / sizeof(*temp)) { |
200 | fprintf(stderr, "Dimensions are too large, or invalid\n"); |
201 | return -2; |
202 | } |
203 | |
204 | frame_size = w*h*3LL/2; |
205 | for( i=0; i<2; i++ ) |
206 | { |
207 | buf[i] = malloc(frame_size); |
208 | plane[i][0] = buf[i]; |
209 | plane[i][1] = plane[i][0] + w*h; |
210 | plane[i][2] = plane[i][1] + w*h/4; |
211 | } |
212 | temp = malloc((2*w+12)*sizeof(*temp)); |
213 | seek = argc<5 ? 0 : atoi(argv[4]); |
214 | fseek(f[seek<0], seek < 0 ? -seek : seek, SEEK_SET); |
215 | |
216 | for( frames=0;; frames++ ) |
217 | { |
218 | uint64_t ssd_one[3]; |
219 | double ssim_one[3]; |
220 | if( fread(buf[0], frame_size, 1, f[0]) != 1) break; |
221 | if( fread(buf[1], frame_size, 1, f[1]) != 1) break; |
222 | for( i=0; i<3; i++ ) |
223 | { |
224 | ssd_one[i] = ssd_plane ( plane[0][i], plane[1][i], w*h>>2*!!i ); |
225 | ssim_one[i] = ssim_plane( plane[0][i], w>>!!i, |
226 | plane[1][i], w>>!!i, |
227 | w>>!!i, h>>!!i, temp, NULL ); |
228 | ssd[i] += ssd_one[i]; |
229 | ssim[i] += ssim_one[i]; |
230 | } |
231 | |
232 | printf("Frame %d | ", frames); |
233 | print_results(ssd_one, ssim_one, 1, w, h); |
234 | printf(" \r"); |
235 | fflush(stdout); |
236 | } |
237 | |
238 | if( !frames ) return 0; |
239 | |
240 | printf("Total %d frames | ", frames); |
241 | print_results(ssd, ssim, frames, w, h); |
242 | printf("\n"); |
243 | |
244 | return 0; |
245 | } |
246 |