### Author Topic: Cosine Transform  (Read 235 times)

#### v

##### Cosine Transform
« on: June 27, 2018, 08:03:36 PM »
Neat program I found:

Code: [Select]
`defint a-zdim shared mx,my,mbl,mbr,mwdim pi as doublepi = 3.141593sw = 800sh = 600dim xx(sw - 1), x(sw - 1)screen _newimage(sw,sh,32)do getmouse line(0, 0)-(sw, sh), _rgb(0,0,0), bf line(sw/2, 0)-(sw/2, sh), _rgb(255,0,0) line(0, sh/4)-(sw, sh/4), _rgb(255,0,0) line(0, 3*sh/4)-(sw, 3*sh/4), _rgb(255,0,0) ox = 0 oy = 3*sh/4 - x(0) for i = 0 to ubound(xx) if xx(i) <> 0 then line(i, sh/4)-step(0,-xx(i)) line(ox, oy)-(i, 3*sh/4 - x(i)) ox = i oy = 3*sh/4 - x(i) next 'if m = 0 then line(mx, sh/4)-(mx, my) 'end if if mbl then do getmouse loop until mbl = 0 xx(mx) =  sh/4 - my for i = 0 to ubound(x) x(i) = x(i) + xx(mx) * cos(2*pi*i*(mx - sw/2)*(sw)) next end if _displayloop until _keyhit=27systemsub getmouse () do mx = _mousex my = _mousey mbl = _mousebutton(1) mbr = _mousebutton(2) mw = mw + _mousewheel loop while _mouseinputend sub`

#### bplus

##### Re: Cosine Transform
« Reply #1 on: June 28, 2018, 11:42:23 AM »
Hi v,

This reminds me of Fourier, something about being able to create a wave formula for any finite set of data ie a curve that will fit any graphed set of data points (as long as there is always one and only one y value for any x value AKA a continuous function), something like that..

The only use for it I've witnessed so far was to create a cartoon alien planet surface for lander practice in Just Basic Lander.bas code.

Has anyone another application idea?
« Last Edit: June 28, 2018, 11:43:30 AM by bplus »
B = B + ...
QB64 x 64 v1.2 2018 0228/86 git b30af92
QB64 v1.2 2018 0228/86 git b30af92
QB64 v1.1 2017 1106/82

#### Petr

##### Re: Cosine Transform
« Reply #2 on: June 28, 2018, 04:55:53 PM »
I think, that this transformation is the basis for JPG and MP3. But as it work, I really do not know. But I have written a head for MP3 if someone wanted to try to disassemble the MP3 samples and then use them in _SNDRAW.
Coding is relax (At least sometimes)

#### codeguy

##### Re: Cosine Transform
« Reply #3 on: June 28, 2018, 08:12:23 PM »
Discrete Cosine Transform is a basis for JPEG compression. I will check this demo out later. I'm sure it will be good.

#### Petr

##### Re: Cosine Transform
« Reply #4 on: June 29, 2018, 11:24:16 AM »
Hi codeguy, one JPG - create / save program is in installation QB64 in programs/samples/n54/big/jpegmake.bas
Coding is relax (At least sometimes)

#### v

##### Re: Cosine Transform
« Reply #5 on: July 01, 2018, 11:19:25 AM »
Here's a demonstration of how the DCT can be used in image compression.  The second image is about 15% the size of the first since it is created by only using 9 out of 64 DCT coefficients quantized to 256 bits in the inverse computation yet is still quite legible.  It's several steps short of the full JPEG algorithm but demonstrates the same idea.

Code: [Select]
`deflng a-zconst n = 8type dct_type        r as double        g as double        b as doubleend typetype q_type        r as _unsigned _byte        g as _unsigned _byte        b as _unsigned _byteend typedim shared pi as doublepi = _piimg1 = _loadimage("img/greenland1.png", 32)w = _width(img1)h = _height(img1)ww = (w\n+1)*nhh = (h\n+1)*ndim dct(ww, hh) as dct_typedim q(ww, hh) as q_typedim sr as double, sg as double, sb as doubledim c as double, cu as double, cv as doubleimg2 = _newimage(w, h, 32)screen _newimage(w, 2*h, 32)_putimage (0,0),img1_source img1_dest img2'forward DCTfor y=0 to hh-1 step nfor x=0 to ww-1 step n        for yy=0 to n-1        for xx=0 to n-1                sr = 0                sg = 0                sb = 0                for y3=0 to n-1                for x3=0 to n-1                        if (x + x3 > w - 1) then px = x + x3 - n else px = x + x3                        if (y + y3 > h - 1) then py = y + x3 - n else py = y + y3                        z = point(px, py)                        r = _red(z)                        g = _green(z)                        b = _blue(z)                        c = cos((2*x3 + 1)*xx*pi/(2*n))*cos((2*y3 + 1)*yy*pi/(2*n))                        sr = sr + r*c                        sg = sg + g*c                        sb = sb + b*c                next                next                if xx = 0 then cu = 1/sqr(2) else cu = 1                if yy = 0 then cv = 1/sqr(2) else cv = 1                dct(x + xx, y + yy).r = sr*cu*cv/(0.5*n)                dct(x + xx, y + yy).g = sg*cu*cv/(0.5*n)                dct(x + xx, y + yy).b = sb*cu*cv/(0.5*n)        next        nextnextnext'quantizationdim minr as double, ming as double, minb as doubledim maxr as double, maxg as double, maxb as doubleminr = 1000000ming = 1000000minb = 1000000maxr = -1000000 maxg = -1000000 maxb = -1000000 for y=0 to hhfor x=0 to ww        if dct(x, y).r < minr then minr = dct(x, y).r        if dct(x, y).g < ming then ming = dct(x, y).g        if dct(x, y).b < minb then minb = dct(x, y).b        if dct(x, y).r > maxr then maxr = dct(x, y).r        if dct(x, y).g > maxg then maxg = dct(x, y).g        if dct(x, y).b > maxb then maxb = dct(x, y).bnextnextfor y=0 to hhfor x=0 to ww        q(x, y).r = 255*(dct(x,y).r - minr)/(maxr - minr)        q(x, y).g = 255*(dct(x,y).g - ming)/(maxg - ming)        q(x, y).b = 255*(dct(x,y).b - minb)/(maxb - minb)nextnext'inverse DCTfor y=0 to hh-1 step nfor x=0 to ww-1 step n        for yy=0 to n-1        for xx=0 to n-1                sr = 0                sg = 0                sb = 0                for y3=0 to 2 'n-1                for x3=0 to 2 'n-1                        c = cos((2*xx + 1)*x3*pi/(2*n))*cos((2*yy + 1)*y3*pi/(2*n))                        if x3 = 0 then cu = 1/sqr(2) else cu = 1                        if y3 = 0 then cv = 1/sqr(2) else cv = 1                        'sr = sr + dct(x + x3, y + y3).r*c*cu*cv                        'sg = sg + dct(x + x3, y + y3).g*c*cu*cv                        'sb = sb + dct(x + x3, y + y3).b*c*cu*cv                        r = q(x + x3, y + y3).r                         g =     q(x + x3, y + y3).g                         b = q(x + x3, y + y3).b                         sr = sr + c*cu*cv*((r/255)*(maxr - minr) + minr)                        sg = sg + c*cu*cv*((g/255)*(maxg - ming) + ming)                        sb = sb + c*cu*cv*((b/255)*(maxb - minb) + minb)                next                next                sr = sr/(0.5*n)                sg = sg/(0.5*n)                sb = sb/(0.5*n)                if (x + xx < w) and (y + yy < h) then pset (x + xx, y + yy), _rgb(sr, sg, sb)        next        nextnextnext_dest 0_putimage (0,h), img2doloop until _keyhit=27system`

#### Petr

##### Re: Cosine Transform
« Reply #6 on: July 01, 2018, 02:15:25 PM »
I say straight away that I will need some time to study and understand the function of your program and to study more details of this compression. I appreciate for your demos.
Coding is relax (At least sometimes)

#### v

##### Re: Cosine Transform
« Reply #7 on: July 02, 2018, 12:14:21 AM »
Hi Petr,

I thought I should explain the program. In this particular case I used N=16 for better clarity.  You take the NxN 2-dimensional DCT of each NxN square in the image shown in the first pair of images.  The NxN DCT produces NxN coefficients from NxN pixels.  You notice that most of the 'information' about each square is contained in the top-left corner corresponding to low frequencies.  I guess the property of photographs is that most of their frequency content is focused in the low frequency region.  High frequencies correspond to minute details that might be important in line art but less so in natural photographs and appear mostly gray in the attached image.

Taking the inverse DCT on the NxN DCT coefficients will produce the original images (minus some quantization error if applicable).  We can do a sort of low pass filtering by removing the high frequencies and then reconstructing the image, for example taking the inverse DCT on only (N/2)x(N/2) of the coefficients.  The reconstructed image remains quite legible but takes up only 25% of storage space.  The JPEG algorithm does something similar, but rather than removing high frequencies altogether, it preforms entropy encoding (i.e. Huffman coding) on them which assigns less bits to the less relevant frequencies.

The DCT has several other applications in image processing other than just compression.  Taking the DCT of an entire image and zeroing out frequency components can be done to filter out undesired parts of the image.  For example, if you take a photograph behind a chain-link fence, its repetitive pattern corresponds to particular frequencies in the DCT domain.  Zeroing them out and reconstructing the image might remove the chain-link fence in the new image.  If there's interest, I can make more DCT image processing demos.

Source for attached image:
Code: [Select]
`deflng a-zconst n = 16type dct_type        r as double        g as double        b as doubleend typetype q_type        r as _unsigned _byte        g as _unsigned _byte        b as _unsigned _byteend typedim shared pi as doublepi = _piimg1 = _loadimage("img/greenland1.png", 32)w = _width(img1)h = _height(img1)ww = (w\n+1)*nhh = (h\n+1)*ndim dct(ww, hh) as dct_typedim q(ww, hh) as q_typedim sr as double, sg as double, sb as doubledim c as double, cu as double, cv as doubleimg2 = _newimage(w, h, 32)img3 = _newimage(w, h, 32)img1_dct = _newimage(w, h, 32)img2_dct = _newimage(w, h, 32)img3_dct = _newimage(w, h, 32)screen _newimage(3*w, 2*h, 32)_putimage (0,0),img1_source img1'forward DCTfor y0=0 to hh-1 step nfor x0=0 to ww-1 step n        for y=0 to n-1        for x=0 to n-1                sr = 0                sg = 0                sb = 0                for v=0 to n-1                for u=0 to n-1                        if (x0 + u > w - 1) then px = x0 + u - n else px = x0 + u                        if (y0 + v > h - 1) then py = y0 + v - n else py = y0 + v                        z = point(px, py)                        r = _red(z)                        g = _green(z)                        b = _blue(z)                        c = cos((2*u + 1)*x*pi/(2*n)) * cos((2*v + 1)*y*pi/(2*n))                        sr = sr + r*c                        sg = sg + g*c                        sb = sb + b*c                next                next                if x = 0 then cu = 1/sqr(2) else cu = 1                if y = 0 then cv = 1/sqr(2) else cv = 1                dct(x0 + x, y0 + y).r = sr*cu*cv/(0.5*n)                dct(x0 + x, y0 + y).g = sg*cu*cv/(0.5*n)                dct(x0 + x, y0 + y).b = sb*cu*cv/(0.5*n)        next        nextnextnext'quantizationdim minr as double, ming as double, minb as doubledim maxr as double, maxg as double, maxb as doubleminr = 1000000ming = 1000000minb = 1000000maxr = -1000000 maxg = -1000000 maxb = -1000000 for y=0 to hhfor x=0 to ww        if dct(x, y).r < minr then minr = dct(x, y).r        if dct(x, y).g < ming then ming = dct(x, y).g        if dct(x, y).b < minb then minb = dct(x, y).b        if dct(x, y).r > maxr then maxr = dct(x, y).r        if dct(x, y).g > maxg then maxg = dct(x, y).g        if dct(x, y).b > maxb then maxb = dct(x, y).bnextnext_dest img1_dctfor y=0 to hhfor x=0 to ww        r = q(x, y).r         g = q(x, y).g         b = q(x, y).b         pset (x, y), _rgb(r, g, b)nextnext_dest img1_dctfor y=0 to hhfor x=0 to ww        q(x, y).r = 255*(dct(x,y).r - minr)/(maxr - minr)        q(x, y).g = 255*(dct(x,y).g - ming)/(maxg - ming)        q(x, y).b = 255*(dct(x,y).b - minb)/(maxb - minb)        r = q(x, y).r        g = q(x, y).g        b = q(x, y).b        pset (x, y), _rgb(r, g, b)nextnext_dest img2_dctfor y0=0 to hh-1 step nfor x0=0 to ww-1 step n        for y=0 to 7 'n-1        for x=0 to 7 'n-1                r = q(x0 + x, y0 + y).r                 g = q(x0 + x, y0 + y).g                 b = q(x0 + x, y0 + y).b                 if (x0 + x < w) and (y0 + y < h) then pset (x0 + x, y0 + y), _rgb(r, g, b)        next        nextnextnext_dest img3_dctfor y0=0 to hh-1 step nfor x0=0 to ww-1 step n        for y=0 to 2 'n-1        for x=0 to 2 'n-1                r = q(x0 + x, y0 + y).r                 g = q(x0 + x, y0 + y).g                 b = q(x0 + x, y0 + y).b                 if (x0 + x < w) and (y0 + y < h) then pset (x0 + x, y0 + y), _rgb(r, g, b)        next        nextnextnext_dest img2'inverse DCTfor y0=0 to hh-1 step nfor x0=0 to ww-1 step n        for y=0 to n-1        for x=0 to n-1                sr = 0                sg = 0                sb = 0                for v=0 to 7 'n-1                for u=0 to 7 'n-1                        c = cos((2*x + 1)*u*pi/(2*n))*cos((2*y + 1)*v*pi/(2*n))                        if u = 0 then cu = 1/sqr(2) else cu = 1                        if v = 0 then cv = 1/sqr(2) else cv = 1                        'sr = sr + dct(x + x3, y + y3).r*c*cu*cv                        'sg = sg + dct(x + x3, y + y3).g*c*cu*cv                        'sb = sb + dct(x + x3, y + y3).b*c*cu*cv                        r = q(x0 + u, y0 + v).r                         g =     q(x0 + u, y0 + v).g                         b = q(x0 + u, y0 + v).b                         sr = sr + c*cu*cv*((r/255)*(maxr - minr) + minr)                        sg = sg + c*cu*cv*((g/255)*(maxg - ming) + ming)                        sb = sb + c*cu*cv*((b/255)*(maxb - minb) + minb)                next                next                sr = sr/(0.5*n)                sg = sg/(0.5*n)                sb = sb/(0.5*n)                if (x0 + x < w) and (y0 + y < h) then pset (x0 + x, y0 + y), _rgb(sr, sg, sb)        next        nextnextnext_dest img3'inverse DCTfor y0=0 to hh-1 step nfor x0=0 to ww-1 step n        for y=0 to n-1        for x=0 to n-1                sr = 0                sg = 0                sb = 0                for v=0 to 2                for u=0 to 2                        c = cos((2*x + 1)*u*pi/(2*n))*cos((2*y + 1)*v*pi/(2*n))                        if u = 0 then cu = 1/sqr(2) else cu = 1                        if v = 0 then cv = 1/sqr(2) else cv = 1                        'sr = sr + dct(x + x3, y + y3).r*c*cu*cv                        'sg = sg + dct(x + x3, y + y3).g*c*cu*cv                        'sb = sb + dct(x + x3, y + y3).b*c*cu*cv                        r = q(x0 + u, y0 + v).r                         g = q(x0 + u, y0 + v).g                         b = q(x0 + u, y0 + v).b                         sr = sr + c*cu*cv*((r/255)*(maxr - minr) + minr)                        sg = sg + c*cu*cv*((g/255)*(maxg - ming) + ming)                        sb = sb + c*cu*cv*((b/255)*(maxb - minb) + minb)                next                next                sr = sr/(0.5*n)                sg = sg/(0.5*n)                sb = sb/(0.5*n)                if (x0 + x < w) and (y0 + y < h) then pset (x0 + x, y0 + y), _rgb(sr, sg, sb)        next        nextnextnext_dest 0_putimage (w,0), img2_putimage (2*w,0), img3_putimage (0,h), img1_dct_putimage (w,h), img2_dct_putimage (2*w,h), img3_dctdoloop until _keyhit=27system`
« Last Edit: July 02, 2018, 11:35:02 AM by odin »

#### v

##### Re: Cosine Transform
« Reply #8 on: July 02, 2018, 04:14:44 AM »
Out of interest, in this image I filtered out the lower frequencies rather than high frequencies and enhanced the otherwise very faint colors.  This highlights the information about the image lost in the above compression schemes.

#### Petr

##### Re: Cosine Transform
« Reply #9 on: July 02, 2018, 01:09:50 PM »
Thank you very much, V. I muss to study :-D
Coding is relax (At least sometimes)