Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

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

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 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

41
typedef 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

53
METHODDEF(void)
54
start_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

179
METHODDEF(void)
180
forward_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

269
METHODDEF(void)
270
forward_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

344
GLOBAL(void)
345
jinit_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

389