## ffmpeg / libavcodec / jfdctint.c @ 0469baf1

History | View | Annotate | Download (11.2 KB)

1 | 28db7fce | Michael Niedermayer | ```
/*
``` |
---|---|---|---|

2 | ```
* jfdctint.c
``` |
||

3 | ```
*
``` |
||

4 | ```
* Copyright (C) 1991-1996, 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 | ```
* forward DCT (Discrete Cosine Transform).
``` |
||

10 | ```
*
``` |
||

11 | ```
* A 2-D DCT can be done by 1-D DCT on each row followed by 1-D DCT
``` |
||

12 | ```
* on each column. Direct algorithms are also available, but they are
``` |
||

13 | ```
* much more complex and seem not to be any faster when reduced to code.
``` |
||

14 | ```
*
``` |
||

15 | ```
* This implementation is based on an algorithm described in
``` |
||

16 | ```
* C. Loeffler, A. Ligtenberg and G. Moschytz, "Practical Fast 1-D DCT
``` |
||

17 | ```
* Algorithms with 11 Multiplications", Proc. Int'l. Conf. on Acoustics,
``` |
||

18 | ```
* Speech, and Signal Processing 1989 (ICASSP '89), pp. 988-991.
``` |
||

19 | ```
* The primary algorithm described there uses 11 multiplies and 29 adds.
``` |
||

20 | ```
* We use their alternate method with 12 multiplies and 32 adds.
``` |
||

21 | ```
* The advantage of this method is that no data path contains more than one
``` |
||

22 | ```
* multiplication; this allows a very simple and accurate implementation in
``` |
||

23 | ```
* scaled fixed-point arithmetic, with a minimal number of shifts.
``` |
||

24 | ```
*/
``` |
||

25 | |||

26 | 983e3246 | Michael Niedermayer | ```
/**
``` |

27 | ```
* @file jfdctint.c
``` |
||

28 | ```
* Independent JPEG Group's slow & accurate dct.
``` |
||

29 | ```
*/
``` |
||

30 | |||

31 | 28db7fce | Michael Niedermayer | #include <stdlib.h> |

32 | #include <stdio.h> |
||

33 | #include "common.h" |
||

34 | #include "dsputil.h" |
||

35 | |||

36 | ```
#define SHIFT_TEMPS
``` |
||

37 | #define DCTSIZE 8 |
||

38 | 004c18ee | Michael Niedermayer | #define BITS_IN_JSAMPLE 8 |

39 | 28db7fce | Michael Niedermayer | ```
#define GLOBAL(x) x
``` |

40 | ```
#define RIGHT_SHIFT(x, n) ((x) >> (n))
``` |
||

41 | 004c18ee | Michael Niedermayer | #define MULTIPLY16C16(var,const) ((var)*(const)) |

42 | 28db7fce | Michael Niedermayer | |

43 | #if 1 //def USE_ACCURATE_ROUNDING |
||

44 | #define DESCALE(x,n) RIGHT_SHIFT((x) + (1 << ((n) - 1)), n) |
||

45 | ```
#else
``` |
||

46 | ```
#define DESCALE(x,n) RIGHT_SHIFT(x, n)
``` |
||

47 | ```
#endif
``` |
||

48 | |||

49 | |||

50 | ```
/*
``` |
||

51 | ```
* This module is specialized to the case DCTSIZE = 8.
``` |
||

52 | ```
*/
``` |
||

53 | |||

54 | #if DCTSIZE != 8 |
||

55 | Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */ |
||

56 | ```
#endif
``` |
||

57 | |||

58 | |||

59 | ```
/*
``` |
||

60 | ```
* The poop on this scaling stuff is as follows:
``` |
||

61 | ```
*
``` |
||

62 | ```
* Each 1-D DCT step produces outputs which are a factor of sqrt(N)
``` |
||

63 | ```
* larger than the true DCT outputs. The final outputs are therefore
``` |
||

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

65 | ```
* a simple right shift at the end of the algorithm. The advantage of
``` |
||

66 | ```
* this arrangement is that we save two multiplications per 1-D DCT,
``` |
||

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

68 | ```
* In the IJG code, this factor of 8 is removed by the quantization step
``` |
||

69 | ```
* (in jcdctmgr.c), NOT in this module.
``` |
||

70 | ```
*
``` |
||

71 | ```
* We have to do addition and subtraction of the integer inputs, which
``` |
||

72 | ```
* is no problem, and multiplication by fractional constants, which is
``` |
||

73 | ```
* a problem to do in integer arithmetic. We multiply all the constants
``` |
||

74 | ```
* by CONST_SCALE and convert them to integer constants (thus retaining
``` |
||

75 | ```
* CONST_BITS bits of precision in the constants). After doing a
``` |
||

76 | ```
* multiplication we have to divide the product by CONST_SCALE, with proper
``` |
||

77 | ```
* rounding, to produce the correct output. This division can be done
``` |
||

78 | ```
* cheaply as a right shift of CONST_BITS bits. We postpone shifting
``` |
||

79 | ```
* as long as possible so that partial sums can be added together with
``` |
||

80 | ```
* full fractional precision.
``` |
||

81 | ```
*
``` |
||

82 | ```
* The outputs of the first pass are scaled up by PASS1_BITS bits so that
``` |
||

83 | ```
* they are represented to better-than-integral precision. These outputs
``` |
||

84 | ```
* require BITS_IN_JSAMPLE + PASS1_BITS + 3 bits; this fits in a 16-bit word
``` |
||

85 | ```
* with the recommended scaling. (For 12-bit sample data, the intermediate
``` |
||

86 | 0c1a9eda | Zdenek Kabelac | ```
* array is int32_t anyway.)
``` |

87 | 28db7fce | Michael Niedermayer | ```
*
``` |

88 | ```
* To avoid overflow of the 32-bit intermediate results in pass 2, we must
``` |
||

89 | ```
* have BITS_IN_JSAMPLE + CONST_BITS + PASS1_BITS <= 26. Error analysis
``` |
||

90 | ```
* shows that the values given below are the most effective.
``` |
||

91 | ```
*/
``` |
||

92 | |||

93 | #if BITS_IN_JSAMPLE == 8 |
||

94 | #define CONST_BITS 13 |
||

95 | 004c18ee | Michael Niedermayer | #define PASS1_BITS 4 /* set this to 2 if 16x16 multiplies are faster */ |

96 | 28db7fce | Michael Niedermayer | ```
#else
``` |

97 | #define CONST_BITS 13 |
||

98 | #define PASS1_BITS 1 /* lose a little precision to avoid overflow */ |
||

99 | ```
#endif
``` |
||

100 | |||

101 | ```
/* Some C compilers fail to reduce "FIX(constant)" at compile time, thus
``` |
||

102 | ```
* causing a lot of useless floating-point operations at run time.
``` |
||

103 | ```
* To get around this we use the following pre-calculated constants.
``` |
||

104 | ```
* If you change CONST_BITS you may want to add appropriate values.
``` |
||

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

106 | ```
*/
``` |
||

107 | |||

108 | #if CONST_BITS == 13 |
||

109 | 0c1a9eda | Zdenek Kabelac | #define FIX_0_298631336 ((int32_t) 2446) /* FIX(0.298631336) */ |

110 | #define FIX_0_390180644 ((int32_t) 3196) /* FIX(0.390180644) */ |
||

111 | #define FIX_0_541196100 ((int32_t) 4433) /* FIX(0.541196100) */ |
||

112 | #define FIX_0_765366865 ((int32_t) 6270) /* FIX(0.765366865) */ |
||

113 | #define FIX_0_899976223 ((int32_t) 7373) /* FIX(0.899976223) */ |
||

114 | #define FIX_1_175875602 ((int32_t) 9633) /* FIX(1.175875602) */ |
||

115 | #define FIX_1_501321110 ((int32_t) 12299) /* FIX(1.501321110) */ |
||

116 | #define FIX_1_847759065 ((int32_t) 15137) /* FIX(1.847759065) */ |
||

117 | #define FIX_1_961570560 ((int32_t) 16069) /* FIX(1.961570560) */ |
||

118 | #define FIX_2_053119869 ((int32_t) 16819) /* FIX(2.053119869) */ |
||

119 | #define FIX_2_562915447 ((int32_t) 20995) /* FIX(2.562915447) */ |
||

120 | #define FIX_3_072711026 ((int32_t) 25172) /* FIX(3.072711026) */ |
||

121 | 28db7fce | Michael Niedermayer | ```
#else
``` |

122 | #define FIX_0_298631336 FIX(0.298631336) |
||

123 | #define FIX_0_390180644 FIX(0.390180644) |
||

124 | #define FIX_0_541196100 FIX(0.541196100) |
||

125 | #define FIX_0_765366865 FIX(0.765366865) |
||

126 | #define FIX_0_899976223 FIX(0.899976223) |
||

127 | #define FIX_1_175875602 FIX(1.175875602) |
||

128 | #define FIX_1_501321110 FIX(1.501321110) |
||

129 | #define FIX_1_847759065 FIX(1.847759065) |
||

130 | #define FIX_1_961570560 FIX(1.961570560) |
||

131 | #define FIX_2_053119869 FIX(2.053119869) |
||

132 | #define FIX_2_562915447 FIX(2.562915447) |
||

133 | #define FIX_3_072711026 FIX(3.072711026) |
||

134 | ```
#endif
``` |
||

135 | |||

136 | |||

137 | 0c1a9eda | Zdenek Kabelac | ```
/* Multiply an int32_t variable by an int32_t constant to yield an int32_t result.
``` |

138 | 28db7fce | Michael Niedermayer | ```
* For 8-bit samples with the recommended scaling, all the variable
``` |

139 | ```
* and constant values involved are no more than 16 bits wide, so a
``` |
||

140 | ```
* 16x16->32 bit multiply can be used instead of a full 32x32 multiply.
``` |
||

141 | ```
* For 12-bit samples, a full 32-bit multiplication will be needed.
``` |
||

142 | ```
*/
``` |
||

143 | |||

144 | 004c18ee | Michael Niedermayer | #if BITS_IN_JSAMPLE == 8 && CONST_BITS<=13 && PASS1_BITS<=2 |

145 | 28db7fce | Michael Niedermayer | #define MULTIPLY(var,const) MULTIPLY16C16(var,const) |

146 | ```
#else
``` |
||

147 | #define MULTIPLY(var,const) ((var) * (const)) |
||

148 | ```
#endif
``` |
||

149 | |||

150 | |||

151 | ```
/*
``` |
||

152 | ```
* Perform the forward DCT on one block of samples.
``` |
||

153 | ```
*/
``` |
||

154 | |||

155 | ```
GLOBAL(void)
``` |
||

156 | ff_jpeg_fdct_islow (DCTELEM * data) |
||

157 | { |
||

158 | 0c1a9eda | Zdenek Kabelac | int32_t tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7; |

159 | int32_t tmp10, tmp11, tmp12, tmp13; |
||

160 | int32_t z1, z2, z3, z4, z5; |
||

161 | 28db7fce | Michael Niedermayer | DCTELEM *dataptr; |

162 | ```
int ctr;
``` |
||

163 | SHIFT_TEMPS |
||

164 | |||

165 | ```
/* Pass 1: process rows. */
``` |
||

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

167 | ```
/* furthermore, we scale the results by 2**PASS1_BITS. */
``` |
||

168 | |||

169 | dataptr = data; |
||

170 | for (ctr = DCTSIZE-1; ctr >= 0; ctr--) { |
||

171 | tmp0 = dataptr[0] + dataptr[7]; |
||

172 | tmp7 = dataptr[0] - dataptr[7]; |
||

173 | tmp1 = dataptr[1] + dataptr[6]; |
||

174 | tmp6 = dataptr[1] - dataptr[6]; |
||

175 | tmp2 = dataptr[2] + dataptr[5]; |
||

176 | tmp5 = dataptr[2] - dataptr[5]; |
||

177 | tmp3 = dataptr[3] + dataptr[4]; |
||

178 | tmp4 = dataptr[3] - dataptr[4]; |
||

179 | |||

180 | ```
/* Even part per LL&M figure 1 --- note that published figure is faulty;
``` |
||

181 | ```
* rotator "sqrt(2)*c1" should be "sqrt(2)*c6".
``` |
||

182 | ```
*/
``` |
||

183 | |||

184 | tmp10 = tmp0 + tmp3; |
||

185 | tmp13 = tmp0 - tmp3; |
||

186 | tmp11 = tmp1 + tmp2; |
||

187 | tmp12 = tmp1 - tmp2; |
||

188 | |||

189 | ```
dataptr[0] = (DCTELEM) ((tmp10 + tmp11) << PASS1_BITS);
``` |
||

190 | ```
dataptr[4] = (DCTELEM) ((tmp10 - tmp11) << PASS1_BITS);
``` |
||

191 | |||

192 | z1 = MULTIPLY(tmp12 + tmp13, FIX_0_541196100); |
||

193 | ```
dataptr[2] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp13, FIX_0_765366865),
``` |
||

194 | CONST_BITS-PASS1_BITS); |
||

195 | ```
dataptr[6] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp12, - FIX_1_847759065),
``` |
||

196 | CONST_BITS-PASS1_BITS); |
||

197 | |||

198 | ```
/* Odd part per figure 8 --- note paper omits factor of sqrt(2).
``` |
||

199 | ```
* cK represents cos(K*pi/16).
``` |
||

200 | ```
* i0..i3 in the paper are tmp4..tmp7 here.
``` |
||

201 | ```
*/
``` |
||

202 | |||

203 | z1 = tmp4 + tmp7; |
||

204 | z2 = tmp5 + tmp6; |
||

205 | z3 = tmp4 + tmp6; |
||

206 | z4 = tmp5 + tmp7; |
||

207 | ```
z5 = MULTIPLY(z3 + z4, FIX_1_175875602); /* sqrt(2) * c3 */
``` |
||

208 | |||

209 | ```
tmp4 = MULTIPLY(tmp4, FIX_0_298631336); /* sqrt(2) * (-c1+c3+c5-c7) */
``` |
||

210 | ```
tmp5 = MULTIPLY(tmp5, FIX_2_053119869); /* sqrt(2) * ( c1+c3-c5+c7) */
``` |
||

211 | ```
tmp6 = MULTIPLY(tmp6, FIX_3_072711026); /* sqrt(2) * ( c1+c3+c5-c7) */
``` |
||

212 | ```
tmp7 = MULTIPLY(tmp7, FIX_1_501321110); /* sqrt(2) * ( c1+c3-c5-c7) */
``` |
||

213 | ```
z1 = MULTIPLY(z1, - FIX_0_899976223); /* sqrt(2) * (c7-c3) */
``` |
||

214 | ```
z2 = MULTIPLY(z2, - FIX_2_562915447); /* sqrt(2) * (-c1-c3) */
``` |
||

215 | ```
z3 = MULTIPLY(z3, - FIX_1_961570560); /* sqrt(2) * (-c3-c5) */
``` |
||

216 | ```
z4 = MULTIPLY(z4, - FIX_0_390180644); /* sqrt(2) * (c5-c3) */
``` |
||

217 | |||

218 | z3 += z5; |
||

219 | z4 += z5; |
||

220 | |||

221 | ```
dataptr[7] = (DCTELEM) DESCALE(tmp4 + z1 + z3, CONST_BITS-PASS1_BITS);
``` |
||

222 | ```
dataptr[5] = (DCTELEM) DESCALE(tmp5 + z2 + z4, CONST_BITS-PASS1_BITS);
``` |
||

223 | ```
dataptr[3] = (DCTELEM) DESCALE(tmp6 + z2 + z3, CONST_BITS-PASS1_BITS);
``` |
||

224 | ```
dataptr[1] = (DCTELEM) DESCALE(tmp7 + z1 + z4, CONST_BITS-PASS1_BITS);
``` |
||

225 | |||

226 | ```
dataptr += DCTSIZE; /* advance pointer to next row */
``` |
||

227 | } |
||

228 | |||

229 | ```
/* Pass 2: process columns.
``` |
||

230 | ```
* We remove the PASS1_BITS scaling, but leave the results scaled up
``` |
||

231 | ```
* by an overall factor of 8.
``` |
||

232 | ```
*/
``` |
||

233 | |||

234 | dataptr = data; |
||

235 | for (ctr = DCTSIZE-1; ctr >= 0; ctr--) { |
||

236 | tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*7]; |
||

237 | tmp7 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*7]; |
||

238 | tmp1 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*6]; |
||

239 | tmp6 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*6]; |
||

240 | tmp2 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*5]; |
||

241 | tmp5 = dataptr[DCTSIZE*2] - dataptr[DCTSIZE*5]; |
||

242 | tmp3 = dataptr[DCTSIZE*3] + dataptr[DCTSIZE*4]; |
||

243 | tmp4 = dataptr[DCTSIZE*3] - dataptr[DCTSIZE*4]; |
||

244 | |||

245 | ```
/* Even part per LL&M figure 1 --- note that published figure is faulty;
``` |
||

246 | ```
* rotator "sqrt(2)*c1" should be "sqrt(2)*c6".
``` |
||

247 | ```
*/
``` |
||

248 | |||

249 | tmp10 = tmp0 + tmp3; |
||

250 | tmp13 = tmp0 - tmp3; |
||

251 | tmp11 = tmp1 + tmp2; |
||

252 | tmp12 = tmp1 - tmp2; |
||

253 | |||

254 | ```
dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp11, PASS1_BITS);
``` |
||

255 | ```
dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp10 - tmp11, PASS1_BITS);
``` |
||

256 | |||

257 | z1 = MULTIPLY(tmp12 + tmp13, FIX_0_541196100); |
||

258 | ```
dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp13, FIX_0_765366865),
``` |
||

259 | CONST_BITS+PASS1_BITS); |
||

260 | ```
dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp12, - FIX_1_847759065),
``` |
||

261 | CONST_BITS+PASS1_BITS); |
||

262 | |||

263 | ```
/* Odd part per figure 8 --- note paper omits factor of sqrt(2).
``` |
||

264 | ```
* cK represents cos(K*pi/16).
``` |
||

265 | ```
* i0..i3 in the paper are tmp4..tmp7 here.
``` |
||

266 | ```
*/
``` |
||

267 | |||

268 | z1 = tmp4 + tmp7; |
||

269 | z2 = tmp5 + tmp6; |
||

270 | z3 = tmp4 + tmp6; |
||

271 | z4 = tmp5 + tmp7; |
||

272 | ```
z5 = MULTIPLY(z3 + z4, FIX_1_175875602); /* sqrt(2) * c3 */
``` |
||

273 | |||

274 | ```
tmp4 = MULTIPLY(tmp4, FIX_0_298631336); /* sqrt(2) * (-c1+c3+c5-c7) */
``` |
||

275 | ```
tmp5 = MULTIPLY(tmp5, FIX_2_053119869); /* sqrt(2) * ( c1+c3-c5+c7) */
``` |
||

276 | ```
tmp6 = MULTIPLY(tmp6, FIX_3_072711026); /* sqrt(2) * ( c1+c3+c5-c7) */
``` |
||

277 | ```
tmp7 = MULTIPLY(tmp7, FIX_1_501321110); /* sqrt(2) * ( c1+c3-c5-c7) */
``` |
||

278 | ```
z1 = MULTIPLY(z1, - FIX_0_899976223); /* sqrt(2) * (c7-c3) */
``` |
||

279 | ```
z2 = MULTIPLY(z2, - FIX_2_562915447); /* sqrt(2) * (-c1-c3) */
``` |
||

280 | ```
z3 = MULTIPLY(z3, - FIX_1_961570560); /* sqrt(2) * (-c3-c5) */
``` |
||

281 | ```
z4 = MULTIPLY(z4, - FIX_0_390180644); /* sqrt(2) * (c5-c3) */
``` |
||

282 | |||

283 | z3 += z5; |
||

284 | z4 += z5; |
||

285 | |||

286 | ```
dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp4 + z1 + z3,
``` |
||

287 | CONST_BITS+PASS1_BITS); |
||

288 | ```
dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp5 + z2 + z4,
``` |
||

289 | CONST_BITS+PASS1_BITS); |
||

290 | ```
dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp6 + z2 + z3,
``` |
||

291 | CONST_BITS+PASS1_BITS); |
||

292 | ```
dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp7 + z1 + z4,
``` |
||

293 | CONST_BITS+PASS1_BITS); |
||

294 | |||

295 | ```
dataptr++; /* advance pointer to next column */
``` |
||

296 | } |
||

297 | } |