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 {
00058     int i;
00059     int out = 0;
00060 
00061     if (a < 1 << 16)
00062         return -log16((1LL << 32) / a);
00063     a <<= 16;
00064 
00065     for (i = 20; i >= 0; i--) {
00066         int64_t b = exp16_table[i];
00067         if (a < (b << 16))
00068             continue;
00069         out |= 1 << i;
00070         a    = ((a / b) << 16) + (((a % b) << 16) + b / 2) / b;
00071     }
00072     return out;
00073 }
00074 
00075 static uint64_t int_sqrt(uint64_t a)
00076 {
00077     uint64_t ret    = 0;
00078     uint64_t ret_sq = 0;
00079     int s;
00080 
00081     for (s = 31; s >= 0; s--) {
00082         uint64_t b = ret_sq + (1ULL << (s * 2)) + (ret << s) * 2;
00083         if (b <= a) {
00084             ret_sq = b;
00085             ret   += 1ULL << s;
00086         }
00087     }
00088     return ret;
00089 }
00090 
00091 int main(int argc, char *argv[])
00092 {
00093     int i, j;
00094     uint64_t sse = 0;
00095     uint64_t dev;
00096     FILE *f[2];
00097     uint8_t buf[2][SIZE];
00098     uint64_t psnr;
00099     int len        = argc < 4 ? 1 : atoi(argv[3]);
00100     int64_t max    = (1 << (8 * len)) - 1;
00101     int shift      = argc < 5 ? 0 : atoi(argv[4]);
00102     int skip_bytes = argc < 6 ? 0 : atoi(argv[5]);
00103     int size0      = 0;
00104     int size1      = 0;
00105     int maxdist    = 0;
00106 
00107     if (argc < 3) {
00108         printf("tiny_psnr <file1> <file2> [<elem size> [<shift> [<skip bytes>]]]\n");
00109         printf("WAV headers are skipped automatically.\n");
00110         return 1;
00111     }
00112 
00113     f[0] = fopen(argv[1], "rb");
00114     f[1] = fopen(argv[2], "rb");
00115     if (!f[0] || !f[1]) {
00116         fprintf(stderr, "Could not open input files.\n");
00117         return 1;
00118     }
00119 
00120     for (i = 0; i < 2; i++) {
00121         uint8_t *p = buf[i];
00122         if (fread(p, 1, 12, f[i]) != 12)
00123             return 1;
00124         if (!memcmp(p, "RIFF", 4) &&
00125             !memcmp(p + 8, "WAVE", 4)) {
00126             if (fread(p, 1, 8, f[i]) != 8)
00127                 return 1;
00128             while (memcmp(p, "data", 4)) {
00129                 int s = p[4] | p[5] << 8 | p[6] << 16 | p[7] << 24;
00130                 fseek(f[i], s, SEEK_CUR);
00131                 if (fread(p, 1, 8, f[i]) != 8)
00132                     return 1;
00133             }
00134         } else {
00135             fseek(f[i], -12, SEEK_CUR);
00136         }
00137     }
00138 
00139     fseek(f[shift < 0], abs(shift), SEEK_CUR);
00140 
00141     fseek(f[0], skip_bytes, SEEK_CUR);
00142     fseek(f[1], skip_bytes, SEEK_CUR);
00143 
00144     for (;;) {
00145         int s0 = fread(buf[0], 1, SIZE, f[0]);
00146         int s1 = fread(buf[1], 1, SIZE, f[1]);
00147 
00148         for (j = 0; j < FFMIN(s0, s1); j++) {
00149             int64_t a = buf[0][j];
00150             int64_t b = buf[1][j];
00151             int dist;
00152             if (len == 2) {
00153                 a = (int16_t)(a | (buf[0][++j] << 8));
00154                 b = (int16_t)(b | (buf[1][  j] << 8));
00155             }
00156             sse += (a - b) * (a - b);
00157             dist = abs(a - b);
00158             if (dist > maxdist)
00159                 maxdist = dist;
00160         }
00161         size0 += s0;
00162         size1 += s1;
00163         if (s0 + s1 <= 0)
00164             break;
00165     }
00166 
00167     i = FFMIN(size0, size1) / len;
00168     if (!i)
00169         i = 1;
00170     dev = int_sqrt(((sse / i) * F * F) + (((sse % i) * F * F) + i / 2) / i);
00171     if (sse)
00172         psnr = ((2 * log16(max << 16) + log16(i) - log16(sse)) *
00173                 284619LL * F + (1LL << 31)) / (1LL << 32);
00174     else
00175         psnr = 1000 * F - 1; // floating point free infinity :)
00176 
00177     printf("stddev:%5d.%02d PSNR:%3d.%02d MAXDIFF:%5d bytes:%9d/%9d\n",
00178            (int)(dev / F), (int)(dev % F),
00179            (int)(psnr / F), (int)(psnr % F),
00180            maxdist, size0, size1);
00181     return 0;
00182 }