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
