| 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 | #define type float |
| 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 | #define type int |
| 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 | #undef type |
| 99 | } |
| 100 | |
| 101 | static float ssim_end4( int sum0[5][4], int sum1[5][4], int width ) |
| 102 | { |
| 103 | float ssim = 0.0; |
| 104 | int i; |
| 105 | |
| 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] ); |
| 111 | return ssim; |
| 112 | } |
| 113 | |
| 114 | float ssim_plane( |
| 115 | pixel *pix1, intptr_t stride1, |
| 116 | pixel *pix2, intptr_t stride2, |
| 117 | int width, int height, void *buf, int *cnt ) |
| 118 | { |
| 119 | int z = 0; |
| 120 | int x, y; |
| 121 | float ssim = 0.0; |
| 122 | int (*sum0)[4] = buf; |
| 123 | int (*sum1)[4] = sum0 + (width >> 2) + 3; |
| 124 | width >>= 2; |
| 125 | height >>= 2; |
| 126 | for( y = 1; y < height; y++ ) |
| 127 | { |
| 128 | for( ; z <= y; z++ ) |
| 129 | { |
| 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] ); |
| 133 | } |
| 134 | for( x = 0; x < width-1; x += 4 ) |
| 135 | ssim += ssim_end4( sum0+x, sum1+x, FFMIN(4,width-x-1) ); |
| 136 | } |
| 137 | // *cnt = (height-1) * (width-1); |
| 138 | return ssim / ((height-1) * (width-1)); |
| 139 | } |
| 140 | |
| 141 | |
| 142 | uint64_t ssd_plane( const uint8_t *pix1, const uint8_t *pix2, int size ) |
| 143 | { |
| 144 | uint64_t ssd = 0; |
| 145 | int i; |
| 146 | for( i=0; i<size; i++ ) |
| 147 | { |
| 148 | int d = pix1[i] - pix2[i]; |
| 149 | ssd += d*d; |
| 150 | } |
| 151 | return ssd; |
| 152 | } |
| 153 | |
| 154 | static double ssd_to_psnr( uint64_t ssd, uint64_t denom ) |
| 155 | { |
| 156 | return -10*log((double)ssd/(denom*255*255))/log(10); |
| 157 | } |
| 158 | |
| 159 | static double ssim_db( double ssim, double weight ) |
| 160 | { |
| 161 | return 10*(log(weight)/log(10)-log(weight-ssim)/log(10)); |
| 162 | } |
| 163 | |
| 164 | static void print_results(uint64_t ssd[3], double ssim[3], int frames, int w, int h) |
| 165 | { |
| 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)", |
| 172 | ssim[0] / frames, |
| 173 | ssim[1] / frames, |
| 174 | ssim[2] / frames, |
| 175 | (ssim[0]*4 + ssim[1] + ssim[2]) / (frames*6), |
| 176 | ssim_db(ssim[0] * 4 + ssim[1] + ssim[2], frames*6)); |
| 177 | } |
| 178 | |
| 179 | int main(int argc, char* argv[]) |
| 180 | { |
| 181 | FILE *f[2]; |
| 182 | uint8_t *buf[2], *plane[2][3]; |
| 183 | int *temp; |
| 184 | uint64_t ssd[3] = {0,0,0}; |
| 185 | double ssim[3] = {0,0,0}; |
| 186 | int frame_size, w, h; |
| 187 | int frames, seek; |
| 188 | int i; |
| 189 | |
| 190 | if( argc<4 || 2 != sscanf(argv[3], "%dx%d", &w, &h) ) |
| 191 | { |
| 192 | printf("tiny_ssim <file1.yuv> <file2.yuv> <width>x<height> [<seek>]\n"); |
| 193 | return -1; |
| 194 | } |
| 195 | |
| 196 | f[0] = fopen(argv[1], "rb"); |
| 197 | f[1] = fopen(argv[2], "rb"); |
| 198 | sscanf(argv[3], "%dx%d", &w, &h); |
| 199 | |
| 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"); |
| 202 | return -2; |
| 203 | } |
| 204 | |
| 205 | frame_size = w*h*3LL/2; |
| 206 | for( i=0; i<2; i++ ) |
| 207 | { |
| 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; |
| 212 | } |
| 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); |
| 216 | |
| 217 | for( frames=0;; frames++ ) |
| 218 | { |
| 219 | uint64_t ssd_one[3]; |
| 220 | double ssim_one[3]; |
| 221 | if( fread(buf[0], frame_size, 1, f[0]) != 1) break; |
| 222 | if( fread(buf[1], frame_size, 1, f[1]) != 1) break; |
| 223 | for( i=0; i<3; i++ ) |
| 224 | { |
| 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, |
| 227 | plane[1][i], w>>!!i, |
| 228 | w>>!!i, h>>!!i, temp, NULL ); |
| 229 | ssd[i] += ssd_one[i]; |
| 230 | ssim[i] += ssim_one[i]; |
| 231 | } |
| 232 | |
| 233 | printf("Frame %d | ", frames); |
| 234 | print_results(ssd_one, ssim_one, 1, w, h); |
| 235 | printf(" \r"); |
| 236 | fflush(stdout); |
| 237 | } |
| 238 | |
| 239 | if( !frames ) return 0; |
| 240 | |
| 241 | printf("Total %d frames | ", frames); |
| 242 | print_results(ssd, ssim, frames, w, h); |
| 243 | printf("\n"); |
| 244 | |
| 245 | return 0; |
| 246 | } |