CentrED/Imaging/JpegLib/imjfdctfst.pas

238 lines
7.7 KiB
Plaintext
Raw Normal View History

2022-05-08 10:47:53 +02:00
unit imjfdctfst;
{ This file contains a fast, not so accurate integer implementation of the
forward DCT (Discrete Cosine Transform).
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 fixed-point math,
accuracy is lost due to imprecise representation of the scaled
quantization values. The smaller the quantization table entry, the less
precise the scaled value, so this implementation does worse with high-
quality-setting files than with low-quality ones. }
{ Original: jfdctfst.c ; Copyright (C) 1994-1996, Thomas G. Lane. }
interface
{$I imjconfig.inc}
uses
imjmorecfg,
imjinclude,
imjpeglib,
imjdct; { Private declarations for DCT subsystem }
{ Perform the forward DCT on one block of samples. }
{GLOBAL}
procedure jpeg_fdct_ifast (var data : array of DCTELEM);
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}
{ Scaling decisions are generally the same as in the LL&M algorithm;
see jfdctint.c for more details. However, we choose to descale
(right shift) multiplication products as soon as they are formed,
rather than carrying additional fractional bits into subsequent additions.
This compromises accuracy slightly, but it lets us save a few shifts.
More importantly, 16-bit arithmetic is then adequate (for 8-bit samples)
everywhere except in the multiplications proper; this saves a good deal
of work on 16-bit-int machines.
Again to save a few shifts, the intermediate results between pass 1 and
pass 2 are not upscaled, but are represented only to integral precision.
A final compromise is to represent the multiplicative constants to only
8 fractional bits, rather than 13. This saves some shifting work on some
machines, and may also reduce the cost of multiplication (since there
are fewer one-bits in the constants). }
const
CONST_BITS = 8;
const
CONST_SCALE = (INT32(1) shl CONST_BITS);
const
FIX_0_382683433 = INT32(Round(CONST_SCALE * 0.382683433)); {98}
FIX_0_541196100 = INT32(Round(CONST_SCALE * 0.541196100)); {139}
FIX_0_707106781 = INT32(Round(CONST_SCALE * 0.707106781)); {181}
FIX_1_306562965 = INT32(Round(CONST_SCALE * 1.306562965)); {334}
{ Descale and correctly round an INT32 value that's scaled by N bits.
We assume RIGHT_SHIFT rounds towards minus infinity, so adding
the fudge factor is correct for either sign of X. }
function DESCALE(x : INT32; n : int) : INT32;
var
shift_temp : INT32;
begin
{ We can gain a little more speed, with a further compromise in accuracy,
by omitting the addition in a descaling shift. This yields an incorrectly
rounded result half the time... }
{$ifndef USE_ACCURATE_ROUNDING}
shift_temp := x;
{$else}
shift_temp := x + (INT32(1) shl (n-1));
{$endif}
{$ifdef RIGHT_SHIFT_IS_UNSIGNED}
if shift_temp < 0 then
Descale := (shift_temp shr n) or ((not INT32(0)) shl (32-n))
else
{$endif}
Descale := (shift_temp shr n);
end;
{ Multiply a DCTELEM variable by an INT32 constant, and immediately
descale to yield a DCTELEM result. }
function MULTIPLY(X : DCTELEM; Y: INT32): DCTELEM;
begin
Multiply := DeScale((X) * (Y), CONST_BITS);
end;
{ Perform the forward DCT on one block of samples. }
{GLOBAL}
procedure jpeg_fdct_ifast (var data : array of DCTELEM);
type
PWorkspace = ^TWorkspace;
TWorkspace = array [0..DCTSIZE2-1] of DCTELEM;
var
tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7 : DCTELEM;
tmp10, tmp11, tmp12, tmp13 : DCTELEM;
z1, z2, z3, z4, z5, z11, z13 : DCTELEM;
dataptr : PWorkspace;
ctr : int;
{SHIFT_TEMPS}
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 := MULTIPLY(tmp12 + tmp13, FIX_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 := MULTIPLY(tmp10 - tmp12, FIX_0_382683433); { c6 }
z2 := MULTIPLY(tmp10, FIX_0_541196100) + z5; { c2-c6 }
z4 := MULTIPLY(tmp12, FIX_1_306562965) + z5; { c2+c6 }
z3 := MULTIPLY(tmp11, FIX_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(DCTELEMPTR(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 := MULTIPLY(tmp12 + tmp13, FIX_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 := MULTIPLY(tmp10 - tmp12, FIX_0_382683433); { c6 }
z2 := MULTIPLY(tmp10, FIX_0_541196100) + z5; { c2-c6 }
z4 := MULTIPLY(tmp12, FIX_1_306562965) + z5; { c2+c6 }
z3 := MULTIPLY(tmp11, FIX_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(DCTELEMPTR(dataptr)); { advance pointer to next column }
end;
end;
end.