Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 1 | /* |
| 2 | * jfwddct.c |
| 3 | * |
| 4 | * Copyright (C) 1991, 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 basic DCT (Discrete Cosine Transform) |
| 9 | * transformation subroutine. |
| 10 | * |
| 11 | * This implementation is based on Appendix A.2 of the book |
| 12 | * "Discrete Cosine Transform---Algorithms, Advantages, Applications" |
| 13 | * by K.R. Rao and P. Yip (Academic Press, Inc, London, 1990). |
| 14 | * It uses scaled fixed-point arithmetic instead of floating point. |
| 15 | */ |
| 16 | |
| 17 | #include "jinclude.h" |
| 18 | |
| 19 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 20 | /* We assume that right shift corresponds to signed division by 2 with |
| 21 | * rounding towards minus infinity. This is correct for typical "arithmetic |
| 22 | * shift" instructions that shift in copies of the sign bit. But some |
| 23 | * C compilers implement >> with an unsigned shift. For these machines you |
| 24 | * must define RIGHT_SHIFT_IS_UNSIGNED. |
| 25 | * RIGHT_SHIFT provides a signed right shift of an INT32 quantity. |
| 26 | * It is only applied with constant shift counts. |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 27 | */ |
| 28 | |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 29 | #ifdef RIGHT_SHIFT_IS_UNSIGNED |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 30 | #define SHIFT_TEMPS INT32 shift_temp; |
| 31 | #define RIGHT_SHIFT(x,shft) \ |
| 32 | ((shift_temp = (x)) < 0 ? \ |
| 33 | (shift_temp >> (shft)) | ((~0) << (32-(shft))) : \ |
| 34 | (shift_temp >> (shft))) |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 35 | #else |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 36 | #define SHIFT_TEMPS |
| 37 | #define RIGHT_SHIFT(x,shft) ((x) >> (shft)) |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 38 | #endif |
| 39 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 40 | |
| 41 | /* The poop on this scaling stuff is as follows: |
| 42 | * |
| 43 | * We have to do addition and subtraction of the integer inputs, which |
| 44 | * is no problem, and multiplication by fractional constants, which is |
| 45 | * a problem to do in integer arithmetic. We multiply all the constants |
| 46 | * by DCT_SCALE and convert them to integer constants (thus retaining |
| 47 | * LG2_DCT_SCALE bits of precision in the constants). After doing a |
| 48 | * multiplication we have to divide the product by DCT_SCALE, with proper |
| 49 | * rounding, to produce the correct output. The division can be implemented |
| 50 | * cheaply as a right shift of LG2_DCT_SCALE bits. The DCT equations also |
| 51 | * specify an additional division by 2 on the final outputs; this can be |
| 52 | * folded into the right-shift by shifting one more bit (see UNFIXH). |
| 53 | * |
| 54 | * If you are planning to recode this in assembler, you might want to set |
| 55 | * LG2_DCT_SCALE to 15. This loses a bit of precision, but then all the |
| 56 | * multiplications are between 16-bit quantities (given 8-bit JSAMPLEs!) |
| 57 | * so you could use a signed 16x16=>32 bit multiply instruction instead of |
| 58 | * full 32x32 multiply. Unfortunately there's no way to describe such a |
| 59 | * multiply portably in C, so we've gone for the extra bit of accuracy here. |
| 60 | */ |
| 61 | |
| 62 | #ifdef EIGHT_BIT_SAMPLES |
| 63 | #define LG2_DCT_SCALE 16 |
| 64 | #else |
| 65 | #define LG2_DCT_SCALE 15 /* lose a little precision to avoid overflow */ |
| 66 | #endif |
| 67 | |
| 68 | #define ONE ((INT32) 1) |
| 69 | |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 70 | #define DCT_SCALE (ONE << LG2_DCT_SCALE) |
| 71 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 72 | /* In some places we shift the inputs left by a couple more bits, */ |
| 73 | /* so that they can be added to fractional results without too much */ |
| 74 | /* loss of precision. */ |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 75 | #define LG2_OVERSCALE 2 |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 76 | #define OVERSCALE (ONE << LG2_OVERSCALE) |
| 77 | #define OVERSHIFT(x) ((x) <<= LG2_OVERSCALE) |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 78 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 79 | /* Scale a fractional constant by DCT_SCALE */ |
| 80 | #define FIX(x) ((INT32) ((x) * DCT_SCALE + 0.5)) |
| 81 | |
| 82 | /* Scale a fractional constant by DCT_SCALE/OVERSCALE */ |
| 83 | /* Such a constant can be multiplied with an overscaled input */ |
| 84 | /* to produce something that's scaled by DCT_SCALE */ |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 85 | #define FIXO(x) ((INT32) ((x) * DCT_SCALE / OVERSCALE + 0.5)) |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 86 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 87 | /* Descale and correctly round a value that's scaled by DCT_SCALE */ |
| 88 | #define UNFIX(x) RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1)), LG2_DCT_SCALE) |
| 89 | |
| 90 | /* Same with an additional division by 2, ie, correctly rounded UNFIX(x/2) */ |
| 91 | #define UNFIXH(x) RIGHT_SHIFT((x) + (ONE << LG2_DCT_SCALE), LG2_DCT_SCALE+1) |
| 92 | |
| 93 | /* Take a value scaled by DCT_SCALE and round to integer scaled by OVERSCALE */ |
| 94 | #define UNFIXO(x) RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1-LG2_OVERSCALE)),\ |
| 95 | LG2_DCT_SCALE-LG2_OVERSCALE) |
| 96 | |
| 97 | /* Here are the constants we need */ |
| 98 | /* SIN_i_j is sine of i*pi/j, scaled by DCT_SCALE */ |
| 99 | /* COS_i_j is cosine of i*pi/j, scaled by DCT_SCALE */ |
| 100 | |
| 101 | #define SIN_1_4 FIX(0.707106781) |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 102 | #define COS_1_4 SIN_1_4 |
| 103 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 104 | #define SIN_1_8 FIX(0.382683432) |
| 105 | #define COS_1_8 FIX(0.923879533) |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 106 | #define SIN_3_8 COS_1_8 |
| 107 | #define COS_3_8 SIN_1_8 |
| 108 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 109 | #define SIN_1_16 FIX(0.195090322) |
| 110 | #define COS_1_16 FIX(0.980785280) |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 111 | #define SIN_7_16 COS_1_16 |
| 112 | #define COS_7_16 SIN_1_16 |
| 113 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 114 | #define SIN_3_16 FIX(0.555570233) |
| 115 | #define COS_3_16 FIX(0.831469612) |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 116 | #define SIN_5_16 COS_3_16 |
| 117 | #define COS_5_16 SIN_3_16 |
| 118 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 119 | /* OSIN_i_j is sine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */ |
| 120 | /* OCOS_i_j is cosine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */ |
| 121 | |
| 122 | #define OSIN_1_4 FIXO(0.707106781) |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 123 | #define OCOS_1_4 OSIN_1_4 |
| 124 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 125 | #define OSIN_1_8 FIXO(0.382683432) |
| 126 | #define OCOS_1_8 FIXO(0.923879533) |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 127 | #define OSIN_3_8 OCOS_1_8 |
| 128 | #define OCOS_3_8 OSIN_1_8 |
| 129 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 130 | #define OSIN_1_16 FIXO(0.195090322) |
| 131 | #define OCOS_1_16 FIXO(0.980785280) |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 132 | #define OSIN_7_16 OCOS_1_16 |
| 133 | #define OCOS_7_16 OSIN_1_16 |
| 134 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 135 | #define OSIN_3_16 FIXO(0.555570233) |
| 136 | #define OCOS_3_16 FIXO(0.831469612) |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 137 | #define OSIN_5_16 OCOS_3_16 |
| 138 | #define OCOS_5_16 OSIN_3_16 |
| 139 | |
| 140 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 141 | /* |
| 142 | * Perform a 1-dimensional DCT. |
| 143 | * Note that this code is specialized to the case DCTSIZE = 8. |
| 144 | */ |
| 145 | |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 146 | INLINE |
| 147 | LOCAL void |
| 148 | fast_dct_8 (DCTELEM *in, int stride) |
| 149 | { |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 150 | /* many tmps have nonoverlapping lifetime -- flashy register colourers |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 151 | * should be able to do this lot very well |
| 152 | */ |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 153 | INT32 in0, in1, in2, in3, in4, in5, in6, in7; |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 154 | INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7; |
| 155 | INT32 tmp10, tmp11, tmp12, tmp13; |
| 156 | INT32 tmp14, tmp15, tmp16, tmp17; |
| 157 | INT32 tmp25, tmp26; |
| 158 | SHIFT_TEMPS |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 159 | |
| 160 | in0 = in[ 0]; |
| 161 | in1 = in[stride ]; |
| 162 | in2 = in[stride*2]; |
| 163 | in3 = in[stride*3]; |
| 164 | in4 = in[stride*4]; |
| 165 | in5 = in[stride*5]; |
| 166 | in6 = in[stride*6]; |
| 167 | in7 = in[stride*7]; |
| 168 | |
| 169 | tmp0 = in7 + in0; |
| 170 | tmp1 = in6 + in1; |
| 171 | tmp2 = in5 + in2; |
| 172 | tmp3 = in4 + in3; |
| 173 | tmp4 = in3 - in4; |
| 174 | tmp5 = in2 - in5; |
| 175 | tmp6 = in1 - in6; |
| 176 | tmp7 = in0 - in7; |
| 177 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 178 | tmp10 = tmp3 + tmp0; |
| 179 | tmp11 = tmp2 + tmp1; |
| 180 | tmp12 = tmp1 - tmp2; |
| 181 | tmp13 = tmp0 - tmp3; |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 182 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 183 | in[ 0] = (DCTELEM) UNFIXH((tmp10 + tmp11) * SIN_1_4); |
| 184 | in[stride*4] = (DCTELEM) UNFIXH((tmp10 - tmp11) * COS_1_4); |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 185 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 186 | in[stride*2] = (DCTELEM) UNFIXH(tmp13*COS_1_8 + tmp12*SIN_1_8); |
| 187 | in[stride*6] = (DCTELEM) UNFIXH(tmp13*SIN_1_8 - tmp12*COS_1_8); |
| 188 | |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 189 | tmp16 = UNFIXO((tmp6 + tmp5) * SIN_1_4); |
| 190 | tmp15 = UNFIXO((tmp6 - tmp5) * COS_1_4); |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 191 | |
| 192 | OVERSHIFT(tmp4); |
| 193 | OVERSHIFT(tmp7); |
| 194 | |
| 195 | /* tmp4, tmp7, tmp15, tmp16 are overscaled by OVERSCALE */ |
| 196 | |
| 197 | tmp14 = tmp4 + tmp15; |
| 198 | tmp25 = tmp4 - tmp15; |
| 199 | tmp26 = tmp7 - tmp16; |
| 200 | tmp17 = tmp7 + tmp16; |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 201 | |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 202 | in[stride ] = (DCTELEM) UNFIXH(tmp17*OCOS_1_16 + tmp14*OSIN_1_16); |
| 203 | in[stride*7] = (DCTELEM) UNFIXH(tmp17*OCOS_7_16 - tmp14*OSIN_7_16); |
| 204 | in[stride*5] = (DCTELEM) UNFIXH(tmp26*OCOS_5_16 + tmp25*OSIN_5_16); |
| 205 | in[stride*3] = (DCTELEM) UNFIXH(tmp26*OCOS_3_16 - tmp25*OSIN_3_16); |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 206 | } |
| 207 | |
| 208 | |
| 209 | /* |
| 210 | * Perform the forward DCT on one block of samples. |
| 211 | * |
Thomas G. Lane | bd543f0 | 1991-12-13 00:00:00 +0000 | [diff] [blame^] | 212 | * A 2-D DCT can be done by 1-D DCT on each row |
| 213 | * followed by 1-D DCT on each column. |
Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame] | 214 | */ |
| 215 | |
| 216 | GLOBAL void |
| 217 | j_fwd_dct (DCTBLOCK data) |
| 218 | { |
| 219 | int i; |
| 220 | |
| 221 | for (i = 0; i < DCTSIZE; i++) |
| 222 | fast_dct_8(data+i*DCTSIZE, 1); |
| 223 | |
| 224 | for (i = 0; i < DCTSIZE; i++) |
| 225 | fast_dct_8(data+i, DCTSIZE); |
| 226 | } |