In [98]:
# Here's what they look like.
figure(figsize=(9,12))
for u in xrange(N):
subplot(4, 2, u+1)
ylim((-1, 1))
title(str(u))
plot(dct[u, :])
plot(dct[u, :],'ro')
Jim Mahoney | Computer Science @ Marlboro College | cs.marlboro.edu | Feb 2014 | MIT License
An explanation and illustration of the math behind the Discrete Cosine Transform and the concepts used in lossy JPEG image compression - low pass filtering.
You can think of a vector - a list of numbers - as coefficients times basis vectors.
\[ f_0 \left[ \begin{array}{c} 1 \\ 0 \end{array} \right] + f_1 \left[ \begin{array}{c} 0 \\ 1 \end{array} \right] \]
Using a different basis, different coefficients can describe the same vector.
\[ G_0 \frac{1}{\sqrt{2}} \left[ \begin{array}{c} 1 \\ 1 \end{array} \right] + G_1 \frac{1}{\sqrt{2}} \left[ \begin{array}{c} 1 \\ -1 \end{array} \right] \]
(The sqrt(2)'s give the basis vectors length 1, i.e. "normalizes" them.)
This transormation f to G is a DCT (Discrete Cosine Transform). For a vector with 2 components, this perhaps isn't all that exciting, but does still transform the original \((f_0, f_1)\) into low and high frequency components \((G_0, G_1)\).
If \(f_0 = 3\) and \(f_1 = 5\), what are the G's?
This transform can be written as a matrix multiplication.
\[ f_0 \left[ \begin{array}{c} 1 \\ 0 \end{array} \right] + f_1 \left[ \begin{array}{c} 0 \\ 1 \end{array} \right] = \left[ \begin{array}{c} f_0 \\ f_1 \end{array} \right] = G_0 \frac{1}{\sqrt{2}} \left[ \begin{array}{c} 1 \\ 1 \end{array} \right] + G_1 \frac{1}{\sqrt{2}} \left[ \begin{array}{c} 1 \\ -1 \end{array} \right] = \frac{1}{\sqrt{2}} \left[ \begin{array}{cc} 1 & 1 \\ 1 & -1 \end{array} \right] \left[ \begin{array}{c} G_0 \\ G_1 \end{array} \right] \]
Moreover, this orthnormal matrix has the interesting and useful property that its transpose is its inverse. That makes the equation easy to invert.
Show that the matrix times its transpose is the identity, and use that to find the G's.
The same idea can be applied to 2D images rather than 1D vectors, by applying the 1D transform to each row and column of the image.
The 2D basis images for N=2 are then the outer products of the 1D basis vectors. From lowest (0,0) to highest (1,1) spatial frequency these basis images are :
basis = (1/sqrt(2) * array([1, 1]), 1/sqrt(2) * array([1, -1]))
for i in [0,1]:
for j in [0,1]:
print "{}, {} :".format(i,j)
print outer(basis[i], basis[j])
print
For an image \(f = \left[ \begin{array}{cc} 5 & 8 \\ 4 & -1 \end{array} \right]\), what are the correspoding four \(G\) coefficients?
JPEG image compression uses the same sort of transform but with 8 coefficients, not 2.
The matrix is defined by this formula :
\[ G_u = \frac{0.5}{\sqrt{2}} f_0 + 0.5 \sum_{x=1}^7 f_x \, cos\left( \frac{\pi}{8} (u + \frac{1}{2}) x \right) \]
(See http://en.wikipedia.org/wiki/Discrete_cosine_transform and http://www.whydomath.org/node/wavlets/dct.html for the details.)
# The 8 x 8 DCT matrix thus looks like this.
N = 8
dct = zeros((8,8))
for x in range(N):
dct[0,x] = 0.5/sqrt(2)
for u in xrange(1,N):
for x in xrange(N):
dct[u,x] = 0.5 * cos((pi/N) * u * (x + 0.5) )
np.set_printoptions(precision=3)
dct
The corresponding eight 1D basis vectors (the matrix rows) oscillate with successively higher spatial frequencies.
# Here's what they look like.
figure(figsize=(9,12))
for u in xrange(N):
subplot(4, 2, u+1)
ylim((-1, 1))
title(str(u))
plot(dct[u, :])
plot(dct[u, :],'ro')
Like the N=2 case, the vectors are orthnormal. In other words, their dot products are 0, and each has length 1. Here are a few illustrative products.
def rowdot(i,j):
return dot(dct[i, :], dct[j, :])
rowdot(0,0), rowdot(3,3), rowdot(0,3), rowdot(1, 7), rowdot(1,5)
This also implies the inverse of this matrix is just its transpose.
dct_transpose = dct.transpose()
dct_transpose
# Is the dot product of dct and its transpose the identity?
maybe_identity = dot(dct, dct_transpose)
# Since there are many nearly zero like 3.2334e-17 in this numerical result,
# the output will look much nicer if we round them all of to (say) 6 places.
roundoff = vectorize(lambda m: round(m, 6))
roundoff(maybe_identity)
To make all this more concrete, let's apply the 2D DCT transform to part of a real image.
Here's one, takenly randomly from the web.
# See http://matplotlib.org/users/image_tutorial.html for the image manipulation syntax.
# The image itself is a small piece from http://www.cordwainer-smith.com/virgil_finlay.htm.
import matplotlib.image as mpimg
img = mpimg.imread('stormplanet112.jpg')
p=plt.imshow(img, origin='lower')
# The image itself contains 3 dimensions: rows, columns, and colors
img.shape
All three of the R,G,B color values in the greyscale image are the same for each pixel.
Let's just look at values from one tiny 8 x 8 block (which is what's used JPEG compression) near his nose.
(The next images use a false color spectrum to display pixel intensity.)
tiny = img[40:48, 40:48, 0] # a tiny 8 x 8 block, in the color=0 (Red) channel
def show_image(img):
plt.imshow(img)
plt.colorbar()
show_image(tiny)
# And here are the numbers.
tiny
Now we define the 2D version of the N=8 DCT described above.
The trick is to apply the 1D DCT to every column, and then also apply it to every row, i.e.
\[ G = {DCT} \cdot f \cdot {DCT}^{T} \]
def doDCT(grid):
return dot(dot(dct, grid), dct_transpose)
def undoDCT(grid):
return dot(dot(dct_transpose, grid), dct)
# test : do DCT, then undo DCT; should get back the same image.
tiny_do_undo = undoDCT(doDCT(tiny))
show_image(tiny_do_undo) # Yup, looks the same.
# And the numbers are the same.
tiny_do_undo
The DCT transform looks like this. Note that most of the intensity is at the top left, in the lowest frequencies.
tinyDCT = doDCT(tiny)
show_image(tinyDCT)
set_printoptions(linewidth=100) # output line width (default is 75)
round6 = vectorize(lambda m: '{:6.1f}'.format(m))
round6(tinyDCT)
The grid positions in that last image correspond to spatial frequencies, with the lowest DC component at the top left, and the highest vertical and horizontal frequency at the bottom right.
These 2D basis functions can be visualized with this image from wikimedia commons.
The details of what I'm doing here don't really match the JPEG transformations: I haven't done the color space transforms, and I haven't handled the DC offsets as the JPEG spec does (which centers the values around 0 explicitly.)
But the concept is visible in the last two pictures: after the DCT, most of the power is in fewer pixels, which are typically near the top left DC part of the grid.
So here's a simple lossy "low pass filter" of the data : let's chop some of the high frequency numbers. One (somewhat arbitrary) choice to to set the frequencies higher than the (1,7) to (7,1) line, to zero.
This is a lossy transormation since we're throwing away information - it can't be undone. But since there are fewer numbers, it's a form of compression.
# First make a copy to work on.
tinyDCT_chopped = tinyDCT.copy()
# Then zero the pieces below the x + y = 8 line.
for x in xrange(N):
for u in xrange(N):
if x + u > 8:
tinyDCT_chopped[x,u] = 0.0
show_image(tinyDCT_chopped)
round6(tinyDCT_chopped)
# Notice all the zeros at the bottom right - those are the chopped high frequences.
# We've essentially done a "low pass filter" on the spacial frequencies.
To see what this did to the original, we just transform it back.
tiny_chopped_float = undoDCT(tinyDCT_chopped)
# Also convert the floats back to uint8, which was the original format
tiny_chopped = vectorize(lambda x: uint8(x))(tiny_chopped_float)
show_image(tiny_chopped)
tiny_chopped
And we have something close to the original back again - even though close to half of the transformed image was set to zero.
The procedue here isn't what happens in JPEG compression, but does illustrate one of the central concepts - throwing away some of higher spatial frequency information after a DCT transform.
In the real JPEG lossy compression algorithm, the steps are
the color space is transformed from R,G,B to Y,Cb,Cr to take advantage of human visual prejudices
the values are shifted so that they center around zero
the values after the DCT are "quantized" (i.e. rounded off) by different amounts at different spots in the grid. (This* is the lossy step, and how lossy depends on the JPEG quality.)
a zigzag (keeping similar frequencies together) pattern turns this to a 1D stream of 64 values
which are then huffman encoded by, typically by a pre-chosen code (part of the JPEG standard)
For all the JPEG details, see http://en.wikipedia.org/wiki/JPEG .