libavcodec/dct-test.c

   1 /*
   2  * (c) 2001 Fabrice Bellard
   3  *     2007 Marc Hoffman <marc.hoffman@analog.com>
   4  *
   5  * This file is part of Libav.
   6  *
   7  * Libav is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU Lesser General Public
   9  * License as published by the Free Software Foundation; either
  10  * version 2.1 of the License, or (at your option) any later version.
  11  *
  12  * Libav is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15  * Lesser General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU Lesser General Public
  18  * License along with Libav; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  20  */
  21
  22 /**
  23  * @file
  24  * DCT test (c) 2001 Fabrice Bellard
  25  * Started from sample code by Juan J. Sierralta P.
  26  */
  27
  28 #include <stdlib.h>
  29 #include <stdio.h>
  30 #include <string.h>
  31 #include <sys/time.h>
  32 #include <unistd.h>
  33 #include <math.h>
  34
  35 #include "libavutil/cpu.h"
  36 #include "libavutil/common.h"
  37 #include "libavutil/lfg.h"
  38
  39 #include "simple_idct.h"
  40 #include "aandcttab.h"
  41 #include "faandct.h"
  42 #include "faanidct.h"
  43 #include "x86/idct_xvid.h"
  44 #include "dctref.h"
  45
  46 #undef printf
  47
  48 void ff_mmx_idct(DCTELEM *data);
  49 void ff_mmxext_idct(DCTELEM *data);
  50
  51 void odivx_idct_c(short *block);
  52
  53 // BFIN
  54 void ff_bfin_idct(DCTELEM *block);
  55 void ff_bfin_fdct(DCTELEM *block);
  56
  57 // ALTIVEC
  58 void fdct_altivec(DCTELEM *block);
  59 //void idct_altivec(DCTELEM *block);?? no routine
  60
  61 // ARM
  62 void ff_j_rev_dct_arm(DCTELEM *data);
  63 void ff_simple_idct_arm(DCTELEM *data);
  64 void ff_simple_idct_armv5te(DCTELEM *data);
  65 void ff_simple_idct_armv6(DCTELEM *data);
  66 void ff_simple_idct_neon(DCTELEM *data);
  67
  68 void ff_simple_idct_axp(DCTELEM *data);
  69
  70 struct algo {
  71     const char *name;
  72     enum { FDCT, IDCT } is_idct;
  73     void (*func)(DCTELEM *block);
  74     void (*ref) (DCTELEM *block);
  75     enum formattag { NO_PERM, MMX_PERM, MMX_SIMPLE_PERM, SCALE_PERM,
  76                      SSE2_PERM, PARTTRANS_PERM } format;
  77     int mm_support;
  78 };
  79
  80 #ifndef FAAN_POSTSCALE
  81 #define FAAN_SCALE SCALE_PERM
  82 #else
  83 #define FAAN_SCALE NO_PERM
  84 #endif
  85
  86 static int cpu_flags;
  87
  88 struct algo algos[] = {
  89     {"REF-DBL",         0, ff_ref_fdct,        ff_ref_fdct, NO_PERM},
  90     {"FAAN",            0, ff_faandct,         ff_ref_fdct, FAAN_SCALE},
  91     {"FAANI",           1, ff_faanidct,        ff_ref_idct, NO_PERM},
  92     {"IJG-AAN-INT",     0, fdct_ifast,         ff_ref_fdct, SCALE_PERM},
  93     {"IJG-LLM-INT",     0, ff_jpeg_fdct_islow, ff_ref_fdct, NO_PERM},
  94     {"REF-DBL",         1, ff_ref_idct,        ff_ref_idct, NO_PERM},
  95     {"INT",             1, j_rev_dct,          ff_ref_idct, MMX_PERM},
  96     {"SIMPLE-C",        1, ff_simple_idct,     ff_ref_idct, NO_PERM},
  97
  98 #if HAVE_MMX
  99     {"MMX",             0, ff_fdct_mmx,        ff_ref_fdct, NO_PERM, AV_CPU_FLAG_MMX},
 100 #if HAVE_MMX2
 101     {"MMX2",            0, ff_fdct_mmx2,       ff_ref_fdct, NO_PERM, AV_CPU_FLAG_MMX2},
 102     {"SSE2",            0, ff_fdct_sse2,       ff_ref_fdct, NO_PERM, AV_CPU_FLAG_SSE2},
 103 #endif
 104
 105 #if CONFIG_GPL
 106     {"LIBMPEG2-MMX",    1, ff_mmx_idct,        ff_ref_idct, MMX_PERM, AV_CPU_FLAG_MMX},
 107     {"LIBMPEG2-MMX2",   1, ff_mmxext_idct,     ff_ref_idct, MMX_PERM, AV_CPU_FLAG_MMX2},
 108 #endif
 109     {"SIMPLE-MMX",      1, ff_simple_idct_mmx, ff_ref_idct, MMX_SIMPLE_PERM, AV_CPU_FLAG_MMX},
 110     {"XVID-MMX",        1, ff_idct_xvid_mmx,   ff_ref_idct, NO_PERM, AV_CPU_FLAG_MMX},
 111     {"XVID-MMX2",       1, ff_idct_xvid_mmx2,  ff_ref_idct, NO_PERM, AV_CPU_FLAG_MMX2},
 112     {"XVID-SSE2",       1, ff_idct_xvid_sse2,  ff_ref_idct, SSE2_PERM, AV_CPU_FLAG_SSE2},
 113 #endif
 114
 115 #if HAVE_ALTIVEC
 116     {"altivecfdct",     0, fdct_altivec,       ff_ref_fdct, NO_PERM, AV_CPU_FLAG_ALTIVEC},
 117 #endif
 118
 119 #if ARCH_BFIN
 120     {"BFINfdct",        0, ff_bfin_fdct,       ff_ref_fdct, NO_PERM},
 121     {"BFINidct",        1, ff_bfin_idct,       ff_ref_idct, NO_PERM},
 122 #endif
 123
 124 #if ARCH_ARM
 125     {"SIMPLE-ARM",      1, ff_simple_idct_arm, ff_ref_idct, NO_PERM },
 126     {"INT-ARM",         1, ff_j_rev_dct_arm,   ff_ref_idct, MMX_PERM },
 127 #if HAVE_ARMV5TE
 128     {"SIMPLE-ARMV5TE",  1, ff_simple_idct_armv5te, ff_ref_idct, NO_PERM },
 129 #endif
 130 #if HAVE_ARMV6
 131     {"SIMPLE-ARMV6",    1, ff_simple_idct_armv6, ff_ref_idct, MMX_PERM },
 132 #endif
 133 #if HAVE_NEON
 134     {"SIMPLE-NEON",     1, ff_simple_idct_neon, ff_ref_idct, PARTTRANS_PERM },
 135 #endif
 136 #endif /* ARCH_ARM */
 137
 138 #if ARCH_ALPHA
 139     {"SIMPLE-ALPHA",    1, ff_simple_idct_axp,  ff_ref_idct, NO_PERM },
 140 #endif
 141
 142     { 0 }
 143 };
 144
 145 #define AANSCALE_BITS 12
 146
 147 uint8_t cropTbl[256 + 2 * MAX_NEG_CROP];
 148
 149 static int64_t gettime(void)
 150 {
 151     struct timeval tv;
 152     gettimeofday(&tv, NULL);
 153     return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
 154 }
 155
 156 #define NB_ITS 20000
 157 #define NB_ITS_SPEED 50000
 158
 159 static short idct_mmx_perm[64];
 160
 161 static short idct_simple_mmx_perm[64] = {
 162     0x00, 0x08, 0x04, 0x09, 0x01, 0x0C, 0x05, 0x0D,
 163     0x10, 0x18, 0x14, 0x19, 0x11, 0x1C, 0x15, 0x1D,
 164     0x20, 0x28, 0x24, 0x29, 0x21, 0x2C, 0x25, 0x2D,
 165     0x12, 0x1A, 0x16, 0x1B, 0x13, 0x1E, 0x17, 0x1F,
 166     0x02, 0x0A, 0x06, 0x0B, 0x03, 0x0E, 0x07, 0x0F,
 167     0x30, 0x38, 0x34, 0x39, 0x31, 0x3C, 0x35, 0x3D,
 168     0x22, 0x2A, 0x26, 0x2B, 0x23, 0x2E, 0x27, 0x2F,
 169     0x32, 0x3A, 0x36, 0x3B, 0x33, 0x3E, 0x37, 0x3F,
 170 };
 171
 172 static const uint8_t idct_sse2_row_perm[8] = { 0, 4, 1, 5, 2, 6, 3, 7 };
 173
 174 static void idct_mmx_init(void)
 175 {
 176     int i;
 177
 178     /* the mmx/mmxext idct uses a reordered input, so we patch scan tables */
 179     for (i = 0; i < 64; i++) {
 180         idct_mmx_perm[i] = (i & 0x38) | ((i & 6) >> 1) | ((i & 1) << 2);
 181     }
 182 }
 183
 184 DECLARE_ALIGNED(16, static DCTELEM, block)[64];
 185 DECLARE_ALIGNED(8,  static DCTELEM, block1)[64];
 186 DECLARE_ALIGNED(8,  static DCTELEM, block_org)[64];
 187
 188 static inline void mmx_emms(void)
 189 {
 190 #if HAVE_MMX
 191     if (cpu_flags & AV_CPU_FLAG_MMX)
 192         __asm__ volatile ("emms\n\t");
 193 #endif
 194 }
 195
 196 static void dct_error(const struct algo *dct, int test)
 197 {
 198     int it, i, scale;
 199     int err_inf, v;
 200     int64_t err2, ti, ti1, it1;
 201     int64_t sysErr[64], sysErrMax = 0;
 202     int maxout = 0;
 203     int blockSumErrMax = 0, blockSumErr;
 204     AVLFG prng;
 205
 206     av_lfg_init(&prng, 1);
 207
 208     err_inf = 0;
 209     err2 = 0;
 210     for (i = 0; i < 64; i++)
 211         sysErr[i] = 0;
 212     for (it = 0; it < NB_ITS; it++) {
 213         for (i = 0; i < 64; i++)
 214             block1[i] = 0;
 215         switch (test) {
 216         case 0:
 217             for (i = 0; i < 64; i++)
 218                 block1[i] = (av_lfg_get(&prng) % 512) - 256;
 219             if (dct->is_idct) {
 220                 ff_ref_fdct(block1);
 221                 for (i = 0; i < 64; i++)
 222                     block1[i] >>= 3;
 223             }
 224             break;
 225         case 1: {
 226                 int num = av_lfg_get(&prng) % 10 + 1;
 227                 for (i = 0; i < num; i++)
 228                     block1[av_lfg_get(&prng) % 64] =
 229                         av_lfg_get(&prng) % 512 - 256;
 230             }
 231             break;
 232         case 2:
 233             block1[0] = av_lfg_get(&prng) % 4096 - 2048;
 234             block1[63] = (block1[0] & 1) ^ 1;
 235             break;
 236         }
 237
 238         for (i = 0; i < 64; i++)
 239             block_org[i] = block1[i];
 240
 241         if (dct->format == MMX_PERM) {
 242             for (i = 0; i < 64; i++)
 243                 block[idct_mmx_perm[i]] = block1[i];
 244         } else if (dct->format == MMX_SIMPLE_PERM) {
 245             for (i = 0; i < 64; i++)
 246                 block[idct_simple_mmx_perm[i]] = block1[i];
 247         } else if (dct->format == SSE2_PERM) {
 248             for (i = 0; i < 64; i++)
 249                 block[(i & 0x38) | idct_sse2_row_perm[i & 7]] = block1[i];
 250         } else if (dct->format == PARTTRANS_PERM) {
 251             for (i = 0; i < 64; i++)
 252                 block[(i & 0x24) | ((i & 3) << 3) | ((i >> 3) & 3)] = block1[i];
 253         } else {
 254             for (i = 0; i < 64; i++)
 255                 block[i] = block1[i];
 256         }
 257
 258         dct->func(block);
 259         mmx_emms();
 260
 261         if (dct->format == SCALE_PERM) {
 262             for (i = 0; i < 64; i++) {
 263                 scale = 8 * (1 << (AANSCALE_BITS + 11)) / ff_aanscales[i];
 264                 block[i] = (block[i] * scale) >> AANSCALE_BITS;
 265             }
 266         }
 267
 268         dct->ref(block1);
 269
 270         blockSumErr = 0;
 271         for (i = 0; i < 64; i++) {
 272             v = abs(block[i] - block1[i]);
 273             if (v > err_inf)
 274                 err_inf = v;
 275             err2 += v * v;
 276             sysErr[i] += block[i] - block1[i];
 277             blockSumErr += v;
 278             if (abs(block[i]) > maxout)
 279                 maxout = abs(block[i]);
 280         }
 281         if (blockSumErrMax < blockSumErr)
 282             blockSumErrMax = blockSumErr;
 283     }
 284     for (i = 0; i < 64; i++)
 285         sysErrMax = FFMAX(sysErrMax, FFABS(sysErr[i]));
 286
 287     for (i = 0; i < 64; i++) {
 288         if (i % 8 == 0)
 289             printf("\n");
 290         printf("%7d ", (int) sysErr[i]);
 291     }
 292     printf("\n");
 293
 294     printf("%s %s: err_inf=%d err2=%0.8f syserr=%0.8f maxout=%d blockSumErr=%d\n",
 295            dct->is_idct ? "IDCT" : "DCT", dct->name, err_inf,
 296            (double) err2 / NB_ITS / 64.0, (double) sysErrMax / NB_ITS,
 297            maxout, blockSumErrMax);
 298
 299     /* speed test */
 300     for (i = 0; i < 64; i++)
 301         block1[i] = 0;
 302
 303     switch (test) {
 304     case 0:
 305         for (i = 0; i < 64; i++)
 306             block1[i] = av_lfg_get(&prng) % 512 - 256;
 307         if (dct->is_idct) {
 308             ff_ref_fdct(block1);
 309             for (i = 0; i < 64; i++)
 310                 block1[i] >>= 3;
 311         }
 312         break;
 313     case 1:
 314     case 2:
 315         block1[0] = av_lfg_get(&prng) % 512 - 256;
 316         block1[1] = av_lfg_get(&prng) % 512 - 256;
 317         block1[2] = av_lfg_get(&prng) % 512 - 256;
 318         block1[3] = av_lfg_get(&prng) % 512 - 256;
 319         break;
 320     }
 321
 322     if (dct->format == MMX_PERM) {
 323         for (i = 0; i < 64; i++)
 324             block[idct_mmx_perm[i]] = block1[i];
 325     } else if (dct->format == MMX_SIMPLE_PERM) {
 326         for (i = 0; i < 64; i++)
 327             block[idct_simple_mmx_perm[i]] = block1[i];
 328     } else {
 329         for (i = 0; i < 64; i++)
 330             block[i] = block1[i];
 331     }
 332
 333     ti = gettime();
 334     it1 = 0;
 335     do {
 336         for (it = 0; it < NB_ITS_SPEED; it++) {
 337             for (i = 0; i < 64; i++)
 338                 block[i] = block1[i];
 339             dct->func(block);
 340         }
 341         it1 += NB_ITS_SPEED;
 342         ti1 = gettime() - ti;
 343     } while (ti1 < 1000000);
 344     mmx_emms();
 345
 346     printf("%s %s: %0.1f kdct/s\n", dct->is_idct ? "IDCT" : "DCT", dct->name,
 347            (double) it1 * 1000.0 / (double) ti1);
 348 }
 349
 350 DECLARE_ALIGNED(8, static uint8_t, img_dest)[64];
 351 DECLARE_ALIGNED(8, static uint8_t, img_dest1)[64];
 352
 353 static void idct248_ref(uint8_t *dest, int linesize, int16_t *block)
 354 {
 355     static int init;
 356     static double c8[8][8];
 357     static double c4[4][4];
 358     double block1[64], block2[64], block3[64];
 359     double s, sum, v;
 360     int i, j, k;
 361
 362     if (!init) {
 363         init = 1;
 364
 365         for (i = 0; i < 8; i++) {
 366             sum = 0;
 367             for (j = 0; j < 8; j++) {
 368                 s = (i == 0) ? sqrt(1.0 / 8.0) : sqrt(1.0 / 4.0);
 369                 c8[i][j] = s * cos(M_PI * i * (j + 0.5) / 8.0);
 370                 sum += c8[i][j] * c8[i][j];
 371             }
 372         }
 373
 374         for (i = 0; i < 4; i++) {
 375             sum = 0;
 376             for (j = 0; j < 4; j++) {
 377                 s = (i == 0) ? sqrt(1.0 / 4.0) : sqrt(1.0 / 2.0);
 378                 c4[i][j] = s * cos(M_PI * i * (j + 0.5) / 4.0);
 379                 sum += c4[i][j] * c4[i][j];
 380             }
 381         }
 382     }
 383
 384     /* butterfly */
 385     s = 0.5 * sqrt(2.0);
 386     for (i = 0; i < 4; i++) {
 387         for (j = 0; j < 8; j++) {
 388             block1[8 * (2 * i) + j] =
 389                 (block[8 * (2 * i) + j] + block[8 * (2 * i + 1) + j]) * s;
 390             block1[8 * (2 * i + 1) + j] =
 391                 (block[8 * (2 * i) + j] - block[8 * (2 * i + 1) + j]) * s;
 392         }
 393     }
 394
 395     /* idct8 on lines */
 396     for (i = 0; i < 8; i++) {
 397         for (j = 0; j < 8; j++) {
 398             sum = 0;
 399             for (k = 0; k < 8; k++)
 400                 sum += c8[k][j] * block1[8 * i + k];
 401             block2[8 * i + j] = sum;
 402         }
 403     }
 404
 405     /* idct4 */
 406     for (i = 0; i < 8; i++) {
 407         for (j = 0; j < 4; j++) {
 408             /* top */
 409             sum = 0;
 410             for (k = 0; k < 4; k++)
 411                 sum += c4[k][j] * block2[8 * (2 * k) + i];
 412             block3[8 * (2 * j) + i] = sum;
 413
 414             /* bottom */
 415             sum = 0;
 416             for (k = 0; k < 4; k++)
 417                 sum += c4[k][j] * block2[8 * (2 * k + 1) + i];
 418             block3[8 * (2 * j + 1) + i] = sum;
 419         }
 420     }
 421
 422     /* clamp and store the result */
 423     for (i = 0; i < 8; i++) {
 424         for (j = 0; j < 8; j++) {
 425             v = block3[8 * i + j];
 426             if      (v < 0)   v = 0;
 427             else if (v > 255) v = 255;
 428             dest[i * linesize + j] = (int) rint(v);
 429         }
 430     }
 431 }
 432
 433 static void idct248_error(const char *name,
 434                           void (*idct248_put)(uint8_t *dest, int line_size,
 435                                               int16_t *block))
 436 {
 437     int it, i, it1, ti, ti1, err_max, v;
 438     AVLFG prng;
 439
 440     av_lfg_init(&prng, 1);
 441
 442     /* just one test to see if code is correct (precision is less
 443        important here) */
 444     err_max = 0;
 445     for (it = 0; it < NB_ITS; it++) {
 446         /* XXX: use forward transform to generate values */
 447         for (i = 0; i < 64; i++)
 448             block1[i] = av_lfg_get(&prng) % 256 - 128;
 449         block1[0] += 1024;
 450
 451         for (i = 0; i < 64; i++)
 452             block[i] = block1[i];
 453         idct248_ref(img_dest1, 8, block);
 454
 455         for (i = 0; i < 64; i++)
 456             block[i] = block1[i];
 457         idct248_put(img_dest, 8, block);
 458
 459         for (i = 0; i < 64; i++) {
 460             v = abs((int) img_dest[i] - (int) img_dest1[i]);
 461             if (v == 255)
 462                 printf("%d %d\n", img_dest[i], img_dest1[i]);
 463             if (v > err_max)
 464                 err_max = v;
 465         }
 466     }
 467     printf("%s %s: err_inf=%d\n", 1 ? "IDCT248" : "DCT248", name, err_max);
 468
 469     ti = gettime();
 470     it1 = 0;
 471     do {
 472         for (it = 0; it < NB_ITS_SPEED; it++) {
 473             for (i = 0; i < 64; i++)
 474                 block[i] = block1[i];
 475             idct248_put(img_dest, 8, block);
 476         }
 477         it1 += NB_ITS_SPEED;
 478         ti1 = gettime() - ti;
 479     } while (ti1 < 1000000);
 480     mmx_emms();
 481
 482     printf("%s %s: %0.1f kdct/s\n", 1 ? "IDCT248" : "DCT248", name,
 483            (double) it1 * 1000.0 / (double) ti1);
 484 }
 485
 486 static void help(void)
 487 {
 488     printf("dct-test [-i] [<test-number>]\n"
 489            "test-number 0 -> test with random matrixes\n"
 490            "            1 -> test with random sparse matrixes\n"
 491            "            2 -> do 3. test from mpeg4 std\n"
 492            "-i          test IDCT implementations\n"
 493            "-4          test IDCT248 implementations\n");
 494 }
 495
 496 int main(int argc, char **argv)
 497 {
 498     int test_idct = 0, test_248_dct = 0;
 499     int c, i;
 500     int test = 1;
 501
 502     cpu_flags = av_get_cpu_flags();
 503
 504     ff_ref_dct_init();
 505     idct_mmx_init();
 506
 507     for (i = 0; i < 256; i++)
 508         cropTbl[i + MAX_NEG_CROP] = i;
 509     for (i = 0; i < MAX_NEG_CROP; i++) {
 510         cropTbl[i] = 0;
 511         cropTbl[i + MAX_NEG_CROP + 256] = 255;
 512     }
 513
 514     for (;;) {
 515         c = getopt(argc, argv, "ih4");
 516         if (c == -1)
 517             break;
 518         switch (c) {
 519         case 'i':
 520             test_idct = 1;
 521             break;
 522         case '4':
 523             test_248_dct = 1;
 524             break;
 525         default:
 526         case 'h':
 527             help();
 528             return 0;
 529         }
 530     }
 531
 532     if (optind < argc)
 533         test = atoi(argv[optind]);
 534
 535     printf("ffmpeg DCT/IDCT test\n");
 536
 537     if (test_248_dct) {
 538         idct248_error("SIMPLE-C", ff_simple_idct248_put);
 539     } else {
 540         for (i = 0; algos[i].name; i++)
 541             if (algos[i].is_idct == test_idct &&
 542                 !(~cpu_flags & algos[i].mm_support)) {
 543                 dct_error(&algos[i], test);
 544             }
 545     }
 546     return 0;
 547 }