1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
|
Unit JFDctFlt;
{$N+}
{ This file contains a floating-point implementation of the
forward DCT (Discrete Cosine Transform).
This implementation should be more accurate than either of the integer
DCT implementations. However, it may not give the same results on all
machines because of differences in roundoff behavior. Speed will depend
on the hardware's floating point capacity.
A 2-D DCT can be done by 1-D DCT on each row followed by 1-D DCT
on each column. Direct algorithms are also available, but they are
much more complex and seem not to be any faster when reduced to code.
This implementation is based on Arai, Agui, and Nakajima's algorithm for
scaled DCT. Their original paper (Trans. IEICE E-71(11):1095) is in
Japanese, but the algorithm is described in the Pennebaker & Mitchell
JPEG textbook (see REFERENCES section in file README). The following code
is based directly on figure 4-8 in P&M.
While an 8-point DCT cannot be done in less than 11 multiplies, it is
possible to arrange the computation so that many of the multiplies are
simple scalings of the final outputs. These multiplies can then be
folded into the multiplications or divisions by the JPEG quantization
table entries. The AA&N method leaves only 5 multiplies and 29 adds
to be done in the DCT itself.
The primary disadvantage of this method is that with a fixed-point
implementation, accuracy is lost due to imprecise representation of the
scaled quantization values. However, that problem does not arise if
we use floating point arithmetic. }
{ Original : jfdctflt.c ; Copyright (C) 1994-1996, Thomas G. Lane. }
interface
{$I jconfig.inc}
uses
jmorecfg,
jinclude,
jpeglib,
jdct; { Private declarations for DCT subsystem }
{ Perform the forward DCT on one block of samples.}
{GLOBAL}
procedure jpeg_fdct_float (var data : array of FAST_FLOAT);
implementation
{ This module is specialized to the case DCTSIZE = 8. }
{$ifndef DCTSIZE_IS_8}
Sorry, this code only copes with 8x8 DCTs. { deliberate syntax err }
{$endif}
{ Perform the forward DCT on one block of samples.}
{GLOBAL}
procedure jpeg_fdct_float (var data : array of FAST_FLOAT);
type
PWorkspace = ^TWorkspace;
TWorkspace = array [0..DCTSIZE2-1] of FAST_FLOAT;
var
tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7 : FAST_FLOAT;
tmp10, tmp11, tmp12, tmp13 : FAST_FLOAT;
z1, z2, z3, z4, z5, z11, z13 : FAST_FLOAT;
dataptr : PWorkspace;
ctr : int;
begin
{ Pass 1: process rows. }
dataptr := PWorkspace(@data);
for ctr := DCTSIZE-1 downto 0 do
begin
tmp0 := dataptr^[0] + dataptr^[7];
tmp7 := dataptr^[0] - dataptr^[7];
tmp1 := dataptr^[1] + dataptr^[6];
tmp6 := dataptr^[1] - dataptr^[6];
tmp2 := dataptr^[2] + dataptr^[5];
tmp5 := dataptr^[2] - dataptr^[5];
tmp3 := dataptr^[3] + dataptr^[4];
tmp4 := dataptr^[3] - dataptr^[4];
{ Even part }
tmp10 := tmp0 + tmp3; { phase 2 }
tmp13 := tmp0 - tmp3;
tmp11 := tmp1 + tmp2;
tmp12 := tmp1 - tmp2;
dataptr^[0] := tmp10 + tmp11; { phase 3 }
dataptr^[4] := tmp10 - tmp11;
z1 := (tmp12 + tmp13) * ({FAST_FLOAT}(0.707106781)); { c4 }
dataptr^[2] := tmp13 + z1; { phase 5 }
dataptr^[6] := tmp13 - z1;
{ Odd part }
tmp10 := tmp4 + tmp5; { phase 2 }
tmp11 := tmp5 + tmp6;
tmp12 := tmp6 + tmp7;
{ The rotator is modified from fig 4-8 to avoid extra negations. }
z5 := (tmp10 - tmp12) * ( {FAST_FLOAT}(0.382683433)); { c6 }
z2 := {FAST_FLOAT}(0.541196100) * tmp10 + z5; { c2-c6 }
z4 := {FAST_FLOAT}(1.306562965) * tmp12 + z5; { c2+c6 }
z3 := tmp11 * {FAST_FLOAT} (0.707106781); { c4 }
z11 := tmp7 + z3; { phase 5 }
z13 := tmp7 - z3;
dataptr^[5] := z13 + z2; { phase 6 }
dataptr^[3] := z13 - z2;
dataptr^[1] := z11 + z4;
dataptr^[7] := z11 - z4;
Inc(FAST_FLOAT_PTR(dataptr), DCTSIZE); { advance pointer to next row }
end;
{ Pass 2: process columns. }
dataptr := PWorkspace(@data);
for ctr := DCTSIZE-1 downto 0 do
begin
tmp0 := dataptr^[DCTSIZE*0] + dataptr^[DCTSIZE*7];
tmp7 := dataptr^[DCTSIZE*0] - dataptr^[DCTSIZE*7];
tmp1 := dataptr^[DCTSIZE*1] + dataptr^[DCTSIZE*6];
tmp6 := dataptr^[DCTSIZE*1] - dataptr^[DCTSIZE*6];
tmp2 := dataptr^[DCTSIZE*2] + dataptr^[DCTSIZE*5];
tmp5 := dataptr^[DCTSIZE*2] - dataptr^[DCTSIZE*5];
tmp3 := dataptr^[DCTSIZE*3] + dataptr^[DCTSIZE*4];
tmp4 := dataptr^[DCTSIZE*3] - dataptr^[DCTSIZE*4];
{ Even part }
tmp10 := tmp0 + tmp3; { phase 2 }
tmp13 := tmp0 - tmp3;
tmp11 := tmp1 + tmp2;
tmp12 := tmp1 - tmp2;
dataptr^[DCTSIZE*0] := tmp10 + tmp11; { phase 3 }
dataptr^[DCTSIZE*4] := tmp10 - tmp11;
z1 := (tmp12 + tmp13) * {FAST_FLOAT} (0.707106781); { c4 }
dataptr^[DCTSIZE*2] := tmp13 + z1; { phase 5 }
dataptr^[DCTSIZE*6] := tmp13 - z1;
{ Odd part }
tmp10 := tmp4 + tmp5; { phase 2 }
tmp11 := tmp5 + tmp6;
tmp12 := tmp6 + tmp7;
{ The rotator is modified from fig 4-8 to avoid extra negations. }
z5 := (tmp10 - tmp12) * {FAST_FLOAT} (0.382683433); { c6 }
z2 := {FAST_FLOAT} (0.541196100) * tmp10 + z5; { c2-c6 }
z4 := {FAST_FLOAT} (1.306562965) * tmp12 + z5; { c2+c6 }
z3 := tmp11 * {FAST_FLOAT} (0.707106781); { c4 }
z11 := tmp7 + z3; { phase 5 }
z13 := tmp7 - z3;
dataptr^[DCTSIZE*5] := z13 + z2; { phase 6 }
dataptr^[DCTSIZE*3] := z13 - z2;
dataptr^[DCTSIZE*1] := z11 + z4;
dataptr^[DCTSIZE*7] := z11 - z4;
Inc(FAST_FLOAT_PTR(dataptr)); { advance pointer to next column }
end;
end;
end.
|
