d21e5ed34075abec57d96f8f162a46d6939363bf
[ffmpeg.git] / libavcodec / simple_idct.c
1 /*
2  * Simple IDCT
3  *
4  * Copyright (c) 2001 Michael Niedermayer <michaelni@gmx.at>
5  *
6  * This file is part of Libav.
7  *
8  * Libav is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU Lesser General Public
10  * License as published by the Free Software Foundation; either
11  * version 2.1 of the License, or (at your option) any later version.
12  *
13  * Libav is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
16  * Lesser General Public License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with Libav; if not, write to the Free Software
20  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21  */
22
23 /**
24  * @file
25  * simpleidct in C.
26  */
27
28 /*
29   based upon some outcommented c code from mpeg2dec (idct_mmx.c
30   written by Aaron Holtzman <aholtzma@ess.engr.uvic.ca>)
31  */
32 #include "avcodec.h"
33 #include "dsputil.h"
34 #include "mathops.h"
35 #include "simple_idct.h"
36
37 #define W1  22725  //cos(i*M_PI/16)*sqrt(2)*(1<<14) + 0.5
38 #define W2  21407  //cos(i*M_PI/16)*sqrt(2)*(1<<14) + 0.5
39 #define W3  19266  //cos(i*M_PI/16)*sqrt(2)*(1<<14) + 0.5
40 #define W4  16383  //cos(i*M_PI/16)*sqrt(2)*(1<<14) + 0.5
41 #define W5  12873  //cos(i*M_PI/16)*sqrt(2)*(1<<14) + 0.5
42 #define W6  8867   //cos(i*M_PI/16)*sqrt(2)*(1<<14) + 0.5
43 #define W7  4520   //cos(i*M_PI/16)*sqrt(2)*(1<<14) + 0.5
44 #define ROW_SHIFT 11
45 #define COL_SHIFT 20 // 6
46
47 static inline void idctRowCondDC (DCTELEM * row)
48 {
49         int a0, a1, a2, a3, b0, b1, b2, b3;
50 #if HAVE_FAST_64BIT
51         uint64_t temp;
52 #else
53         uint32_t temp;
54 #endif
55
56 #if HAVE_FAST_64BIT
57 #if HAVE_BIGENDIAN
58 #define ROW0_MASK 0xffff000000000000LL
59 #else
60 #define ROW0_MASK 0xffffLL
61 #endif
62         if(sizeof(DCTELEM)==2){
63             if ( ((((uint64_t *)row)[0] & ~ROW0_MASK) |
64                   ((uint64_t *)row)[1]) == 0) {
65                 temp = (row[0] << 3) & 0xffff;
66                 temp += temp << 16;
67                 temp += temp << 32;
68                 ((uint64_t *)row)[0] = temp;
69                 ((uint64_t *)row)[1] = temp;
70                 return;
71             }
72         }else{
73             if (!(row[1]|row[2]|row[3]|row[4]|row[5]|row[6]|row[7])) {
74                 row[0]=row[1]=row[2]=row[3]=row[4]=row[5]=row[6]=row[7]= row[0] << 3;
75                 return;
76             }
77         }
78 #else
79         if(sizeof(DCTELEM)==2){
80             if (!(((uint32_t*)row)[1] |
81                   ((uint32_t*)row)[2] |
82                   ((uint32_t*)row)[3] |
83                   row[1])) {
84                 temp = (row[0] << 3) & 0xffff;
85                 temp += temp << 16;
86                 ((uint32_t*)row)[0]=((uint32_t*)row)[1] =
87                 ((uint32_t*)row)[2]=((uint32_t*)row)[3] = temp;
88                 return;
89             }
90         }else{
91             if (!(row[1]|row[2]|row[3]|row[4]|row[5]|row[6]|row[7])) {
92                 row[0]=row[1]=row[2]=row[3]=row[4]=row[5]=row[6]=row[7]= row[0] << 3;
93                 return;
94             }
95         }
96 #endif
97
98         a0 = (W4 * row[0]) + (1 << (ROW_SHIFT - 1));
99         a1 = a0;
100         a2 = a0;
101         a3 = a0;
102
103         /* no need to optimize : gcc does it */
104         a0 += W2 * row[2];
105         a1 += W6 * row[2];
106         a2 -= W6 * row[2];
107         a3 -= W2 * row[2];
108
109         b0 = MUL16(W1, row[1]);
110         MAC16(b0, W3, row[3]);
111         b1 = MUL16(W3, row[1]);
112         MAC16(b1, -W7, row[3]);
113         b2 = MUL16(W5, row[1]);
114         MAC16(b2, -W1, row[3]);
115         b3 = MUL16(W7, row[1]);
116         MAC16(b3, -W5, row[3]);
117
118 #if HAVE_FAST_64BIT
119         temp = ((uint64_t*)row)[1];
120 #else
121         temp = ((uint32_t*)row)[2] | ((uint32_t*)row)[3];
122 #endif
123         if (temp != 0) {
124             a0 += W4*row[4] + W6*row[6];
125             a1 += - W4*row[4] - W2*row[6];
126             a2 += - W4*row[4] + W2*row[6];
127             a3 += W4*row[4] - W6*row[6];
128
129             MAC16(b0, W5, row[5]);
130             MAC16(b0, W7, row[7]);
131
132             MAC16(b1, -W1, row[5]);
133             MAC16(b1, -W5, row[7]);
134
135             MAC16(b2, W7, row[5]);
136             MAC16(b2, W3, row[7]);
137
138             MAC16(b3, W3, row[5]);
139             MAC16(b3, -W1, row[7]);
140         }
141
142         row[0] = (a0 + b0) >> ROW_SHIFT;
143         row[7] = (a0 - b0) >> ROW_SHIFT;
144         row[1] = (a1 + b1) >> ROW_SHIFT;
145         row[6] = (a1 - b1) >> ROW_SHIFT;
146         row[2] = (a2 + b2) >> ROW_SHIFT;
147         row[5] = (a2 - b2) >> ROW_SHIFT;
148         row[3] = (a3 + b3) >> ROW_SHIFT;
149         row[4] = (a3 - b3) >> ROW_SHIFT;
150 }
151
152 static inline void idctSparseColPut (uint8_t *dest, int line_size,
153                                      DCTELEM * col)
154 {
155         int a0, a1, a2, a3, b0, b1, b2, b3;
156         uint8_t *cm = ff_cropTbl + MAX_NEG_CROP;
157
158         /* XXX: I did that only to give same values as previous code */
159         a0 = W4 * (col[8*0] + ((1<<(COL_SHIFT-1))/W4));
160         a1 = a0;
161         a2 = a0;
162         a3 = a0;
163
164         a0 +=  + W2*col[8*2];
165         a1 +=  + W6*col[8*2];
166         a2 +=  - W6*col[8*2];
167         a3 +=  - W2*col[8*2];
168
169         b0 = MUL16(W1, col[8*1]);
170         b1 = MUL16(W3, col[8*1]);
171         b2 = MUL16(W5, col[8*1]);
172         b3 = MUL16(W7, col[8*1]);
173
174         MAC16(b0, + W3, col[8*3]);
175         MAC16(b1, - W7, col[8*3]);
176         MAC16(b2, - W1, col[8*3]);
177         MAC16(b3, - W5, col[8*3]);
178
179         if(col[8*4]){
180             a0 += + W4*col[8*4];
181             a1 += - W4*col[8*4];
182             a2 += - W4*col[8*4];
183             a3 += + W4*col[8*4];
184         }
185
186         if (col[8*5]) {
187             MAC16(b0, + W5, col[8*5]);
188             MAC16(b1, - W1, col[8*5]);
189             MAC16(b2, + W7, col[8*5]);
190             MAC16(b3, + W3, col[8*5]);
191         }
192
193         if(col[8*6]){
194             a0 += + W6*col[8*6];
195             a1 += - W2*col[8*6];
196             a2 += + W2*col[8*6];
197             a3 += - W6*col[8*6];
198         }
199
200         if (col[8*7]) {
201             MAC16(b0, + W7, col[8*7]);
202             MAC16(b1, - W5, col[8*7]);
203             MAC16(b2, + W3, col[8*7]);
204             MAC16(b3, - W1, col[8*7]);
205         }
206
207         dest[0] = cm[(a0 + b0) >> COL_SHIFT];
208         dest += line_size;
209         dest[0] = cm[(a1 + b1) >> COL_SHIFT];
210         dest += line_size;
211         dest[0] = cm[(a2 + b2) >> COL_SHIFT];
212         dest += line_size;
213         dest[0] = cm[(a3 + b3) >> COL_SHIFT];
214         dest += line_size;
215         dest[0] = cm[(a3 - b3) >> COL_SHIFT];
216         dest += line_size;
217         dest[0] = cm[(a2 - b2) >> COL_SHIFT];
218         dest += line_size;
219         dest[0] = cm[(a1 - b1) >> COL_SHIFT];
220         dest += line_size;
221         dest[0] = cm[(a0 - b0) >> COL_SHIFT];
222 }
223
224 static inline void idctSparseColAdd (uint8_t *dest, int line_size,
225                                      DCTELEM * col)
226 {
227         int a0, a1, a2, a3, b0, b1, b2, b3;
228         uint8_t *cm = ff_cropTbl + MAX_NEG_CROP;
229
230         /* XXX: I did that only to give same values as previous code */
231         a0 = W4 * (col[8*0] + ((1<<(COL_SHIFT-1))/W4));
232         a1 = a0;
233         a2 = a0;
234         a3 = a0;
235
236         a0 +=  + W2*col[8*2];
237         a1 +=  + W6*col[8*2];
238         a2 +=  - W6*col[8*2];
239         a3 +=  - W2*col[8*2];
240
241         b0 = MUL16(W1, col[8*1]);
242         b1 = MUL16(W3, col[8*1]);
243         b2 = MUL16(W5, col[8*1]);
244         b3 = MUL16(W7, col[8*1]);
245
246         MAC16(b0, + W3, col[8*3]);
247         MAC16(b1, - W7, col[8*3]);
248         MAC16(b2, - W1, col[8*3]);
249         MAC16(b3, - W5, col[8*3]);
250
251         if(col[8*4]){
252             a0 += + W4*col[8*4];
253             a1 += - W4*col[8*4];
254             a2 += - W4*col[8*4];
255             a3 += + W4*col[8*4];
256         }
257
258         if (col[8*5]) {
259             MAC16(b0, + W5, col[8*5]);
260             MAC16(b1, - W1, col[8*5]);
261             MAC16(b2, + W7, col[8*5]);
262             MAC16(b3, + W3, col[8*5]);
263         }
264
265         if(col[8*6]){
266             a0 += + W6*col[8*6];
267             a1 += - W2*col[8*6];
268             a2 += + W2*col[8*6];
269             a3 += - W6*col[8*6];
270         }
271
272         if (col[8*7]) {
273             MAC16(b0, + W7, col[8*7]);
274             MAC16(b1, - W5, col[8*7]);
275             MAC16(b2, + W3, col[8*7]);
276             MAC16(b3, - W1, col[8*7]);
277         }
278
279         dest[0] = cm[dest[0] + ((a0 + b0) >> COL_SHIFT)];
280         dest += line_size;
281         dest[0] = cm[dest[0] + ((a1 + b1) >> COL_SHIFT)];
282         dest += line_size;
283         dest[0] = cm[dest[0] + ((a2 + b2) >> COL_SHIFT)];
284         dest += line_size;
285         dest[0] = cm[dest[0] + ((a3 + b3) >> COL_SHIFT)];
286         dest += line_size;
287         dest[0] = cm[dest[0] + ((a3 - b3) >> COL_SHIFT)];
288         dest += line_size;
289         dest[0] = cm[dest[0] + ((a2 - b2) >> COL_SHIFT)];
290         dest += line_size;
291         dest[0] = cm[dest[0] + ((a1 - b1) >> COL_SHIFT)];
292         dest += line_size;
293         dest[0] = cm[dest[0] + ((a0 - b0) >> COL_SHIFT)];
294 }
295
296 static inline void idctSparseCol (DCTELEM * col)
297 {
298         int a0, a1, a2, a3, b0, b1, b2, b3;
299
300         /* XXX: I did that only to give same values as previous code */
301         a0 = W4 * (col[8*0] + ((1<<(COL_SHIFT-1))/W4));
302         a1 = a0;
303         a2 = a0;
304         a3 = a0;
305
306         a0 +=  + W2*col[8*2];
307         a1 +=  + W6*col[8*2];
308         a2 +=  - W6*col[8*2];
309         a3 +=  - W2*col[8*2];
310
311         b0 = MUL16(W1, col[8*1]);
312         b1 = MUL16(W3, col[8*1]);
313         b2 = MUL16(W5, col[8*1]);
314         b3 = MUL16(W7, col[8*1]);
315
316         MAC16(b0, + W3, col[8*3]);
317         MAC16(b1, - W7, col[8*3]);
318         MAC16(b2, - W1, col[8*3]);
319         MAC16(b3, - W5, col[8*3]);
320
321         if(col[8*4]){
322             a0 += + W4*col[8*4];
323             a1 += - W4*col[8*4];
324             a2 += - W4*col[8*4];
325             a3 += + W4*col[8*4];
326         }
327
328         if (col[8*5]) {
329             MAC16(b0, + W5, col[8*5]);
330             MAC16(b1, - W1, col[8*5]);
331             MAC16(b2, + W7, col[8*5]);
332             MAC16(b3, + W3, col[8*5]);
333         }
334
335         if(col[8*6]){
336             a0 += + W6*col[8*6];
337             a1 += - W2*col[8*6];
338             a2 += + W2*col[8*6];
339             a3 += - W6*col[8*6];
340         }
341
342         if (col[8*7]) {
343             MAC16(b0, + W7, col[8*7]);
344             MAC16(b1, - W5, col[8*7]);
345             MAC16(b2, + W3, col[8*7]);
346             MAC16(b3, - W1, col[8*7]);
347         }
348
349         col[0 ] = ((a0 + b0) >> COL_SHIFT);
350         col[8 ] = ((a1 + b1) >> COL_SHIFT);
351         col[16] = ((a2 + b2) >> COL_SHIFT);
352         col[24] = ((a3 + b3) >> COL_SHIFT);
353         col[32] = ((a3 - b3) >> COL_SHIFT);
354         col[40] = ((a2 - b2) >> COL_SHIFT);
355         col[48] = ((a1 - b1) >> COL_SHIFT);
356         col[56] = ((a0 - b0) >> COL_SHIFT);
357 }
358
359 void ff_simple_idct_put(uint8_t *dest, int line_size, DCTELEM *block)
360 {
361     int i;
362     for(i=0; i<8; i++)
363         idctRowCondDC(block + i*8);
364
365     for(i=0; i<8; i++)
366         idctSparseColPut(dest + i, line_size, block + i);
367 }
368
369 void ff_simple_idct_add(uint8_t *dest, int line_size, DCTELEM *block)
370 {
371     int i;
372     for(i=0; i<8; i++)
373         idctRowCondDC(block + i*8);
374
375     for(i=0; i<8; i++)
376         idctSparseColAdd(dest + i, line_size, block + i);
377 }
378
379 void ff_simple_idct(DCTELEM *block)
380 {
381     int i;
382     for(i=0; i<8; i++)
383         idctRowCondDC(block + i*8);
384
385     for(i=0; i<8; i++)
386         idctSparseCol(block + i);
387 }
388
389 /* 2x4x8 idct */
390
391 #define CN_SHIFT 12
392 #define C_FIX(x) ((int)((x) * (1 << CN_SHIFT) + 0.5))
393 #define C1 C_FIX(0.6532814824)
394 #define C2 C_FIX(0.2705980501)
395
396 /* row idct is multiple by 16 * sqrt(2.0), col idct4 is normalized,
397    and the butterfly must be multiplied by 0.5 * sqrt(2.0) */
398 #define C_SHIFT (4+1+12)
399
400 static inline void idct4col_put(uint8_t *dest, int line_size, const DCTELEM *col)
401 {
402     int c0, c1, c2, c3, a0, a1, a2, a3;
403     const uint8_t *cm = ff_cropTbl + MAX_NEG_CROP;
404
405     a0 = col[8*0];
406     a1 = col[8*2];
407     a2 = col[8*4];
408     a3 = col[8*6];
409     c0 = ((a0 + a2) << (CN_SHIFT - 1)) + (1 << (C_SHIFT - 1));
410     c2 = ((a0 - a2) << (CN_SHIFT - 1)) + (1 << (C_SHIFT - 1));
411     c1 = a1 * C1 + a3 * C2;
412     c3 = a1 * C2 - a3 * C1;
413     dest[0] = cm[(c0 + c1) >> C_SHIFT];
414     dest += line_size;
415     dest[0] = cm[(c2 + c3) >> C_SHIFT];
416     dest += line_size;
417     dest[0] = cm[(c2 - c3) >> C_SHIFT];
418     dest += line_size;
419     dest[0] = cm[(c0 - c1) >> C_SHIFT];
420 }
421
422 #define BF(k) \
423 {\
424     int a0, a1;\
425     a0 = ptr[k];\
426     a1 = ptr[8 + k];\
427     ptr[k] = a0 + a1;\
428     ptr[8 + k] = a0 - a1;\
429 }
430
431 /* only used by DV codec. The input must be interlaced. 128 is added
432    to the pixels before clamping to avoid systematic error
433    (1024*sqrt(2)) offset would be needed otherwise. */
434 /* XXX: I think a 1.0/sqrt(2) normalization should be needed to
435    compensate the extra butterfly stage - I don't have the full DV
436    specification */
437 void ff_simple_idct248_put(uint8_t *dest, int line_size, DCTELEM *block)
438 {
439     int i;
440     DCTELEM *ptr;
441
442     /* butterfly */
443     ptr = block;
444     for(i=0;i<4;i++) {
445         BF(0);
446         BF(1);
447         BF(2);
448         BF(3);
449         BF(4);
450         BF(5);
451         BF(6);
452         BF(7);
453         ptr += 2 * 8;
454     }
455
456     /* IDCT8 on each line */
457     for(i=0; i<8; i++) {
458         idctRowCondDC(block + i*8);
459     }
460
461     /* IDCT4 and store */
462     for(i=0;i<8;i++) {
463         idct4col_put(dest + i, 2 * line_size, block + i);
464         idct4col_put(dest + line_size + i, 2 * line_size, block + 8 + i);
465     }
466 }
467
468 /* 8x4 & 4x8 WMV2 IDCT */
469 #undef CN_SHIFT
470 #undef C_SHIFT
471 #undef C_FIX
472 #undef C1
473 #undef C2
474 #define CN_SHIFT 12
475 #define C_FIX(x) ((int)((x) * 1.414213562 * (1 << CN_SHIFT) + 0.5))
476 #define C1 C_FIX(0.6532814824)
477 #define C2 C_FIX(0.2705980501)
478 #define C3 C_FIX(0.5)
479 #define C_SHIFT (4+1+12)
480 static inline void idct4col_add(uint8_t *dest, int line_size, const DCTELEM *col)
481 {
482     int c0, c1, c2, c3, a0, a1, a2, a3;
483     const uint8_t *cm = ff_cropTbl + MAX_NEG_CROP;
484
485     a0 = col[8*0];
486     a1 = col[8*1];
487     a2 = col[8*2];
488     a3 = col[8*3];
489     c0 = (a0 + a2)*C3 + (1 << (C_SHIFT - 1));
490     c2 = (a0 - a2)*C3 + (1 << (C_SHIFT - 1));
491     c1 = a1 * C1 + a3 * C2;
492     c3 = a1 * C2 - a3 * C1;
493     dest[0] = cm[dest[0] + ((c0 + c1) >> C_SHIFT)];
494     dest += line_size;
495     dest[0] = cm[dest[0] + ((c2 + c3) >> C_SHIFT)];
496     dest += line_size;
497     dest[0] = cm[dest[0] + ((c2 - c3) >> C_SHIFT)];
498     dest += line_size;
499     dest[0] = cm[dest[0] + ((c0 - c1) >> C_SHIFT)];
500 }
501
502 #define RN_SHIFT 15
503 #define R_FIX(x) ((int)((x) * 1.414213562 * (1 << RN_SHIFT) + 0.5))
504 #define R1 R_FIX(0.6532814824)
505 #define R2 R_FIX(0.2705980501)
506 #define R3 R_FIX(0.5)
507 #define R_SHIFT 11
508 static inline void idct4row(DCTELEM *row)
509 {
510     int c0, c1, c2, c3, a0, a1, a2, a3;
511     //const uint8_t *cm = ff_cropTbl + MAX_NEG_CROP;
512
513     a0 = row[0];
514     a1 = row[1];
515     a2 = row[2];
516     a3 = row[3];
517     c0 = (a0 + a2)*R3 + (1 << (R_SHIFT - 1));
518     c2 = (a0 - a2)*R3 + (1 << (R_SHIFT - 1));
519     c1 = a1 * R1 + a3 * R2;
520     c3 = a1 * R2 - a3 * R1;
521     row[0]= (c0 + c1) >> R_SHIFT;
522     row[1]= (c2 + c3) >> R_SHIFT;
523     row[2]= (c2 - c3) >> R_SHIFT;
524     row[3]= (c0 - c1) >> R_SHIFT;
525 }
526
527 void ff_simple_idct84_add(uint8_t *dest, int line_size, DCTELEM *block)
528 {
529     int i;
530
531     /* IDCT8 on each line */
532     for(i=0; i<4; i++) {
533         idctRowCondDC(block + i*8);
534     }
535
536     /* IDCT4 and store */
537     for(i=0;i<8;i++) {
538         idct4col_add(dest + i, line_size, block + i);
539     }
540 }
541
542 void ff_simple_idct48_add(uint8_t *dest, int line_size, DCTELEM *block)
543 {
544     int i;
545
546     /* IDCT4 on each line */
547     for(i=0; i<8; i++) {
548         idct4row(block + i*8);
549     }
550
551     /* IDCT8 and store */
552     for(i=0; i<4; i++){
553         idctSparseColAdd(dest + i, line_size, block + i);
554     }
555 }
556
557 void ff_simple_idct44_add(uint8_t *dest, int line_size, DCTELEM *block)
558 {
559     int i;
560
561     /* IDCT4 on each line */
562     for(i=0; i<4; i++) {
563         idct4row(block + i*8);
564     }
565
566     /* IDCT4 and store */
567     for(i=0; i<4; i++){
568         idct4col_add(dest + i, line_size, block + i);
569     }
570 }