Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | RSS feed

  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.  
  36. uint64_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()
  62. static 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()
  76. static 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.  
  95. static 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.  
  111. static 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.  
  121. static 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.  
  128. static double get_f64l(uint8_t *p)
  129. {
  130.     return av_int2double(AV_RL64(p));
  131. }
  132.  
  133. static 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.  
  266. int 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\n", argv[3]);
  291.                 return 1;
  292.             }
  293.         }
  294.     }
  295.  
  296.     if (argc < 3) {
  297.         printf("tiny_psnr <file1> <file2> [<elem size> [<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.