1 /*
2 * idct.c
3 * Copyright (C) 2000-2002 Michel Lespinasse <walken@zoy.org>
4 * Copyright (C) 1999-2000 Aaron Holtzman <aholtzma@ess.engr.uvic.ca>
5 *
6 * Portions of this code are from the MPEG software simulation group
7 * idct implementation. This code will be replaced with a new
8 * implementation soon.
9 *
10 * This file is part of mpeg2dec, a free MPEG-2 video stream decoder.
11 * See http://libmpeg2.sourceforge.net/ for updates.
12 *
13 * mpeg2dec is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License as published by
15 * the Free Software Foundation; either version 2 of the License, or
16 * (at your option) any later version.
17 *
18 * mpeg2dec is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU General Public License for more details.
22 *
23 * You should have received a copy of the GNU General Public License
24 * along with this program; if not, write to the Free Software
25 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
26 */
27
28 /**********************************************************/
29 /* inverse two dimensional DCT, Chen-Wang algorithm */
30 /* (cf. IEEE ASSP-32, pp. 803-816, Aug. 1984) */
31 /* 32-bit integer arithmetic (8 bit coefficients) */
32 /* 11 mults, 29 adds per DCT */
33 /* sE, 18.8.91 */
34 /**********************************************************/
35 /* coefficients extended to 12 bit for IEEE1180-1990 */
36 /* compliance sE, 2.1.94 */
37 /**********************************************************/
38
39 /* this code assumes >> to be a two's-complement arithmetic */
40 /* right shift: (-2)>>1 == -1 , (-3)>>1 == -2 */
41
42 #include "config.h"
43
44 #include <stdio.h>
45 #include <stdlib.h>
46 #include <inttypes.h>
47
48 #include "mpeg2_internal.h"
49 #include <xine/xineutils.h>
50
51 #define W1 2841 /* 2048*sqrt (2)*cos (1*pi/16) */
52 #define W2 2676 /* 2048*sqrt (2)*cos (2*pi/16) */
53 #define W3 2408 /* 2048*sqrt (2)*cos (3*pi/16) */
54 #define W5 1609 /* 2048*sqrt (2)*cos (5*pi/16) */
55 #define W6 1108 /* 2048*sqrt (2)*cos (6*pi/16) */
56 #define W7 565 /* 2048*sqrt (2)*cos (7*pi/16) */
57
58 /* idct main entry points */
59 void (* mpeg2_idct_copy) (int16_t * block, uint8_t * dest, int stride);
60 void (* mpeg2_idct_add) (int16_t * block, uint8_t * dest, int stride);
61 void (* mpeg2_idct) (int16_t * block);
62 void (* mpeg2_zero_block) (int16_t * block);
63
64 static uint8_t clip_lut[1024];
65 #define CLIP(i) ((clip_lut+384)[ (i)])
66
67 /* row (horizontal) IDCT
68 *
69 * 7 pi 1
70 * dst[k] = sum c[l] * src[l] * cos ( -- * ( k + - ) * l )
71 * l=0 8 2
72 *
73 * where: c[0] = 128
74 * c[1..7] = 128*sqrt (2)
75 */
76
idct_row(int16_t * block)77 static inline void idct_row (int16_t * block)
78 {
79 int x0, x1, x2, x3, x4, x5, x6, x7, x8;
80
81 x1 = block[4] << 11;
82 x2 = block[6];
83 x3 = block[2];
84 x4 = block[1];
85 x5 = block[7];
86 x6 = block[5];
87 x7 = block[3];
88
89 /* shortcut */
90 if (! (x1 | x2 | x3 | x4 | x5 | x6 | x7 )) {
91 block[0] = block[1] = block[2] = block[3] = block[4] =
92 block[5] = block[6] = block[7] = block[0]<<3;
93 return;
94 }
95
96 x0 = (block[0] << 11) + 128; /* for proper rounding in the fourth stage */
97
98 /* first stage */
99 x8 = W7 * (x4 + x5);
100 x4 = x8 + (W1 - W7) * x4;
101 x5 = x8 - (W1 + W7) * x5;
102 x8 = W3 * (x6 + x7);
103 x6 = x8 - (W3 - W5) * x6;
104 x7 = x8 - (W3 + W5) * x7;
105
106 /* second stage */
107 x8 = x0 + x1;
108 x0 -= x1;
109 x1 = W6 * (x3 + x2);
110 x2 = x1 - (W2 + W6) * x2;
111 x3 = x1 + (W2 - W6) * x3;
112 x1 = x4 + x6;
113 x4 -= x6;
114 x6 = x5 + x7;
115 x5 -= x7;
116
117 /* third stage */
118 x7 = x8 + x3;
119 x8 -= x3;
120 x3 = x0 + x2;
121 x0 -= x2;
122 x2 = (181 * (x4 + x5) + 128) >> 8;
123 x4 = (181 * (x4 - x5) + 128) >> 8;
124
125 /* fourth stage */
126 block[0] = (x7 + x1) >> 8;
127 block[1] = (x3 + x2) >> 8;
128 block[2] = (x0 + x4) >> 8;
129 block[3] = (x8 + x6) >> 8;
130 block[4] = (x8 - x6) >> 8;
131 block[5] = (x0 - x4) >> 8;
132 block[6] = (x3 - x2) >> 8;
133 block[7] = (x7 - x1) >> 8;
134 }
135
136 /* column (vertical) IDCT
137 *
138 * 7 pi 1
139 * dst[8*k] = sum c[l] * src[8*l] * cos ( -- * ( k + - ) * l )
140 * l=0 8 2
141 *
142 * where: c[0] = 1/1024
143 * c[1..7] = (1/1024)*sqrt (2)
144 */
145
idct_col(int16_t * block)146 static inline void idct_col (int16_t *block)
147 {
148 int x0, x1, x2, x3, x4, x5, x6, x7, x8;
149
150 /* shortcut */
151 x1 = block [8*4] << 8;
152 x2 = block [8*6];
153 x3 = block [8*2];
154 x4 = block [8*1];
155 x5 = block [8*7];
156 x6 = block [8*5];
157 x7 = block [8*3];
158
159 #if 0
160 if (! (x1 | x2 | x3 | x4 | x5 | x6 | x7 )) {
161 block[8*0] = block[8*1] = block[8*2] = block[8*3] = block[8*4] =
162 block[8*5] = block[8*6] = block[8*7] = (block[8*0] + 32) >> 6;
163 return;
164 }
165 #endif
166
167 x0 = (block[8*0] << 8) + 8192;
168
169 /* first stage */
170 x8 = W7 * (x4 + x5) + 4;
171 x4 = (x8 + (W1 - W7) * x4) >> 3;
172 x5 = (x8 - (W1 + W7) * x5) >> 3;
173 x8 = W3 * (x6 + x7) + 4;
174 x6 = (x8 - (W3 - W5) * x6) >> 3;
175 x7 = (x8 - (W3 + W5) * x7) >> 3;
176
177 /* second stage */
178 x8 = x0 + x1;
179 x0 -= x1;
180 x1 = W6 * (x3 + x2) + 4;
181 x2 = (x1 - (W2 + W6) * x2) >> 3;
182 x3 = (x1 + (W2 - W6) * x3) >> 3;
183 x1 = x4 + x6;
184 x4 -= x6;
185 x6 = x5 + x7;
186 x5 -= x7;
187
188 /* third stage */
189 x7 = x8 + x3;
190 x8 -= x3;
191 x3 = x0 + x2;
192 x0 -= x2;
193 x2 = (181 * (x4 + x5) + 128) >> 8;
194 x4 = (181 * (x4 - x5) + 128) >> 8;
195
196 /* fourth stage */
197 block[8*0] = (x7 + x1) >> 14;
198 block[8*1] = (x3 + x2) >> 14;
199 block[8*2] = (x0 + x4) >> 14;
200 block[8*3] = (x8 + x6) >> 14;
201 block[8*4] = (x8 - x6) >> 14;
202 block[8*5] = (x0 - x4) >> 14;
203 block[8*6] = (x3 - x2) >> 14;
204 block[8*7] = (x7 - x1) >> 14;
205 }
206
mpeg2_idct_copy_c(int16_t * block,uint8_t * dest,int stride)207 static void mpeg2_idct_copy_c (int16_t * block, uint8_t * dest, int stride)
208 {
209 int i;
210
211 for (i = 0; i < 8; i++)
212 idct_row (block + 8 * i);
213
214 for (i = 0; i < 8; i++)
215 idct_col (block + i);
216
217 i = 8;
218 do {
219 dest[0] = CLIP (block[0]);
220 dest[1] = CLIP (block[1]);
221 dest[2] = CLIP (block[2]);
222 dest[3] = CLIP (block[3]);
223 dest[4] = CLIP (block[4]);
224 dest[5] = CLIP (block[5]);
225 dest[6] = CLIP (block[6]);
226 dest[7] = CLIP (block[7]);
227
228 block[0] = 0; block[1] = 0; block[2] = 0; block[3] = 0;
229 block[4] = 0; block[5] = 0; block[6] = 0; block[7] = 0;
230
231 dest += stride;
232 block += 8;
233 } while (--i);
234 }
235
mpeg2_idct_add_c(int16_t * block,uint8_t * dest,int stride)236 static void mpeg2_idct_add_c (int16_t * block, uint8_t * dest, int stride)
237 {
238 int i;
239
240 for (i = 0; i < 8; i++)
241 idct_row (block + 8 * i);
242
243 for (i = 0; i < 8; i++)
244 idct_col (block + i);
245
246 i = 8;
247 do {
248 dest[0] = CLIP (block[0] + dest[0]);
249 dest[1] = CLIP (block[1] + dest[1]);
250 dest[2] = CLIP (block[2] + dest[2]);
251 dest[3] = CLIP (block[3] + dest[3]);
252 dest[4] = CLIP (block[4] + dest[4]);
253 dest[5] = CLIP (block[5] + dest[5]);
254 dest[6] = CLIP (block[6] + dest[6]);
255 dest[7] = CLIP (block[7] + dest[7]);
256
257 block[0] = 0; block[1] = 0; block[2] = 0; block[3] = 0;
258 block[4] = 0; block[5] = 0; block[6] = 0; block[7] = 0;
259
260 dest += stride;
261 block += 8;
262 } while (--i);
263 }
264
mpeg2_idct_c(int16_t * block)265 static void mpeg2_idct_c (int16_t * block)
266 {
267 int i;
268
269 for (i = 0; i < 8; i++)
270 idct_row (block + 8 * i);
271
272 for (i = 0; i < 8; i++)
273 idct_col (block + i);
274 }
275
mpeg2_zero_block_c(int16_t * wblock)276 static void mpeg2_zero_block_c (int16_t * wblock)
277 {
278 memset( wblock, 0, sizeof(int16_t) * 64 );
279 }
280
mpeg2_idct_init(uint32_t mm_accel)281 void mpeg2_idct_init (uint32_t mm_accel)
282 {
283 mpeg2_zero_block = mpeg2_zero_block_c;
284
285 #if defined(ARCH_X86)
286 if (mm_accel & MM_ACCEL_X86_MMXEXT) {
287 #ifdef LOG
288 fprintf (stderr, "Using MMXEXT for IDCT transform\n");
289 #endif
290 mpeg2_idct_copy = mpeg2_idct_copy_mmxext;
291 mpeg2_idct_add = mpeg2_idct_add_mmxext;
292 mpeg2_idct = mpeg2_idct_mmxext;
293 mpeg2_zero_block = mpeg2_zero_block_mmx;
294 mpeg2_idct_mmx_init ();
295 } else if (mm_accel & MM_ACCEL_X86_MMX) {
296 #ifdef LOG
297 fprintf (stderr, "Using MMX for IDCT transform\n");
298 #endif
299 mpeg2_idct_copy = mpeg2_idct_copy_mmx;
300 mpeg2_idct_add = mpeg2_idct_add_mmx;
301 mpeg2_idct = mpeg2_idct_mmx;
302 mpeg2_zero_block = mpeg2_zero_block_mmx;
303 mpeg2_idct_mmx_init ();
304 } else
305 #endif
306 #if defined (ARCH_PPC) && defined (ENABLE_ALTIVEC)
307 if (mm_accel & MM_ACCEL_PPC_ALTIVEC) {
308 #ifdef LOG
309 fprintf (stderr, "Using altivec for IDCT transform\n");
310 #endif
311 mpeg2_idct_copy = mpeg2_idct_copy_altivec;
312 mpeg2_idct_add = mpeg2_idct_add_altivec;
313 mpeg2_idct_altivec_init ();
314 mpeg2_idct = mpeg2_idct_c;
315 } else
316 #endif
317 #ifdef LIBMPEG2_MLIB
318 if (mm_accel & MM_ACCEL_MLIB) {
319 char * env_var;
320
321 env_var = getenv ("MLIB_NON_IEEE");
322
323 mpeg2_idct = mpeg2_idct_mlib;
324 if (env_var == NULL) {
325 #ifdef LOG
326 fprintf (stderr, "Using mlib for IDCT transform\n");
327 #endif
328 mpeg2_idct_add = mpeg2_idct_add_mlib;
329 } else {
330 fprintf (stderr, "Using non-IEEE mlib for IDCT transform\n");
331 mpeg2_idct_add = mpeg2_idct_add_mlib_non_ieee;
332 }
333 mpeg2_idct_copy = mpeg2_idct_copy_mlib_non_ieee;
334 } else
335 #endif
336 {
337 int i;
338
339 #ifdef LOG
340 fprintf (stderr, "No accelerated IDCT transform found\n");
341 #endif
342 mpeg2_idct_copy = mpeg2_idct_copy_c;
343 mpeg2_idct_add = mpeg2_idct_add_c;
344 mpeg2_idct = mpeg2_idct_c;
345 for (i = -384; i < 640; i++)
346 clip_lut[i+384] = (i < 0) ? 0 : ((i > 255) ? 255 : i);
347 }
348 }
349