00001 /*
00002  * jcdctmgr.c
00003  *
00004  * Copyright (C) 1994-1996, Thomas G. Lane.
00005  * This file is part of the Independent JPEG Group's software.
00006  * For conditions of distribution and use, see the accompanying README file.
00007  *
00008  * This file contains the forward-DCT management logic.
00009  * This code selects a particular DCT implementation to be used,
00010  * and it performs related housekeeping chores including coefficient
00011  * quantization.
00012  */
00013 
00014 #define JPEG_INTERNALS
00015 #include "jinclude.h"
00016 #include "jpeglib.h"
00017 #include "jdct.h"               /* Private declarations for DCT subsystem */
00018 
00019 
00020 /* Private subobject for this module */
00021 
00022 typedef struct {
00023   struct jpeg_forward_dct pub;  /* public fields */
00024 
00025   /* Pointer to the DCT routine actually in use */
00026   forward_DCT_method_ptr do_dct;
00027 
00028   /* The actual post-DCT divisors --- not identical to the quant table
00029    * entries, because of scaling (especially for an unnormalized DCT).
00030    * Each table is given in normal array order.
00031    */
00032   DCTELEM * divisors[NUM_QUANT_TBLS];
00033 
00034 #ifdef DCT_FLOAT_SUPPORTED
00035   /* Same as above for the floating-point case. */
00036   float_DCT_method_ptr do_float_dct;
00037   FAST_FLOAT * float_divisors[NUM_QUANT_TBLS];
00038 #endif
00039 } my_fdct_controller;
00040 
00041 typedef my_fdct_controller * my_fdct_ptr;
00042 
00043 
00044 /*
00045  * Initialize for a processing pass.
00046  * Verify that all referenced Q-tables are present, and set up
00047  * the divisor table for each one.
00048  * In the current implementation, DCT of all components is done during
00049  * the first pass, even if only some components will be output in the
00050  * first scan.  Hence all components should be examined here.
00051  */
00052 
00053 METHODDEF(void)
00054 start_pass_fdctmgr (j_compress_ptr cinfo)
00055 {
00056   my_fdct_ptr fdct = (my_fdct_ptr) cinfo->fdct;
00057   int ci, qtblno, i;
00058   jpeg_component_info *compptr;
00059   JQUANT_TBL * qtbl;
00060   DCTELEM * dtbl;
00061 
00062   for (ci = 0, compptr = cinfo->comp_info; ci < cinfo->num_components;
00063        ci++, compptr++) {
00064     qtblno = compptr->quant_tbl_no;
00065     /* Make sure specified quantization table is present */
00066     if (qtblno < 0 || qtblno >= NUM_QUANT_TBLS ||
00067         cinfo->quant_tbl_ptrs[qtblno] == NULL)
00068       ERREXIT1(cinfo, JERR_NO_QUANT_TABLE, qtblno);
00069     qtbl = cinfo->quant_tbl_ptrs[qtblno];
00070     /* Compute divisors for this quant table */
00071     /* We may do this more than once for same table, but it's not a big deal */
00072     switch (cinfo->dct_method) {
00073 #ifdef DCT_ISLOW_SUPPORTED
00074     case JDCT_ISLOW:
00075       /* For LL&M IDCT method, divisors are equal to raw quantization
00076        * coefficients multiplied by 8 (to counteract scaling).
00077        */
00078       if (fdct->divisors[qtblno] == NULL) {
00079         fdct->divisors[qtblno] = (DCTELEM *)
00080           (*cinfo->mem->alloc_small) ((j_common_ptr) cinfo, JPOOL_IMAGE,
00081                                       DCTSIZE2 * SIZEOF(DCTELEM));
00082       }
00083       dtbl = fdct->divisors[qtblno];
00084       for (i = 0; i < DCTSIZE2; i++) {
00085         dtbl[i] = ((DCTELEM) qtbl->quantval[i]) << 3;
00086       }
00087       break;
00088 #endif
00089 #ifdef DCT_IFAST_SUPPORTED
00090     case JDCT_IFAST:
00091       {
00092         /* For AA&N IDCT method, divisors are equal to quantization
00093          * coefficients scaled by scalefactor[row]*scalefactor[col], where
00094          *   scalefactor[0] = 1
00095          *   scalefactor[k] = cos(k*PI/16) * sqrt(2)    for k=1..7
00096          * We apply a further scale factor of 8.
00097          */
00098 #define CONST_BITS 14
00099         static const INT16 aanscales[DCTSIZE2] = {
00100           /* precomputed values scaled up by 14 bits */
00101           16384, 22725, 21407, 19266, 16384, 12873,  8867,  4520,
00102           22725, 31521, 29692, 26722, 22725, 17855, 12299,  6270,
00103           21407, 29692, 27969, 25172, 21407, 16819, 11585,  5906,
00104           19266, 26722, 25172, 22654, 19266, 15137, 10426,  5315,
00105           16384, 22725, 21407, 19266, 16384, 12873,  8867,  4520,
00106           12873, 17855, 16819, 15137, 12873, 10114,  6967,  3552,
00107            8867, 12299, 11585, 10426,  8867,  6967,  4799,  2446,
00108            4520,  6270,  5906,  5315,  4520,  3552,  2446,  1247
00109         };
00110         SHIFT_TEMPS
00111 
00112         if (fdct->divisors[qtblno] == NULL) {
00113           fdct->divisors[qtblno] = (DCTELEM *)
00114             (*cinfo->mem->alloc_small) ((j_common_ptr) cinfo, JPOOL_IMAGE,
00115                                         DCTSIZE2 * SIZEOF(DCTELEM));
00116         }
00117         dtbl = fdct->divisors[qtblno];
00118         for (i = 0; i < DCTSIZE2; i++) {
00119           dtbl[i] = (DCTELEM)
00120             DESCALE(MULTIPLY16V16((INT32) qtbl->quantval[i],
00121                                   (INT32) aanscales[i]),
00122                     CONST_BITS-3);
00123         }
00124       }
00125       break;
00126 #endif
00127 #ifdef DCT_FLOAT_SUPPORTED
00128     case JDCT_FLOAT:
00129       {
00130         /* For float AA&N IDCT method, divisors are equal to quantization
00131          * coefficients scaled by scalefactor[row]*scalefactor[col], where
00132          *   scalefactor[0] = 1
00133          *   scalefactor[k] = cos(k*PI/16) * sqrt(2)    for k=1..7
00134          * We apply a further scale factor of 8.
00135          * What's actually stored is 1/divisor so that the inner loop can
00136          * use a multiplication rather than a division.
00137          */
00138         FAST_FLOAT * fdtbl;
00139         int row, col;
00140         static const double aanscalefactor[DCTSIZE] = {
00141           1.0, 1.387039845, 1.306562965, 1.175875602,
00142           1.0, 0.785694958, 0.541196100, 0.275899379
00143         };
00144 
00145         if (fdct->float_divisors[qtblno] == NULL) {
00146           fdct->float_divisors[qtblno] = (FAST_FLOAT *)
00147             (*cinfo->mem->alloc_small) ((j_common_ptr) cinfo, JPOOL_IMAGE,
00148                                         DCTSIZE2 * SIZEOF(FAST_FLOAT));
00149         }
00150         fdtbl = fdct->float_divisors[qtblno];
00151         i = 0;
00152         for (row = 0; row < DCTSIZE; row++) {
00153           for (col = 0; col < DCTSIZE; col++) {
00154             fdtbl[i] = (FAST_FLOAT)
00155               (1.0 / (((double) qtbl->quantval[i] *
00156                        aanscalefactor[row] * aanscalefactor[col] * 8.0)));
00157             i++;
00158           }
00159         }
00160       }
00161       break;
00162 #endif
00163     default:
00164       ERREXIT(cinfo, JERR_NOT_COMPILED);
00165       break;
00166     }
00167   }
00168 }
00169 
00170 
00171 /*
00172  * Perform forward DCT on one or more blocks of a component.
00173  *
00174  * The input samples are taken from the sample_data[] array starting at
00175  * position start_row/start_col, and moving to the right for any additional
00176  * blocks. The quantized coefficients are returned in coef_blocks[].
00177  */
00178 
00179 METHODDEF(void)
00180 forward_DCT (j_compress_ptr cinfo, jpeg_component_info * compptr,
00181              JSAMPARRAY sample_data, JBLOCKROW coef_blocks,
00182              JDIMENSION start_row, JDIMENSION start_col,
00183              JDIMENSION num_blocks)
00184 /* This version is used for integer DCT implementations. */
00185 {
00186   /* This routine is heavily used, so it's worth coding it tightly. */
00187   my_fdct_ptr fdct = (my_fdct_ptr) cinfo->fdct;
00188   forward_DCT_method_ptr do_dct = fdct->do_dct;
00189   DCTELEM * divisors = fdct->divisors[compptr->quant_tbl_no];
00190   DCTELEM workspace[DCTSIZE2];  /* work area for FDCT subroutine */
00191   JDIMENSION bi;
00192 
00193   sample_data += start_row;     /* fold in the vertical offset once */
00194 
00195   for (bi = 0; bi < num_blocks; bi++, start_col += DCTSIZE) {
00196     /* Load data into workspace, applying unsigned->signed conversion */
00197     { register DCTELEM *workspaceptr;
00198       register JSAMPROW elemptr;
00199       register int elemr;
00200 
00201       workspaceptr = workspace;
00202       for (elemr = 0; elemr < DCTSIZE; elemr++) {
00203         elemptr = sample_data[elemr] + start_col;
00204 #if DCTSIZE == 8                /* unroll the inner loop */
00205         *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
00206         *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
00207         *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
00208         *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
00209         *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
00210         *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
00211         *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
00212         *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
00213 #else
00214         { register int elemc;
00215           for (elemc = DCTSIZE; elemc > 0; elemc--) {
00216             *workspaceptr++ = GETJSAMPLE(*elemptr++) - CENTERJSAMPLE;
00217           }
00218         }
00219 #endif
00220       }
00221     }
00222 
00223     /* Perform the DCT */
00224     (*do_dct) (workspace);
00225 
00226     /* Quantize/descale the coefficients, and store into coef_blocks[] */
00227     { register DCTELEM temp, qval;
00228       register int i;
00229       register JCOEFPTR output_ptr = coef_blocks[bi];
00230 
00231       for (i = 0; i < DCTSIZE2; i++) {
00232         qval = divisors[i];
00233         temp = workspace[i];
00234         /* Divide the coefficient value by qval, ensuring proper rounding.
00235          * Since C does not specify the direction of rounding for negative
00236          * quotients, we have to force the dividend positive for portability.
00237          *
00238          * In most files, at least half of the output values will be zero
00239          * (at default quantization settings, more like three-quarters...)
00240          * so we should ensure that this case is fast.  On many machines,
00241          * a comparison is enough cheaper than a divide to make a special test
00242          * a win.  Since both inputs will be nonnegative, we need only test
00243          * for a < b to discover whether a/b is 0.
00244          * If your machine's division is fast enough, define FAST_DIVIDE.
00245          */
00246 #ifdef FAST_DIVIDE
00247 #define DIVIDE_BY(a,b)  a /= b
00248 #else
00249 #define DIVIDE_BY(a,b)  if (a >= b) a /= b; else a = 0
00250 #endif
00251         if (temp < 0) {
00252           temp = -temp;
00253           temp += qval>>1;      /* for rounding */
00254           DIVIDE_BY(temp, qval);
00255           temp = -temp;
00256         } else {
00257           temp += qval>>1;      /* for rounding */
00258           DIVIDE_BY(temp, qval);
00259         }
00260         output_ptr[i] = (JCOEF) temp;
00261       }
00262     }
00263   }
00264 }
00265 
00266 
00267 #ifdef DCT_FLOAT_SUPPORTED
00268 
00269 METHODDEF(void)
00270 forward_DCT_float (j_compress_ptr cinfo, jpeg_component_info * compptr,
00271                    JSAMPARRAY sample_data, JBLOCKROW coef_blocks,
00272                    JDIMENSION start_row, JDIMENSION start_col,
00273                    JDIMENSION num_blocks)
00274 /* This version is used for floating-point DCT implementations. */
00275 {
00276   /* This routine is heavily used, so it's worth coding it tightly. */
00277   my_fdct_ptr fdct = (my_fdct_ptr) cinfo->fdct;
00278   float_DCT_method_ptr do_dct = fdct->do_float_dct;
00279   FAST_FLOAT * divisors = fdct->float_divisors[compptr->quant_tbl_no];
00280   FAST_FLOAT workspace[DCTSIZE2]; /* work area for FDCT subroutine */
00281   JDIMENSION bi;
00282 
00283   sample_data += start_row;     /* fold in the vertical offset once */
00284 
00285   for (bi = 0; bi < num_blocks; bi++, start_col += DCTSIZE) {
00286     /* Load data into workspace, applying unsigned->signed conversion */
00287     { register FAST_FLOAT *workspaceptr;
00288       register JSAMPROW elemptr;
00289       register int elemr;
00290 
00291       workspaceptr = workspace;
00292       for (elemr = 0; elemr < DCTSIZE; elemr++) {
00293         elemptr = sample_data[elemr] + start_col;
00294 #if DCTSIZE == 8                /* unroll the inner loop */
00295         *workspaceptr++ = (FAST_FLOAT)(GETJSAMPLE(*elemptr++) - CENTERJSAMPLE);
00296         *workspaceptr++ = (FAST_FLOAT)(GETJSAMPLE(*elemptr++) - CENTERJSAMPLE);
00297         *workspaceptr++ = (FAST_FLOAT)(GETJSAMPLE(*elemptr++) - CENTERJSAMPLE);
00298         *workspaceptr++ = (FAST_FLOAT)(GETJSAMPLE(*elemptr++) - CENTERJSAMPLE);
00299         *workspaceptr++ = (FAST_FLOAT)(GETJSAMPLE(*elemptr++) - CENTERJSAMPLE);
00300         *workspaceptr++ = (FAST_FLOAT)(GETJSAMPLE(*elemptr++) - CENTERJSAMPLE);
00301         *workspaceptr++ = (FAST_FLOAT)(GETJSAMPLE(*elemptr++) - CENTERJSAMPLE);
00302         *workspaceptr++ = (FAST_FLOAT)(GETJSAMPLE(*elemptr++) - CENTERJSAMPLE);
00303 #else
00304         { register int elemc;
00305           for (elemc = DCTSIZE; elemc > 0; elemc--) {
00306             *workspaceptr++ = (FAST_FLOAT)
00307               (GETJSAMPLE(*elemptr++) - CENTERJSAMPLE);
00308           }
00309         }
00310 #endif
00311       }
00312     }
00313 
00314     /* Perform the DCT */
00315     (*do_dct) (workspace);
00316 
00317     /* Quantize/descale the coefficients, and store into coef_blocks[] */
00318     { register FAST_FLOAT temp;
00319       register int i;
00320       register JCOEFPTR output_ptr = coef_blocks[bi];
00321 
00322       for (i = 0; i < DCTSIZE2; i++) {
00323         /* Apply the quantization and scaling factor */
00324         temp = workspace[i] * divisors[i];
00325         /* Round to nearest integer.
00326          * Since C does not specify the direction of rounding for negative
00327          * quotients, we have to force the dividend positive for portability.
00328          * The maximum coefficient size is +-16K (for 12-bit data), so this
00329          * code should work for either 16-bit or 32-bit ints.
00330          */
00331         output_ptr[i] = (JCOEF) ((int) (temp + (FAST_FLOAT) 16384.5) - 16384);
00332       }
00333     }
00334   }
00335 }
00336 
00337 #endif /* DCT_FLOAT_SUPPORTED */
00338 
00339 
00340 /*
00341  * Initialize FDCT manager.
00342  */
00343 
00344 GLOBAL(void)
00345 jinit_forward_dct (j_compress_ptr cinfo)
00346 {
00347   my_fdct_ptr fdct;
00348   int i;
00349 
00350   fdct = (my_fdct_ptr)
00351     (*cinfo->mem->alloc_small) ((j_common_ptr) cinfo, JPOOL_IMAGE,
00352                                 SIZEOF(my_fdct_controller));
00353   cinfo->fdct = (struct jpeg_forward_dct *) fdct;
00354   fdct->pub.start_pass = start_pass_fdctmgr;
00355 
00356   switch (cinfo->dct_method) {
00357 #ifdef DCT_ISLOW_SUPPORTED
00358   case JDCT_ISLOW:
00359     fdct->pub.forward_DCT = forward_DCT;
00360     fdct->do_dct = jpeg_fdct_islow;
00361     break;
00362 #endif
00363 #ifdef DCT_IFAST_SUPPORTED
00364   case JDCT_IFAST:
00365     fdct->pub.forward_DCT = forward_DCT;
00366     fdct->do_dct = jpeg_fdct_ifast;
00367     break;
00368 #endif
00369 #ifdef DCT_FLOAT_SUPPORTED
00370   case JDCT_FLOAT:
00371     fdct->pub.forward_DCT = forward_DCT_float;
00372     fdct->do_float_dct = jpeg_fdct_float;
00373     break;
00374 #endif
00375   default:
00376     ERREXIT(cinfo, JERR_NOT_COMPILED);
00377     break;
00378   }
00379 
00380   /* Mark divisor tables unallocated */
00381   for (i = 0; i < NUM_QUANT_TBLS; i++) {
00382     fdct->divisors[i] = NULL;
00383 #ifdef DCT_FLOAT_SUPPORTED
00384     fdct->float_divisors[i] = NULL;
00385 #endif
00386   }
00387 }