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