00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
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
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;
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 }