jfdctfst.c

00001 /*
00002  * jfdctfst.c
00003  *
00004  * Copyright (C) 1994-1996, Thomas G. Lane.
00005  * Modified 2003-2015 by Guido Vollbeding.
00006  * This file is part of the Independent JPEG Group's software.
00007  * For conditions of distribution and use, see the accompanying README file.
00008  *
00009  * This file contains a fast, not so accurate integer implementation of the
00010  * forward DCT (Discrete Cosine Transform).
00011  *
00012  * A 2-D DCT can be done by 1-D DCT on each row followed by 1-D DCT
00013  * on each column.  Direct algorithms are also available, but they are
00014  * much more complex and seem not to be any faster when reduced to code.
00015  *
00016  * This implementation is based on Arai, Agui, and Nakajima's algorithm for
00017  * scaled DCT.  Their original paper (Trans. IEICE E-71(11):1095) is in
00018  * Japanese, but the algorithm is described in the Pennebaker & Mitchell
00019  * JPEG textbook (see REFERENCES section in file README).  The following code
00020  * is based directly on figure 4-8 in P&M.
00021  * While an 8-point DCT cannot be done in less than 11 multiplies, it is
00022  * possible to arrange the computation so that many of the multiplies are
00023  * simple scalings of the final outputs.  These multiplies can then be
00024  * folded into the multiplications or divisions by the JPEG quantization
00025  * table entries.  The AA&N method leaves only 5 multiplies and 29 adds
00026  * to be done in the DCT itself.
00027  * The primary disadvantage of this method is that with fixed-point math,
00028  * accuracy is lost due to imprecise representation of the scaled
00029  * quantization values.  The smaller the quantization table entry, the less
00030  * precise the scaled value, so this implementation does worse with high-
00031  * quality-setting files than with low-quality ones.
00032  */
00033 
00034 #define JPEG_INTERNALS
00035 #include "jinclude.h"
00036 #include "jpeglib.h"
00037 #include "jdct.h"       /* Private declarations for DCT subsystem */
00038 
00039 #ifdef DCT_IFAST_SUPPORTED
00040 
00041 
00042 /*
00043  * This module is specialized to the case DCTSIZE = 8.
00044  */
00045 
00046 #if DCTSIZE != 8
00047   Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */
00048 #endif
00049 
00050 
00051 /* Scaling decisions are generally the same as in the LL&M algorithm;
00052  * see jfdctint.c for more details.  However, we choose to descale
00053  * (right shift) multiplication products as soon as they are formed,
00054  * rather than carrying additional fractional bits into subsequent additions.
00055  * This compromises accuracy slightly, but it lets us save a few shifts.
00056  * More importantly, 16-bit arithmetic is then adequate (for 8-bit samples)
00057  * everywhere except in the multiplications proper; this saves a good deal
00058  * of work on 16-bit-int machines.
00059  *
00060  * Again to save a few shifts, the intermediate results between pass 1 and
00061  * pass 2 are not upscaled, but are represented only to integral precision.
00062  *
00063  * A final compromise is to represent the multiplicative constants to only
00064  * 8 fractional bits, rather than 13.  This saves some shifting work on some
00065  * machines, and may also reduce the cost of multiplication (since there
00066  * are fewer one-bits in the constants).
00067  */
00068 
00069 #define CONST_BITS  8
00070 
00071 
00072 /* Some C compilers fail to reduce "FIX(constant)" at compile time, thus
00073  * causing a lot of useless floating-point operations at run time.
00074  * To get around this we use the following pre-calculated constants.
00075  * If you change CONST_BITS you may want to add appropriate values.
00076  * (With a reasonable C compiler, you can just rely on the FIX() macro...)
00077  */
00078 
00079 #if CONST_BITS == 8
00080 #define FIX_0_382683433  ((INT32)   98)     /* FIX(0.382683433) */
00081 #define FIX_0_541196100  ((INT32)  139)     /* FIX(0.541196100) */
00082 #define FIX_0_707106781  ((INT32)  181)     /* FIX(0.707106781) */
00083 #define FIX_1_306562965  ((INT32)  334)     /* FIX(1.306562965) */
00084 #else
00085 #define FIX_0_382683433  FIX(0.382683433)
00086 #define FIX_0_541196100  FIX(0.541196100)
00087 #define FIX_0_707106781  FIX(0.707106781)
00088 #define FIX_1_306562965  FIX(1.306562965)
00089 #endif
00090 
00091 
00092 /* We can gain a little more speed, with a further compromise in accuracy,
00093  * by omitting the addition in a descaling shift.  This yields an incorrectly
00094  * rounded result half the time...
00095  */
00096 
00097 #ifndef USE_ACCURATE_ROUNDING
00098 #undef DESCALE
00099 #define DESCALE(x,n)  RIGHT_SHIFT(x, n)
00100 #endif
00101 
00102 
00103 /* Multiply a DCTELEM variable by an INT32 constant, and immediately
00104  * descale to yield a DCTELEM result.
00105  */
00106 
00107 #define MULTIPLY(var,const)  ((DCTELEM) DESCALE((var) * (const), CONST_BITS))
00108 
00109 
00110 /*
00111  * Perform the forward DCT on one block of samples.
00112  *
00113  * cK represents cos(K*pi/16).
00114  */
00115 
00116 GLOBAL(void)
00117 jpeg_fdct_ifast (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
00118 {
00119   DCTELEM tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
00120   DCTELEM tmp10, tmp11, tmp12, tmp13;
00121   DCTELEM z1, z2, z3, z4, z5, z11, z13;
00122   DCTELEM *dataptr;
00123   JSAMPROW elemptr;
00124   int ctr;
00125   SHIFT_TEMPS
00126 
00127   /* Pass 1: process rows. */
00128 
00129   dataptr = data;
00130   for (ctr = 0; ctr < DCTSIZE; ctr++) {
00131     elemptr = sample_data[ctr] + start_col;
00132 
00133     /* Load data into workspace */
00134     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[7]);
00135     tmp7 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[7]);
00136     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[6]);
00137     tmp6 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[6]);
00138     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[5]);
00139     tmp5 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[5]);
00140     tmp3 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[4]);
00141     tmp4 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[4]);
00142 
00143     /* Even part */
00144 
00145     tmp10 = tmp0 + tmp3;    /* phase 2 */
00146     tmp13 = tmp0 - tmp3;
00147     tmp11 = tmp1 + tmp2;
00148     tmp12 = tmp1 - tmp2;
00149 
00150     /* Apply unsigned->signed conversion. */
00151     dataptr[0] = tmp10 + tmp11 - 8 * CENTERJSAMPLE; /* phase 3 */
00152     dataptr[4] = tmp10 - tmp11;
00153 
00154     z1 = MULTIPLY(tmp12 + tmp13, FIX_0_707106781); /* c4 */
00155     dataptr[2] = tmp13 + z1;    /* phase 5 */
00156     dataptr[6] = tmp13 - z1;
00157 
00158     /* Odd part */
00159 
00160     tmp10 = tmp4 + tmp5;    /* phase 2 */
00161     tmp11 = tmp5 + tmp6;
00162     tmp12 = tmp6 + tmp7;
00163 
00164     /* The rotator is modified from fig 4-8 to avoid extra negations. */
00165     z5 = MULTIPLY(tmp10 - tmp12, FIX_0_382683433); /* c6 */
00166     z2 = MULTIPLY(tmp10, FIX_0_541196100) + z5; /* c2-c6 */
00167     z4 = MULTIPLY(tmp12, FIX_1_306562965) + z5; /* c2+c6 */
00168     z3 = MULTIPLY(tmp11, FIX_0_707106781); /* c4 */
00169 
00170     z11 = tmp7 + z3;        /* phase 5 */
00171     z13 = tmp7 - z3;
00172 
00173     dataptr[5] = z13 + z2;  /* phase 6 */
00174     dataptr[3] = z13 - z2;
00175     dataptr[1] = z11 + z4;
00176     dataptr[7] = z11 - z4;
00177 
00178     dataptr += DCTSIZE;     /* advance pointer to next row */
00179   }
00180 
00181   /* Pass 2: process columns. */
00182 
00183   dataptr = data;
00184   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
00185     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*7];
00186     tmp7 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*7];
00187     tmp1 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*6];
00188     tmp6 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*6];
00189     tmp2 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*5];
00190     tmp5 = dataptr[DCTSIZE*2] - dataptr[DCTSIZE*5];
00191     tmp3 = dataptr[DCTSIZE*3] + dataptr[DCTSIZE*4];
00192     tmp4 = dataptr[DCTSIZE*3] - dataptr[DCTSIZE*4];
00193 
00194     /* Even part */
00195 
00196     tmp10 = tmp0 + tmp3;    /* phase 2 */
00197     tmp13 = tmp0 - tmp3;
00198     tmp11 = tmp1 + tmp2;
00199     tmp12 = tmp1 - tmp2;
00200 
00201     dataptr[DCTSIZE*0] = tmp10 + tmp11; /* phase 3 */
00202     dataptr[DCTSIZE*4] = tmp10 - tmp11;
00203 
00204     z1 = MULTIPLY(tmp12 + tmp13, FIX_0_707106781); /* c4 */
00205     dataptr[DCTSIZE*2] = tmp13 + z1; /* phase 5 */
00206     dataptr[DCTSIZE*6] = tmp13 - z1;
00207 
00208     /* Odd part */
00209 
00210     tmp10 = tmp4 + tmp5;    /* phase 2 */
00211     tmp11 = tmp5 + tmp6;
00212     tmp12 = tmp6 + tmp7;
00213 
00214     /* The rotator is modified from fig 4-8 to avoid extra negations. */
00215     z5 = MULTIPLY(tmp10 - tmp12, FIX_0_382683433); /* c6 */
00216     z2 = MULTIPLY(tmp10, FIX_0_541196100) + z5; /* c2-c6 */
00217     z4 = MULTIPLY(tmp12, FIX_1_306562965) + z5; /* c2+c6 */
00218     z3 = MULTIPLY(tmp11, FIX_0_707106781); /* c4 */
00219 
00220     z11 = tmp7 + z3;        /* phase 5 */
00221     z13 = tmp7 - z3;
00222 
00223     dataptr[DCTSIZE*5] = z13 + z2; /* phase 6 */
00224     dataptr[DCTSIZE*3] = z13 - z2;
00225     dataptr[DCTSIZE*1] = z11 + z4;
00226     dataptr[DCTSIZE*7] = z11 - z4;
00227 
00228     dataptr++;          /* advance pointer to next column */
00229   }
00230 }
00231 
00232 #endif /* DCT_IFAST_SUPPORTED */