Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Annotate
Revision Log

jidctint.c @ 21695

Visit:

Revision 21695, 31.3 KB checked in by rbasch, 20 years ago (diff)
This commit was generated by cvs2svn to compensate for changes in r21694, which included commits to RCS files with non-trunk default branches.

Line
1	/*
2	* jidctint.c
3	*
4	* Copyright (C) 1991-1998, Thomas G. Lane.
5	* This file is part of the Independent JPEG Group's software.
6	* For conditions of distribution and use, see the accompanying README file.
7	*
8	* This file contains a slow-but-accurate integer implementation of the
9	* inverse DCT (Discrete Cosine Transform). In the IJG code, this routine
10	* must also perform dequantization of the input coefficients.
11	*
12	* A 2-D IDCT can be done by 1-D IDCT on each column followed by 1-D IDCT
13	* on each row (or vice versa, but it's more convenient to emit a row at
14	* a time). Direct algorithms are also available, but they are much more
15	* complex and seem not to be any faster when reduced to code.
16	*
17	* This implementation is based on an algorithm described in
18	* C. Loeffler, A. Ligtenberg and G. Moschytz, "Practical Fast 1-D DCT
19	* Algorithms with 11 Multiplications", Proc. Int'l. Conf. on Acoustics,
20	* Speech, and Signal Processing 1989 (ICASSP '89), pp. 988-991.
21	* The primary algorithm described there uses 11 multiplies and 29 adds.
22	* We use their alternate method with 12 multiplies and 32 adds.
23	* The advantage of this method is that no data path contains more than one
24	* multiplication; this allows a very simple and accurate implementation in
25	* scaled fixed-point arithmetic, with a minimal number of shifts.
26	*/
27
28	#define JPEG_INTERNALS
29	#include "jinclude.h"
30	#include "jpeglib.h"
31	#include "jdct.h" /* Private declarations for DCT subsystem */
32
33	#ifdef DCT_ISLOW_SUPPORTED
34
35
36	/*
37	* This module is specialized to the case DCTSIZE = 8.
38	*/
39
40	#if DCTSIZE != 8
41	Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */
42	#endif
43
44
45	/*
46	* The poop on this scaling stuff is as follows:
47	*
48	* Each 1-D IDCT step produces outputs which are a factor of sqrt(N)
49	* larger than the true IDCT outputs. The final outputs are therefore
50	* a factor of N larger than desired; since N=8 this can be cured by
51	* a simple right shift at the end of the algorithm. The advantage of
52	* this arrangement is that we save two multiplications per 1-D IDCT,
53	* because the y0 and y4 inputs need not be divided by sqrt(N).
54	*
55	* We have to do addition and subtraction of the integer inputs, which
56	* is no problem, and multiplication by fractional constants, which is
57	* a problem to do in integer arithmetic. We multiply all the constants
58	* by CONST_SCALE and convert them to integer constants (thus retaining
59	* CONST_BITS bits of precision in the constants). After doing a
60	* multiplication we have to divide the product by CONST_SCALE, with proper
61	* rounding, to produce the correct output. This division can be done
62	* cheaply as a right shift of CONST_BITS bits. We postpone shifting
63	* as long as possible so that partial sums can be added together with
64	* full fractional precision.
65	*
66	* The outputs of the first pass are scaled up by PASS1_BITS bits so that
67	* they are represented to better-than-integral precision. These outputs
68	* require BITS_IN_JSAMPLE + PASS1_BITS + 3 bits; this fits in a 16-bit word
69	* with the recommended scaling. (To scale up 12-bit sample data further, an
70	* intermediate INT32 array would be needed.)
71	*
72	* To avoid overflow of the 32-bit intermediate results in pass 2, we must
73	* have BITS_IN_JSAMPLE + CONST_BITS + PASS1_BITS <= 26. Error analysis
74	* shows that the values given below are the most effective.
75	*/
76
77	#if BITS_IN_JSAMPLE == 8
78	#define CONST_BITS 13
79	#define PASS1_BITS 2
80	#else
81	#define CONST_BITS 13
82	#define PASS1_BITS 1 /* lose a little precision to avoid overflow */
83	#endif
84
85	/* Some C compilers fail to reduce "FIX(constant)" at compile time, thus
86	* causing a lot of useless floating-point operations at run time.
87	* To get around this we use the following pre-calculated constants.
88	* If you change CONST_BITS you may want to add appropriate values.
89	* (With a reasonable C compiler, you can just rely on the FIX() macro...)
90	*/
91
92	#if CONST_BITS == 13
93	#define FIX_0_298631336 ((INT32) 2446) /* FIX(0.298631336) */
94	#define FIX_0_390180644 ((INT32) 3196) /* FIX(0.390180644) */
95	#define FIX_0_541196100 ((INT32) 4433) /* FIX(0.541196100) */
96	#define FIX_0_765366865 ((INT32) 6270) /* FIX(0.765366865) */
97	#define FIX_0_899976223 ((INT32) 7373) /* FIX(0.899976223) */
98	#define FIX_1_175875602 ((INT32) 9633) /* FIX(1.175875602) */
99	#define FIX_1_501321110 ((INT32) 12299) /* FIX(1.501321110) */
100	#define FIX_1_847759065 ((INT32) 15137) /* FIX(1.847759065) */
101	#define FIX_1_961570560 ((INT32) 16069) /* FIX(1.961570560) */
102	#define FIX_2_053119869 ((INT32) 16819) /* FIX(2.053119869) */
103	#define FIX_2_562915447 ((INT32) 20995) /* FIX(2.562915447) */
104	#define FIX_3_072711026 ((INT32) 25172) /* FIX(3.072711026) */
105	#else
106	#define FIX_0_298631336 FIX(0.298631336)
107	#define FIX_0_390180644 FIX(0.390180644)
108	#define FIX_0_541196100 FIX(0.541196100)
109	#define FIX_0_765366865 FIX(0.765366865)
110	#define FIX_0_899976223 FIX(0.899976223)
111	#define FIX_1_175875602 FIX(1.175875602)
112	#define FIX_1_501321110 FIX(1.501321110)
113	#define FIX_1_847759065 FIX(1.847759065)
114	#define FIX_1_961570560 FIX(1.961570560)
115	#define FIX_2_053119869 FIX(2.053119869)
116	#define FIX_2_562915447 FIX(2.562915447)
117	#define FIX_3_072711026 FIX(3.072711026)
118	#endif
119
120
121	/* Multiply an INT32 variable by an INT32 constant to yield an INT32 result.
122	* For 8-bit samples with the recommended scaling, all the variable
123	* and constant values involved are no more than 16 bits wide, so a
124	* 16x16->32 bit multiply can be used instead of a full 32x32 multiply.
125	* For 12-bit samples, a full 32-bit multiplication will be needed.
126	*/
127
128	#if BITS_IN_JSAMPLE == 8
129	#define MULTIPLY(var,const) MULTIPLY16C16(var,const)
130	#else
131	#define MULTIPLY(var,const) ((var) * (const))
132	#endif
133
134
135	/* Dequantize a coefficient by multiplying it by the multiplier-table
136	* entry; produce an int result. In this module, both inputs and result
137	* are 16 bits or less, so either int or short multiply will work.
138	*/
139
140	#define DEQUANTIZE(coef,quantval) (((ISLOW_MULT_TYPE) (coef)) * (quantval))
141
142
143	/*
144	* Perform dequantization and inverse DCT on one block of coefficients.
145	*/
146
147	GLOBAL(void)
148	jpeg_idct_islow (j_decompress_ptr cinfo, jpeg_component_info * compptr,
149	JCOEFPTR coef_block,
150	JSAMPARRAY output_buf, JDIMENSION output_col)
151	{
152	INT32 tmp0, tmp1, tmp2, tmp3;
153	INT32 tmp10, tmp11, tmp12, tmp13;
154	INT32 z1, z2, z3, z4, z5;
155	JCOEFPTR inptr;
156	ISLOW_MULT_TYPE * quantptr;
157	int * wsptr;
158	JSAMPROW outptr;
159	JSAMPLE *range_limit = IDCT_range_limit(cinfo);
160	int ctr;
161	int workspace[DCTSIZE2]; /* buffers data between passes */
162	SHIFT_TEMPS
163
164	/* Pass 1: process columns from input, store into work array. */
165	/* Note results are scaled up by sqrt(8) compared to a true IDCT; */
166	/* furthermore, we scale the results by 2*PASS1_BITS. /
167
168	inptr = coef_block;
169	quantptr = (ISLOW_MULT_TYPE *) compptr->dct_table;
170	wsptr = workspace;
171	for (ctr = DCTSIZE; ctr > 0; ctr--) {
172	/* Due to quantization, we will usually find that many of the input
173	* coefficients are zero, especially the AC terms. We can exploit this
174	* by short-circuiting the IDCT calculation for any column in which all
175	* the AC terms are zero. In that case each output is equal to the
176	* DC coefficient (with scale factor as needed).
177	* With typical images and quantization tables, half or more of the
178	* column DCT calculations can be simplified this way.
179	*/
180
181	if (inptr[DCTSIZE1] == 0 && inptr[DCTSIZE2] == 0 &&
182	inptr[DCTSIZE3] == 0 && inptr[DCTSIZE4] == 0 &&
183	inptr[DCTSIZE5] == 0 && inptr[DCTSIZE6] == 0 &&
184	inptr[DCTSIZE*7] == 0) {
185	/* AC terms all zero */
186	int dcval = DEQUANTIZE(inptr[DCTSIZE0], quantptr[DCTSIZE0]) << PASS1_BITS;
187
188	wsptr[DCTSIZE*0] = dcval;
189	wsptr[DCTSIZE*1] = dcval;
190	wsptr[DCTSIZE*2] = dcval;
191	wsptr[DCTSIZE*3] = dcval;
192	wsptr[DCTSIZE*4] = dcval;
193	wsptr[DCTSIZE*5] = dcval;
194	wsptr[DCTSIZE*6] = dcval;
195	wsptr[DCTSIZE*7] = dcval;
196
197	inptr++; /* advance pointers to next column */
198	quantptr++;
199	wsptr++;
200	continue;
201	}
202
203	/* Even part: reverse the even part of the forward DCT. */
204	/* The rotator is sqrt(2)c(-6). /
205
206	z2 = DEQUANTIZE(inptr[DCTSIZE2], quantptr[DCTSIZE2]);
207	z3 = DEQUANTIZE(inptr[DCTSIZE6], quantptr[DCTSIZE6]);
208
209	z1 = MULTIPLY(z2 + z3, FIX_0_541196100);
210	tmp2 = z1 + MULTIPLY(z3, - FIX_1_847759065);
211	tmp3 = z1 + MULTIPLY(z2, FIX_0_765366865);
212
213	z2 = DEQUANTIZE(inptr[DCTSIZE0], quantptr[DCTSIZE0]);
214	z3 = DEQUANTIZE(inptr[DCTSIZE4], quantptr[DCTSIZE4]);
215
216	tmp0 = (z2 + z3) << CONST_BITS;
217	tmp1 = (z2 - z3) << CONST_BITS;
218
219	tmp10 = tmp0 + tmp3;
220	tmp13 = tmp0 - tmp3;
221	tmp11 = tmp1 + tmp2;
222	tmp12 = tmp1 - tmp2;
223
224	/* Odd part per figure 8; the matrix is unitary and hence its
225	* transpose is its inverse. i0..i3 are y7,y5,y3,y1 respectively.
226	*/
227
228	tmp0 = DEQUANTIZE(inptr[DCTSIZE7], quantptr[DCTSIZE7]);
229	tmp1 = DEQUANTIZE(inptr[DCTSIZE5], quantptr[DCTSIZE5]);
230	tmp2 = DEQUANTIZE(inptr[DCTSIZE3], quantptr[DCTSIZE3]);
231	tmp3 = DEQUANTIZE(inptr[DCTSIZE1], quantptr[DCTSIZE1]);
232
233	z1 = tmp0 + tmp3;
234	z2 = tmp1 + tmp2;
235	z3 = tmp0 + tmp2;
236	z4 = tmp1 + tmp3;
237	z5 = MULTIPLY(z3 + z4, FIX_1_175875602); /* sqrt(2) * c3 */
238
239	tmp0 = MULTIPLY(tmp0, FIX_0_298631336); /* sqrt(2) * (-c1+c3+c5-c7) */
240	tmp1 = MULTIPLY(tmp1, FIX_2_053119869); /* sqrt(2) * ( c1+c3-c5+c7) */
241	tmp2 = MULTIPLY(tmp2, FIX_3_072711026); /* sqrt(2) * ( c1+c3+c5-c7) */
242	tmp3 = MULTIPLY(tmp3, FIX_1_501321110); /* sqrt(2) * ( c1+c3-c5-c7) */
243	z1 = MULTIPLY(z1, - FIX_0_899976223); /* sqrt(2) * (c7-c3) */
244	z2 = MULTIPLY(z2, - FIX_2_562915447); /* sqrt(2) * (-c1-c3) */
245	z3 = MULTIPLY(z3, - FIX_1_961570560); /* sqrt(2) * (-c3-c5) */
246	z4 = MULTIPLY(z4, - FIX_0_390180644); /* sqrt(2) * (c5-c3) */
247
248	z3 += z5;
249	z4 += z5;
250
251	tmp0 += z1 + z3;
252	tmp1 += z2 + z4;
253	tmp2 += z2 + z3;
254	tmp3 += z1 + z4;
255
256	/* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
257
258	wsptr[DCTSIZE*0] = (int) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
259	wsptr[DCTSIZE*7] = (int) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
260	wsptr[DCTSIZE*1] = (int) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
261	wsptr[DCTSIZE*6] = (int) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
262	wsptr[DCTSIZE*2] = (int) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
263	wsptr[DCTSIZE*5] = (int) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
264	wsptr[DCTSIZE*3] = (int) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
265	wsptr[DCTSIZE*4] = (int) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);
266
267	inptr++; /* advance pointers to next column */
268	quantptr++;
269	wsptr++;
270	}
271
272	/* Pass 2: process rows from work array, store into output array. */
273	/* Note that we must descale the results by a factor of 8 == 2*3, /
274	/* and also undo the PASS1_BITS scaling. */
275
276	wsptr = workspace;
277	for (ctr = 0; ctr < DCTSIZE; ctr++) {
278	outptr = output_buf[ctr] + output_col;
279	/* Rows of zeroes can be exploited in the same way as we did with columns.
280	* However, the column calculation has created many nonzero AC terms, so
281	* the simplification applies less often (typically 5% to 10% of the time).
282	* On machines with very fast multiplication, it's possible that the
283	* test takes more time than it's worth. In that case this section
284	* may be commented out.
285	*/
286
287	#ifndef NO_ZERO_ROW_TEST
288	if (wsptr[1] == 0 && wsptr[2] == 0 && wsptr[3] == 0 && wsptr[4] == 0 &&
289	wsptr[5] == 0 && wsptr[6] == 0 && wsptr[7] == 0) {
290	/* AC terms all zero */
291	JSAMPLE dcval = range_limit[(int) DESCALE((INT32) wsptr[0], PASS1_BITS+3)
292	& RANGE_MASK];
293
294	outptr[0] = dcval;
295	outptr[1] = dcval;
296	outptr[2] = dcval;
297	outptr[3] = dcval;
298	outptr[4] = dcval;
299	outptr[5] = dcval;
300	outptr[6] = dcval;
301	outptr[7] = dcval;
302
303	wsptr += DCTSIZE; /* advance pointer to next row */
304	continue;
305	}
306	#endif
307
308	/* Even part: reverse the even part of the forward DCT. */
309	/* The rotator is sqrt(2)c(-6). /
310
311	z2 = (INT32) wsptr[2];
312	z3 = (INT32) wsptr[6];
313
314	z1 = MULTIPLY(z2 + z3, FIX_0_541196100);
315	tmp2 = z1 + MULTIPLY(z3, - FIX_1_847759065);
316	tmp3 = z1 + MULTIPLY(z2, FIX_0_765366865);
317
318	tmp0 = ((INT32) wsptr[0] + (INT32) wsptr[4]) << CONST_BITS;
319	tmp1 = ((INT32) wsptr[0] - (INT32) wsptr[4]) << CONST_BITS;
320
321	tmp10 = tmp0 + tmp3;
322	tmp13 = tmp0 - tmp3;
323	tmp11 = tmp1 + tmp2;
324	tmp12 = tmp1 - tmp2;
325
326	/* Odd part per figure 8; the matrix is unitary and hence its
327	* transpose is its inverse. i0..i3 are y7,y5,y3,y1 respectively.
328	*/
329
330	tmp0 = (INT32) wsptr[7];
331	tmp1 = (INT32) wsptr[5];
332	tmp2 = (INT32) wsptr[3];
333	tmp3 = (INT32) wsptr[1];
334
335	z1 = tmp0 + tmp3;
336	z2 = tmp1 + tmp2;
337	z3 = tmp0 + tmp2;
338	z4 = tmp1 + tmp3;
339	z5 = MULTIPLY(z3 + z4, FIX_1_175875602); /* sqrt(2) * c3 */
340
341	tmp0 = MULTIPLY(tmp0, FIX_0_298631336); /* sqrt(2) * (-c1+c3+c5-c7) */
342	tmp1 = MULTIPLY(tmp1, FIX_2_053119869); /* sqrt(2) * ( c1+c3-c5+c7) */
343	tmp2 = MULTIPLY(tmp2, FIX_3_072711026); /* sqrt(2) * ( c1+c3+c5-c7) */
344	tmp3 = MULTIPLY(tmp3, FIX_1_501321110); /* sqrt(2) * ( c1+c3-c5-c7) */
345	z1 = MULTIPLY(z1, - FIX_0_899976223); /* sqrt(2) * (c7-c3) */
346	z2 = MULTIPLY(z2, - FIX_2_562915447); /* sqrt(2) * (-c1-c3) */
347	z3 = MULTIPLY(z3, - FIX_1_961570560); /* sqrt(2) * (-c3-c5) */
348	z4 = MULTIPLY(z4, - FIX_0_390180644); /* sqrt(2) * (c5-c3) */
349
350	z3 += z5;
351	z4 += z5;
352
353	tmp0 += z1 + z3;
354	tmp1 += z2 + z4;
355	tmp2 += z2 + z3;
356	tmp3 += z1 + z4;
357
358	/* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
359
360	outptr[0] = range_limit[(int) DESCALE(tmp10 + tmp3,
361	CONST_BITS+PASS1_BITS+3)
362	& RANGE_MASK];
363	outptr[7] = range_limit[(int) DESCALE(tmp10 - tmp3,
364	CONST_BITS+PASS1_BITS+3)
365	& RANGE_MASK];
366	outptr[1] = range_limit[(int) DESCALE(tmp11 + tmp2,
367	CONST_BITS+PASS1_BITS+3)
368	& RANGE_MASK];
369	outptr[6] = range_limit[(int) DESCALE(tmp11 - tmp2,
370	CONST_BITS+PASS1_BITS+3)
371	& RANGE_MASK];
372	outptr[2] = range_limit[(int) DESCALE(tmp12 + tmp1,
373	CONST_BITS+PASS1_BITS+3)
374	& RANGE_MASK];
375	outptr[5] = range_limit[(int) DESCALE(tmp12 - tmp1,
376	CONST_BITS+PASS1_BITS+3)
377	& RANGE_MASK];
378	outptr[3] = range_limit[(int) DESCALE(tmp13 + tmp0,
379	CONST_BITS+PASS1_BITS+3)
380	& RANGE_MASK];
381	outptr[4] = range_limit[(int) DESCALE(tmp13 - tmp0,
382	CONST_BITS+PASS1_BITS+3)
383	& RANGE_MASK];
384
385	wsptr += DCTSIZE; /* advance pointer to next row */
386	}
387	}
388
389
390	#ifdef HAVE_SSE2_INTEL_MNEMONICS
391
392	/*
393	* Intel SSE2 optimized Inverse Discrete Cosine Transform
394	*
395	*
396	* Copyright (c) 2001-2002 Intel Corporation
397	* All Rights Reserved
398	*
399	*
400	* Authors:
401	* Danilov G.
402	*
403	*
404	*-----------------------------------------------------------------------------
405	*
406	* References:
407	* K.R. Rao and P. Yip
408	* Discrete Cosine Transform.
409	* Algorithms, Advantages, Applications.
410	* Academic Press, Inc, London, 1990.
411	* JPEG Group's software.
412	* This implementation is based on Appendix A.2 of the book (R&Y) ...
413	*
414	*-----------------------------------------------------------------------------
415	*/
416
417	typedef unsigned char Ipp8u;
418	typedef unsigned short Ipp16u;
419	typedef unsigned int Ipp32u;
420
421	typedef signed char Ipp8s;
422	typedef signed short Ipp16s;
423	typedef signed int Ipp32s;
424
425	#define BITS_INV_ACC 4
426	#define SHIFT_INV_ROW 16 - BITS_INV_ACC
427	#define SHIFT_INV_COL 1 + BITS_INV_ACC
428
429	#define RND_INV_ROW 1024 * (6 - BITS_INV_ACC) /* 1 << (SHIFT_INV_ROW-1) */
430	#define RND_INV_COL = 16 * (BITS_INV_ACC - 3) /* 1 << (SHIFT_INV_COL-1) */
431	#define RND_INV_CORR = RND_INV_COL - 1 /* correction -1.0 and round */
432
433	#define c_inv_corr_0 -1024 * (6 - BITS_INV_ACC) + 65536 /* -0.5 + (16.0 or 32.0) */
434	#define c_inv_corr_1 1877 * (6 - BITS_INV_ACC) /* 0.9167 */
435	#define c_inv_corr_2 1236 * (6 - BITS_INV_ACC) /* 0.6035 */
436	#define c_inv_corr_3 680 * (6 - BITS_INV_ACC) /* 0.3322 */
437	#define c_inv_corr_4 0 * (6 - BITS_INV_ACC) /* 0.0 */
438	#define c_inv_corr_5 -569 * (6 - BITS_INV_ACC) /* -0.278 */
439	#define c_inv_corr_6 -512 * (6 - BITS_INV_ACC) /* -0.25 */
440	#define c_inv_corr_7 -651 * (6 - BITS_INV_ACC) /* -0.3176 */
441
442	#define RND_INV_ROW_0 RND_INV_ROW + c_inv_corr_0
443	#define RND_INV_ROW_1 RND_INV_ROW + c_inv_corr_1
444	#define RND_INV_ROW_2 RND_INV_ROW + c_inv_corr_2
445	#define RND_INV_ROW_3 RND_INV_ROW + c_inv_corr_3
446	#define RND_INV_ROW_4 RND_INV_ROW + c_inv_corr_4
447	#define RND_INV_ROW_5 RND_INV_ROW + c_inv_corr_5
448	#define RND_INV_ROW_6 RND_INV_ROW + c_inv_corr_6
449	#define RND_INV_ROW_7 RND_INV_ROW + c_inv_corr_7
450
451	/* Table for rows 0,4 - constants are multiplied on cos_4_16 */
452
453	__declspec(align(16)) short tab_i_04[] = {
454	16384, 21407, 16384, 8867,
455	-16384, 21407, 16384, -8867,
456	16384, -8867, 16384, -21407,
457	16384, 8867, -16384, -21407,
458	22725, 19266, 19266, -4520,
459	4520, 19266, 19266, -22725,
460	12873, -22725, 4520, -12873,
461	12873, 4520, -22725, -12873};
462
463	/* Table for rows 1,7 - constants are multiplied on cos_1_16 */
464
465	__declspec(align(16)) short tab_i_17[] = {
466	22725, 29692, 22725, 12299,
467	-22725, 29692, 22725, -12299,
468	22725, -12299, 22725, -29692,
469	22725, 12299, -22725, -29692,
470	31521, 26722, 26722, -6270,
471	6270, 26722, 26722, -31521,
472	17855, -31521, 6270, -17855,
473	17855, 6270, -31521, -17855};
474
475	/* Table for rows 2,6 - constants are multiplied on cos_2_16 */
476
477	__declspec(align(16)) short tab_i_26[] = {
478	21407, 27969, 21407, 11585,
479	-21407, 27969, 21407, -11585,
480	21407, -11585, 21407, -27969,
481	21407, 11585, -21407, -27969,
482	29692, 25172, 25172, -5906,
483	5906, 25172, 25172, -29692,
484	16819, -29692, 5906, -16819,
485	16819, 5906, -29692, -16819};
486
487	/* Table for rows 3,5 - constants are multiplied on cos_3_16 */
488
489	__declspec(align(16)) short tab_i_35[] = {
490	19266, 25172, 19266, 10426,
491	-19266, 25172, 19266, -10426,
492	19266, -10426, 19266, -25172,
493	19266, 10426, -19266, -25172,
494	26722, 22654, 22654, -5315,
495	5315, 22654, 22654, -26722,
496	15137, -26722, 5315, -15137,
497	15137, 5315, -26722, -15137};
498
499	__declspec(align(16)) long round_i_0[] = {RND_INV_ROW_0,RND_INV_ROW_0,
500	RND_INV_ROW_0,RND_INV_ROW_0};
501	__declspec(align(16)) long round_i_1[] = {RND_INV_ROW_1,RND_INV_ROW_1,
502	RND_INV_ROW_1,RND_INV_ROW_1};
503	__declspec(align(16)) long round_i_2[] = {RND_INV_ROW_2,RND_INV_ROW_2,
504	RND_INV_ROW_2,RND_INV_ROW_2};
505	__declspec(align(16)) long round_i_3[] = {RND_INV_ROW_3,RND_INV_ROW_3,
506	RND_INV_ROW_3,RND_INV_ROW_3};
507	__declspec(align(16)) long round_i_4[] = {RND_INV_ROW_4,RND_INV_ROW_4,
508	RND_INV_ROW_4,RND_INV_ROW_4};
509	__declspec(align(16)) long round_i_5[] = {RND_INV_ROW_5,RND_INV_ROW_5,
510	RND_INV_ROW_5,RND_INV_ROW_5};
511	__declspec(align(16)) long round_i_6[] = {RND_INV_ROW_6,RND_INV_ROW_6,
512	RND_INV_ROW_6,RND_INV_ROW_6};
513	__declspec(align(16)) long round_i_7[] = {RND_INV_ROW_7,RND_INV_ROW_7,
514	RND_INV_ROW_7,RND_INV_ROW_7};
515
516	__declspec(align(16)) short tg_1_16[] = {
517	13036, 13036, 13036, 13036, /* tg * (2<<16) + 0.5 */
518	13036, 13036, 13036, 13036};
519	__declspec(align(16)) short tg_2_16[] = {
520	27146, 27146, 27146, 27146, /* tg * (2<<16) + 0.5 */
521	27146, 27146, 27146, 27146};
522	__declspec(align(16)) short tg_3_16[] = {
523	-21746, -21746, -21746, -21746, /* tg * (2<<16) + 0.5 */
524	-21746, -21746, -21746, -21746};
525	__declspec(align(16)) short cos_4_16[] = {
526	-19195, -19195, -19195, -19195, /* cos * (2<<16) + 0.5 */
527	-19195, -19195, -19195, -19195};
528
529	/*
530	* In this implementation the outputs of the iDCT-1D are multiplied
531	* for rows 0,4 - on cos_4_16,
532	* for rows 1,7 - on cos_1_16,
533	* for rows 2,6 - on cos_2_16,
534	* for rows 3,5 - on cos_3_16
535	* and are shifted to the left for rise of accuracy
536	*
537	* For used constants
538	* FIX(float_const) = (short) (float_const * (1<<15) + 0.5)
539	*
540	*-----------------------------------------------------------------------------
541	*
542	* On the first stage the calculation is executed at once for two rows.
543	* The permutation for each output row is done on second stage
544	* t7 t6 t5 t4 t3 t2 t1 t0 -> t4 t5 t6 t7 t3 t2 t1 t0
545	*
546	*-----------------------------------------------------------------------------
547	*/
548
549	#define DCT_8_INV_ROW_2R(TABLE, ROUND1, ROUND2) __asm { \
550	__asm pshuflw xmm1, xmm0, 10001000b \
551	__asm pshuflw xmm0, xmm0, 11011101b \
552	__asm pshufhw xmm1, xmm1, 10001000b \
553	__asm pshufhw xmm0, xmm0, 11011101b \
554	__asm movdqa xmm2, XMMWORD PTR [TABLE] \
555	__asm pmaddwd xmm2, xmm1 \
556	__asm movdqa xmm3, XMMWORD PTR [TABLE + 32] \
557	__asm pmaddwd xmm3, xmm0 \
558	__asm pmaddwd xmm1, XMMWORD PTR [TABLE + 16] \
559	__asm pmaddwd xmm0, XMMWORD PTR [TABLE + 48] \
560	__asm pshuflw xmm5, xmm4, 10001000b \
561	__asm pshuflw xmm4, xmm4, 11011101b \
562	__asm pshufhw xmm5, xmm5, 10001000b \
563	__asm pshufhw xmm4, xmm4, 11011101b \
564	__asm movdqa xmm6, XMMWORD PTR [TABLE] \
565	__asm pmaddwd xmm6, xmm5 \
566	__asm movdqa xmm7, XMMWORD PTR [TABLE + 32] \
567	__asm pmaddwd xmm7, xmm4 \
568	__asm pmaddwd xmm5, XMMWORD PTR [TABLE + 16] \
569	__asm pmaddwd xmm4, XMMWORD PTR [TABLE + 48] \
570	__asm pshufd xmm1, xmm1, 01001110b \
571	__asm pshufd xmm0, xmm0, 01001110b \
572	__asm paddd xmm2, XMMWORD PTR [ROUND1] \
573	__asm paddd xmm3, xmm0 \
574	__asm paddd xmm1, xmm2 \
575	__asm pshufd xmm5, xmm5, 01001110b \
576	__asm pshufd xmm4, xmm4, 01001110b \
577	__asm movdqa xmm2, xmm1 \
578	__asm psubd xmm2, xmm3 \
579	__asm psrad xmm2, SHIFT_INV_ROW \
580	__asm paddd xmm1, xmm3 \
581	__asm psrad xmm1, SHIFT_INV_ROW \
582	__asm packssdw xmm1, xmm2 \
583	__asm paddd xmm6, XMMWORD PTR [ROUND2] \
584	__asm paddd xmm7, xmm4 \
585	__asm paddd xmm5, xmm6 \
586	__asm movdqa xmm6, xmm5 \
587	__asm psubd xmm6, xmm7 \
588	__asm psrad xmm6, SHIFT_INV_ROW \
589	__asm paddd xmm5, xmm7 \
590	__asm psrad xmm5, SHIFT_INV_ROW \
591	__asm packssdw xmm5, xmm6 \
592	}
593
594	/*
595	*
596	* The second stage - inverse DCTs of columns
597	*
598	* The inputs are multiplied
599	* for rows 0,4 - on cos_4_16,
600	* for rows 1,7 - on cos_1_16,
601	* for rows 2,6 - on cos_2_16,
602	* for rows 3,5 - on cos_3_16
603	* and are shifted to the left for rise of accuracy
604	*/
605
606	#define DCT_8_INV_COL_8R(INP, OUTP) __asm { \
607	__asm movdqa xmm0, [INP + 5*16] \
608	__asm movdqa xmm1, XMMWORD PTR tg_3_16 \
609	__asm movdqa xmm2, xmm0 \
610	__asm movdqa xmm3, [INP + 3*16] \
611	__asm pmulhw xmm0, xmm1 \
612	__asm movdqa xmm4, [INP + 7*16] \
613	__asm pmulhw xmm1, xmm3 \
614	__asm movdqa xmm5, XMMWORD PTR tg_1_16 \
615	__asm movdqa xmm6, xmm4 \
616	__asm pmulhw xmm4, xmm5 \
617	__asm paddsw xmm0, xmm2 \
618	__asm pmulhw xmm5, [INP + 1*16] \
619	__asm paddsw xmm1, xmm3 \
620	__asm movdqa xmm7, [INP + 6*16] \
621	__asm paddsw xmm0, xmm3 \
622	__asm movdqa xmm3, XMMWORD PTR tg_2_16 \
623	__asm psubsw xmm2, xmm1 \
624	__asm pmulhw xmm7, xmm3 \
625	__asm movdqa xmm1, xmm0 \
626	__asm pmulhw xmm3, [INP + 2*16] \
627	__asm psubsw xmm5, xmm6 \
628	__asm paddsw xmm4, [INP + 1*16] \
629	__asm paddsw xmm0, xmm4 \
630	__asm psubsw xmm4, xmm1 \
631	__asm pshufhw xmm0, xmm0, 00011011b \
632	__asm paddsw xmm7, [INP + 2*16] \
633	__asm movdqa xmm6, xmm5 \
634	__asm psubsw xmm3, [INP + 6*16] \
635	__asm psubsw xmm5, xmm2 \
636	__asm paddsw xmm6, xmm2 \
637	__asm movdqa [OUTP + 7*16], xmm0 \
638	__asm movdqa xmm1, xmm4 \
639	__asm movdqa xmm2, XMMWORD PTR cos_4_16 \
640	__asm paddsw xmm4, xmm5 \
641	__asm movdqa xmm0, XMMWORD PTR cos_4_16 \
642	__asm pmulhw xmm2, xmm4 \
643	__asm pshufhw xmm6, xmm6, 00011011b \
644	__asm movdqa [OUTP + 3*16], xmm6 \
645	__asm psubsw xmm1, xmm5 \
646	__asm movdqa xmm6, [INP + 0*16] \
647	__asm pmulhw xmm0, xmm1 \
648	__asm movdqa xmm5, [INP + 4*16] \
649	__asm paddsw xmm4, xmm2 \
650	__asm paddsw xmm5, xmm6 \
651	__asm psubsw xmm6, [INP + 4*16] \
652	__asm paddsw xmm0, xmm1 \
653	__asm pshufhw xmm4, xmm4, 00011011b \
654	__asm movdqa xmm2, xmm5 \
655	__asm paddsw xmm5, xmm7 \
656	__asm movdqa xmm1, xmm6 \
657	__asm psubsw xmm2, xmm7 \
658	__asm movdqa xmm7, [OUTP + 7*16] \
659	__asm paddsw xmm6, xmm3 \
660	__asm pshufhw xmm5, xmm5, 00011011b \
661	__asm paddsw xmm7, xmm5 \
662	__asm psubsw xmm1, xmm3 \
663	__asm pshufhw xmm6, xmm6, 00011011b \
664	__asm movdqa xmm3, xmm6 \
665	__asm paddsw xmm6, xmm4 \
666	__asm pshufhw xmm2, xmm2, 00011011b \
667	__asm psraw xmm7, SHIFT_INV_COL \
668	__asm movdqa [OUTP + 0*16], xmm7 \
669	__asm movdqa xmm7, xmm1 \
670	__asm paddsw xmm1, xmm0 \
671	__asm psraw xmm6, SHIFT_INV_COL \
672	__asm movdqa [OUTP + 1*16], xmm6 \
673	__asm pshufhw xmm1, xmm1, 00011011b \
674	__asm movdqa xmm6, [OUTP + 3*16] \
675	__asm psubsw xmm7, xmm0 \
676	__asm psraw xmm1, SHIFT_INV_COL \
677	__asm movdqa [OUTP + 2*16], xmm1 \
678	__asm psubsw xmm5, [OUTP + 7*16] \
679	__asm paddsw xmm6, xmm2 \
680	__asm psubsw xmm2, [OUTP + 3*16] \
681	__asm psubsw xmm3, xmm4 \
682	__asm psraw xmm7, SHIFT_INV_COL \
683	__asm pshufhw xmm7, xmm7, 00011011b \
684	__asm movdqa [OUTP + 5*16], xmm7 \
685	__asm psraw xmm5, SHIFT_INV_COL \
686	__asm movdqa [OUTP + 7*16], xmm5 \
687	__asm psraw xmm6, SHIFT_INV_COL \
688	__asm movdqa [OUTP + 3*16], xmm6 \
689	__asm psraw xmm2, SHIFT_INV_COL \
690	__asm movdqa [OUTP + 4*16], xmm2 \
691	__asm psraw xmm3, SHIFT_INV_COL \
692	__asm movdqa [OUTP + 6*16], xmm3 \
693	}
694
695	/*
696	*
697	* Name: dct_8x8_inv_16s
698	* Purpose: Inverse Discrete Cosine Transform 8x8 with
699	* 2D buffer of short int data
700	* Context:
701	* void dct_8x8_inv_16s ( short src, short dst )
702	* Parameters:
703	* src - Pointer to the source buffer
704	* dst - Pointer to the destination buffer
705	*
706	*/
707
708	GLOBAL(void)
709	dct_8x8_inv_16s ( short src, short dst ) {
710
711	__asm {
712
713	mov ecx, src
714	mov edx, dst
715
716	movdqa xmm0, [ecx+0*16]
717	movdqa xmm4, [ecx+4*16]
718	DCT_8_INV_ROW_2R(tab_i_04, round_i_0, round_i_4)
719	movdqa [edx+0*16], xmm1
720	movdqa [edx+4*16], xmm5
721
722	movdqa xmm0, [ecx+1*16]
723	movdqa xmm4, [ecx+7*16]
724	DCT_8_INV_ROW_2R(tab_i_17, round_i_1, round_i_7)
725	movdqa [edx+1*16], xmm1
726	movdqa [edx+7*16], xmm5
727
728	movdqa xmm0, [ecx+3*16]
729	movdqa xmm4, [ecx+5*16]
730	DCT_8_INV_ROW_2R(tab_i_35, round_i_3, round_i_5);
731	movdqa [edx+3*16], xmm1
732	movdqa [edx+5*16], xmm5
733
734	movdqa xmm0, [ecx+2*16]
735	movdqa xmm4, [ecx+6*16]
736	DCT_8_INV_ROW_2R(tab_i_26, round_i_2, round_i_6);
737	movdqa [edx+2*16], xmm1
738	movdqa [edx+6*16], xmm5
739
740	DCT_8_INV_COL_8R(edx+0, edx+0);
741	}
742	}
743
744
745	/*
746	* Name:
747	* ownpj_QuantInv_8x8_16s
748	*
749	* Purpose:
750	* Dequantize 8x8 block of DCT coefficients
751	*
752	* Context:
753	* void ownpj_QuantInv_8x8_16s
754	* Ipp16s* pSrc,
755	* Ipp16s* pDst,
756	* const Ipp16u* pQTbl)*
757	*
758	*/
759
760	GLOBAL(void)
761	ownpj_QuantInv_8x8_16s(short * pSrc, short * pDst, const unsigned short * pQTbl)
762	{
763	__asm {
764
765	push ebx
766	push ecx
767	push edx
768	push esi
769	push edi
770
771	mov esi, pSrc
772	mov edi, pDst
773	mov edx, pQTbl
774	mov ecx, 4
775	mov ebx, 32
776
777	again:
778
779	movq mm0, QWORD PTR [esi+0]
780	movq mm1, QWORD PTR [esi+8]
781	movq mm2, QWORD PTR [esi+16]
782	movq mm3, QWORD PTR [esi+24]
783
784	prefetcht0 [esi+ebx] ; fetch next cache line
785
786	pmullw mm0, QWORD PTR [edx+0]
787	pmullw mm1, QWORD PTR [edx+8]
788	pmullw mm2, QWORD PTR [edx+16]
789	pmullw mm3, QWORD PTR [edx+24]
790
791	movq QWORD PTR [edi+0], mm0
792	movq QWORD PTR [edi+8], mm1
793	movq QWORD PTR [edi+16], mm2
794	movq QWORD PTR [edi+24], mm3
795
796	add esi, ebx
797	add edi, ebx
798	add edx, ebx
799	dec ecx
800	jnz again
801
802	emms
803
804	pop edi
805	pop esi
806	pop edx
807	pop ecx
808	pop ebx
809	}
810	}
811
812
813	/*
814	* Name:
815	* ownpj_Add128_8x8_16s8u
816	*
817	* Purpose:
818	* signed to unsigned conversion (level shift)
819	* for 8x8 block of DCT coefficients
820	*
821	* Context:
822	* void ownpj_Add128_8x8_16s8u
823	* const Ipp16s* pSrc,
824	* Ipp8u* pDst,
825	* int DstStep);
826	*
827	*/
828
829	__declspec(align(16)) long const_128[]= {0x00800080, 0x00800080, 0x00800080, 0x00800080};
830
831	GLOBAL(void)
832	ownpj_Add128_8x8_16s8u(const short * pSrc, unsigned char * pDst, int DstStep)
833	{
834	__asm {
835	push eax
836	push ebx
837	push ecx
838	push edx
839	push esi
840	push edi
841
842	mov esi, pSrc
843	mov edi, pDst
844	mov edx, DstStep
845	mov ecx, 2
846	mov ebx, edx
847	mov eax, edx
848	sal ebx, 1
849	add eax, ebx
850	movdqa xmm7, XMMWORD PTR const_128
851
852	again:
853
854	movdqa xmm0, XMMWORD PTR [esi+0] ; line 0
855	movdqa xmm1, XMMWORD PTR [esi+16] ; line 1
856	movdqa xmm2, XMMWORD PTR [esi+32] ; line 2
857	movdqa xmm3, XMMWORD PTR [esi+48] ; line 3
858
859	paddw xmm0, xmm7
860	paddw xmm1, xmm7
861	paddw xmm2, xmm7
862	paddw xmm3, xmm7
863
864	packuswb xmm0, xmm1
865	packuswb xmm2, xmm3
866
867	movq QWORD PTR [edi], xmm0 ;0*DstStep
868	movq QWORD PTR [edi+ebx], xmm2 ;2*DstStep
869
870	psrldq xmm0, 8
871	psrldq xmm2, 8
872
873	movq QWORD PTR [edi+edx], xmm0 ;1*DstStep
874	movq QWORD PTR [edi+eax], xmm2 ;3*DstStep
875
876	add edi, ebx
877	add esi, 64
878	add edi, ebx
879	dec ecx
880	jnz again
881
882	pop edi
883	pop esi
884	pop edx
885	pop ecx
886	pop ebx
887	pop eax
888	}
889	}
890
891
892	/*
893	* Name:
894	* ippiDCTQuantInv8x8LS_JPEG_16s8u_C1R
895	*
896	* Purpose:
897	* Inverse DCT transform, de-quantization and level shift
898	*
899	* Parameters:
900	* pSrc - pointer to source
901	* pDst - pointer to output array
902	* DstStep - line offset for output data
903	* pEncoderQuantTable - pointer to Quantization table
904	*
905	*/
906
907	GLOBAL(void)
908	ippiDCTQuantInv8x8LS_JPEG_16s8u_C1R(
909	short * pSrc,
910	unsigned char * pDst,
911	int DstStep,
912	const unsigned short * pQuantInvTable)
913	{
914
915	__declspec(align(16)) Ipp8u buf[DCTSIZE2*sizeof(Ipp16s)];
916	Ipp16s * workbuf = (Ipp16s *)buf;
917
918	ownpj_QuantInv_8x8_16s(pSrc,workbuf,pQuantInvTable);
919	dct_8x8_inv_16s(workbuf,workbuf);
920	ownpj_Add128_8x8_16s8u(workbuf,pDst,DstStep);
921
922	}
923
924	GLOBAL(void)
925	jpeg_idct_islow_sse2 (
926	j_decompress_ptr cinfo,
927	jpeg_component_info * compptr,
928	JCOEFPTR coef_block,
929	JSAMPARRAY output_buf,
930	JDIMENSION output_col)
931	{
932	int ctr;
933	JCOEFPTR inptr;
934	Ipp16u* quantptr;
935	Ipp8u* wsptr;
936	__declspec(align(16)) Ipp8u workspace[DCTSIZE2];
937	JSAMPROW outptr;
938
939	inptr = coef_block;
940	quantptr = (Ipp16u*)compptr->dct_table;
941	wsptr = workspace;
942
943	ippiDCTQuantInv8x8LS_JPEG_16s8u_C1R(inptr, workspace, 8, quantptr);
944
945	for(ctr = 0; ctr < DCTSIZE; ctr++)
946	{
947	outptr = output_buf[ctr] + output_col;
948
949	outptr[0] = wsptr[0];
950	outptr[1] = wsptr[1];
951	outptr[2] = wsptr[2];
952	outptr[3] = wsptr[3];
953	outptr[4] = wsptr[4];
954	outptr[5] = wsptr[5];
955	outptr[6] = wsptr[6];
956	outptr[7] = wsptr[7];
957
958	wsptr += DCTSIZE;
959	}
960	}
961	#endif /* HAVE_SSE2_INTEL_MNEMONICS */
962
963	#endif /* DCT_ISLOW_SUPPORTED */

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: trunk/third/firefox/jpeg/jidctint.c @ 21695

Download in other formats: