+static
[fjl] / idct.c
1 #include <math.h>
2
3 #include "idct.h"
4
5 void idct_reference(const int16_t* input, const uint32_t* quant_table, uint8_t* output)
6 {
7         double temp[DCTSIZE2];
8
9         for (unsigned y = 0; y < 8; ++y) {
10                 for (unsigned x = 0; x < 8; ++x) {
11                         double acc = 0.0;
12                         for (unsigned u = 0; u < 8; ++u) {
13                                 for (unsigned v = 0; v < 8; ++v) {
14                                         double c_u = (u == 0) ? 1/sqrt(2.0) : 1.0;
15                                         double c_v = (v == 0) ? 1/sqrt(2.0) : 1.0;
16                                         acc += c_u * c_v
17                                                 * input[u * DCTSIZE + v] * quant_table[u * DCTSIZE + v]
18                                                 * cos((2 * x + 1) * v * M_PI / 16.0)
19                                                 * cos((2 * y + 1) * u * M_PI / 16.0);
20                                 }
21                         }
22                         temp[y * DCTSIZE + x] = 0.25 * acc;
23                 }
24         }
25
26         for (unsigned y = 0; y < 8; ++y) {
27                 for (unsigned x = 0; x < 8; ++x) {
28                         double val = temp[y * DCTSIZE + x];
29                         if (val < 0.0) {
30                                 output[y * DCTSIZE + x] = 0;
31                         } else if (val >= 255.0) {
32                                 output[y * DCTSIZE + x] = 255;
33                         } else {
34                                 output[y * DCTSIZE + x] = (uint8_t)(val + 0.5);
35                         }
36                 }
37         }
38 }
39
40 // AA&N (Arai, Agui and Nakajima) floating-point IDCT.
41 // This IDCT is based on the same DCT that libjpeg uses -- in fact, exactly the
42 // same figure from the same book ("JPEG: Still Image Data Compression Standard",
43 // page 52, figure 4-8). However, it is coded from scratch, and uses the
44 // transposition method for converting DCT -> IDCT suggested in the book.
45 // (libjpeg seems to use some other method that yields similar, but not
46 // the same, code.) 
47
48 // As this is generally meant as a reference and not useful code (we expect
49 // a SIMD fixed-point algorithm to be used in most cases), it has not been
50 // attempted significantly optimized. We assume the compiler will be smart
51 // enough to do all the variable propagation for us anyway.
52
53 // Scale factors; 1.0 / (sqrt(2.0) * cos(k * M_PI / 16.0)), except for the first which is 1.
54 static const double scalefac[] = {
55         1.0, 0.7209598220069479, 0.765366864730180, 0.8504300947672564,
56         1.0, 1.2727585805728336, 1.847759065022573, 3.6245097854115502
57 };
58
59 // 1D 8-point DCT.
60 static inline void idct1d_float(double y0, double y1, double y2, double y3, double y4, double y5, double y6, double y7, double *x)
61 {
62         // constants
63         static const double a1 = 0.7071067811865474;   // sqrt(2)
64         static const double a2 = 0.5411961001461971;   // cos(3/8 pi) * sqrt(2)
65         static const double a3 = a1;
66         static const double a4 = 1.3065629648763766;   // cos(pi/8) * sqrt(2)
67         static const double a5 = 0.5 * (a4 - a2);
68
69         // phase 1
70         const double p1_0 = y0 * scalefac[0];
71         const double p1_1 = y4 * scalefac[4];
72         const double p1_2 = y2 * scalefac[2];
73         const double p1_3 = y6 * scalefac[6];
74         const double p1_4 = y5 * scalefac[5];
75         const double p1_5 = y1 * scalefac[1];
76         const double p1_6 = y7 * scalefac[7];
77         const double p1_7 = y3 * scalefac[3];
78
79         // phase 2
80         const double p2_0 = p1_0;
81         const double p2_1 = p1_1;
82         const double p2_2 = p1_2;
83         const double p2_3 = p1_3;
84         const double p2_4 = p1_4 - p1_7;
85         const double p2_5 = p1_5 + p1_6;
86         const double p2_6 = p1_5 - p1_6;
87         const double p2_7 = p1_4 + p1_7;
88
89         // phase 3
90         const double p3_0 = p2_0;
91         const double p3_1 = p2_1;
92         const double p3_2 = p2_2 - p2_3;
93         const double p3_3 = p2_2 + p2_3;
94         const double p3_4 = p2_4;
95         const double p3_5 = p2_5 - p2_7;
96         const double p3_6 = p2_6;
97         const double p3_7 = p2_5 + p2_7;
98         
99         // phase 4
100         const double p4_0 = p3_0;
101         const double p4_1 = p3_1;
102         const double p4_2 = a1 * p3_2;
103         const double p4_3 = p3_3;
104         const double p4_4 = p3_4 * -a2 + (p3_4 + p3_6) * -a5;
105         const double p4_5 = a3 * p3_5;
106         const double p4_6 = p3_6 * a4 + (p3_4 + p3_6) * -a5;
107         const double p4_7 = p3_7;
108
109         // phase 5
110         const double p5_0 = p4_0 + p4_1;
111         const double p5_1 = p4_0 - p4_1;
112         const double p5_2 = p4_2;
113         const double p5_3 = p4_2 + p4_3;
114         const double p5_4 = p4_4;
115         const double p5_5 = p4_5;
116         const double p5_6 = p4_6;
117         const double p5_7 = p4_7;
118
119         // phase 6
120         const double p6_0 = p5_0 + p5_3;
121         const double p6_1 = p5_1 + p5_2;
122         const double p6_2 = p5_1 - p5_2;
123         const double p6_3 = p5_0 - p5_3;
124         const double p6_4 = -p5_4;
125         const double p6_5 = p5_5 - p5_4;
126         const double p6_6 = p5_5 + p5_6;
127         const double p6_7 = p5_6 + p5_7;
128
129         // phase 7
130         x[0] = p6_0 + p6_7;
131         x[1] = p6_1 + p6_6;
132         x[2] = p6_2 + p6_5;
133         x[3] = p6_3 + p6_4;
134         x[4] = p6_3 - p6_4;
135         x[5] = p6_2 - p6_5;
136         x[6] = p6_1 - p6_6;
137         x[7] = p6_0 - p6_7;
138 }
139
140 void idct_float(const int16_t* input, const uint32_t* quant_table, uint8_t* output)
141 {
142         double temp[DCTSIZE2];
143
144         // IDCT columns.
145         for (unsigned x = 0; x < DCTSIZE; ++x) {
146                 idct1d_float(input[DCTSIZE * 0 + x] * (int32_t)quant_table[DCTSIZE * 0 + x],
147                              input[DCTSIZE * 1 + x] * (int32_t)quant_table[DCTSIZE * 1 + x],
148                              input[DCTSIZE * 2 + x] * (int32_t)quant_table[DCTSIZE * 2 + x],
149                              input[DCTSIZE * 3 + x] * (int32_t)quant_table[DCTSIZE * 3 + x],
150                              input[DCTSIZE * 4 + x] * (int32_t)quant_table[DCTSIZE * 4 + x],
151                              input[DCTSIZE * 5 + x] * (int32_t)quant_table[DCTSIZE * 5 + x],
152                              input[DCTSIZE * 6 + x] * (int32_t)quant_table[DCTSIZE * 6 + x],
153                              input[DCTSIZE * 7 + x] * (int32_t)quant_table[DCTSIZE * 7 + x],
154                              temp + x * DCTSIZE);
155         }
156         
157         // IDCT rows.
158         for (unsigned y = 0; y < DCTSIZE; ++y) {
159                 double temp2[DCTSIZE];
160                 idct1d_float(temp[DCTSIZE * 0 + y],
161                              temp[DCTSIZE * 1 + y],
162                              temp[DCTSIZE * 2 + y],
163                              temp[DCTSIZE * 3 + y],
164                              temp[DCTSIZE * 4 + y],
165                              temp[DCTSIZE * 5 + y],
166                              temp[DCTSIZE * 6 + y],
167                              temp[DCTSIZE * 7 + y],
168                              temp2);
169                 for (unsigned x = 0; x < DCTSIZE; ++x) {
170                         double val = (1.0/8.0) * temp2[x];
171                         if (val < 0.0) {
172                                 output[y * DCTSIZE + x] = 0;
173                         } else if (val >= 255.0) {
174                                 output[y * DCTSIZE + x] = 255;
175                         } else {
176                                 output[y * DCTSIZE + x] = (uint8_t)(val + 0.5);
177                         }
178                 }
179         }
180 }