Coverage Report

Created: 2023-12-08 06:53

/src/freeimage-svn/FreeImage/trunk/Source/LibJPEG/jidctflt.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
 * jidctflt.c
3
 *
4
 * Copyright (C) 1994-1998, Thomas G. Lane.
5
 * Modified 2010-2017 by Guido Vollbeding.
6
 * This file is part of the Independent JPEG Group's software.
7
 * For conditions of distribution and use, see the accompanying README file.
8
 *
9
 * This file contains a floating-point implementation of the
10
 * inverse DCT (Discrete Cosine Transform).  In the IJG code, this routine
11
 * must also perform dequantization of the input coefficients.
12
 *
13
 * This implementation should be more accurate than either of the integer
14
 * IDCT implementations.  However, it may not give the same results on all
15
 * machines because of differences in roundoff behavior.  Speed will depend
16
 * on the hardware's floating point capacity.
17
 *
18
 * A 2-D IDCT can be done by 1-D IDCT on each column followed by 1-D IDCT
19
 * on each row (or vice versa, but it's more convenient to emit a row at
20
 * a time).  Direct algorithms are also available, but they are much more
21
 * complex and seem not to be any faster when reduced to code.
22
 *
23
 * This implementation is based on Arai, Agui, and Nakajima's algorithm for
24
 * scaled DCT.  Their original paper (Trans. IEICE E-71(11):1095) is in
25
 * Japanese, but the algorithm is described in the Pennebaker & Mitchell
26
 * JPEG textbook (see REFERENCES section in file README).  The following code
27
 * is based directly on figure 4-8 in P&M.
28
 * While an 8-point DCT cannot be done in less than 11 multiplies, it is
29
 * possible to arrange the computation so that many of the multiplies are
30
 * simple scalings of the final outputs.  These multiplies can then be
31
 * folded into the multiplications or divisions by the JPEG quantization
32
 * table entries.  The AA&N method leaves only 5 multiplies and 29 adds
33
 * to be done in the DCT itself.
34
 * The primary disadvantage of this method is that with a fixed-point
35
 * implementation, accuracy is lost due to imprecise representation of the
36
 * scaled quantization values.  However, that problem does not arise if
37
 * we use floating point arithmetic.
38
 */
39
40
#define JPEG_INTERNALS
41
#include "jinclude.h"
42
#include "jpeglib.h"
43
#include "jdct.h"   /* Private declarations for DCT subsystem */
44
45
#ifdef DCT_FLOAT_SUPPORTED
46
47
48
/*
49
 * This module is specialized to the case DCTSIZE = 8.
50
 */
51
52
#if DCTSIZE != 8
53
  Sorry, this code only copes with 8x8 DCT blocks. /* deliberate syntax err */
54
#endif
55
56
57
/* Dequantize a coefficient by multiplying it by the multiplier-table
58
 * entry; produce a float result.
59
 */
60
61
0
#define DEQUANTIZE(coef,quantval)  (((FAST_FLOAT) (coef)) * (quantval))
62
63
64
/*
65
 * Perform dequantization and inverse DCT on one block of coefficients.
66
 *
67
 * cK represents cos(K*pi/16).
68
 */
69
70
GLOBAL(void)
71
jpeg_idct_float (j_decompress_ptr cinfo, jpeg_component_info * compptr,
72
     JCOEFPTR coef_block,
73
     JSAMPARRAY output_buf, JDIMENSION output_col)
74
0
{
75
0
  FAST_FLOAT tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
76
0
  FAST_FLOAT tmp10, tmp11, tmp12, tmp13;
77
0
  FAST_FLOAT z5, z10, z11, z12, z13;
78
0
  JCOEFPTR inptr;
79
0
  FLOAT_MULT_TYPE * quantptr;
80
0
  FAST_FLOAT * wsptr;
81
0
  JSAMPROW outptr;
82
0
  JSAMPLE *range_limit = IDCT_range_limit(cinfo);
83
0
  int ctr;
84
0
  FAST_FLOAT workspace[DCTSIZE2]; /* buffers data between passes */
85
86
  /* Pass 1: process columns from input, store into work array. */
87
88
0
  inptr = coef_block;
89
0
  quantptr = (FLOAT_MULT_TYPE *) compptr->dct_table;
90
0
  wsptr = workspace;
91
0
  for (ctr = DCTSIZE; ctr > 0; ctr--) {
92
    /* Due to quantization, we will usually find that many of the input
93
     * coefficients are zero, especially the AC terms.  We can exploit this
94
     * by short-circuiting the IDCT calculation for any column in which all
95
     * the AC terms are zero.  In that case each output is equal to the
96
     * DC coefficient (with scale factor as needed).
97
     * With typical images and quantization tables, half or more of the
98
     * column DCT calculations can be simplified this way.
99
     */
100
101
0
    if (inptr[DCTSIZE*1] == 0 && inptr[DCTSIZE*2] == 0 &&
102
0
  inptr[DCTSIZE*3] == 0 && inptr[DCTSIZE*4] == 0 &&
103
0
  inptr[DCTSIZE*5] == 0 && inptr[DCTSIZE*6] == 0 &&
104
0
  inptr[DCTSIZE*7] == 0) {
105
      /* AC terms all zero */
106
0
      FAST_FLOAT dcval = DEQUANTIZE(inptr[DCTSIZE*0], quantptr[DCTSIZE*0]);
107
108
0
      wsptr[DCTSIZE*0] = dcval;
109
0
      wsptr[DCTSIZE*1] = dcval;
110
0
      wsptr[DCTSIZE*2] = dcval;
111
0
      wsptr[DCTSIZE*3] = dcval;
112
0
      wsptr[DCTSIZE*4] = dcval;
113
0
      wsptr[DCTSIZE*5] = dcval;
114
0
      wsptr[DCTSIZE*6] = dcval;
115
0
      wsptr[DCTSIZE*7] = dcval;
116
117
0
      inptr++;      /* advance pointers to next column */
118
0
      quantptr++;
119
0
      wsptr++;
120
0
      continue;
121
0
    }
122
123
    /* Even part */
124
125
0
    tmp0 = DEQUANTIZE(inptr[DCTSIZE*0], quantptr[DCTSIZE*0]);
126
0
    tmp1 = DEQUANTIZE(inptr[DCTSIZE*2], quantptr[DCTSIZE*2]);
127
0
    tmp2 = DEQUANTIZE(inptr[DCTSIZE*4], quantptr[DCTSIZE*4]);
128
0
    tmp3 = DEQUANTIZE(inptr[DCTSIZE*6], quantptr[DCTSIZE*6]);
129
130
0
    tmp10 = tmp0 + tmp2;  /* phase 3 */
131
0
    tmp11 = tmp0 - tmp2;
132
133
0
    tmp13 = tmp1 + tmp3;  /* phases 5-3 */
134
0
    tmp12 = (tmp1 - tmp3) * ((FAST_FLOAT) 1.414213562) - tmp13; /* 2*c4 */
135
136
0
    tmp0 = tmp10 + tmp13; /* phase 2 */
137
0
    tmp3 = tmp10 - tmp13;
138
0
    tmp1 = tmp11 + tmp12;
139
0
    tmp2 = tmp11 - tmp12;
140
141
    /* Odd part */
142
143
0
    tmp4 = DEQUANTIZE(inptr[DCTSIZE*1], quantptr[DCTSIZE*1]);
144
0
    tmp5 = DEQUANTIZE(inptr[DCTSIZE*3], quantptr[DCTSIZE*3]);
145
0
    tmp6 = DEQUANTIZE(inptr[DCTSIZE*5], quantptr[DCTSIZE*5]);
146
0
    tmp7 = DEQUANTIZE(inptr[DCTSIZE*7], quantptr[DCTSIZE*7]);
147
148
0
    z13 = tmp6 + tmp5;    /* phase 6 */
149
0
    z10 = tmp6 - tmp5;
150
0
    z11 = tmp4 + tmp7;
151
0
    z12 = tmp4 - tmp7;
152
153
0
    tmp7 = z11 + z13;   /* phase 5 */
154
0
    tmp11 = (z11 - z13) * ((FAST_FLOAT) 1.414213562); /* 2*c4 */
155
156
0
    z5 = (z10 + z12) * ((FAST_FLOAT) 1.847759065); /* 2*c2 */
157
0
    tmp10 = z5 - z12 * ((FAST_FLOAT) 1.082392200); /* 2*(c2-c6) */
158
0
    tmp12 = z5 - z10 * ((FAST_FLOAT) 2.613125930); /* 2*(c2+c6) */
159
160
0
    tmp6 = tmp12 - tmp7;  /* phase 2 */
161
0
    tmp5 = tmp11 - tmp6;
162
0
    tmp4 = tmp10 - tmp5;
163
164
0
    wsptr[DCTSIZE*0] = tmp0 + tmp7;
165
0
    wsptr[DCTSIZE*7] = tmp0 - tmp7;
166
0
    wsptr[DCTSIZE*1] = tmp1 + tmp6;
167
0
    wsptr[DCTSIZE*6] = tmp1 - tmp6;
168
0
    wsptr[DCTSIZE*2] = tmp2 + tmp5;
169
0
    wsptr[DCTSIZE*5] = tmp2 - tmp5;
170
0
    wsptr[DCTSIZE*3] = tmp3 + tmp4;
171
0
    wsptr[DCTSIZE*4] = tmp3 - tmp4;
172
173
0
    inptr++;      /* advance pointers to next column */
174
0
    quantptr++;
175
0
    wsptr++;
176
0
  }
177
178
  /* Pass 2: process rows from work array, store into output array. */
179
180
0
  wsptr = workspace;
181
0
  for (ctr = 0; ctr < DCTSIZE; ctr++) {
182
0
    outptr = output_buf[ctr] + output_col;
183
    /* Rows of zeroes can be exploited in the same way as we did with columns.
184
     * However, the column calculation has created many nonzero AC terms, so
185
     * the simplification applies less often (typically 5% to 10% of the time).
186
     * And testing floats for zero is relatively expensive, so we don't bother.
187
     */
188
189
    /* Even part */
190
191
    /* Prepare range-limit and float->int conversion */
192
0
    z5 = wsptr[0] + (((FAST_FLOAT) RANGE_CENTER) + ((FAST_FLOAT) 0.5));
193
0
    tmp10 = z5 + wsptr[4];
194
0
    tmp11 = z5 - wsptr[4];
195
196
0
    tmp13 = wsptr[2] + wsptr[6];
197
0
    tmp12 = (wsptr[2] - wsptr[6]) *
198
0
        ((FAST_FLOAT) 1.414213562) - tmp13; /* 2*c4 */
199
200
0
    tmp0 = tmp10 + tmp13;
201
0
    tmp3 = tmp10 - tmp13;
202
0
    tmp1 = tmp11 + tmp12;
203
0
    tmp2 = tmp11 - tmp12;
204
205
    /* Odd part */
206
207
0
    z13 = wsptr[5] + wsptr[3];
208
0
    z10 = wsptr[5] - wsptr[3];
209
0
    z11 = wsptr[1] + wsptr[7];
210
0
    z12 = wsptr[1] - wsptr[7];
211
212
0
    tmp7 = z11 + z13;   /* phase 5 */
213
0
    tmp11 = (z11 - z13) * ((FAST_FLOAT) 1.414213562); /* 2*c4 */
214
215
0
    z5 = (z10 + z12) * ((FAST_FLOAT) 1.847759065); /* 2*c2 */
216
0
    tmp10 = z5 - z12 * ((FAST_FLOAT) 1.082392200); /* 2*(c2-c6) */
217
0
    tmp12 = z5 - z10 * ((FAST_FLOAT) 2.613125930); /* 2*(c2+c6) */
218
219
0
    tmp6 = tmp12 - tmp7;  /* phase 2 */
220
0
    tmp5 = tmp11 - tmp6;
221
0
    tmp4 = tmp10 - tmp5;
222
223
    /* Final output stage: float->int conversion and range-limit */
224
225
0
    outptr[0] = range_limit[(int) (tmp0 + tmp7) & RANGE_MASK];
226
0
    outptr[7] = range_limit[(int) (tmp0 - tmp7) & RANGE_MASK];
227
0
    outptr[1] = range_limit[(int) (tmp1 + tmp6) & RANGE_MASK];
228
0
    outptr[6] = range_limit[(int) (tmp1 - tmp6) & RANGE_MASK];
229
0
    outptr[2] = range_limit[(int) (tmp2 + tmp5) & RANGE_MASK];
230
0
    outptr[5] = range_limit[(int) (tmp2 - tmp5) & RANGE_MASK];
231
0
    outptr[3] = range_limit[(int) (tmp3 + tmp4) & RANGE_MASK];
232
0
    outptr[4] = range_limit[(int) (tmp3 - tmp4) & RANGE_MASK];
233
234
0
    wsptr += DCTSIZE;   /* advance pointer to next row */
235
0
  }
236
0
}
237
238
#endif /* DCT_FLOAT_SUPPORTED */