Libav
|
00001 /* 00002 * Copyright (c) 2003 Michael Niedermayer <michaelni@gmx.at> 00003 * 00004 * This file is part of FFmpeg. 00005 * 00006 * FFmpeg is free software; you can redistribute it and/or 00007 * modify it under the terms of the GNU Lesser General Public 00008 * License as published by the Free Software Foundation; either 00009 * version 2.1 of the License, or (at your option) any later version. 00010 * 00011 * FFmpeg is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00014 * Lesser General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU Lesser General Public 00017 * License along with FFmpeg; if not, write to the Free Software 00018 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00019 */ 00020 00021 #include <stdio.h> 00022 #include <stdlib.h> 00023 #include <inttypes.h> 00024 #include <assert.h> 00025 00026 #define FFMIN(a,b) ((a) > (b) ? (b) : (a)) 00027 #define F 100 00028 #define SIZE 2048 00029 00030 uint64_t exp16_table[21]={ 00031 65537, 00032 65538, 00033 65540, 00034 65544, 00035 65552, 00036 65568, 00037 65600, 00038 65664, 00039 65793, 00040 66050, 00041 66568, 00042 67616, 00043 69763, 00044 74262, 00045 84150, 00046 108051, 00047 178145, 00048 484249, 00049 3578144, 00050 195360063, 00051 582360139072LL, 00052 }; 00053 00054 #if 0 00055 // 16.16 fixpoint exp() 00056 static unsigned int exp16(unsigned int a){ 00057 int i; 00058 int out= 1<<16; 00059 00060 for(i=19;i>=0;i--){ 00061 if(a&(1<<i)) 00062 out= (out*exp16_table[i] + (1<<15))>>16; 00063 } 00064 00065 return out; 00066 } 00067 #endif 00068 00069 // 16.16 fixpoint log() 00070 static int64_t log16(uint64_t a){ 00071 int i; 00072 int out=0; 00073 00074 if(a < 1<<16) 00075 return -log16((1LL<<32) / a); 00076 a<<=16; 00077 00078 for(i=20;i>=0;i--){ 00079 int64_t b= exp16_table[i]; 00080 if(a<(b<<16)) continue; 00081 out |= 1<<i; 00082 a = ((a/b)<<16) + (((a%b)<<16) + b/2)/b; 00083 } 00084 return out; 00085 } 00086 00087 static uint64_t int_sqrt(uint64_t a) 00088 { 00089 uint64_t ret=0; 00090 int s; 00091 uint64_t ret_sq=0; 00092 00093 for(s=31; s>=0; s--){ 00094 uint64_t b= ret_sq + (1ULL<<(s*2)) + (ret<<s)*2; 00095 if(b<=a){ 00096 ret_sq=b; 00097 ret+= 1ULL<<s; 00098 } 00099 } 00100 return ret; 00101 } 00102 00103 int main(int argc,char* argv[]){ 00104 int i, j; 00105 uint64_t sse=0; 00106 uint64_t dev; 00107 FILE *f[2]; 00108 uint8_t buf[2][SIZE]; 00109 uint64_t psnr; 00110 int len= argc<4 ? 1 : atoi(argv[3]); 00111 int64_t max= (1<<(8*len))-1; 00112 int shift= argc<5 ? 0 : atoi(argv[4]); 00113 int skip_bytes = argc<6 ? 0 : atoi(argv[5]); 00114 int size0=0; 00115 int size1=0; 00116 00117 if(argc<3){ 00118 printf("tiny_psnr <file1> <file2> [<elem size> [<shift> [<skip bytes>]]]\n"); 00119 printf("For WAV files use the following:\n"); 00120 printf("./tiny_psnr file1.wav file2.wav 2 0 44 to skip the header.\n"); 00121 return -1; 00122 } 00123 00124 f[0]= fopen(argv[1], "rb"); 00125 f[1]= fopen(argv[2], "rb"); 00126 if(!f[0] || !f[1]){ 00127 fprintf(stderr, "Could not open input files.\n"); 00128 return -1; 00129 } 00130 fseek(f[shift<0], shift < 0 ? -shift : shift, SEEK_SET); 00131 00132 fseek(f[0],skip_bytes,SEEK_CUR); 00133 fseek(f[1],skip_bytes,SEEK_CUR); 00134 00135 for(;;){ 00136 int s0= fread(buf[0], 1, SIZE, f[0]); 00137 int s1= fread(buf[1], 1, SIZE, f[1]); 00138 00139 for(j=0; j<FFMIN(s0,s1); j++){ 00140 int64_t a= buf[0][j]; 00141 int64_t b= buf[1][j]; 00142 if(len==2){ 00143 a= (int16_t)(a | (buf[0][++j]<<8)); 00144 b= (int16_t)(b | (buf[1][ j]<<8)); 00145 } 00146 sse += (a-b) * (a-b); 00147 } 00148 size0 += s0; 00149 size1 += s1; 00150 if(s0+s1<=0) 00151 break; 00152 } 00153 00154 i= FFMIN(size0,size1)/len; 00155 if(!i) i=1; 00156 dev= int_sqrt( ((sse/i)*F*F) + (((sse%i)*F*F) + i/2)/i ); 00157 if(sse) 00158 psnr= ((2*log16(max<<16) + log16(i) - log16(sse))*284619LL*F + (1LL<<31)) / (1LL<<32); 00159 else 00160 psnr= 1000*F-1; //floating point free infinity :) 00161 00162 printf("stddev:%5d.%02d PSNR:%3d.%02d bytes:%9d/%9d\n", 00163 (int)(dev/F), (int)(dev%F), 00164 (int)(psnr/F), (int)(psnr%F), 00165 size0, size1); 00166 return 0; 00167 } 00168 00169