Cosine Transform

[_vince]

Neat program I found:

Code: QB64: [Select]

1. defint a-z
2.  
3. dim shared mx,my,mbl,mbr,mw
4. dim pi as double
5. pi = 3.141593
6.  
7. sw = 800
8. sh = 600
9.  
10. dim xx(sw - 1), x(sw - 1)
11.  
12. screen _newimage(sw,sh,32)
13.  
14. do
15.         getmouse
16.  
17.         line(0, 0)-(sw, sh), _rgb(0,0,0), bf
18.  
19.         line(sw/2, 0)-(sw/2, sh), _rgb(255,0,0)
20.         line(0, sh/4)-(sw, sh/4), _rgb(255,0,0)
21.         line(0, 3*sh/4)-(sw, 3*sh/4), _rgb(255,0,0)
22.  
23.         ox = 0
24.         oy = 3*sh/4 - x(0)
25.         for i = 0 to ubound(xx)
26.                 if xx(i) <> 0 then line(i, sh/4)-step(0,-xx(i))
27.  
28.                 line(ox, oy)-(i, 3*sh/4 - x(i))
29.                 ox = i
30.                 oy = 3*sh/4 - x(i)
31.         next
32.  
33.         'if m = 0 then
34.                 line(mx, sh/4)-(mx, my)
35.         'end if
36.  
37.         if mbl then
38.                 do
39.                         getmouse
40.                 loop until mbl = 0
41.  
42.                 xx(mx) =  sh/4 - my
43.  
44.                 for i = 0 to ubound(x)
45.                         x(i) = x(i) + xx(mx) * cos(2*pi*i*(mx - sw/2)*(sw))
46.                 next
47.         end if
48.  
49.         _display
50. loop until _keyhit=27
51. system
52.  
53.  
54. sub getmouse ()
55.         do
56.                 mx = _mousex
57.                 my = _mousey
58.                 mbl = _mousebutton(1)
59.                 mbr = _mousebutton(2)
60.                 mw = mw + _mousewheel
61.         loop while _mouseinput
62. end sub
63.  
64.  

[bplus]

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?

[Petr]

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.

[codeguy]

Discrete Cosine Transform is a basis for JPEG compression. I will check this demo out later. I'm sure it will be good.

[Petr]

Hi codeguy, one JPG - create / save program is in installation QB64 in programs/samples/n54/big/jpegmake.bas   

[_vince]

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: QB64: [Select]

1. deflng a-z
2.  
3. const n = 8
4.  
5. type dct_type
6.         r as double
7.         g as double
8.         b as double
9. end type
10.  
11. type q_type
12.         r as _unsigned _byte
13.         g as _unsigned _byte
14.         b as _unsigned _byte
15. end type
16.  
17. dim shared pi as double
18. pi = _pi
19.  
20. img1 = _loadimage("img/greenland1.png", 32)
21.  
22. w = _width(img1)
23. h = _height(img1)
24.  
25. ww = (w\n+1)*n
26. hh = (h\n+1)*n
27.  
28.  
29. dim dct(ww, hh) as dct_type
30. dim q(ww, hh) as q_type
31.  
32. dim sr as double, sg as double, sb as double
33. dim c as double, cu as double, cv as double
34.  
35. img2 = _newimage(w, h, 32)
36.  
37. screen _newimage(w, 2*h, 32)
38. _putimage (0,0),img1
39.  
40. _source img1
41. _dest img2
42.  
43. 'forward DCT
44. for y=0 to hh-1 step n
45. for x=0 to ww-1 step n
46.         for yy=0 to n-1
47.         for xx=0 to n-1
48.                 sr = 0
49.                 sg = 0
50.                 sb = 0
51.  
52.                 for y3=0 to n-1
53.                 for x3=0 to n-1
54.                         if (x + x3 > w - 1) then px = x + x3 - n else px = x + x3
55.                         if (y + y3 > h - 1) then py = y + x3 - n else py = y + y3
56.  
57.                         z = point(px, py)
58.                         r = _red(z)
59.                         g = _green(z)
60.                         b = _blue(z)
61.  
62.                         c = cos((2*x3 + 1)*xx*pi/(2*n))*cos((2*y3 + 1)*yy*pi/(2*n))
63.  
64.                         sr = sr + r*c
65.                         sg = sg + g*c
66.                         sb = sb + b*c
67.                 next
68.                 next
69.  
70.                 if xx = 0 then cu = 1/sqr(2) else cu = 1
71.                 if yy = 0 then cv = 1/sqr(2) else cv = 1
72.  
73.                 dct(x + xx, y + yy).r = sr*cu*cv/(0.5*n)
74.                 dct(x + xx, y + yy).g = sg*cu*cv/(0.5*n)
75.                 dct(x + xx, y + yy).b = sb*cu*cv/(0.5*n)
76.         next
77.         next
78. next
79. next
80.  
81. 'quantization
82. dim minr as double, ming as double, minb as double
83. dim maxr as double, maxg as double, maxb as double
84.  
85. minr = 1000000
86. ming = 1000000
87. minb = 1000000
88.  
89. maxr = -1000000
90. maxg = -1000000
91. maxb = -1000000
92.  
93. for y=0 to hh
94. for x=0 to ww
95.         if dct(x, y).r < minr then minr = dct(x, y).r
96.         if dct(x, y).g < ming then ming = dct(x, y).g
97.         if dct(x, y).b < minb then minb = dct(x, y).b
98.  
99.         if dct(x, y).r > maxr then maxr = dct(x, y).r
100.         if dct(x, y).g > maxg then maxg = dct(x, y).g
101.         if dct(x, y).b > maxb then maxb = dct(x, y).b
102. next
103. next
104.  
105. for y=0 to hh
106. for x=0 to ww
107.         q(x, y).r = 255*(dct(x,y).r - minr)/(maxr - minr)
108.         q(x, y).g = 255*(dct(x,y).g - ming)/(maxg - ming)
109.         q(x, y).b = 255*(dct(x,y).b - minb)/(maxb - minb)
110. next
111. next
112.  
113. 'inverse DCT
114. for y=0 to hh-1 step n
115. for x=0 to ww-1 step n
116.         for yy=0 to n-1
117.         for xx=0 to n-1
118.                 sr = 0
119.                 sg = 0
120.                 sb = 0
121.  
122.                 for y3=0 to 2 'n-1
123.                 for x3=0 to 2 'n-1
124.                         c = cos((2*xx + 1)*x3*pi/(2*n))*cos((2*yy + 1)*y3*pi/(2*n))
125.  
126.                         if x3 = 0 then cu = 1/sqr(2) else cu = 1
127.                         if y3 = 0 then cv = 1/sqr(2) else cv = 1
128.  
129.                         'sr = sr + dct(x + x3, y + y3).r*c*cu*cv
130.                         'sg = sg + dct(x + x3, y + y3).g*c*cu*cv
131.                         'sb = sb + dct(x + x3, y + y3).b*c*cu*cv
132.  
133.                         r = q(x + x3, y + y3).r
134.                         g =     q(x + x3, y + y3).g
135.                         b = q(x + x3, y + y3).b
136.  
137.                         sr = sr + c*cu*cv*((r/255)*(maxr - minr) + minr)
138.                         sg = sg + c*cu*cv*((g/255)*(maxg - ming) + ming)
139.                         sb = sb + c*cu*cv*((b/255)*(maxb - minb) + minb)
140.                 next
141.                 next
142.  
143.                 sr = sr/(0.5*n)
144.                 sg = sg/(0.5*n)
145.                 sb = sb/(0.5*n)
146.  
147.                 if (x + xx < w) and (y + yy < h) then pset (x + xx, y + yy), _rgb(sr, sg, sb)
148.         next
149.         next
150. next
151. next
152.  
153. _dest 0
154. _putimage (0,h), img2
155.  
156. do
157. loop until _keyhit=27
158. system
159.  

[Petr]

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.

[_vince]

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: QB64: [Select]

1. deflng a-z
2.  
3. const n = 16
4.  
5. type dct_type
6.         r as double
7.         g as double
8.         b as double
9. end type
10.  
11. type q_type
12.         r as _unsigned _byte
13.         g as _unsigned _byte
14.         b as _unsigned _byte
15. end type
16.  
17. dim shared pi as double
18. pi = _pi
19.  
20. img1 = _loadimage("img/greenland1.png", 32)
21.  
22. w = _width(img1)
23. h = _height(img1)
24.  
25. ww = (w\n+1)*n
26. hh = (h\n+1)*n
27.  
28. dim dct(ww, hh) as dct_type
29. dim q(ww, hh) as q_type
30.  
31. dim sr as double, sg as double, sb as double
32. dim c as double, cu as double, cv as double
33.  
34. img2 = _newimage(w, h, 32)
35. img3 = _newimage(w, h, 32)
36.  
37. img1_dct = _newimage(w, h, 32)
38. img2_dct = _newimage(w, h, 32)
39. img3_dct = _newimage(w, h, 32)
40.  
41. screen _newimage(3*w, 2*h, 32)
42. _putimage (0,0),img1
43.  
44. _source img1
45.  
46. 'forward DCT
47. for y0=0 to hh-1 step n
48. for x0=0 to ww-1 step n
49.         for y=0 to n-1
50.         for x=0 to n-1
51.                 sr = 0
52.                 sg = 0
53.                 sb = 0
54.  
55.                 for v=0 to n-1
56.                 for u=0 to n-1
57.                         if (x0 + u > w - 1) then px = x0 + u - n else px = x0 + u
58.                         if (y0 + v > h - 1) then py = y0 + v - n else py = y0 + v
59.  
60.                         z = point(px, py)
61.                         r = _red(z)
62.                         g = _green(z)
63.                         b = _blue(z)
64.  
65.                         c = cos((2*u + 1)*x*pi/(2*n)) * cos((2*v + 1)*y*pi/(2*n))
66.  
67.                         sr = sr + r*c
68.                         sg = sg + g*c
69.                         sb = sb + b*c
70.                 next
71.                 next
72.  
73.                 if x = 0 then cu = 1/sqr(2) else cu = 1
74.                 if y = 0 then cv = 1/sqr(2) else cv = 1
75.  
76.                 dct(x0 + x, y0 + y).r = sr*cu*cv/(0.5*n)
77.                 dct(x0 + x, y0 + y).g = sg*cu*cv/(0.5*n)
78.                 dct(x0 + x, y0 + y).b = sb*cu*cv/(0.5*n)
79.         next
80.         next
81. next
82. next
83.  
84. 'quantization
85. dim minr as double, ming as double, minb as double
86. dim maxr as double, maxg as double, maxb as double
87.  
88. minr = 1000000
89. ming = 1000000
90. minb = 1000000
91.  
92. maxr = -1000000
93. maxg = -1000000
94. maxb = -1000000
95.  
96. for y=0 to hh
97. for x=0 to ww
98.         if dct(x, y).r < minr then minr = dct(x, y).r
99.         if dct(x, y).g < ming then ming = dct(x, y).g
100.         if dct(x, y).b < minb then minb = dct(x, y).b
101.  
102.         if dct(x, y).r > maxr then maxr = dct(x, y).r
103.         if dct(x, y).g > maxg then maxg = dct(x, y).g
104.         if dct(x, y).b > maxb then maxb = dct(x, y).b
105. next
106. next
107.  
108. _dest img1_dct
109. for y=0 to hh
110. for x=0 to ww
111.         r = q(x, y).r
112.         g = q(x, y).g
113.         b = q(x, y).b
114.  
115.         pset (x, y), _rgb(r, g, b)
116. next
117. next
118.  
119.  
120. _dest img1_dct
121. for y=0 to hh
122. for x=0 to ww
123.         q(x, y).r = 255*(dct(x,y).r - minr)/(maxr - minr)
124.         q(x, y).g = 255*(dct(x,y).g - ming)/(maxg - ming)
125.         q(x, y).b = 255*(dct(x,y).b - minb)/(maxb - minb)
126.  
127.         r = q(x, y).r
128.         g = q(x, y).g
129.         b = q(x, y).b
130.  
131.         pset (x, y), _rgb(r, g, b)
132. next
133. next
134.  
135.  
136. _dest img2_dct
137. for y0=0 to hh-1 step n
138. for x0=0 to ww-1 step n
139.         for y=0 to 7 'n-1
140.         for x=0 to 7 'n-1
141.                 r = q(x0 + x, y0 + y).r
142.                 g = q(x0 + x, y0 + y).g
143.                 b = q(x0 + x, y0 + y).b
144.  
145.                 if (x0 + x < w) and (y0 + y < h) then pset (x0 + x, y0 + y), _rgb(r, g, b)
146.         next
147.         next
148. next
149. next
150.  
151. _dest img3_dct
152. for y0=0 to hh-1 step n
153. for x0=0 to ww-1 step n
154.         for y=0 to 2 'n-1
155.         for x=0 to 2 'n-1
156.                 r = q(x0 + x, y0 + y).r
157.                 g = q(x0 + x, y0 + y).g
158.                 b = q(x0 + x, y0 + y).b
159.  
160.                 if (x0 + x < w) and (y0 + y < h) then pset (x0 + x, y0 + y), _rgb(r, g, b)
161.         next
162.         next
163. next
164. next
165.  
166. _dest img2
167. 'inverse DCT
168. for y0=0 to hh-1 step n
169. for x0=0 to ww-1 step n
170.         for y=0 to n-1
171.         for x=0 to n-1
172.                 sr = 0
173.                 sg = 0
174.                 sb = 0
175.  
176.                 for v=0 to 7 'n-1
177.                 for u=0 to 7 'n-1
178.                         c = cos((2*x + 1)*u*pi/(2*n))*cos((2*y + 1)*v*pi/(2*n))
179.  
180.                         if u = 0 then cu = 1/sqr(2) else cu = 1
181.                         if v = 0 then cv = 1/sqr(2) else cv = 1
182.  
183.                         'sr = sr + dct(x + x3, y + y3).r*c*cu*cv
184.                         'sg = sg + dct(x + x3, y + y3).g*c*cu*cv
185.                         'sb = sb + dct(x + x3, y + y3).b*c*cu*cv
186.  
187.                         r = q(x0 + u, y0 + v).r
188.                         g =     q(x0 + u, y0 + v).g
189.                         b = q(x0 + u, y0 + v).b
190.  
191.                         sr = sr + c*cu*cv*((r/255)*(maxr - minr) + minr)
192.                         sg = sg + c*cu*cv*((g/255)*(maxg - ming) + ming)
193.                         sb = sb + c*cu*cv*((b/255)*(maxb - minb) + minb)
194.                 next
195.                 next
196.  
197.                 sr = sr/(0.5*n)
198.                 sg = sg/(0.5*n)
199.                 sb = sb/(0.5*n)
200.  
201.                 if (x0 + x < w) and (y0 + y < h) then pset (x0 + x, y0 + y), _rgb(sr, sg, sb)
202.         next
203.         next
204. next
205. next
206.  
207. _dest img3
208. 'inverse DCT
209. for y0=0 to hh-1 step n
210. for x0=0 to ww-1 step n
211.         for y=0 to n-1
212.         for x=0 to n-1
213.                 sr = 0
214.                 sg = 0
215.                 sb = 0
216.  
217.                 for v=0 to 2
218.                 for u=0 to 2
219.                         c = cos((2*x + 1)*u*pi/(2*n))*cos((2*y + 1)*v*pi/(2*n))
220.  
221.                         if u = 0 then cu = 1/sqr(2) else cu = 1
222.                         if v = 0 then cv = 1/sqr(2) else cv = 1
223.  
224.                         'sr = sr + dct(x + x3, y + y3).r*c*cu*cv
225.                         'sg = sg + dct(x + x3, y + y3).g*c*cu*cv
226.                         'sb = sb + dct(x + x3, y + y3).b*c*cu*cv
227.  
228.                         r = q(x0 + u, y0 + v).r
229.                         g = q(x0 + u, y0 + v).g
230.                         b = q(x0 + u, y0 + v).b
231.  
232.                         sr = sr + c*cu*cv*((r/255)*(maxr - minr) + minr)
233.                         sg = sg + c*cu*cv*((g/255)*(maxg - ming) + ming)
234.                         sb = sb + c*cu*cv*((b/255)*(maxb - minb) + minb)
235.                 next
236.                 next
237.  
238.                 sr = sr/(0.5*n)
239.                 sg = sg/(0.5*n)
240.                 sb = sb/(0.5*n)
241.  
242.                 if (x0 + x < w) and (y0 + y < h) then pset (x0 + x, y0 + y), _rgb(sr, sg, sb)
243.         next
244.         next
245. next
246. next
247.  
248.  
249. _dest 0
250. _putimage (w,0), img2
251. _putimage (2*w,0), img3
252. _putimage (0,h), img1_dct
253. _putimage (w,h), img2_dct
254. _putimage (2*w,h), img3_dct
255.  
256. do
257. loop until _keyhit=27
258. system
259.  

[_vince]

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]

Thank you very much, V. I muss to study :-D