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