diff options
Diffstat (limited to 'packages/pasjpeg/src/jfdctfst.pas')
| -rw-r--r-- | packages/pasjpeg/src/jfdctfst.pas | 237 |
1 files changed, 237 insertions, 0 deletions
diff --git a/packages/pasjpeg/src/jfdctfst.pas b/packages/pasjpeg/src/jfdctfst.pas new file mode 100644 index 0000000000..faf4121bc7 --- /dev/null +++ b/packages/pasjpeg/src/jfdctfst.pas @@ -0,0 +1,237 @@ +Unit JFDctFst; + +{ 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 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_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. |