Libav 0.7.1
tests/tiny_psnr.c
Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 2003 Michael Niedermayer <michaelni@gmx.at>
00003  *
00004  * This file is part of Libav.
00005  *
00006  * Libav 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  * Libav 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 Libav; 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 <string.h>
00024 #include <inttypes.h>
00025 #include <assert.h>
00026 
00027 #define FFMIN(a,b) ((a) > (b) ? (b) : (a))
00028 #define F 100
00029 #define SIZE 2048
00030 
00031 uint64_t exp16_table[21]={
00032      65537,
00033      65538,
00034      65540,
00035      65544,
00036      65552,
00037      65568,
00038      65600,
00039      65664,
00040      65793,
00041      66050,
00042      66568,
00043      67616,
00044      69763,
00045      74262,
00046      84150,
00047     108051,
00048     178145,
00049     484249,
00050    3578144,
00051  195360063,
00052  582360139072LL,
00053 };
00054 
00055 // 16.16 fixpoint log()
00056 static int64_t log16(uint64_t a){
00057     int i;
00058     int out=0;
00059 
00060     if(a < 1<<16)
00061         return -log16((1LL<<32) / a);
00062     a<<=16;
00063 
00064     for(i=20;i>=0;i--){
00065         int64_t b= exp16_table[i];
00066         if(a<(b<<16)) continue;
00067         out |= 1<<i;
00068         a = ((a/b)<<16) + (((a%b)<<16) + b/2)/b;
00069     }
00070     return out;
00071 }
00072 
00073 static uint64_t int_sqrt(uint64_t a)
00074 {
00075     uint64_t ret=0;
00076     int s;
00077     uint64_t ret_sq=0;
00078 
00079     for(s=31; s>=0; s--){
00080         uint64_t b= ret_sq + (1ULL<<(s*2)) + (ret<<s)*2;
00081         if(b<=a){
00082             ret_sq=b;
00083             ret+= 1ULL<<s;
00084         }
00085     }
00086     return ret;
00087 }
00088 
00089 int main(int argc,char* argv[]){
00090     int i, j;
00091     uint64_t sse=0;
00092     uint64_t dev;
00093     FILE *f[2];
00094     uint8_t buf[2][SIZE];
00095     uint64_t psnr;
00096     int len= argc<4 ? 1 : atoi(argv[3]);
00097     int64_t max= (1<<(8*len))-1;
00098     int shift= argc<5 ? 0 : atoi(argv[4]);
00099     int skip_bytes = argc<6 ? 0 : atoi(argv[5]);
00100     int size0=0;
00101     int size1=0;
00102     int maxdist = 0;
00103 
00104     if(argc<3){
00105         printf("tiny_psnr <file1> <file2> [<elem size> [<shift> [<skip bytes>]]]\n");
00106         printf("WAV headers are skipped automatically.\n");
00107         return 1;
00108     }
00109 
00110     f[0]= fopen(argv[1], "rb");
00111     f[1]= fopen(argv[2], "rb");
00112     if(!f[0] || !f[1]){
00113         fprintf(stderr, "Could not open input files.\n");
00114         return 1;
00115     }
00116 
00117     for (i = 0; i < 2; i++) {
00118         uint8_t *p = buf[i];
00119         if (fread(p, 1, 12, f[i]) != 12)
00120             return 1;
00121         if (!memcmp(p,   "RIFF", 4) &&
00122             !memcmp(p+8, "WAVE", 4)) {
00123             if (fread(p, 1, 8, f[i]) != 8)
00124                 return 1;
00125             while (memcmp(p, "data", 4)) {
00126                 int s = p[4] | p[5]<<8 | p[6]<<16 | p[7]<<24;
00127                 fseek(f[i], s, SEEK_CUR);
00128                 if (fread(p, 1, 8, f[i]) != 8)
00129                     return 1;
00130             }
00131         } else {
00132             fseek(f[i], -12, SEEK_CUR);
00133         }
00134     }
00135 
00136     fseek(f[shift<0], abs(shift), SEEK_CUR);
00137 
00138     fseek(f[0],skip_bytes,SEEK_CUR);
00139     fseek(f[1],skip_bytes,SEEK_CUR);
00140 
00141     for(;;){
00142         int s0= fread(buf[0], 1, SIZE, f[0]);
00143         int s1= fread(buf[1], 1, SIZE, f[1]);
00144 
00145         for(j=0; j<FFMIN(s0,s1); j++){
00146             int64_t a= buf[0][j];
00147             int64_t b= buf[1][j];
00148             int dist;
00149             if(len==2){
00150                 a= (int16_t)(a | (buf[0][++j]<<8));
00151                 b= (int16_t)(b | (buf[1][  j]<<8));
00152             }
00153             sse += (a-b) * (a-b);
00154             dist = abs(a-b);
00155             if (dist > maxdist) maxdist = dist;
00156         }
00157         size0 += s0;
00158         size1 += s1;
00159         if(s0+s1<=0)
00160             break;
00161     }
00162 
00163     i= FFMIN(size0,size1)/len;
00164     if(!i) i=1;
00165     dev= int_sqrt( ((sse/i)*F*F) + (((sse%i)*F*F) + i/2)/i );
00166     if(sse)
00167         psnr= ((2*log16(max<<16) + log16(i) - log16(sse))*284619LL*F + (1LL<<31)) / (1LL<<32);
00168     else
00169         psnr= 1000*F-1; //floating point free infinity :)
00170 
00171     printf("stddev:%5d.%02d PSNR:%3d.%02d MAXDIFF:%5d bytes:%9d/%9d\n",
00172         (int)(dev/F), (int)(dev%F),
00173         (int)(psnr/F), (int)(psnr%F),
00174         maxdist,
00175         size0, size1);
00176     return 0;
00177 }
00178 
00179