1 /*
2 * jcdctmgr.c
3 *
4 * Copyright (C) 1994, 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 the forward-DCT management logic.
9 * This code selects a particular DCT implementation to be used,
10 * and it performs related housekeeping chores including coefficient
11 * quantization.
12 */
13
14 #define JPEG_INTERNALS
15 #include "jinclude.h"
16 #include "jpeglib.h"
17 #include "jdct.h" /* Private declarations for DCT subsystem */
18
19
20 /* Private subobject for this module */
21
22 typedef struct {
23 struct jpeg_forward_dct pub; /* public fields */
24
25 /* Pointer to the DCT routine actually in use */
26 forward_DCT_method_ptr do_dct;
27
28 /* The actual post-DCT divisors --- not identical to the quant table
29 * entries, because of scaling (especially for an unnormalized DCT).
30 * Each table is given in zigzag order.
31 */
32 DCTELEM * divisors[NUM_QUANT_TBLS];
33
34 #ifdef DCT_FLOAT_SUPPORTED
35 /* Same as above for the floating-point case. */
36 float_DCT_method_ptr do_float_dct;
37 FAST_FLOAT * float_divisors[NUM_QUANT_TBLS];
38 #endif
39 } my_fdct_controller;
40
41 typedef my_fdct_controller * my_fdct_ptr;
42
43
44 /* ZAG[i] is the natural-order position of the i'th element of zigzag order. */
45
46 static const int ZAG[DCTSIZE2] = {
47 0, 1, 8, 16, 9, 2, 3, 10,
48 17, 24, 32, 25, 18, 11, 4, 5,
49 12, 19, 26, 33, 40, 48, 41, 34,
50 27, 20, 13, 6, 7, 14, 21, 28,
51 35, 42, 49, 56, 57, 50, 43, 36,
52 29, 22, 15, 23, 30, 37, 44, 51,
53 58, 59, 52, 45, 38, 31, 39, 46,
54 53, 60, 61, 54, 47, 55, 62, 63
55 };
56
57
58 /*
59 * Initialize for a processing pass.
60 * Verify that all referenced Q-tables are present, and set up
61 * the divisor table for each one.
62 * In the current implementation, DCT of all components is done during
63 * the first pass, even if only some components will be output in the
64 * first scan. Hence all components should be examined here.
65 */
66
67 METHODDEF void
start_pass_fdctmgr(j_compress_ptr cinfo)68 start_pass_fdctmgr (j_compress_ptr cinfo)
69 {
70 my_fdct_ptr fdct = (my_fdct_ptr) cinfo->fdct;
71 int ci, qtblno, i;
72 jpeg_component_info *compptr;
73 JQUANT_TBL * qtbl;
74 DCTELEM * dtbl;
75
76 for (ci = 0, compptr = cinfo->comp_info; ci < cinfo->num_components;
77 ci++, compptr++) {
78 qtblno = compptr->quant_tbl_no;
79 /* Make sure specified quantization table is present */
80 if (qtblno < 0 || qtblno >= NUM_QUANT_TBLS ||
81 cinfo->quant_tbl_ptrs[qtblno] == NULL)
82 ERREXIT1(cinfo, JERR_NO_QUANT_TABLE, qtblno);
83 qtbl = cinfo->quant_tbl_ptrs[qtblno];
84 /* Compute divisors for this quant table */
85 /* We may do this more than once for same table, but it's not a big deal */
86 switch (cinfo->dct_method) {
87 #ifdef DCT_ISLOW_SUPPORTED
88 case JDCT_ISLOW:
89 /* For LL&M IDCT method, divisors are equal to raw quantization
90 * coefficients multiplied by 8 (to counteract scaling).
91 */
92 if (fdct->divisors[qtblno] == NULL) {
93 fdct->divisors[qtblno] = (DCTELEM *)
94 (*cinfo->mem->alloc_small) ((j_common_ptr) cinfo, JPOOL_IMAGE,
95 DCTSIZE2 * SIZEOF(DCTELEM));
96 }
97 dtbl = fdct->divisors[qtblno];
98 for (i = 0; i < DCTSIZE2; i++) {
99 dtbl[i] = ((DCTELEM) qtbl->quantval[i]) << 3;
100 }
101 break;
102 #endif
103 #ifdef DCT_IFAST_SUPPORTED
104 case JDCT_IFAST:
105 {
106 /* For AA&N IDCT method, divisors are equal to quantization
107 * coefficients scaled by scalefactor[row]*scalefactor[col], where
108 * scalefactor[0] = 1
109 * scalefactor[k] = cos(k*PI/16) * sqrt(2) for k=1..7
110 * We apply a further scale factor of 8.
111 */
112 #define CONST_BITS 14
113 static const INT16 aanscales[DCTSIZE2] = {
114 /* precomputed values scaled up by 14 bits: in natural order */
115 16384, 22725, 21407, 19266, 16384, 12873, 8867, 4520,
116 22725, 31521, 29692, 26722, 22725, 17855, 12299, 6270,
117 21407, 29692, 27969, 25172, 21407, 16819, 11585, 5906,
118 19266, 26722, 25172, 22654, 19266, 15137, 10426, 5315,
119 16384, 22725, 21407, 19266, 16384, 12873, 8867, 4520,
120 12873, 17855, 16819, 15137, 12873, 10114, 6967, 3552,
121 8867, 12299, 11585, 10426, 8867, 6967, 4799, 2446,
122 4520, 6270, 5906, 5315, 4520, 3552, 2446, 1247
123 };
124 SHIFT_TEMPS
125
126 if (fdct->divisors[qtblno] == NULL) {
127 fdct->divisors[qtblno] = (DCTELEM *)
128 (*cinfo->mem->alloc_small) ((j_common_ptr) cinfo, JPOOL_IMAGE,
129 DCTSIZE2 * SIZEOF(DCTELEM));
130 }
131 dtbl = fdct->divisors[qtblno];
132 for (i = 0; i < DCTSIZE2; i++) {
133 dtbl[i] = (DCTELEM)
134 DESCALE(MULTIPLY16V16((INT32) qtbl->quantval[i],
135 (INT32) aanscales[ZAG[i]]),
136 CONST_BITS-3);
137 }
138 }
139 break;
140 #endif
141 #ifdef DCT_FLOAT_SUPPORTED
142 case JDCT_FLOAT:
143 {
144 /* For float AA&N IDCT method, divisors are equal to quantization
145 * coefficients scaled by scalefactor[row]*scalefactor[col], where
146 * scalefactor[0] = 1
147 * scalefactor[k] = cos(k*PI/16) * sqrt(2) for k=1..7
148 * We apply a further scale factor of 8.
149 * What's actually stored is 1/divisor so that the inner loop can
150 * use a multiplication rather than a division.
151 */
152 FAST_FLOAT * fdtbl;
153 int row, col;
154 static const double aanscalefactor[DCTSIZE] = {
155 1.0, 1.387039845, 1.306562965, 1.175875602,
156 1.0, 0.785694958, 0.541196100, 0.275899379
157 };
158
159 if (fdct->float_divisors[qtblno] == NULL) {
160 fdct->float_divisors[qtblno] = (FAST_FLOAT *)
161 (*cinfo->mem->alloc_small) ((j_common_ptr) cinfo, JPOOL_IMAGE,
162 DCTSIZE2 * SIZEOF(FAST_FLOAT));
163 }
164 fdtbl = fdct->float_divisors[qtblno];
165 for (i = 0; i < DCTSIZE2; i++) {
166 row = ZAG[i] >> 3;
167 col = ZAG[i] & 7;
168 fdtbl[i] = (FAST_FLOAT)
169 (1.0 / (((double) qtbl->quantval[i] *
170 aanscalefactor[row] * aanscalefactor[col] * 8.0)));
171 }
172 }
173 break;
174 #endif
175 default:
176 ERREXIT(cinfo, JERR_NOT_COMPILED);
177 break;
178 }
179 }
180 }
181
182
183 /*
184 * Perform forward DCT on one or more blocks of a component.
185 *
186 * The input samples are taken from the sample_data[] array starting at
187 * position start_row/start_col, and moving to the right for any additional
188 * blocks. The quantized, zigzagged coefficients are returned in coef_blocks[].
189 */
190
191 METHODDEF void
forward_DCT(j_compress_ptr cinfo,jpeg_component_info * compptr,JSAMPARRAY sample_data,JBLOCKROW coef_blocks,JDIMENSION start_row,JDIMENSION start_col,JDIMENSION num_blocks)192 forward_DCT (j_compress_ptr cinfo, jpeg_component_info * compptr,
193 JSAMPARRAY sample_data, JBLOCKROW coef_blocks,
194 JDIMENSION start_row, JDIMENSION start_col,
195 JDIMENSION num_blocks)
196 /* This version is used for integer DCT implementations. */
197 {
198 /* This routine is heavily used, so it's worth coding it tightly. */
199 my_fdct_ptr fdct = (my_fdct_ptr) cinfo->fdct;
200 forward_DCT_method_ptr do_dct = fdct->do_dct;
201 DCTELEM * divisors = fdct->divisors[compptr->quant_tbl_no];
202 DCTELEM workspace[DCTSIZE2]; /* work area for FDCT subroutine */
203 JDIMENSION bi;
204
205 sample_data += start_row; /* fold in the vertical offset once */
206
207 for (bi = 0; bi < num_blocks; bi++, start_col += DCTSIZE) {
208 /* Load data into workspace, applying unsigned->signed conversion */
209 { register DCTELEM *workspaceptr;
210 register JSAMPROW elemptr;
211 register int elemr;
212
213 workspaceptr = workspace;
214 for (elemr = 0; elemr < DCTSIZE; elemr++) {
215 elemptr = sample_data[elemr] + start_col;
216 #if DCTSIZE == 8 /* unroll the inner loop */
217 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
218 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
219 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
220 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
221 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
222 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
223 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
224 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
225 #else
226 { register int elemc;
227 for (elemc = DCTSIZE; elemc > 0; elemc--) {
228 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
229 }
230 }
231 #endif
232 }
233 }
234
235 /* Perform the DCT */
236 (*do_dct) (workspace);
237
238 /* Quantize/descale the coefficients, and store into coef_blocks[] */
239 { register DCTELEM temp, qval;
240 register int i;
241 register JCOEFPTR output_ptr = coef_blocks[bi];
242
243 for (i = 0; i < DCTSIZE2; i++) {
244 qval = divisors[i];
245 temp = workspace[ZAG[i]];
246 /* Divide the coefficient value by qval, ensuring proper rounding.
247 * Since C does not specify the direction of rounding for negative
248 * quotients, we have to force the dividend positive for portability.
249 *
250 * In most files, at least half of the output values will be zero
251 * (at default quantization settings, more like three-quarters...)
252 * so we should ensure that this case is fast. On many machines,
253 * a comparison is enough cheaper than a divide to make a special test
254 * a win. Since both inputs will be nonnegative, we need only test
255 * for a < b to discover whether a/b is 0.
256 * If your machine's division is fast enough, define FAST_DIVIDE.
257 */
258 #ifdef FAST_DIVIDE
259 #define DIVIDE_BY(a,b) a /= b
260 #else
261 #define DIVIDE_BY(a,b) if (a >= b) a /= b; else a = 0
262 #endif
263 if (temp < 0) {
264 temp = -temp;
265 temp += qval>>1; /* for rounding */
266 DIVIDE_BY(temp, qval);
267 temp = -temp;
268 } else {
269 temp += qval>>1; /* for rounding */
270 DIVIDE_BY(temp, qval);
271 }
272 output_ptr[i] = (JCOEF) temp;
273 }
274 }
275 }
276 }
277
278
279 #ifdef DCT_FLOAT_SUPPORTED
280
281 METHODDEF void
forward_DCT_float(j_compress_ptr cinfo,jpeg_component_info * compptr,JSAMPARRAY sample_data,JBLOCKROW coef_blocks,JDIMENSION start_row,JDIMENSION start_col,JDIMENSION num_blocks)282 forward_DCT_float (j_compress_ptr cinfo, jpeg_component_info * compptr,
283 JSAMPARRAY sample_data, JBLOCKROW coef_blocks,
284 JDIMENSION start_row, JDIMENSION start_col,
285 JDIMENSION num_blocks)
286 /* This version is used for floating-point DCT implementations. */
287 {
288 /* This routine is heavily used, so it's worth coding it tightly. */
289 my_fdct_ptr fdct = (my_fdct_ptr) cinfo->fdct;
290 float_DCT_method_ptr do_dct = fdct->do_float_dct;
291 FAST_FLOAT * divisors = fdct->float_divisors[compptr->quant_tbl_no];
292 FAST_FLOAT workspace[DCTSIZE2]; /* work area for FDCT subroutine */
293 JDIMENSION bi;
294
295 sample_data += start_row; /* fold in the vertical offset once */
296
297 for (bi = 0; bi < num_blocks; bi++, start_col += DCTSIZE) {
298 /* Load data into workspace, applying unsigned->signed conversion */
299 { register FAST_FLOAT *workspaceptr;
300 register JSAMPROW elemptr;
301 register int elemr;
302
303 workspaceptr = workspace;
304 for (elemr = 0; elemr < DCTSIZE; elemr++) {
305 elemptr = sample_data[elemr] + start_col;
306 #if DCTSIZE == 8 /* unroll the inner loop */
307 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
308 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
309 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
310 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
311 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
312 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
313 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
314 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
315 #else
316 { register int elemc;
317 for (elemc = DCTSIZE; elemc > 0; elemc--) {
318 *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
319 }
320 }
321 #endif
322 }
323 }
324
325 /* Perform the DCT */
326 (*do_dct) (workspace);
327
328 /* Quantize/descale the coefficients, and store into coef_blocks[] */
329 { register FAST_FLOAT temp;
330 register int i;
331 register JCOEFPTR output_ptr = coef_blocks[bi];
332
333 for (i = 0; i < DCTSIZE2; i++) {
334 /* Apply the quantization and scaling factor */
335 temp = workspace[ZAG[i]] * divisors[i];
336 /* Round to nearest integer.
337 * Since C does not specify the direction of rounding for negative
338 * quotients, we have to force the dividend positive for portability.
339 * The maximum coefficient size is +-16K (for 12-bit data), so this
340 * code should work for either 16-bit or 32-bit ints.
341 */
342 output_ptr[i] = (JCOEF) ((int) (temp + (FAST_FLOAT) 16384.5) - 16384);
343 }
344 }
345 }
346 }
347
348 #endif /* DCT_FLOAT_SUPPORTED */
349
350
351 /*
352 * Initialize FDCT manager.
353 */
354
355 GLOBAL void
jinit_forward_dct(j_compress_ptr cinfo)356 jinit_forward_dct (j_compress_ptr cinfo)
357 {
358 my_fdct_ptr fdct;
359 int i;
360
361 fdct = (my_fdct_ptr)
362 (*cinfo->mem->alloc_small) ((j_common_ptr) cinfo, JPOOL_IMAGE,
363 SIZEOF(my_fdct_controller));
364 cinfo->fdct = (struct jpeg_forward_dct *) fdct;
365 fdct->pub.start_pass = start_pass_fdctmgr;
366
367 switch (cinfo->dct_method) {
368 #ifdef DCT_ISLOW_SUPPORTED
369 case JDCT_ISLOW:
370 fdct->pub.forward_DCT = forward_DCT;
371 fdct->do_dct = jpeg_fdct_islow;
372 break;
373 #endif
374 #ifdef DCT_IFAST_SUPPORTED
375 case JDCT_IFAST:
376 fdct->pub.forward_DCT = forward_DCT;
377 fdct->do_dct = jpeg_fdct_ifast;
378 break;
379 #endif
380 #ifdef DCT_FLOAT_SUPPORTED
381 case JDCT_FLOAT:
382 fdct->pub.forward_DCT = forward_DCT_float;
383 fdct->do_float_dct = jpeg_fdct_float;
384 break;
385 #endif
386 default:
387 ERREXIT(cinfo, JERR_NOT_COMPILED);
388 break;
389 }
390
391 /* Mark divisor tables unallocated */
392 for (i = 0; i < NUM_QUANT_TBLS; i++) {
393 fdct->divisors[i] = NULL;
394 #ifdef DCT_FLOAT_SUPPORTED
395 fdct->float_divisors[i] = NULL;
396 #endif
397 }
398 }
399