Author Topic: Cosine Transform  (Read 286 times)

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

Code: [Select]
defint a-z

dim shared mx,my,mbl,mbr,mw
dim pi as double
pi = 3.141593

sw = 800
sh = 600

dim 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

_display
loop until _keyhit=27
system


sub getmouse ()
do
mx = _mousex
my = _mousey
mbl = _mousebutton(1)
mbr = _mousebutton(2)
mw = mw + _mousewheel
loop while _mouseinput
end sub


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 20180228/86 git 6fde149
QB64 v1.2 [dev build]_d84bb00

Offline Petr

  • I am instructed.
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.

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.

Offline Petr

  • I am instructed.
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   

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-z

const n = 8

type dct_type
        r as double
        g as double
        b as double
end type

type q_type
        r as _unsigned _byte
        g as _unsigned _byte
        b as _unsigned _byte
end type

dim shared pi as double
pi = _pi

img1 = _loadimage("img/greenland1.png", 32)

w = _width(img1)
h = _height(img1)

ww = (w\n+1)*n
hh = (h\n+1)*n


dim dct(ww, hh) as dct_type
dim q(ww, hh) as q_type

dim sr as double, sg as double, sb as double
dim c as double, cu as double, cv as double

img2 = _newimage(w, h, 32)

screen _newimage(w, 2*h, 32)
_putimage (0,0),img1

_source img1
_dest img2

'forward DCT
for y=0 to hh-1 step n
for 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
        next
next
next

'quantization
dim minr as double, ming as double, minb as double
dim maxr as double, maxg as double, maxb as double

minr = 1000000
ming = 1000000
minb = 1000000

maxr = -1000000
maxg = -1000000
maxb = -1000000

for y=0 to hh
for 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).b
next
next

for y=0 to hh
for 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)
next
next

'inverse DCT
for y=0 to hh-1 step n
for 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
        next
next
next

_dest 0
_putimage (0,h), img2

do
loop until _keyhit=27
system

Offline Petr

  • I am instructed.
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.

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-z

const n = 16

type dct_type
        r as double
        g as double
        b as double
end type

type q_type
        r as _unsigned _byte
        g as _unsigned _byte
        b as _unsigned _byte
end type

dim shared pi as double
pi = _pi

img1 = _loadimage("img/greenland1.png", 32)

w = _width(img1)
h = _height(img1)

ww = (w\n+1)*n
hh = (h\n+1)*n

dim dct(ww, hh) as dct_type
dim q(ww, hh) as q_type

dim sr as double, sg as double, sb as double
dim c as double, cu as double, cv as double

img2 = _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 DCT
for y0=0 to hh-1 step n
for 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
        next
next
next

'quantization
dim minr as double, ming as double, minb as double
dim maxr as double, maxg as double, maxb as double

minr = 1000000
ming = 1000000
minb = 1000000

maxr = -1000000
maxg = -1000000
maxb = -1000000

for y=0 to hh
for 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).b
next
next

_dest img1_dct
for y=0 to hh
for 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)
next
next


_dest img1_dct
for y=0 to hh
for 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)
next
next


_dest img2_dct
for y0=0 to hh-1 step n
for 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
        next
next
next

_dest img3_dct
for y0=0 to hh-1 step n
for 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
        next
next
next

_dest img2
'inverse DCT
for y0=0 to hh-1 step n
for 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
        next
next
next

_dest img3
'inverse DCT
for y0=0 to hh-1 step n
for 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
        next
next
next


_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_dct

do
loop until _keyhit=27
system
« Last Edit: July 02, 2018, 11:35:02 AM by odin »

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.

Offline Petr

  • I am instructed.
Re: Cosine Transform
« Reply #9 on: July 02, 2018, 01:09:50 PM »
Thank you very much, V. I muss to study :-D