1/* 2 * Copyright (c) 2003 Michael Niedermayer <michaelni@gmx.at> 3 * 4 * This file is part of FFmpeg. 5 * 6 * FFmpeg 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 * FFmpeg 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 FFmpeg; 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 <math.h> 26#include <float.h> 27#include <limits.h> 28 29#include "libavutil/intfloat.h" 30#include "libavutil/intreadwrite.h" 31 32#define FFMIN(a, b) ((a) > (b) ? (b) : (a)) 33#define F 100 34#define SIZE 2048 35 36uint64_t exp16_table[21] = { 37 65537, 38 65538, 39 65540, 40 65544, 41 65552, 42 65568, 43 65600, 44 65664, 45 65793, 46 66050, 47 66568, 48 67616, 49 69763, 50 74262, 51 84150, 52 108051, 53 178145, 54 484249, 55 3578144, 56 195360063, 57 582360139072LL, 58}; 59 60#if 0 61// 16.16 fixpoint exp() 62static unsigned int exp16(unsigned int a){ 63 int i; 64 int out= 1<<16; 65 66 for(i=19;i>=0;i--){ 67 if(a&(1<<i)) 68 out= (out*exp16_table[i] + (1<<15))>>16; 69 } 70 71 return out; 72} 73#endif 74 75// 16.16 fixpoint log() 76static int64_t log16(uint64_t a) 77{ 78 int i; 79 int out = 0; 80 81 if (a < 1 << 16) 82 return -log16((1LL << 32) / a); 83 a <<= 16; 84 85 for (i = 20; i >= 0; i--) { 86 int64_t b = exp16_table[i]; 87 if (a < (b << 16)) 88 continue; 89 out |= 1 << i; 90 a = ((a / b) << 16) + (((a % b) << 16) + b / 2) / b; 91 } 92 return out; 93} 94 95static uint64_t int_sqrt(uint64_t a) 96{ 97 uint64_t ret = 0; 98 uint64_t ret_sq = 0; 99 int s; 100 101 for (s = 31; s >= 0; s--) { 102 uint64_t b = ret_sq + (1ULL << (s * 2)) + (ret << s) * 2; 103 if (b <= a) { 104 ret_sq = b; 105 ret += 1ULL << s; 106 } 107 } 108 return ret; 109} 110 111static int16_t get_s16l(uint8_t *p) 112{ 113 union { 114 uint16_t u; 115 int16_t s; 116 } v; 117 v.u = p[0] | p[1] << 8; 118 return v.s; 119} 120 121static float get_f32l(uint8_t *p) 122{ 123 union av_intfloat32 v; 124 v.i = p[0] | p[1] << 8 | p[2] << 16 | p[3] << 24; 125 return v.f; 126} 127 128static double get_f64l(uint8_t *p) 129{ 130 return av_int2double(AV_RL64(p)); 131} 132 133static int run_psnr(FILE *f[2], int len, int shift, int skip_bytes) 134{ 135 int i, j; 136 uint64_t sse = 0; 137 double sse_d = 0.0; 138 uint8_t buf[2][SIZE]; 139 int64_t max = (1LL << (8 * len)) - 1; 140 int size0 = 0; 141 int size1 = 0; 142 uint64_t maxdist = 0; 143 double maxdist_d = 0.0; 144 int noseek; 145 146 noseek = fseek(f[0], 0, SEEK_SET) || 147 fseek(f[1], 0, SEEK_SET); 148 149 if (!noseek) { 150 for (i = 0; i < 2; i++) { 151 uint8_t *p = buf[i]; 152 if (fread(p, 1, 12, f[i]) != 12) 153 return 1; 154 if (!memcmp(p, "RIFF", 4) && 155 !memcmp(p + 8, "WAVE", 4)) { 156 if (fread(p, 1, 8, f[i]) != 8) 157 return 1; 158 while (memcmp(p, "data", 4)) { 159 int s = p[4] | p[5] << 8 | p[6] << 16 | p[7] << 24; 160 fseek(f[i], s, SEEK_CUR); 161 if (fread(p, 1, 8, f[i]) != 8) 162 return 1; 163 } 164 } else { 165 fseek(f[i], -12, SEEK_CUR); 166 } 167 } 168 169 fseek(f[shift < 0], abs(shift), SEEK_CUR); 170 171 fseek(f[0], skip_bytes, SEEK_CUR); 172 fseek(f[1], skip_bytes, SEEK_CUR); 173 } 174 175 for (;;) { 176 int s0 = fread(buf[0], 1, SIZE, f[0]); 177 int s1 = fread(buf[1], 1, SIZE, f[1]); 178 179 for (j = 0; j < FFMIN(s0, s1); j += len) { 180 switch (len) { 181 case 1: 182 case 2: { 183 int64_t a = buf[0][j]; 184 int64_t b = buf[1][j]; 185 int dist; 186 if (len == 2) { 187 a = get_s16l(buf[0] + j); 188 b = get_s16l(buf[1] + j); 189 } else { 190 a = buf[0][j]; 191 b = buf[1][j]; 192 } 193 sse += (a - b) * (a - b); 194 dist = abs(a - b); 195 if (dist > maxdist) 196 maxdist = dist; 197 break; 198 } 199 case 4: 200 case 8: { 201 double dist, a, b; 202 if (len == 8) { 203 a = get_f64l(buf[0] + j); 204 b = get_f64l(buf[1] + j); 205 } else { 206 a = get_f32l(buf[0] + j); 207 b = get_f32l(buf[1] + j); 208 } 209 dist = fabs(a - b); 210 sse_d += (a - b) * (a - b); 211 if (dist > maxdist_d) 212 maxdist_d = dist; 213 break; 214 } 215 } 216 } 217 size0 += s0; 218 size1 += s1; 219 if (s0 + s1 <= 0) 220 break; 221 } 222 223 i = FFMIN(size0, size1) / len; 224 if (!i) 225 i = 1; 226 switch (len) { 227 case 1: 228 case 2: { 229 uint64_t psnr; 230 uint64_t dev = int_sqrt(((sse / i) * F * F) + (((sse % i) * F * F) + i / 2) / i); 231 if (sse) 232 psnr = ((2 * log16(max << 16) + log16(i) - log16(sse)) * 233 284619LL * F + (1LL << 31)) / (1LL << 32); 234 else 235 psnr = 1000 * F - 1; // floating point free infinity :) 236 237 printf("stddev:%5d.%02d PSNR:%3d.%02d MAXDIFF:%5"PRIu64" bytes:%9d/%9d\n", 238 (int)(dev / F), (int)(dev % F), 239 (int)(psnr / F), (int)(psnr % F), 240 maxdist, size0, size1); 241 return psnr; 242 } 243 case 4: 244 case 8: { 245 char psnr_str[64]; 246 double psnr = INT_MAX; 247 double dev = sqrt(sse_d / i); 248 uint64_t scale = (len == 4) ? (1ULL << 24) : (1ULL << 32); 249 250 if (sse_d) { 251 psnr = 2 * log(DBL_MAX) - log(i / sse_d); 252 snprintf(psnr_str, sizeof(psnr_str), "%5.02f", psnr); 253 } else 254 snprintf(psnr_str, sizeof(psnr_str), "inf"); 255 256 maxdist = maxdist_d * scale; 257 258 printf("stddev:%10.2f PSNR:%s MAXDIFF:%10"PRIu64" bytes:%9d/%9d\n", 259 dev * scale, psnr_str, maxdist, size0, size1); 260 return psnr; 261 } 262 } 263 return -1; 264} 265 266int main(int argc, char *argv[]) 267{ 268 FILE *f[2]; 269 int len = 1; 270 int shift_first= argc < 5 ? 0 : atoi(argv[4]); 271 int skip_bytes = argc < 6 ? 0 : atoi(argv[5]); 272 int shift_last = shift_first + (argc < 7 ? 0 : atoi(argv[6])); 273 int shift; 274 int max_psnr = -1; 275 int max_psnr_shift = 0; 276 277 if (argc > 3) { 278 if (!strcmp(argv[3], "u8")) { 279 len = 1; 280 } else if (!strcmp(argv[3], "s16")) { 281 len = 2; 282 } else if (!strcmp(argv[3], "f32")) { 283 len = 4; 284 } else if (!strcmp(argv[3], "f64")) { 285 len = 8; 286 } else { 287 char *end; 288 len = strtol(argv[3], &end, 0); 289 if (*end || len < 1 || len > 2) { 290 fprintf(stderr, "Unsupported sample format: %s\nSupported: u8, s16, f32, f64\n", argv[3]); 291 return 1; 292 } 293 } 294 } 295 296 if (argc < 3) { 297 printf("tiny_psnr <file1> <file2> [<elem size>|u8|s16|f32|f64 [<shift> [<skip bytes> [<shift search range>]]]]\n"); 298 printf("WAV headers are skipped automatically.\n"); 299 return 1; 300 } 301 302 f[0] = fopen(argv[1], "rb"); 303 f[1] = fopen(argv[2], "rb"); 304 if (!f[0] || !f[1]) { 305 fprintf(stderr, "Could not open input files.\n"); 306 return 1; 307 } 308 309 for (shift = shift_first; shift <= shift_last; shift++) { 310 int psnr = run_psnr(f, len, shift, skip_bytes); 311 if (psnr > max_psnr || (shift < 0 && psnr == max_psnr)) { 312 max_psnr = psnr; 313 max_psnr_shift = shift; 314 } 315 } 316 if (shift_last > shift_first) 317 printf("Best PSNR is %3d.%02d for shift %i\n", (int)(max_psnr / F), (int)(max_psnr % F), max_psnr_shift); 318 return 0; 319} 320