Merge remote-tracking branch 'qatar/master'
[ffmpeg.git] / 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 FFmpeg.
6  *
7  * FFmpeg 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  * FFmpeg 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 FFmpeg; 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     void (*func)(DCTELEM *block);
73     enum formattag { NO_PERM, MMX_PERM, MMX_SIMPLE_PERM, SCALE_PERM,
74                      SSE2_PERM, PARTTRANS_PERM } format;
75     int mm_support;
76     int nonspec;
77 };
78
79 #ifndef FAAN_POSTSCALE
80 #define FAAN_SCALE SCALE_PERM
81 #else
82 #define FAAN_SCALE NO_PERM
83 #endif
84
85 static int cpu_flags;
86
87 static const struct algo fdct_tab[] = {
88     { "REF-DBL",        ff_ref_fdct,           NO_PERM    },
89     { "FAAN",           ff_faandct,            FAAN_SCALE },
90     { "IJG-AAN-INT",    fdct_ifast,            SCALE_PERM },
91     { "IJG-LLM-INT",    ff_jpeg_fdct_islow,    NO_PERM    },
92
93 #if HAVE_MMX
94     { "MMX",            ff_fdct_mmx,           NO_PERM,   AV_CPU_FLAG_MMX     },
95     { "MMX2",           ff_fdct_mmx2,          NO_PERM,   AV_CPU_FLAG_MMX2    },
96     { "SSE2",           ff_fdct_sse2,          NO_PERM,   AV_CPU_FLAG_SSE2    },
97 #endif
98
99 #if HAVE_ALTIVEC
100     { "altivecfdct",    fdct_altivec,          NO_PERM,   AV_CPU_FLAG_ALTIVEC },
101 #endif
102
103 #if ARCH_BFIN
104     { "BFINfdct",       ff_bfin_fdct,          NO_PERM  },
105 #endif
106
107     { 0 }
108 };
109
110 static const struct algo idct_tab[] = {
111     { "FAANI",          ff_faanidct,           NO_PERM  },
112     { "REF-DBL",        ff_ref_idct,           NO_PERM  },
113     { "INT",            j_rev_dct,             MMX_PERM },
114     { "SIMPLE-C",       ff_simple_idct_8,      NO_PERM  },
115
116 #if HAVE_MMX
117 #if CONFIG_GPL
118     { "LIBMPEG2-MMX",   ff_mmx_idct,           MMX_PERM,  AV_CPU_FLAG_MMX,  1 },
119     { "LIBMPEG2-MMX2",  ff_mmxext_idct,        MMX_PERM,  AV_CPU_FLAG_MMX2, 1 },
120 #endif
121     { "SIMPLE-MMX",     ff_simple_idct_mmx,  MMX_SIMPLE_PERM, AV_CPU_FLAG_MMX },
122     { "XVID-MMX",       ff_idct_xvid_mmx,      NO_PERM,   AV_CPU_FLAG_MMX,  1 },
123     { "XVID-MMX2",      ff_idct_xvid_mmx2,     NO_PERM,   AV_CPU_FLAG_MMX2, 1 },
124     { "XVID-SSE2",      ff_idct_xvid_sse2,     SSE2_PERM, AV_CPU_FLAG_SSE2, 1 },
125 #endif
126
127 #if ARCH_BFIN
128     { "BFINidct",       ff_bfin_idct,          NO_PERM  },
129 #endif
130
131 #if ARCH_ARM
132     { "SIMPLE-ARM",     ff_simple_idct_arm,    NO_PERM  },
133     { "INT-ARM",        ff_j_rev_dct_arm,      MMX_PERM },
134 #endif
135 #if HAVE_ARMV5TE
136     { "SIMPLE-ARMV5TE", ff_simple_idct_armv5te,NO_PERM  },
137 #endif
138 #if HAVE_ARMV6
139     { "SIMPLE-ARMV6",   ff_simple_idct_armv6,  MMX_PERM },
140 #endif
141 #if HAVE_NEON
142     { "SIMPLE-NEON",    ff_simple_idct_neon,   PARTTRANS_PERM },
143 #endif
144
145 #if ARCH_ALPHA
146     { "SIMPLE-ALPHA",   ff_simple_idct_axp,    NO_PERM },
147 #endif
148
149     { 0 }
150 };
151
152 #define AANSCALE_BITS 12
153
154 uint8_t cropTbl[256 + 2 * MAX_NEG_CROP];
155
156 static int64_t gettime(void)
157 {
158     struct timeval tv;
159     gettimeofday(&tv, NULL);
160     return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
161 }
162
163 #define NB_ITS 20000
164 #define NB_ITS_SPEED 50000
165
166 static short idct_mmx_perm[64];
167
168 static short idct_simple_mmx_perm[64] = {
169     0x00, 0x08, 0x04, 0x09, 0x01, 0x0C, 0x05, 0x0D,
170     0x10, 0x18, 0x14, 0x19, 0x11, 0x1C, 0x15, 0x1D,
171     0x20, 0x28, 0x24, 0x29, 0x21, 0x2C, 0x25, 0x2D,
172     0x12, 0x1A, 0x16, 0x1B, 0x13, 0x1E, 0x17, 0x1F,
173     0x02, 0x0A, 0x06, 0x0B, 0x03, 0x0E, 0x07, 0x0F,
174     0x30, 0x38, 0x34, 0x39, 0x31, 0x3C, 0x35, 0x3D,
175     0x22, 0x2A, 0x26, 0x2B, 0x23, 0x2E, 0x27, 0x2F,
176     0x32, 0x3A, 0x36, 0x3B, 0x33, 0x3E, 0x37, 0x3F,
177 };
178
179 static const uint8_t idct_sse2_row_perm[8] = { 0, 4, 1, 5, 2, 6, 3, 7 };
180
181 static void idct_mmx_init(void)
182 {
183     int i;
184
185     /* the mmx/mmxext idct uses a reordered input, so we patch scan tables */
186     for (i = 0; i < 64; i++) {
187         idct_mmx_perm[i] = (i & 0x38) | ((i & 6) >> 1) | ((i & 1) << 2);
188     }
189 }
190
191 DECLARE_ALIGNED(16, static DCTELEM, block)[64];
192 DECLARE_ALIGNED(8,  static DCTELEM, block1)[64];
193
194 static inline void mmx_emms(void)
195 {
196 #if HAVE_MMX
197     if (cpu_flags & AV_CPU_FLAG_MMX)
198         __asm__ volatile ("emms\n\t");
199 #endif
200 }
201
202
203 static int dct_error(const struct algo *dct, int test, int is_idct, int speed, const int bits)
204 {
205     void (*ref)(DCTELEM *block) = is_idct ? ff_ref_idct : ff_ref_fdct;
206     int it, i, scale;
207     int err_inf, v;
208     int64_t err2, ti, ti1, it1, err_sum = 0;
209     int64_t sysErr[64], sysErrMax = 0;
210     int maxout = 0;
211     int blockSumErrMax = 0, blockSumErr;
212     AVLFG prng;
213     const int vals=1<<bits;
214     double omse, ome;
215     int spec_err;
216
217     av_lfg_init(&prng, 1);
218
219     err_inf = 0;
220     err2 = 0;
221     for (i = 0; i < 64; i++)
222         sysErr[i] = 0;
223     for (it = 0; it < NB_ITS; it++) {
224         for (i = 0; i < 64; i++)
225             block1[i] = 0;
226         switch (test) {
227         case 0:
228             for (i = 0; i < 64; i++)
229                 block1[i] = (av_lfg_get(&prng) % (2*vals)) -vals;
230             if (is_idct) {
231                 ff_ref_fdct(block1);
232                 for (i = 0; i < 64; i++)
233                     block1[i] >>= 3;
234             }
235             break;
236         case 1: {
237                 int num = av_lfg_get(&prng) % 10 + 1;
238                 for (i = 0; i < num; i++)
239                     block1[av_lfg_get(&prng) % 64] = av_lfg_get(&prng) % (2*vals) -vals;
240             }
241             break;
242         case 2:
243             block1[0] = av_lfg_get(&prng) % (16*vals) - (8*vals);
244             block1[63] = (block1[0] & 1) ^ 1;
245             break;
246         }
247
248         if (dct->format == MMX_PERM) {
249             for (i = 0; i < 64; i++)
250                 block[idct_mmx_perm[i]] = block1[i];
251         } else if (dct->format == MMX_SIMPLE_PERM) {
252             for (i = 0; i < 64; i++)
253                 block[idct_simple_mmx_perm[i]] = block1[i];
254         } else if (dct->format == SSE2_PERM) {
255             for (i = 0; i < 64; i++)
256                 block[(i & 0x38) | idct_sse2_row_perm[i & 7]] = block1[i];
257         } else if (dct->format == PARTTRANS_PERM) {
258             for (i = 0; i < 64; i++)
259                 block[(i & 0x24) | ((i & 3) << 3) | ((i >> 3) & 3)] = block1[i];
260         } else {
261             for (i = 0; i < 64; i++)
262                 block[i] = block1[i];
263         }
264
265         dct->func(block);
266         mmx_emms();
267
268         if (dct->format == SCALE_PERM) {
269             for (i = 0; i < 64; i++) {
270                 scale = 8 * (1 << (AANSCALE_BITS + 11)) / ff_aanscales[i];
271                 block[i] = (block[i] * scale) >> AANSCALE_BITS;
272             }
273         }
274
275         ref(block1);
276
277         blockSumErr = 0;
278         for (i = 0; i < 64; i++) {
279             int err = block[i] - block1[i];
280             err_sum += err;
281             v = abs(err);
282             if (v > err_inf)
283                 err_inf = v;
284             err2 += v * v;
285             sysErr[i] += block[i] - block1[i];
286             blockSumErr += v;
287             if (abs(block[i]) > maxout)
288                 maxout = abs(block[i]);
289         }
290         if (blockSumErrMax < blockSumErr)
291             blockSumErrMax = blockSumErr;
292     }
293     for (i = 0; i < 64; i++)
294         sysErrMax = FFMAX(sysErrMax, FFABS(sysErr[i]));
295
296     for (i = 0; i < 64; i++) {
297         if (i % 8 == 0)
298             printf("\n");
299         printf("%7d ", (int) sysErr[i]);
300     }
301     printf("\n");
302
303     omse = (double) err2 / NB_ITS / 64;
304     ome  = (double) err_sum / NB_ITS / 64;
305
306     spec_err = is_idct && (err_inf > 1 || omse > 0.02 || fabs(ome) > 0.0015);
307
308     printf("%s %s: ppe=%d omse=%0.8f ome=%0.8f syserr=%0.8f maxout=%d blockSumErr=%d\n",
309            is_idct ? "IDCT" : "DCT", dct->name, err_inf,
310            omse, ome, (double) sysErrMax / NB_ITS,
311            maxout, blockSumErrMax);
312
313     if (spec_err && !dct->nonspec)
314         return 1;
315
316     if (!speed)
317         return 0;
318
319     /* speed test */
320     for (i = 0; i < 64; i++)
321         block1[i] = 0;
322
323     switch (test) {
324     case 0:
325         for (i = 0; i < 64; i++)
326             block1[i] = av_lfg_get(&prng) % (2*vals) -vals;
327         if (is_idct) {
328             ff_ref_fdct(block1);
329             for (i = 0; i < 64; i++)
330                 block1[i] >>= 3;
331         }
332         break;
333     case 1:
334     case 2:
335         block1[0] = av_lfg_get(&prng) % (2*vals) -vals;
336         block1[1] = av_lfg_get(&prng) % (2*vals) -vals;
337         block1[2] = av_lfg_get(&prng) % (2*vals) -vals;
338         block1[3] = av_lfg_get(&prng) % (2*vals) -vals;
339         break;
340     }
341
342     if (dct->format == MMX_PERM) {
343         for (i = 0; i < 64; i++)
344             block[idct_mmx_perm[i]] = block1[i];
345     } else if (dct->format == MMX_SIMPLE_PERM) {
346         for (i = 0; i < 64; i++)
347             block[idct_simple_mmx_perm[i]] = block1[i];
348     } else {
349         for (i = 0; i < 64; i++)
350             block[i] = block1[i];
351     }
352
353     ti = gettime();
354     it1 = 0;
355     do {
356         for (it = 0; it < NB_ITS_SPEED; it++) {
357             for (i = 0; i < 64; i++)
358                 block[i] = block1[i];
359             dct->func(block);
360         }
361         it1 += NB_ITS_SPEED;
362         ti1 = gettime() - ti;
363     } while (ti1 < 1000000);
364     mmx_emms();
365
366     printf("%s %s: %0.1f kdct/s\n", is_idct ? "IDCT" : "DCT", dct->name,
367            (double) it1 * 1000.0 / (double) ti1);
368
369     return 0;
370 }
371
372 DECLARE_ALIGNED(8, static uint8_t, img_dest)[64];
373 DECLARE_ALIGNED(8, static uint8_t, img_dest1)[64];
374
375 static void idct248_ref(uint8_t *dest, int linesize, int16_t *block)
376 {
377     static int init;
378     static double c8[8][8];
379     static double c4[4][4];
380     double block1[64], block2[64], block3[64];
381     double s, sum, v;
382     int i, j, k;
383
384     if (!init) {
385         init = 1;
386
387         for (i = 0; i < 8; i++) {
388             sum = 0;
389             for (j = 0; j < 8; j++) {
390                 s = (i == 0) ? sqrt(1.0 / 8.0) : sqrt(1.0 / 4.0);
391                 c8[i][j] = s * cos(M_PI * i * (j + 0.5) / 8.0);
392                 sum += c8[i][j] * c8[i][j];
393             }
394         }
395
396         for (i = 0; i < 4; i++) {
397             sum = 0;
398             for (j = 0; j < 4; j++) {
399                 s = (i == 0) ? sqrt(1.0 / 4.0) : sqrt(1.0 / 2.0);
400                 c4[i][j] = s * cos(M_PI * i * (j + 0.5) / 4.0);
401                 sum += c4[i][j] * c4[i][j];
402             }
403         }
404     }
405
406     /* butterfly */
407     s = 0.5 * sqrt(2.0);
408     for (i = 0; i < 4; i++) {
409         for (j = 0; j < 8; j++) {
410             block1[8 * (2 * i) + j] =
411                 (block[8 * (2 * i) + j] + block[8 * (2 * i + 1) + j]) * s;
412             block1[8 * (2 * i + 1) + j] =
413                 (block[8 * (2 * i) + j] - block[8 * (2 * i + 1) + j]) * s;
414         }
415     }
416
417     /* idct8 on lines */
418     for (i = 0; i < 8; i++) {
419         for (j = 0; j < 8; j++) {
420             sum = 0;
421             for (k = 0; k < 8; k++)
422                 sum += c8[k][j] * block1[8 * i + k];
423             block2[8 * i + j] = sum;
424         }
425     }
426
427     /* idct4 */
428     for (i = 0; i < 8; i++) {
429         for (j = 0; j < 4; j++) {
430             /* top */
431             sum = 0;
432             for (k = 0; k < 4; k++)
433                 sum += c4[k][j] * block2[8 * (2 * k) + i];
434             block3[8 * (2 * j) + i] = sum;
435
436             /* bottom */
437             sum = 0;
438             for (k = 0; k < 4; k++)
439                 sum += c4[k][j] * block2[8 * (2 * k + 1) + i];
440             block3[8 * (2 * j + 1) + i] = sum;
441         }
442     }
443
444     /* clamp and store the result */
445     for (i = 0; i < 8; i++) {
446         for (j = 0; j < 8; j++) {
447             v = block3[8 * i + j];
448             if      (v < 0)   v = 0;
449             else if (v > 255) v = 255;
450             dest[i * linesize + j] = (int) rint(v);
451         }
452     }
453 }
454
455 static void idct248_error(const char *name,
456                           void (*idct248_put)(uint8_t *dest, int line_size,
457                                               int16_t *block),
458                           int speed)
459 {
460     int it, i, it1, ti, ti1, err_max, v;
461     AVLFG prng;
462
463     av_lfg_init(&prng, 1);
464
465     /* just one test to see if code is correct (precision is less
466        important here) */
467     err_max = 0;
468     for (it = 0; it < NB_ITS; it++) {
469         /* XXX: use forward transform to generate values */
470         for (i = 0; i < 64; i++)
471             block1[i] = av_lfg_get(&prng) % 256 - 128;
472         block1[0] += 1024;
473
474         for (i = 0; i < 64; i++)
475             block[i] = block1[i];
476         idct248_ref(img_dest1, 8, block);
477
478         for (i = 0; i < 64; i++)
479             block[i] = block1[i];
480         idct248_put(img_dest, 8, block);
481
482         for (i = 0; i < 64; i++) {
483             v = abs((int) img_dest[i] - (int) img_dest1[i]);
484             if (v == 255)
485                 printf("%d %d\n", img_dest[i], img_dest1[i]);
486             if (v > err_max)
487                 err_max = v;
488         }
489 #if 0
490         printf("ref=\n");
491         for(i=0;i<8;i++) {
492             int j;
493             for(j=0;j<8;j++) {
494                 printf(" %3d", img_dest1[i*8+j]);
495             }
496             printf("\n");
497         }
498
499         printf("out=\n");
500         for(i=0;i<8;i++) {
501             int j;
502             for(j=0;j<8;j++) {
503                 printf(" %3d", img_dest[i*8+j]);
504             }
505             printf("\n");
506         }
507 #endif
508     }
509     printf("%s %s: err_inf=%d\n", 1 ? "IDCT248" : "DCT248", name, err_max);
510
511     if (!speed)
512         return;
513
514     ti = gettime();
515     it1 = 0;
516     do {
517         for (it = 0; it < NB_ITS_SPEED; it++) {
518             for (i = 0; i < 64; i++)
519                 block[i] = block1[i];
520             idct248_put(img_dest, 8, block);
521         }
522         it1 += NB_ITS_SPEED;
523         ti1 = gettime() - ti;
524     } while (ti1 < 1000000);
525     mmx_emms();
526
527     printf("%s %s: %0.1f kdct/s\n", 1 ? "IDCT248" : "DCT248", name,
528            (double) it1 * 1000.0 / (double) ti1);
529 }
530
531 static void help(void)
532 {
533     printf("dct-test [-i] [<test-number>]\n"
534            "test-number 0 -> test with random matrixes\n"
535            "            1 -> test with random sparse matrixes\n"
536            "            2 -> do 3. test from mpeg4 std\n"
537            "-i          test IDCT implementations\n"
538            "-4          test IDCT248 implementations\n"
539            "-t          speed test\n");
540 }
541
542 int main(int argc, char **argv)
543 {
544     int test_idct = 0, test_248_dct = 0;
545     int c, i;
546     int test = 1;
547     int speed = 0;
548     int err = 0;
549     int bits=8;
550
551     cpu_flags = av_get_cpu_flags();
552
553     ff_ref_dct_init();
554     idct_mmx_init();
555
556     for (i = 0; i < 256; i++)
557         cropTbl[i + MAX_NEG_CROP] = i;
558     for (i = 0; i < MAX_NEG_CROP; i++) {
559         cropTbl[i] = 0;
560         cropTbl[i + MAX_NEG_CROP + 256] = 255;
561     }
562
563     for (;;) {
564         c = getopt(argc, argv, "ih4t");
565         if (c == -1)
566             break;
567         switch (c) {
568         case 'i':
569             test_idct = 1;
570             break;
571         case '4':
572             test_248_dct = 1;
573             break;
574         case 't':
575             speed = 1;
576             break;
577         default:
578         case 'h':
579             help();
580             return 0;
581         }
582     }
583
584     if (optind < argc)
585         test = atoi(argv[optind]);
586     if(optind+1 < argc) bits= atoi(argv[optind+1]);
587
588     printf("ffmpeg DCT/IDCT test\n");
589
590     if (test_248_dct) {
591         idct248_error("SIMPLE-C", ff_simple_idct248_put, speed);
592     } else {
593         const struct algo *algos = test_idct ? idct_tab : fdct_tab;
594         for (i = 0; algos[i].name; i++)
595             if (!(~cpu_flags & algos[i].mm_support)) {
596                 err |= dct_error(&algos[i], test, test_idct, speed, bits);
597             }
598     }
599
600     return err;
601 }