signed 16bit support
[libav.git] / tests / tiny_psnr.c
CommitLineData
67cbe681
MN
1/*
2 * Copyright (c) 2003 Michael Niedermayer <michaelni@gmx.at>
3 *
4 * This library is free software; you can redistribute it and/or
5 * modify it under the terms of the GNU Lesser General Public
6 * License as published by the Free Software Foundation; either
7 * version 2 of the License, or (at your option) any later version.
8 *
9 * This library 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 GNU
12 * Lesser General Public License for more details.
13 *
14 * You should have received a copy of the GNU Lesser General Public
15 * License along with this library; if not, write to the Free Software
16 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 *
18 */
19
20#include <stdio.h>
21#include <inttypes.h>
9e2a16e1 22#include <assert.h>
67cbe681
MN
23
24#define F 100
25#define SIZE 2048
26
0c90161f 27uint64_t exp16_table[21]={
9e2a16e1
MN
28 65537,
29 65538,
30 65540,
31 65544,
32 65552,
33 65568,
34 65600,
35 65664,
36 65793,
37 66050,
38 66568,
39 67616,
40 69763,
41 74262,
42 84150,
43 108051,
44 178145,
45 484249,
46 3578144,
47 195360063,
0c90161f 48 582360139072LL,
9e2a16e1
MN
49};
50#if 1
51// 16.16 fixpoint exp()
52static unsigned int exp16(unsigned int a){
53 int i;
54 int out= 1<<16;
55
56 for(i=19;i>=0;i--){
57 if(a&(1<<i))
58 out= (out*exp16_table[i] + (1<<15))>>16;
59 }
60
61 return out;
62}
63// 16.16 fixpoint log()
0c90161f 64static int64_t log16(uint64_t a){
9e2a16e1
MN
65 int i;
66 int out=0;
0c90161f
MN
67
68 if(a < 1<<16)
69 return -log16((1LL<<32) / a);
9e2a16e1
MN
70 a<<=16;
71
0c90161f 72 for(i=20;i>=0;i--){
8176bd1a
MN
73 int64_t b= exp16_table[i];
74 if(a<(b<<16)) continue;
9e2a16e1 75 out |= 1<<i;
8176bd1a 76 a = ((a/b)<<16) + (((a%b)<<16) + b/2)/b;
9e2a16e1
MN
77 }
78 return out;
79}
80
81#endif
67cbe681
MN
82static uint64_t int_sqrt(uint64_t a)
83{
84 uint64_t ret=0;
85 int s;
86 uint64_t ret_sq=0;
87
88 for(s=31; s>=0; s--){
89 uint64_t b= ret_sq + (1ULL<<(s*2)) + (ret<<s)*2;
90 if(b<=a){
91 ret_sq=b;
92 ret+= 1ULL<<s;
93 }
94 }
95 return ret;
96}
97
98int main(int argc,char* argv[]){
99 int i, j;
100 uint64_t sse=0;
101 uint64_t dev;
102 FILE *f[2];
103 uint8_t buf[2][SIZE];
9e2a16e1 104 uint64_t psnr;
0c90161f
MN
105 int len= argc<4 ? 1 : 2;
106 int64_t max= (1<<(8*len))-1;
67cbe681 107
0c90161f 108 if(argc<3){
67cbe681
MN
109 printf("tiny_psnr <file1> <file2>\n");
110 return -1;
111 }
112
0d6d0cf9
MN
113 f[0]= fopen(argv[1], "rb");
114 f[1]= fopen(argv[2], "rb");
67cbe681
MN
115
116 for(i=0;;){
117 if( fread(buf[0], SIZE, 1, f[0]) != 1) break;
118 if( fread(buf[1], SIZE, 1, f[1]) != 1) break;
119
120 for(j=0; j<SIZE; i++,j++){
0c90161f
MN
121 int64_t a= buf[0][j];
122 int64_t b= buf[1][j];
123 if(len==2){
124 a= (int16_t)(a | (buf[0][++j]<<8));
125 b= (int16_t)(b | (buf[1][ j]<<8));
126 }
67cbe681
MN
127 sse += (a-b) * (a-b);
128 }
129 }
130
eeaa742c 131 if(!i) i=1;
0c90161f 132 dev= int_sqrt( ((sse/i)*F*F) + (((sse%i)*F*F) + i/2)/i );
9e2a16e1 133 if(sse)
0c90161f 134 psnr= ((2*log16(max<<16) + log16(i) - log16(sse))*284619LL*F + (1<<31)) / (1LL<<32);
9e2a16e1
MN
135 else
136 psnr= 100*F-1; //floating point free infinity :)
67cbe681 137
9e2a16e1
MN
138 printf("stddev:%3d.%02d PSNR:%2d.%02d bytes:%d\n",
139 (int)(dev/F), (int)(dev%F),
140 (int)(psnr/F), (int)(psnr%F),
141 i);
67cbe681
MN
142 return 0;
143}
9e2a16e1
MN
144
145