ffmpeg / libavcodec / jfdctint.c @ 849f1035
History  View  Annotate  Download (15.8 KB)
1 
/*


2 
* jfdctint.c

3 
*

4 
* This file is part of the Independent JPEG Group's software.

5 
*

6 
* The authors make NO WARRANTY or representation, either express or implied,

7 
* with respect to this software, its quality, accuracy, merchantability, or

8 
* fitness for a particular purpose. This software is provided "AS IS", and

9 
* you, its user, assume the entire risk as to its quality and accuracy.

10 
*

11 
* This software is copyright (C) 19911996, Thomas G. Lane.

12 
* All Rights Reserved except as specified below.

13 
*

14 
* Permission is hereby granted to use, copy, modify, and distribute this

15 
* software (or portions thereof) for any purpose, without fee, subject to

16 
* these conditions:

17 
* (1) If any part of the source code for this software is distributed, then

18 
* this README file must be included, with this copyright and nowarranty

19 
* notice unaltered; and any additions, deletions, or changes to the original

20 
* files must be clearly indicated in accompanying documentation.

21 
* (2) If only executable code is distributed, then the accompanying

22 
* documentation must state that "this software is based in part on the work

23 
* of the Independent JPEG Group".

24 
* (3) Permission for use of this software is granted only if the user accepts

25 
* full responsibility for any undesirable consequences; the authors accept

26 
* NO LIABILITY for damages of any kind.

27 
*

28 
* These conditions apply to any software derived from or based on the IJG

29 
* code, not just to the unmodified library. If you use our work, you ought

30 
* to acknowledge us.

31 
*

32 
* Permission is NOT granted for the use of any IJG author's name or company

33 
* name in advertising or publicity relating to this software or products

34 
* derived from it. This software may be referred to only as "the Independent

35 
* JPEG Group's software".

36 
*

37 
* We specifically permit and encourage the use of this software as the basis

38 
* of commercial products, provided that all warranty or liability claims are

39 
* assumed by the product vendor.

40 
*

41 
* This file contains a slowbutaccurate integer implementation of the

42 
* forward DCT (Discrete Cosine Transform).

43 
*

44 
* A 2D DCT can be done by 1D DCT on each row followed by 1D DCT

45 
* on each column. Direct algorithms are also available, but they are

46 
* much more complex and seem not to be any faster when reduced to code.

47 
*

48 
* This implementation is based on an algorithm described in

49 
* C. Loeffler, A. Ligtenberg and G. Moschytz, "Practical Fast 1D DCT

50 
* Algorithms with 11 Multiplications", Proc. Int'l. Conf. on Acoustics,

51 
* Speech, and Signal Processing 1989 (ICASSP '89), pp. 988991.

52 
* The primary algorithm described there uses 11 multiplies and 29 adds.

53 
* We use their alternate method with 12 multiplies and 32 adds.

54 
* The advantage of this method is that no data path contains more than one

55 
* multiplication; this allows a very simple and accurate implementation in

56 
* scaled fixedpoint arithmetic, with a minimal number of shifts.

57 
*/

58  
59 
/**

60 
* @file jfdctint.c

61 
* Independent JPEG Group's slow & accurate dct.

62 
*/

63  
64 
#include <stdlib.h> 
65 
#include <stdio.h> 
66 
#include "common.h" 
67 
#include "dsputil.h" 
68  
69 
#define SHIFT_TEMPS

70 
#define DCTSIZE 8 
71 
#define BITS_IN_JSAMPLE 8 
72 
#define GLOBAL(x) x

73 
#define RIGHT_SHIFT(x, n) ((x) >> (n))

74 
#define MULTIPLY16C16(var,const) ((var)*(const)) 
75  
76 
#if 1 //def USE_ACCURATE_ROUNDING 
77 
#define DESCALE(x,n) RIGHT_SHIFT((x) + (1 << ((n)  1)), n) 
78 
#else

79 
#define DESCALE(x,n) RIGHT_SHIFT(x, n)

80 
#endif

81  
82  
83 
/*

84 
* This module is specialized to the case DCTSIZE = 8.

85 
*/

86  
87 
#if DCTSIZE != 8 
88 
Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */ 
89 
#endif

90  
91  
92 
/*

93 
* The poop on this scaling stuff is as follows:

94 
*

95 
* Each 1D DCT step produces outputs which are a factor of sqrt(N)

96 
* larger than the true DCT outputs. The final outputs are therefore

97 
* a factor of N larger than desired; since N=8 this can be cured by

98 
* a simple right shift at the end of the algorithm. The advantage of

99 
* this arrangement is that we save two multiplications per 1D DCT,

100 
* because the y0 and y4 outputs need not be divided by sqrt(N).

101 
* In the IJG code, this factor of 8 is removed by the quantization step

102 
* (in jcdctmgr.c), NOT in this module.

103 
*

104 
* We have to do addition and subtraction of the integer inputs, which

105 
* is no problem, and multiplication by fractional constants, which is

106 
* a problem to do in integer arithmetic. We multiply all the constants

107 
* by CONST_SCALE and convert them to integer constants (thus retaining

108 
* CONST_BITS bits of precision in the constants). After doing a

109 
* multiplication we have to divide the product by CONST_SCALE, with proper

110 
* rounding, to produce the correct output. This division can be done

111 
* cheaply as a right shift of CONST_BITS bits. We postpone shifting

112 
* as long as possible so that partial sums can be added together with

113 
* full fractional precision.

114 
*

115 
* The outputs of the first pass are scaled up by PASS1_BITS bits so that

116 
* they are represented to betterthanintegral precision. These outputs

117 
* require BITS_IN_JSAMPLE + PASS1_BITS + 3 bits; this fits in a 16bit word

118 
* with the recommended scaling. (For 12bit sample data, the intermediate

119 
* array is int32_t anyway.)

120 
*

121 
* To avoid overflow of the 32bit intermediate results in pass 2, we must

122 
* have BITS_IN_JSAMPLE + CONST_BITS + PASS1_BITS <= 26. Error analysis

123 
* shows that the values given below are the most effective.

124 
*/

125  
126 
#if BITS_IN_JSAMPLE == 8 
127 
#define CONST_BITS 13 
128 
#define PASS1_BITS 4 /* set this to 2 if 16x16 multiplies are faster */ 
129 
#else

130 
#define CONST_BITS 13 
131 
#define PASS1_BITS 1 /* lose a little precision to avoid overflow */ 
132 
#endif

133  
134 
/* Some C compilers fail to reduce "FIX(constant)" at compile time, thus

135 
* causing a lot of useless floatingpoint operations at run time.

136 
* To get around this we use the following precalculated constants.

137 
* If you change CONST_BITS you may want to add appropriate values.

138 
* (With a reasonable C compiler, you can just rely on the FIX() macro...)

139 
*/

140  
141 
#if CONST_BITS == 13 
142 
#define FIX_0_298631336 ((int32_t) 2446) /* FIX(0.298631336) */ 
143 
#define FIX_0_390180644 ((int32_t) 3196) /* FIX(0.390180644) */ 
144 
#define FIX_0_541196100 ((int32_t) 4433) /* FIX(0.541196100) */ 
145 
#define FIX_0_765366865 ((int32_t) 6270) /* FIX(0.765366865) */ 
146 
#define FIX_0_899976223 ((int32_t) 7373) /* FIX(0.899976223) */ 
147 
#define FIX_1_175875602 ((int32_t) 9633) /* FIX(1.175875602) */ 
148 
#define FIX_1_501321110 ((int32_t) 12299) /* FIX(1.501321110) */ 
149 
#define FIX_1_847759065 ((int32_t) 15137) /* FIX(1.847759065) */ 
150 
#define FIX_1_961570560 ((int32_t) 16069) /* FIX(1.961570560) */ 
151 
#define FIX_2_053119869 ((int32_t) 16819) /* FIX(2.053119869) */ 
152 
#define FIX_2_562915447 ((int32_t) 20995) /* FIX(2.562915447) */ 
153 
#define FIX_3_072711026 ((int32_t) 25172) /* FIX(3.072711026) */ 
154 
#else

155 
#define FIX_0_298631336 FIX(0.298631336) 
156 
#define FIX_0_390180644 FIX(0.390180644) 
157 
#define FIX_0_541196100 FIX(0.541196100) 
158 
#define FIX_0_765366865 FIX(0.765366865) 
159 
#define FIX_0_899976223 FIX(0.899976223) 
160 
#define FIX_1_175875602 FIX(1.175875602) 
161 
#define FIX_1_501321110 FIX(1.501321110) 
162 
#define FIX_1_847759065 FIX(1.847759065) 
163 
#define FIX_1_961570560 FIX(1.961570560) 
164 
#define FIX_2_053119869 FIX(2.053119869) 
165 
#define FIX_2_562915447 FIX(2.562915447) 
166 
#define FIX_3_072711026 FIX(3.072711026) 
167 
#endif

168  
169  
170 
/* Multiply an int32_t variable by an int32_t constant to yield an int32_t result.

171 
* For 8bit samples with the recommended scaling, all the variable

172 
* and constant values involved are no more than 16 bits wide, so a

173 
* 16x16>32 bit multiply can be used instead of a full 32x32 multiply.

174 
* For 12bit samples, a full 32bit multiplication will be needed.

175 
*/

176  
177 
#if BITS_IN_JSAMPLE == 8 && CONST_BITS<=13 && PASS1_BITS<=2 
178 
#define MULTIPLY(var,const) MULTIPLY16C16(var,const) 
179 
#else

180 
#define MULTIPLY(var,const) ((var) * (const)) 
181 
#endif

182  
183  
184 
static av_always_inline void row_fdct(DCTELEM * data){ 
185 
int_fast32_t tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7; 
186 
int_fast32_t tmp10, tmp11, tmp12, tmp13; 
187 
int_fast32_t z1, z2, z3, z4, z5; 
188 
DCTELEM *dataptr; 
189 
int ctr;

190 
SHIFT_TEMPS 
191  
192 
/* Pass 1: process rows. */

193 
/* Note results are scaled up by sqrt(8) compared to a true DCT; */

194 
/* furthermore, we scale the results by 2**PASS1_BITS. */

195  
196 
dataptr = data; 
197 
for (ctr = DCTSIZE1; ctr >= 0; ctr) { 
198 
tmp0 = dataptr[0] + dataptr[7]; 
199 
tmp7 = dataptr[0]  dataptr[7]; 
200 
tmp1 = dataptr[1] + dataptr[6]; 
201 
tmp6 = dataptr[1]  dataptr[6]; 
202 
tmp2 = dataptr[2] + dataptr[5]; 
203 
tmp5 = dataptr[2]  dataptr[5]; 
204 
tmp3 = dataptr[3] + dataptr[4]; 
205 
tmp4 = dataptr[3]  dataptr[4]; 
206  
207 
/* Even part per LL&M figure 1  note that published figure is faulty;

208 
* rotator "sqrt(2)*c1" should be "sqrt(2)*c6".

209 
*/

210  
211 
tmp10 = tmp0 + tmp3; 
212 
tmp13 = tmp0  tmp3; 
213 
tmp11 = tmp1 + tmp2; 
214 
tmp12 = tmp1  tmp2; 
215  
216 
dataptr[0] = (DCTELEM) ((tmp10 + tmp11) << PASS1_BITS);

217 
dataptr[4] = (DCTELEM) ((tmp10  tmp11) << PASS1_BITS);

218  
219 
z1 = MULTIPLY(tmp12 + tmp13, FIX_0_541196100); 
220 
dataptr[2] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp13, FIX_0_765366865),

221 
CONST_BITSPASS1_BITS); 
222 
dataptr[6] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp12,  FIX_1_847759065),

223 
CONST_BITSPASS1_BITS); 
224  
225 
/* Odd part per figure 8  note paper omits factor of sqrt(2).

226 
* cK represents cos(K*pi/16).

227 
* i0..i3 in the paper are tmp4..tmp7 here.

228 
*/

229  
230 
z1 = tmp4 + tmp7; 
231 
z2 = tmp5 + tmp6; 
232 
z3 = tmp4 + tmp6; 
233 
z4 = tmp5 + tmp7; 
234 
z5 = MULTIPLY(z3 + z4, FIX_1_175875602); /* sqrt(2) * c3 */

235  
236 
tmp4 = MULTIPLY(tmp4, FIX_0_298631336); /* sqrt(2) * (c1+c3+c5c7) */

237 
tmp5 = MULTIPLY(tmp5, FIX_2_053119869); /* sqrt(2) * ( c1+c3c5+c7) */

238 
tmp6 = MULTIPLY(tmp6, FIX_3_072711026); /* sqrt(2) * ( c1+c3+c5c7) */

239 
tmp7 = MULTIPLY(tmp7, FIX_1_501321110); /* sqrt(2) * ( c1+c3c5c7) */

240 
z1 = MULTIPLY(z1,  FIX_0_899976223); /* sqrt(2) * (c7c3) */

241 
z2 = MULTIPLY(z2,  FIX_2_562915447); /* sqrt(2) * (c1c3) */

242 
z3 = MULTIPLY(z3,  FIX_1_961570560); /* sqrt(2) * (c3c5) */

243 
z4 = MULTIPLY(z4,  FIX_0_390180644); /* sqrt(2) * (c5c3) */

244  
245 
z3 += z5; 
246 
z4 += z5; 
247  
248 
dataptr[7] = (DCTELEM) DESCALE(tmp4 + z1 + z3, CONST_BITSPASS1_BITS);

249 
dataptr[5] = (DCTELEM) DESCALE(tmp5 + z2 + z4, CONST_BITSPASS1_BITS);

250 
dataptr[3] = (DCTELEM) DESCALE(tmp6 + z2 + z3, CONST_BITSPASS1_BITS);

251 
dataptr[1] = (DCTELEM) DESCALE(tmp7 + z1 + z4, CONST_BITSPASS1_BITS);

252  
253 
dataptr += DCTSIZE; /* advance pointer to next row */

254 
} 
255 
} 
256  
257 
/*

258 
* Perform the forward DCT on one block of samples.

259 
*/

260  
261 
GLOBAL(void)

262 
ff_jpeg_fdct_islow (DCTELEM * data) 
263 
{ 
264 
int_fast32_t tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7; 
265 
int_fast32_t tmp10, tmp11, tmp12, tmp13; 
266 
int_fast32_t z1, z2, z3, z4, z5; 
267 
DCTELEM *dataptr; 
268 
int ctr;

269 
SHIFT_TEMPS 
270  
271 
row_fdct(data); 
272  
273 
/* Pass 2: process columns.

274 
* We remove the PASS1_BITS scaling, but leave the results scaled up

275 
* by an overall factor of 8.

276 
*/

277  
278 
dataptr = data; 
279 
for (ctr = DCTSIZE1; ctr >= 0; ctr) { 
280 
tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*7]; 
281 
tmp7 = dataptr[DCTSIZE*0]  dataptr[DCTSIZE*7]; 
282 
tmp1 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*6]; 
283 
tmp6 = dataptr[DCTSIZE*1]  dataptr[DCTSIZE*6]; 
284 
tmp2 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*5]; 
285 
tmp5 = dataptr[DCTSIZE*2]  dataptr[DCTSIZE*5]; 
286 
tmp3 = dataptr[DCTSIZE*3] + dataptr[DCTSIZE*4]; 
287 
tmp4 = dataptr[DCTSIZE*3]  dataptr[DCTSIZE*4]; 
288  
289 
/* Even part per LL&M figure 1  note that published figure is faulty;

290 
* rotator "sqrt(2)*c1" should be "sqrt(2)*c6".

291 
*/

292  
293 
tmp10 = tmp0 + tmp3; 
294 
tmp13 = tmp0  tmp3; 
295 
tmp11 = tmp1 + tmp2; 
296 
tmp12 = tmp1  tmp2; 
297  
298 
dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp11, PASS1_BITS);

299 
dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp10  tmp11, PASS1_BITS);

300  
301 
z1 = MULTIPLY(tmp12 + tmp13, FIX_0_541196100); 
302 
dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp13, FIX_0_765366865),

303 
CONST_BITS+PASS1_BITS); 
304 
dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp12,  FIX_1_847759065),

305 
CONST_BITS+PASS1_BITS); 
306  
307 
/* Odd part per figure 8  note paper omits factor of sqrt(2).

308 
* cK represents cos(K*pi/16).

309 
* i0..i3 in the paper are tmp4..tmp7 here.

310 
*/

311  
312 
z1 = tmp4 + tmp7; 
313 
z2 = tmp5 + tmp6; 
314 
z3 = tmp4 + tmp6; 
315 
z4 = tmp5 + tmp7; 
316 
z5 = MULTIPLY(z3 + z4, FIX_1_175875602); /* sqrt(2) * c3 */

317  
318 
tmp4 = MULTIPLY(tmp4, FIX_0_298631336); /* sqrt(2) * (c1+c3+c5c7) */

319 
tmp5 = MULTIPLY(tmp5, FIX_2_053119869); /* sqrt(2) * ( c1+c3c5+c7) */

320 
tmp6 = MULTIPLY(tmp6, FIX_3_072711026); /* sqrt(2) * ( c1+c3+c5c7) */

321 
tmp7 = MULTIPLY(tmp7, FIX_1_501321110); /* sqrt(2) * ( c1+c3c5c7) */

322 
z1 = MULTIPLY(z1,  FIX_0_899976223); /* sqrt(2) * (c7c3) */

323 
z2 = MULTIPLY(z2,  FIX_2_562915447); /* sqrt(2) * (c1c3) */

324 
z3 = MULTIPLY(z3,  FIX_1_961570560); /* sqrt(2) * (c3c5) */

325 
z4 = MULTIPLY(z4,  FIX_0_390180644); /* sqrt(2) * (c5c3) */

326  
327 
z3 += z5; 
328 
z4 += z5; 
329  
330 
dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp4 + z1 + z3,

331 
CONST_BITS+PASS1_BITS); 
332 
dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp5 + z2 + z4,

333 
CONST_BITS+PASS1_BITS); 
334 
dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp6 + z2 + z3,

335 
CONST_BITS+PASS1_BITS); 
336 
dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp7 + z1 + z4,

337 
CONST_BITS+PASS1_BITS); 
338  
339 
dataptr++; /* advance pointer to next column */

340 
} 
341 
} 
342  
343 
/*

344 
* The secret of DCT248 is really simple  you do the usual 1DCT

345 
* on the rows and then, instead of doing even and odd, part on the colums

346 
* you do even part two times.

347 
*/

348 
GLOBAL(void)

349 
ff_fdct248_islow (DCTELEM * data) 
350 
{ 
351 
int_fast32_t tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7; 
352 
int_fast32_t tmp10, tmp11, tmp12, tmp13; 
353 
int_fast32_t z1; 
354 
DCTELEM *dataptr; 
355 
int ctr;

356 
SHIFT_TEMPS 
357  
358 
row_fdct(data); 
359  
360 
/* Pass 2: process columns.

361 
* We remove the PASS1_BITS scaling, but leave the results scaled up

362 
* by an overall factor of 8.

363 
*/

364  
365 
dataptr = data; 
366 
for (ctr = DCTSIZE1; ctr >= 0; ctr) { 
367 
tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*1]; 
368 
tmp1 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*3]; 
369 
tmp2 = dataptr[DCTSIZE*4] + dataptr[DCTSIZE*5]; 
370 
tmp3 = dataptr[DCTSIZE*6] + dataptr[DCTSIZE*7]; 
371 
tmp4 = dataptr[DCTSIZE*0]  dataptr[DCTSIZE*1]; 
372 
tmp5 = dataptr[DCTSIZE*2]  dataptr[DCTSIZE*3]; 
373 
tmp6 = dataptr[DCTSIZE*4]  dataptr[DCTSIZE*5]; 
374 
tmp7 = dataptr[DCTSIZE*6]  dataptr[DCTSIZE*7]; 
375  
376 
tmp10 = tmp0 + tmp3; 
377 
tmp11 = tmp1 + tmp2; 
378 
tmp12 = tmp1  tmp2; 
379 
tmp13 = tmp0  tmp3; 
380  
381 
dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp11, PASS1_BITS);

382 
dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp10  tmp11, PASS1_BITS);

383  
384 
z1 = MULTIPLY(tmp12 + tmp13, FIX_0_541196100); 
385 
dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp13, FIX_0_765366865),

386 
CONST_BITS+PASS1_BITS); 
387 
dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp12,  FIX_1_847759065),

388 
CONST_BITS+PASS1_BITS); 
389  
390 
tmp10 = tmp4 + tmp7; 
391 
tmp11 = tmp5 + tmp6; 
392 
tmp12 = tmp5  tmp6; 
393 
tmp13 = tmp4  tmp7; 
394  
395 
dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp10 + tmp11, PASS1_BITS);

396 
dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp10  tmp11, PASS1_BITS);

397  
398 
z1 = MULTIPLY(tmp12 + tmp13, FIX_0_541196100); 
399 
dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp13, FIX_0_765366865),

400 
CONST_BITS+PASS1_BITS); 
401 
dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp12,  FIX_1_847759065),

402 
CONST_BITS+PASS1_BITS); 
403  
404 
dataptr++; /* advance pointer to next column */

405 
} 
406 
} 