2 * Copyright (c) 2003 Michael Niedermayer <michaelni@gmx.at>
4 * This file is part of Libav.
6 * Libav is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2.1 of the License, or (at your option) any later version.
11 * Libav is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with Libav; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
27 #include "libavutil/intfloat.h"
29 #define FFMIN(a, b) ((a) > (b) ? (b) : (a))
33 uint64_t exp16_table
[21] = {
57 // 16.16 fixpoint log()
58 static int64_t log16(uint64_t a
)
64 return -log16((1LL << 32) / a
);
67 for (i
= 20; i
>= 0; i
--) {
68 int64_t b
= exp16_table
[i
];
72 a
= ((a
/ b
) << 16) + (((a
% b
) << 16) + b
/ 2) / b
;
77 static uint64_t int_sqrt(uint64_t a
)
83 for (s
= 31; s
>= 0; s
--) {
84 uint64_t b
= ret_sq
+ (1ULL << (s
* 2)) + (ret
<< s
) * 2;
93 static int16_t get_s16l(uint8_t *p
)
99 v
.u
= p
[0] | p
[1] << 8;
103 static float get_f32l(uint8_t *p
)
105 union av_intfloat32 v
;
106 v
.i
= p
[0] | p
[1] << 8 | p
[2] << 16 | p
[3] << 24;
110 int main(int argc
, char *argv
[])
116 uint8_t buf
[2][SIZE
];
120 int shift
= argc
< 5 ? 0 : atoi(argv
[4]);
121 int skip_bytes
= argc
< 6 ? 0 : atoi(argv
[5]);
127 printf("tiny_psnr <file1> <file2> [<elem size> [<shift> [<skip bytes>]]]\n");
128 printf("WAV headers are skipped automatically.\n");
133 if (!strcmp(argv
[3], "u8")) {
135 } else if (!strcmp(argv
[3], "s16")) {
137 } else if (!strcmp(argv
[3], "f32")) {
141 len
= strtol(argv
[3], &end
, 0);
142 if (*end
|| len
< 1 || len
> 2) {
143 fprintf(stderr
, "Unsupported sample format: %s\n", argv
[3]);
149 max
= (1LL << (8 * len
)) - 1;
151 f
[0] = fopen(argv
[1], "rb");
152 f
[1] = fopen(argv
[2], "rb");
153 if (!f
[0] || !f
[1]) {
154 fprintf(stderr
, "Could not open input files.\n");
158 for (i
= 0; i
< 2; i
++) {
160 if (fread(p
, 1, 12, f
[i
]) != 12)
162 if (!memcmp(p
, "RIFF", 4) &&
163 !memcmp(p
+ 8, "WAVE", 4)) {
164 if (fread(p
, 1, 8, f
[i
]) != 8)
166 while (memcmp(p
, "data", 4)) {
167 int s
= p
[4] | p
[5] << 8 | p
[6] << 16 | p
[7] << 24;
168 fseek(f
[i
], s
, SEEK_CUR
);
169 if (fread(p
, 1, 8, f
[i
]) != 8)
173 fseek(f
[i
], -12, SEEK_CUR
);
177 fseek(f
[shift
< 0], abs(shift
), SEEK_CUR
);
179 fseek(f
[0], skip_bytes
, SEEK_CUR
);
180 fseek(f
[1], skip_bytes
, SEEK_CUR
);
183 int s0
= fread(buf
[0], 1, SIZE
, f
[0]);
184 int s1
= fread(buf
[1], 1, SIZE
, f
[1]);
186 for (j
= 0; j
< FFMIN(s0
, s1
); j
+= len
) {
187 int64_t a
= buf
[0][j
];
188 int64_t b
= buf
[1][j
];
191 a
= get_s16l(buf
[0] + j
);
192 b
= get_s16l(buf
[1] + j
);
193 } else if (len
== 4) {
194 a
= get_f32l(buf
[0] + j
) * (1 << 24);
195 b
= get_f32l(buf
[1] + j
) * (1 << 24);
200 sse
+= (a
- b
) * (a
- b
);
211 i
= FFMIN(size0
, size1
) / len
;
214 dev
= int_sqrt(((sse
/ i
) * F
* F
) + (((sse
% i
) * F
* F
) + i
/ 2) / i
);
216 psnr
= ((2 * log16(max
<< 16) + log16(i
) - log16(sse
)) *
217 284619LL * F
+ (1LL << 31)) / (1LL << 32);
219 psnr
= 1000 * F
- 1; // floating point free infinity :)
221 printf("stddev:%5d.%02d PSNR:%3d.%02d MAXDIFF:%5d bytes:%9d/%9d\n",
222 (int)(dev
/ F
), (int)(dev
% F
),
223 (int)(psnr
/ F
), (int)(psnr
% F
),
224 maxdist
, size0
, size1
);