1/* 2 * Copyright (c) 2003 Michael Niedermayer <michaelni@gmx.at> 3 * 4 * This file is part of Libav. 5 * 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. 10 * 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. 15 * 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 19 */ 20 21#include <stdio.h> 22#include <stdlib.h> 23#include <string.h> 24#include <inttypes.h> 25#include <assert.h> 26 27#define FFMIN(a, b) ((a) > (b) ? (b) : (a)) 28#define F 100 29#define SIZE 2048 30 31uint64_t exp16_table[21] = { 32 65537, 33 65538, 34 65540, 35 65544, 36 65552, 37 65568, 38 65600, 39 65664, 40 65793, 41 66050, 42 66568, 43 67616, 44 69763, 45 74262, 46 84150, 47 108051, 48 178145, 49 484249, 50 3578144, 51 195360063, 52 582360139072LL, 53}; 54 55// 16.16 fixpoint log() 56static int64_t log16(uint64_t a) 57{ 58 int i; 59 int out = 0; 60 61 if (a < 1 << 16) 62 return -log16((1LL << 32) / a); 63 a <<= 16; 64 65 for (i = 20; i >= 0; i--) { 66 int64_t b = exp16_table[i]; 67 if (a < (b << 16)) 68 continue; 69 out |= 1 << i; 70 a = ((a / b) << 16) + (((a % b) << 16) + b / 2) / b; 71 } 72 return out; 73} 74 75static uint64_t int_sqrt(uint64_t a) 76{ 77 uint64_t ret = 0; 78 uint64_t ret_sq = 0; 79 int s; 80 81 for (s = 31; s >= 0; s--) { 82 uint64_t b = ret_sq + (1ULL << (s * 2)) + (ret << s) * 2; 83 if (b <= a) { 84 ret_sq = b; 85 ret += 1ULL << s; 86 } 87 } 88 return ret; 89} 90 91int main(int argc, char *argv[]) 92{ 93 int i, j; 94 uint64_t sse = 0; 95 uint64_t dev; 96 FILE *f[2]; 97 uint8_t buf[2][SIZE]; 98 uint64_t psnr; 99 int len = argc < 4 ? 1 : atoi(argv[3]); 100 int64_t max = (1 << (8 * len)) - 1; 101 int shift = argc < 5 ? 0 : atoi(argv[4]); 102 int skip_bytes = argc < 6 ? 0 : atoi(argv[5]); 103 int size0 = 0; 104 int size1 = 0; 105 int maxdist = 0; 106 107 if (argc < 3) { 108 printf("tiny_psnr <file1> <file2> [<elem size> [<shift> [<skip bytes>]]]\n"); 109 printf("WAV headers are skipped automatically.\n"); 110 return 1; 111 } 112 113 f[0] = fopen(argv[1], "rb"); 114 f[1] = fopen(argv[2], "rb"); 115 if (!f[0] || !f[1]) { 116 fprintf(stderr, "Could not open input files.\n"); 117 return 1; 118 } 119 120 for (i = 0; i < 2; i++) { 121 uint8_t *p = buf[i]; 122 if (fread(p, 1, 12, f[i]) != 12) 123 return 1; 124 if (!memcmp(p, "RIFF", 4) && 125 !memcmp(p + 8, "WAVE", 4)) { 126 if (fread(p, 1, 8, f[i]) != 8) 127 return 1; 128 while (memcmp(p, "data", 4)) { 129 int s = p[4] | p[5] << 8 | p[6] << 16 | p[7] << 24; 130 fseek(f[i], s, SEEK_CUR); 131 if (fread(p, 1, 8, f[i]) != 8) 132 return 1; 133 } 134 } else { 135 fseek(f[i], -12, SEEK_CUR); 136 } 137 } 138 139 fseek(f[shift < 0], abs(shift), SEEK_CUR); 140 141 fseek(f[0], skip_bytes, SEEK_CUR); 142 fseek(f[1], skip_bytes, SEEK_CUR); 143 144 for (;;) { 145 int s0 = fread(buf[0], 1, SIZE, f[0]); 146 int s1 = fread(buf[1], 1, SIZE, f[1]); 147 148 for (j = 0; j < FFMIN(s0, s1); j++) { 149 int64_t a = buf[0][j]; 150 int64_t b = buf[1][j]; 151 int dist; 152 if (len == 2) { 153 a = (int16_t)(a | (buf[0][++j] << 8)); 154 b = (int16_t)(b | (buf[1][ j] << 8)); 155 } 156 sse += (a - b) * (a - b); 157 dist = abs(a - b); 158 if (dist > maxdist) 159 maxdist = dist; 160 } 161 size0 += s0; 162 size1 += s1; 163 if (s0 + s1 <= 0) 164 break; 165 } 166 167 i = FFMIN(size0, size1) / len; 168 if (!i) 169 i = 1; 170 dev = int_sqrt(((sse / i) * F * F) + (((sse % i) * F * F) + i / 2) / i); 171 if (sse) 172 psnr = ((2 * log16(max << 16) + log16(i) - log16(sse)) * 173 284619LL * F + (1LL << 31)) / (1LL << 32); 174 else 175 psnr = 1000 * F - 1; // floating point free infinity :) 176 177 printf("stddev:%5d.%02d PSNR:%3d.%02d MAXDIFF:%5d bytes:%9d/%9d\n", 178 (int)(dev / F), (int)(dev % F), 179 (int)(psnr / F), (int)(psnr % F), 180 maxdist, size0, size1); 181 return 0; 182} 183