Home | History | Annotate | Download | only in jpeg-6b
      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