2 * Copyright (c) 2003-2013 Loren Merritt
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.
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.
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
20 * Computes the Structural Similarity Metric between two rawYV12 video files.
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.
26 * To improve speed, this implementation uses the standard approximation of
27 * overlapped 8x8 block sums, rather than the original gaussian weights.
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))
41 #define PIXEL_MAX ((1 << BIT_DEPTH)-1)
42 typedef uint8_t pixel
;
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
,
53 for( z
= 0; z
< 2; z
++ )
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
++ )
59 int a
= pix1
[x
+y
*stride1
];
60 int b
= pix2
[x
+y
*stride2
];
76 static float ssim_end1( int s1
, int s2
, int ss
, int s12
)
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. */
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;
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);
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
));
101 static float ssim_end4( int sum0
[5][4], int sum1
[5][4], int width
)
106 for( i
= 0; i
< width
; i
++ )
107 ssim
+= ssim_end1( sum0
[i
][0] + sum0
[i
+1][0] + sum1
[i
][0] + sum1
[i
+1][0],
108 sum0
[i
][1] + sum0
[i
+1][1] + sum1
[i
][1] + sum1
[i
+1][1],
109 sum0
[i
][2] + sum0
[i
+1][2] + sum1
[i
][2] + sum1
[i
+1][2],
110 sum0
[i
][3] + sum0
[i
+1][3] + sum1
[i
][3] + sum1
[i
+1][3] );
115 pixel
*pix1
, intptr_t stride1
,
116 pixel
*pix2
, intptr_t stride2
,
117 int width
, int height
, void *buf
, int *cnt
)
122 int (*sum0
)[4] = buf
;
123 int (*sum1
)[4] = sum0
+ (width
>> 2) + 3;
126 for( y
= 1; y
< height
; y
++ )
130 FFSWAP( void*, sum0
, sum1
);
131 for( x
= 0; x
< width
; x
+=2 )
132 ssim_4x4x2_core( &pix1
[4*(x
+z
*stride1
)], stride1
, &pix2
[4*(x
+z
*stride2
)], stride2
, &sum0
[x
] );
134 for( x
= 0; x
< width
-1; x
+= 4 )
135 ssim
+= ssim_end4( sum0
+x
, sum1
+x
, FFMIN(4,width
-x
-1) );
137 // *cnt = (height-1) * (width-1);
138 return ssim
/ ((height
-1) * (width
-1));
142 uint64_t ssd_plane( const uint8_t *pix1
, const uint8_t *pix2
, int size
)
146 for( i
=0; i
<size
; i
++ )
148 int d
= pix1
[i
] - pix2
[i
];
154 static double ssd_to_psnr( uint64_t ssd
, uint64_t denom
)
156 return -10*log((double)ssd
/(denom
*255*255))/log(10);
159 static double ssim_db( double ssim
, double weight
)
161 return 10*(log(weight
)/log(10)-log(weight
-ssim
)/log(10));
164 static void print_results(uint64_t ssd
[3], double ssim
[3], int frames
, int w
, int h
)
166 printf( "PSNR Y:%.3f U:%.3f V:%.3f All:%.3f | ",
167 ssd_to_psnr( ssd
[0], (uint64_t)frames
*w
*h
),
168 ssd_to_psnr( ssd
[1], (uint64_t)frames
*w
*h
/4 ),
169 ssd_to_psnr( ssd
[2], (uint64_t)frames
*w
*h
/4 ),
170 ssd_to_psnr( ssd
[0] + ssd
[1] + ssd
[2], (uint64_t)frames
*w
*h
*3/2 ) );
171 printf( "SSIM Y:%.5f U:%.5f V:%.5f All:%.5f (%.5f)",
175 (ssim
[0]*4 + ssim
[1] + ssim
[2]) / (frames
*6),
176 ssim_db(ssim
[0] * 4 + ssim
[1] + ssim
[2], frames
*6));
179 int main(int argc
, char* argv
[])
182 uint8_t *buf
[2], *plane
[2][3];
184 uint64_t ssd
[3] = {0,0,0};
185 double ssim
[3] = {0,0,0};
186 int frame_size
, w
, h
;
190 if( argc
<4 || 2 != sscanf(argv
[3], "%dx%d", &w
, &h
) )
192 printf("tiny_ssim <file1.yuv> <file2.yuv> <width>x<height> [<seek>]\n");
196 f
[0] = fopen(argv
[1], "rb");
197 f
[1] = fopen(argv
[2], "rb");
198 sscanf(argv
[3], "%dx%d", &w
, &h
);
200 if (w
<=0 || h
<=0 || w
*(int64_t)h
>= INT_MAX
/3 || 2LL*w
+12 >= INT_MAX
/ sizeof(*temp
)) {
201 fprintf(stderr
, "Dimensions are too large, or invalid\n");
205 frame_size
= w
*h
*3LL/2;
208 buf
[i
] = malloc(frame_size
);
209 plane
[i
][0] = buf
[i
];
210 plane
[i
][1] = plane
[i
][0] + w
*h
;
211 plane
[i
][2] = plane
[i
][1] + w
*h
/4;
213 temp
= malloc((2*w
+12)*sizeof(*temp
));
214 seek
= argc
<5 ? 0 : atoi(argv
[4]);
215 fseek(f
[seek
<0], seek
< 0 ? -seek
: seek
, SEEK_SET
);
217 for( frames
=0;; frames
++ )
221 if( fread(buf
[0], frame_size
, 1, f
[0]) != 1) break;
222 if( fread(buf
[1], frame_size
, 1, f
[1]) != 1) break;
225 ssd_one
[i
] = ssd_plane ( plane
[0][i
], plane
[1][i
], w
*h
>>2*!!i
);
226 ssim_one
[i
] = ssim_plane( plane
[0][i
], w
>>!!i
,
228 w
>>!!i
, h
>>!!i
, temp
, NULL
);
229 ssd
[i
] += ssd_one
[i
];
230 ssim
[i
] += ssim_one
[i
];
233 printf("Frame %d | ", frames
);
234 print_results(ssd_one
, ssim_one
, 1, w
, h
);
239 if( !frames
) return 0;
241 printf("Total %d frames | ", frames
);
242 print_results(ssd
, ssim
, frames
, w
, h
);