src/gdcmjpeg/jcdctmgr.c

   1 /*
   2  * jcdctmgr.c
   3  *
   4  * Copyright (C) 1994-1998, 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 "jlossy.h"    /* Private declarations for lossy codec */
  18 #include "jdct.h"    /* Private declarations for DCT subsystem */
  19
  20
  21 /* Private subobject for this module */
  22
  23 typedef struct {
  24   /* Pointer to the DCT routine actually in use */
  25   forward_DCT_method_ptr do_dct;
  26
  27   /* The actual post-DCT divisors --- not identical to the quant table
  28    * entries, because of scaling (especially for an unnormalized DCT).
  29    * Each table is given in normal array order.
  30    */
  31   DCTELEM * divisors[NUM_QUANT_TBLS];
  32
  33 #ifdef DCT_FLOAT_SUPPORTED
  34   /* Same as above for the floating-point case. */
  35   float_DCT_method_ptr do_float_dct;
  36   FAST_FLOAT * float_divisors[NUM_QUANT_TBLS];
  37 #endif
  38 } fdct_controller;
  39
  40 typedef fdct_controller * fdct_ptr;
  41
  42
  43 /*
  44  * Initialize for a processing pass.
  45  * Verify that all referenced Q-tables are present, and set up
  46  * the divisor table for each one.
  47  * In the current implementation, DCT of all components is done during
  48  * the first pass, even if only some components will be output in the
  49  * first scan.  Hence all components should be examined here.
  50  */
  51
  52 METHODDEF(void)
  53 start_pass_fdctmgr (j_compress_ptr cinfo)
  54 {
  55   j_lossy_c_ptr lossyc = (j_lossy_c_ptr) cinfo->codec;
  56   fdct_ptr fdct = (fdct_ptr) lossyc->fdct_private;
  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   j_lossy_c_ptr lossyc = (j_lossy_c_ptr) cinfo->codec;
 188   fdct_ptr fdct = (fdct_ptr) lossyc->fdct_private;
 189   forward_DCT_method_ptr do_dct = fdct->do_dct;
 190   DCTELEM * divisors = fdct->divisors[compptr->quant_tbl_no];
 191   DCTELEM workspace[DCTSIZE2];  /* work area for FDCT subroutine */
 192   JDIMENSION bi;
 193
 194   sample_data += start_row;  /* fold in the vertical offset once */
 195
 196   for (bi = 0; bi < num_blocks; bi++, start_col += DCTSIZE) {
 197     /* Load data into workspace, applying unsigned->signed conversion */
 198     { register DCTELEM *workspaceptr;
 199       register JSAMPROW elemptr;
 200       register int elemr;
 201
 202       workspaceptr = workspace;
 203       for (elemr = 0; elemr < DCTSIZE; elemr++) {
 204   elemptr = sample_data[elemr] + start_col;
 205 #if DCTSIZE == 8    /* unroll the inner loop */
 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   *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
 214 #else
 215   { register int elemc;
 216     for (elemc = DCTSIZE; elemc > 0; elemc--) {
 217       *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
 218     }
 219   }
 220 #endif
 221       }
 222     }
 223
 224     /* Perform the DCT */
 225     (*do_dct) (workspace);
 226
 227     /* Quantize/descale the coefficients, and store into coef_blocks[] */
 228     { register DCTELEM temp, qval;
 229       register int i;
 230       register JCOEFPTR output_ptr = coef_blocks[bi];
 231
 232       for (i = 0; i < DCTSIZE2; i++) {
 233   qval = divisors[i];
 234   temp = workspace[i];
 235   /* Divide the coefficient value by qval, ensuring proper rounding.
 236    * Since C does not specify the direction of rounding for negative
 237    * quotients, we have to force the dividend positive for portability.
 238    *
 239    * In most files, at least half of the output values will be zero
 240    * (at default quantization settings, more like three-quarters...)
 241    * so we should ensure that this case is fast.  On many machines,
 242    * a comparison is enough cheaper than a divide to make a special test
 243    * a win.  Since both inputs will be nonnegative, we need only test
 244    * for a < b to discover whether a/b is 0.
 245    * If your machine's division is fast enough, define FAST_DIVIDE.
 246    */
 247 #ifdef FAST_DIVIDE
 248 #define DIVIDE_BY(a,b)  a /= b
 249 #else
 250 #define DIVIDE_BY(a,b)  if (a >= b) a /= b; else a = 0
 251 #endif
 252   if (temp < 0) {
 253     temp = -temp;
 254     temp += qval>>1;  /* for rounding */
 255     DIVIDE_BY(temp, qval);
 256     temp = -temp;
 257   } else {
 258     temp += qval>>1;  /* for rounding */
 259     DIVIDE_BY(temp, qval);
 260   }
 261   output_ptr[i] = (JCOEF) temp;
 262       }
 263     }
 264   }
 265 }
 266
 267
 268 #ifdef DCT_FLOAT_SUPPORTED
 269
 270 METHODDEF(void)
 271 forward_DCT_float (j_compress_ptr cinfo, jpeg_component_info * compptr,
 272        JSAMPARRAY sample_data, JBLOCKROW coef_blocks,
 273        JDIMENSION start_row, JDIMENSION start_col,
 274        JDIMENSION num_blocks)
 275 /* This version is used for floating-point DCT implementations. */
 276 {
 277   /* This routine is heavily used, so it's worth coding it tightly. */
 278   j_lossy_c_ptr lossyc = (j_lossy_c_ptr) cinfo->codec;
 279   fdct_ptr fdct = (fdct_ptr) lossyc->fdct_private;
 280   float_DCT_method_ptr do_dct = fdct->do_float_dct;
 281   FAST_FLOAT * divisors = fdct->float_divisors[compptr->quant_tbl_no];
 282   FAST_FLOAT workspace[DCTSIZE2]; /* work area for FDCT subroutine */
 283   JDIMENSION bi;
 284
 285   sample_data += start_row;  /* fold in the vertical offset once */
 286
 287   for (bi = 0; bi < num_blocks; bi++, start_col += DCTSIZE) {
 288     /* Load data into workspace, applying unsigned->signed conversion */
 289     { register FAST_FLOAT *workspaceptr;
 290       register JSAMPROW elemptr;
 291       register int elemr;
 292
 293       workspaceptr = workspace;
 294       for (elemr = 0; elemr < DCTSIZE; elemr++) {
 295   elemptr = sample_data[elemr] + start_col;
 296 #if DCTSIZE == 8    /* unroll the inner loop */
 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   *workspaceptr++ = (FAST_FLOAT)(GETJSAMPLE(*elemptr++) - CENTERJSAMPLE);
 304   *workspaceptr++ = (FAST_FLOAT)(GETJSAMPLE(*elemptr++) - CENTERJSAMPLE);
 305 #else
 306   { register int elemc;
 307     for (elemc = DCTSIZE; elemc > 0; elemc--) {
 308       *workspaceptr++ = (FAST_FLOAT)
 309         (GETJSAMPLE(*elemptr++) - CENTERJSAMPLE);
 310     }
 311   }
 312 #endif
 313       }
 314     }
 315
 316     /* Perform the DCT */
 317     (*do_dct) (workspace);
 318
 319     /* Quantize/descale the coefficients, and store into coef_blocks[] */
 320     { register FAST_FLOAT temp;
 321       register int i;
 322       register JCOEFPTR output_ptr = coef_blocks[bi];
 323
 324       for (i = 0; i < DCTSIZE2; i++) {
 325   /* Apply the quantization and scaling factor */
 326   temp = workspace[i] * divisors[i];
 327   /* Round to nearest integer.
 328    * Since C does not specify the direction of rounding for negative
 329    * quotients, we have to force the dividend positive for portability.
 330    * The maximum coefficient size is +-16K (for 12-bit data), so this
 331    * code should work for either 16-bit or 32-bit ints.
 332    */
 333   output_ptr[i] = (JCOEF) ((int) (temp + (FAST_FLOAT) 16384.5) - 16384);
 334       }
 335     }
 336   }
 337 }
 338
 339 #endif /* DCT_FLOAT_SUPPORTED */
 340
 341
 342 /*
 343  * Initialize FDCT manager.
 344  */
 345
 346 GLOBAL(void)
 347 jinit_forward_dct (j_compress_ptr cinfo)
 348 {
 349   j_lossy_c_ptr lossyc = (j_lossy_c_ptr) cinfo->codec;
 350   fdct_ptr fdct;
 351   int i;
 352
 353   fdct = (fdct_ptr)
 354     (*cinfo->mem->alloc_small) ((j_common_ptr) cinfo, JPOOL_IMAGE,
 355         SIZEOF(fdct_controller));
 356   lossyc->fdct_private = (struct jpeg_forward_dct *) fdct;
 357   lossyc->fdct_start_pass = start_pass_fdctmgr;
 358
 359   switch (cinfo->dct_method) {
 360 #ifdef DCT_ISLOW_SUPPORTED
 361   case JDCT_ISLOW:
 362     lossyc->fdct_forward_DCT = forward_DCT;
 363     fdct->do_dct = jpeg_fdct_islow;
 364     break;
 365 #endif
 366 #ifdef DCT_IFAST_SUPPORTED
 367   case JDCT_IFAST:
 368     lossyc->fdct_forward_DCT = forward_DCT;
 369     fdct->do_dct = jpeg_fdct_ifast;
 370     break;
 371 #endif
 372 #ifdef DCT_FLOAT_SUPPORTED
 373   case JDCT_FLOAT:
 374     lossyc->fdct_forward_DCT = forward_DCT_float;
 375     fdct->do_float_dct = jpeg_fdct_float;
 376     break;
 377 #endif
 378   default:
 379     ERREXIT(cinfo, JERR_NOT_COMPILED);
 380     break;
 381   }
 382
 383   /* Mark divisor tables unallocated */
 384   for (i = 0; i < NUM_QUANT_TBLS; i++) {
 385     fdct->divisors[i] = NULL;
 386 #ifdef DCT_FLOAT_SUPPORTED
 387     fdct->float_divisors[i] = NULL;
 388 #endif
 389   }
 390 }