flat assembler
Message board for the users of flat assembler. Index > High Level Languages > understanding jpg decompression Goto page 1, 2  Next
Author
 Thread  vivik

Joined: 29 Oct 2016
Posts: 671
vivik
Here is the simplest demonstration of it: https://www.youtube.com/watch?v=Q2aEzeMDHMA

I want to understand the code that is used in libjpeg-turbo for decompressing jpg, specifically the code for inverse discrete cosine transform.

There are multiple versions of this code, but I found "float" version to be the simplest: https://github.com/libjpeg-turbo/libjpeg-turbo/blob/master/jidctflt.c

Quote:
* 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.ijg). 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.

I guess I need to find either that Arai, Agui, and Nakajima's paper, or that Pennebaker & Mitchell book somewhere to understand that code.

Other versions are "slow integer" (jidctint.c simd/jidctint-sse2.asm) and "fast integer" (jfdctfst.c simd/jidctfst-sse2.asm). I'd like to understand all of them. 04 May 2018, 19:27  donn

Joined: 05 Mar 2010
Posts: 214
donn
I'm less familiar with DCT, but I think it's similar to the Discrete Fourier Transform (and its inverse). You can represent a signal (such as sound) in the time domain, or spatially if you are dealing with a single two-dimensional image. The Fourier Transform or DCT I think, provides a means to represent that signal in terms of frequencies (the concept of frequency content).

Frequency, by its nature, involves repetitions, so compression can take advantage of these 'redundancies.' The transform involves coefficients which are these building blocks which can be 'decompressed' back to the spatial or time domain (decompressing a JPEG). Fourier thought you could break any signal (sequence of numbers) down into a sum of sines, which some mathematicians thought was impossible. Where the DCT uses cosines (I think in this same way), DFT uses sines and cosines, or analogous complex numbers using Euler's identities.

I have H. Joseph Weaver's "Applications of Discrete and Continuous Fourier Analysis" and it's on my to-do to implement a Cooley-Tukey FFT down the road. There are online calculators for some of this stuff, so if you feed in a series of numbers, you can verify the output. I've got Discrete Convolution/Faltung down so far (verified here), but this frequency stuff is very interesting.

Mathworks DCT link, showing the principal functions and descriptions of the terms. 04 May 2018, 20:12  pabloreda Joined: 24 Jan 2007 Posts: 105 Location: Argentina pabloreda I code in r4 a jpeg decompress, if you can read a forth like lang, this is the code https://github.com/phreda4/reda4/blob/master/r4/Lib/loadjpg.txt I use fixed point , the code in this file decompress a file for test and can be use like library. Is a translation of a c code. 04 May 2018, 23:40
donn

Joined: 05 Mar 2010
Posts: 214
donn
Wow, r4 is very mature, seems like it has many capabilities.

So the jpeg lib can leverage r4asm to generate fasm-style asm? 05 May 2018, 00:28  pabloreda Joined: 24 Jan 2007 Posts: 105 Location: Argentina pabloreda yes, hit f5 in the editor, and this generate the source code in r4asm/ folder and compile with fasm, the code generator is not very well, I need improve this! 05 May 2018, 03:21
vivik

Joined: 29 Oct 2016
Posts: 671
vivik
@donn
I know the main principle of compression, I need to understand the actual algorithm used in libjpeg. Why this one is used? Why this one is fastest? What are alternatives? How those numbers (0.707106781 0.382683433 1.414213562) are computed?

@pabloreda
I can't read that, unfortunately. Looks like algorithm is slightly different. 05 May 2018, 05:33  Siekmanski

Joined: 18 Jun 2006
Posts: 12
Location: The Netherlands
Siekmanski
sin( PI / 360 * 90 ) = 0.707106781
sin( PI / 360 * 45 ) = 0.382683433
sqrt(2.0) = 1.414213562 05 May 2018, 06:58  vivik

Joined: 29 Oct 2016
Posts: 671
vivik
wow 05 May 2018, 07:29  Siekmanski

Joined: 18 Jun 2006
Posts: 12
Location: The Netherlands
Siekmanski
If you don't want to use sqrt , you can calculate sqrt(2) as,

sin( PI / 360 * 90 ) * 2 = 1.414213562 05 May 2018, 08:56  vivik

Joined: 29 Oct 2016
Posts: 671
vivik
I found this "flowgraph" here:

https://github.com/prtsh/aan_dct

https://web.archive.org/web/20100613221721/http://www.ece.ucdavis.edu/cerl/ReliableJPEG/Cung/dct.html

If you go from left to right, it's "DCT", compression. If you go from right to left, it's decompression. I don't understand what those black dots, white dots and arrows mean.

a1 = 0.707
a2 = 0.541
a3 = 0.707 //same as a1
a4 = 1.307
a5 = 0.383

 Description: Filesize: 21.89 KB Viewed: 17499 Time(s) Description: Filesize: 7.04 KB Viewed: 17743 Time(s) Last edited by vivik on 15 May 2018, 17:02; edited 4 times in total 05 May 2018, 12:30  vivik

Joined: 29 Oct 2016
Posts: 671
vivik

 and  mention Arai-Agui-Nakajima (AAN) algorithm, Loeffler-Ligtenberg-Moschytz (LLM), and Feig and Winograd algorithm.

MPEG   uses incompatible with jpeg implementation of dct, which is fast but less accurate.

 describes dct that is as fast as MPEG decompression

   mention dct that doesn't need any multiplication whatsoever, it uses binary shifts only. I guess this is incompatible with jpg?


http://media.taricorp.net/fast_dct.pdf
Using Streaming SIMD Extensions in a Fast DCT Algorithm for MPEG Encoding


https://www.cs.cmu.edu/~barbic/cs-740/ap922.pdf
A Fast Precise Implementation of 8x8 Discrete Cosine Transform Using the Streaming SIMD Extensions and MMX™ Instructions


https://en.wikipedia.org/wiki/JPEG_XR
"The DCT, the frequency transformation used by JPEG, is slightly lossy because of roundoff error. JPEG XR uses a type of integer transform employing a lifting scheme. The required transform, called the Photo Core Transform (PCT), resembles a 4 × 4 DCT but is lossless (exactly invertible). In fact, it is a particular realization of a larger family of binary-friendly multiplier-less transforms called the binDCT."


http://www.sfu.ca/~jiel/bindct.html


https://www.vcodex.com/h264avc-4x4-transform-and-quantization/
"This approximation is chosen to minimise the complexity of implementing the transform (multiplication by Cf requires only additions and binary shifts) whilst maintaining good compression performance." 05 May 2018, 13:03  vivik

Joined: 29 Oct 2016
Posts: 671
vivik
This binary dct is very intriguing. Can you use it with jpg, or maybe jpg-like format? Decompression speed is just as important for me as compression ratio.

Hm, this paper answers this question: http://thanglong.ece.jhu.edu/Tran/Pub/binDCT-VLSI.pdf

Same speed, mostly same quality. Can save electricity with proper hardware support. 05 May 2018, 13:13  donn

Joined: 05 Mar 2010
Posts: 214
donn
Hi, on the DCT side, I looked into the Matlab equations more and was able to compute and verify some results. Have to catch a train in a couple mins, so this is a bit rushed. I can elaborate points later if needed.

You can start breaking the two dimensional equations down by zeroing out some components and increasing the input size incrementally. I didn't find many online 2D calculators to verify results with, but did find a 1D calculator as one of the first search engine results.

I took a few simple input examples, such as:
{(3,0),(4,0)} or even {(3,0),(0,0)}

each with rowCount and columnCount of 2.

In the transform result matrix at index x,y (each base-0), you would multiply an input (such as one of the 3's above at coords 0,0) by cos((PI(2*rowIndex+1)*x)/(2*rowCount)) * cos((PI(2*columnIndex+1)*x)/(2*columnCount))

You then add each of the other input values, with the x,y values constant (since they are only going into this one destination slot in the result matrix.

Once you have the sums, you multiply by two other scalars, one for the destination x coordinate, the other for the destination matrix y coordinate. The value for each is calculated as:

When the destination x or y is 0, the scalar value is 1/(sqrt(rowCount)) (replacing rowCount with colCount when computing the y equivalent). When the destination index is not 0, the value is similar, but 1/(sqrt(2/rowCount)).

So then you just multiply these, the sum, and these two scalar values.

In the case above where you rowcount and colCount are both 2, the two values would be approximately 0.707106781186547 and 0.707106781186547 at index 0,0.

The inverse transform works on the matrix that was computed to retrive the original using a pretty similar method, but re-arranging when the scalars are multiplied.

So for the inputs used above (the first one), dct matrix result should be

((3.5,3.5),(-0.5,-0.5))

Gotta catch this train. r4's fasm output may have some readable asm of the dct, as mentioned earlier. 09 May 2018, 20:34  alexfru

Joined: 23 Mar 2014
Posts: 80
alexfru
vivik wrote:
If you go from left to right, it's "DCT", compression. If you go from right to left, it's decompression. I don't understand what those black dots, white dots and arrows mean.

https://en.wikipedia.org/wiki/Butterfly_diagram 10 May 2018, 04:25  vivik

Joined: 29 Oct 2016
Posts: 671
vivik
Found yet another variant of dct. Unfortunately, it's behind the paywall.

https://jpeg.org/jpegxs/

>with very low latency and very low complexity.

https://arxiv.org/abs/1606.07414

>Multiplierless 16-point DCT Approximation for Low-complexity Image and Video Coding

>requiring only 44 additions---the lowest arithmetic cost in literature.

https://www.dpreview.com/news/7572625808/the-new-jpeg-xs-image-format-was-built-for-streaming-4k-and-vr-content

>The new JPEG XS image format was built for streaming 4K and VR content

>What’s interesting is that JPEG isn’t trying to shrink the file size with JPEG XS. In fact, quite the opposite. Whereas the JPEG standard has a compression ratio of about 10:1, JPEG XS comes out to a 6:1 ratio. 14 May 2018, 10:30  vivik

Joined: 29 Oct 2016
Posts: 671
vivik
I don't understand anything in that butterfly diagram...

here is code for forward dct, from jfdctflt.c

Code:
```    tmp0 = dataptr + dataptr;
tmp7 = dataptr - dataptr;
tmp1 = dataptr + dataptr;
tmp6 = dataptr - dataptr;
tmp2 = dataptr + dataptr;
tmp5 = dataptr - dataptr;
tmp3 = dataptr + dataptr;
tmp4 = dataptr - dataptr;

/* Even part */

tmp10 = tmp0 + tmp3;        /* phase 2 */
tmp13 = tmp0 - tmp3;
tmp11 = tmp1 + tmp2;
tmp12 = tmp1 - tmp2;

dataptr = tmp10 + tmp11; /* phase 3 */
dataptr = tmp10 - tmp11;

z1 = (tmp12 + tmp13) * 0.707106781; /* c4 */
dataptr = tmp13 + z1;    /* phase 5 */
dataptr = 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) * 0.382683433; /* c6 */
z2 = 0.541196100 * tmp10 + z5; /* c2-c6 */
z4 = 1.306562965 * tmp12 + z5; /* c2+c6 */
z3 = tmp11 * 0.707106781; /* c4 */

z11 = tmp7 + z3;            /* phase 5 */
z13 = tmp7 - z3;

dataptr = z13 + z2;      /* phase 6 */
dataptr = z13 - z2;
dataptr = z11 + z4;
dataptr = z11 - z4;    ```

here is inverse dct, from jidctflt.c

Code:
```    /* Even part */

tmp0 = inptr;
tmp1 = inptr;
tmp2 = inptr;
tmp3 = inptr;

tmp10 = tmp0 + tmp2;        /* phase 3 */
tmp11 = tmp0 - tmp2;

tmp13 = tmp1 + tmp3;        /* phases 5-3 */
tmp12 = (tmp1 - tmp3) * 1.414213562 - tmp13; /* 2*c4 */

tmp0 = tmp10 + tmp13;       /* phase 2 */
tmp3 = tmp10 - tmp13;
tmp1 = tmp11 + tmp12;
tmp2 = tmp11 - tmp12;

/* Odd part */

tmp4 = inptr;
tmp5 = inptr;
tmp6 = inptr;
tmp7 = inptr;

z13 = tmp6 + tmp5;          /* phase 6 */
z10 = tmp6 - tmp5;
z11 = tmp4 + tmp7;
z12 = tmp4 - tmp7;

tmp7 = z11 + z13;           /* phase 5 */
tmp11 = (z11 - z13) * 1.414213562; /* 2*c4 */

z5 = (z10 + z12) * 1.847759065; /* 2*c2 */
tmp10 = z5 - z12 * 1.082392200; /* 2*(c2-c6) */
tmp12 = z5 - z10 * 2.613125930; /* 2*(c2+c6) */

tmp6 = tmp12 - tmp7;        /* phase 2 */
tmp5 = tmp11 - tmp6;
tmp4 = tmp10 - tmp5;

wsptr = tmp0 + tmp7;
wsptr = tmp0 - tmp7;
wsptr = tmp1 + tmp6;
wsptr = tmp1 - tmp6;
wsptr = tmp2 + tmp5;
wsptr = tmp2 - tmp5;
wsptr = tmp3 + tmp4;
wsptr = tmp3 - tmp4;    ```

I can't make sense out of this. I removed some unrelated things like quantization, and turned 2d transformation to 1d, 8 numbers on input and 8 numbers on output. I need to read those papers. 15 May 2018, 14:30  vivik

Joined: 29 Oct 2016
Posts: 671
vivik
Simplified it a bit more, and confirmed that this actually works:

Code:
```#include <stdio.h>

int main() {

float arr = {0, 1, 2, 3, 3, 3, 2, 0};

printf("%f %f %f %f %f %f %f %f\n\n\n",
arr, arr, arr, arr,
arr, arr, arr, arr);

{
float t0 = arr + arr;
float t7 = arr - arr;
float t1 = arr + arr;
float t6 = arr - arr;
float t2 = arr + arr;
float t5 = arr - arr;
float t3 = arr + arr;
float t4 = arr - arr;

/* Even part */

float t10 = t0 + t3;        /* phase 2 */
float t13 = t0 - t3;
float t11 = t1 + t2;
float t12 = t1 - t2;

arr = t10 + t11; /* phase 3 */
arr = t10 - t11;

float z1 = (t12 + t13) * 0.707106781; /* c4 */
arr = t13 + z1;    /* phase 5 */
arr = t13 - z1;

/* Odd part */

t10 = t4 + t5;        /* phase 2 */
t11 = t5 + t6;
t12 = t6 + t7;

/* The rotator is modified from fig 4-8 to avoid extra negations. */
float z5 = (t10 - t12) * 0.382683433; /* c6 */
float z2 = 0.541196100 * t10 + z5; /* c2-c6 */
float z4 = 1.306562965 * t12 + z5; /* c2+c6 */
float z3 = t11 * 0.707106781; /* c4 */

float z11 = t7 + z3;            /* phase 5 */
float z13 = t7 - z3;

arr = z13 + z2;      /* phase 6 */
arr = z13 - z2;
arr = z11 + z4;
arr = z11 - z4;
}

printf("%f %f %f %f %f %f %f %f\n\n\n",
arr, arr, arr, arr,
arr, arr, arr, arr);

{
/* Even part */

float t0 = arr;
float t1 = arr;
float t2 = arr;
float t3 = arr;

float t10 = t0 + t2;        /* phase 3 */
float t11 = t0 - t2;

float t13 = t1 + t3;        /* phases 5-3 */
float t12 = (t1 - t3) * 1.414213562 - t13; /* 2*c4 */

t0 = t10 + t13;       /* phase 2 */
t3 = t10 - t13;
t1 = t11 + t12;
t2 = t11 - t12;

/* Odd part */

float t4 = arr;
float t5 = arr;
float t6 = arr;
float t7 = arr;

float z13 = t6 + t5;          /* phase 6 */
float z10 = t6 - t5;
float z11 = t4 + t7;
float z12 = t4 - t7;

t7 = z11 + z13;           /* phase 5 */
t11 = (z11 - z13) * 1.414213562; /* 2*c4 */

float z5 = (z10 + z12) * 1.847759065; /* 2*c2 */
t10 = z5 - z12 * 1.082392200; /* 2*(c2-c6) */
t12 = z5 - z10 * 2.613125930; /* 2*(c2+c6) */

t6 = t12 - t7;        /* phase 2 */
t5 = t11 - t6;
t4 = t10 - t5;

arr = t0 + t7;
arr = t0 - t7;
arr = t1 + t6;
arr = t1 - t6;
arr = t2 + t5;
arr = t2 - t5;
arr = t3 + t4;
arr = t3 - t4;
}

printf("%f %f %f %f %f %f %f %f\n\n\n",
arr, arr, arr, arr,
arr, arr, arr, arr);

{
arr = arr/8;
arr = arr/8;
arr = arr/8;
arr = arr/8;
arr = arr/8;
arr = arr/8;
arr = arr/8;
arr = arr/8;
}

printf("%f %f %f %f %f %f %f %f\n\n\n",
arr, arr, arr, arr,
arr, arr, arr, arr);

getchar();
}    ```

Output is:

Code:
```0.000000 1.000000 2.000000 3.000000 3.000000 3.000000 2.000000 0.000000

14.000000 -2.720777 -11.656855 1.955410 -2.000000 0.873017 -0.343146 -0.107651

0.000000 8.000000 16.000000 24.000000 24.000000 24.000000 16.000000 0.000000

0.000000 1.000000 2.000000 3.000000 3.000000 3.000000 2.000000 0.000000    ``` 15 May 2018, 16:43  vivik

Joined: 29 Oct 2016
Posts: 671
vivik
Okay, I can read butterfly diagrams now. Lines are additions, lines with arrows are substractions (i guess because for substraction the direction is important, unlike for addition), boxes with a1 a2 a3 a4 a5 are multiplications by this constant number. Black and white circles don't mean anything, they probably indicate direction, doesn't matter.

I can follow most of the code, except this part in inverse dct:

Code:
```                float z5 = (z10 + z12) * 1.847759065; /* 2*c2 */
t10 = z5 - z12 * 1.082392200; /* 2*(c2-c6) */
t12 = z5 - z10 * 2.613125930; /* 2*(c2+c6) */    ```

It looks like code comments use those weird c2 c4 c6 constants, and I know only a1 a2 a3 a4 a5. Most of those constants are just divisions, turned into multiplications, (instead of dividing by a1, they multiply by 1/a1 instead. Same result, but slightly faster). Yet this 1.082392200 constant is 1 / 0.923, and there is no such number.

Code:
```a1 and a3
1.414213562   0.707

a2
1.8477        0.541

?????
1.082392200   0.923 ????
it should be a4 = 1.307

a5
2.613125930   0.3826    ``` 16 May 2018, 06:00  vivik

Joined: 29 Oct 2016
Posts: 671
vivik
wait, I think I found what c2 c4 c6 mean.

https://unix4lyfe.org/dct-1d/

>For brevity, let's define:

>c0 = sqrt(1/2)*sqrt(2/N)
>cj = cos(pi * j / 16) * sqrt(2/N)

then, let's check that 2*(c2-c6) == 1.082392200

c2 = cos(pi * 2 / 16) * sqrt(2 / 8)
c6 = cos(pi * 6 / 16) * sqrt(2 / 8)

well, i just typed this into google:

2 * ( (cos(pi * 2 / 16) * sqrt(2 / 8)) - (cos(pi * 6 / 16) * sqrt(2 / 8)) )

and got 0.54119610014

...that's exactly twice as small as what I expected. How odd. I guess without sqrt(2/8) it will be exact...

Hm, this is what @Siekmanski was talking about... 16 May 2018, 15:12  vivik

Joined: 29 Oct 2016
Posts: 671
vivik
0.707 = 1.414 / 2

0.707 = 1 / 1.414

0.707 = c4 = cos(pi * 4 / 16) = cos(45 grad) = cos(45 degrees) = sqrt(2)/2

wow

How would you prove that 1 / cos(45 degrees) = cos(45 degrees) / 2 ? 16 May 2018, 16:46  Display posts from previous: All Posts1 Day7 Days2 Weeks1 Month3 Months6 Months1 Year Oldest FirstNewest First  Jump to: Select a forum Official----------------AssemblyPeripheria General----------------MainTutorials and ExamplesDOSWindowsLinuxUnixMenuetOS Specific----------------MacroinstructionsOS ConstructionIDE DevelopmentProjects and IdeasNon-x86 architecturesHigh Level LanguagesProgramming Language DesignCompiler Internals Other----------------FeedbackHeapTest Area
Goto page 1, 2  Next

Forum Rules:
 You cannot post new topics in this forumYou cannot reply to topics in this forumYou cannot edit your posts in this forumYou cannot delete your posts in this forumYou cannot vote in polls in this forumYou cannot attach files in this forumYou can download files in this forum