Lesson 17: Image Compression — Lossless and Lossy¶
Storing an image as raw pixels costs 1 byte per channel per pixel, no matter what the image contains. Compression exploits various types of redundancy in images --- i.e., the fact that real images are not random. This lesson covers both lossless compression (the decoded image is bit-for-bit identical to the original — e.g., PNG) and lossy compression (the decoded image is only an approximation — e.g., JPEG).
import numpy as np
import cv2
import heapq
from collections import Counter
import matplotlib.pyplot as plt
Lossless compression¶
Run-length encoding (RLE)¶
The simplest lossless scheme: instead of storing every pixel, store (value, run length) pairs for consecutive runs of identical pixels. This works great on images with large flat regions (synthetic graphics, scanned text, masks), but it does nothing useful on real photographs.
def rle_encode(array):
flat = array.ravel()
change_points = np.where(np.diff(flat) != 0)[0] + 1
starts = np.concatenate([[0], change_points])
ends = np.concatenate([change_points, [len(flat)]])
return [(flat[s], e - s) for s, e in zip(starts, ends)]
def rle_decode(runs, shape):
flat = np.concatenate([np.full(length, value) for value, length in runs])
return flat.reshape(shape).astype(np.uint8)
flat_regions_img = np.zeros((80, 80), dtype=np.uint8)
cv2.rectangle(flat_regions_img, (10, 10), (70, 70), 200, -1)
sp = flat_regions_img.shape
plt.imshow(flat_regions_img, cmap='gray')
plt.title(f'Image with flat regions ({sp[1]}x{sp[0]}) = {sp[0]*sp[1]} bytes')
plt.axis('off')
plt.show()
runs = rle_encode(flat_regions_img)
decoded = rle_decode(runs, flat_regions_img.shape)
print('exact reconstruction?', np.array_equal(decoded, flat_regions_img))
raw_bytes = flat_regions_img.size
rle_bytes = len(runs) * 3 # roughly: 1 byte value + 2 bytes run length, per run
print(f'raw size: {raw_bytes} bytes')
print(f'RLE size: ~{rle_bytes} bytes ({raw_bytes / rle_bytes:.1f}x smaller)')
Huffman coding: fewer bits for common values¶
A more sophisticated lossless method is Huffman coding, which takes advantage of the fact that some pixel values (e.g., a common background gray) occur far more often than others. It builds a variable-length code where frequent values get short codes and rare values get longer ones — which is provably optimal among prefix codes.
Shannon's entropy gives the theoretical floor on average bits/pixel for any code based only on the value distribution (ignoring spatial structure):
$$H = -\sum_i p_i \log_2 p_i$$
where $p_i$ is the probability that a certain pixel value occurs in the image, and the summation is over all possible pixel values.
rng = np.random.default_rng(0)
photo_like = np.zeros((100, 100), dtype=np.uint8)
cv2.rectangle(photo_like, (10, 10), (90, 90), 200, -1)
cv2.circle(photo_like, (50, 50), 25, 120, -1)
photo_like = np.clip(photo_like.astype(np.float64) + rng.normal(0, 5, photo_like.shape), 0, 255).astype(np.uint8)
sp = photo_like.shape
counts = Counter(photo_like.ravel().tolist())
total = sum(counts.values())
probs = {value: count / total for value, count in counts.items()}
entropy = -sum(p * np.log2(p) for p in probs.values())
def build_huffman_codes(probs):
heap = [[p, [symbol, '']] for symbol, p in probs.items()]
heapq.heapify(heap)
while len(heap) > 1:
lo = heapq.heappop(heap)
hi = heapq.heappop(heap)
for pair in lo[1:]:
pair[1] = '0' + pair[1]
for pair in hi[1:]:
pair[1] = '1' + pair[1]
heapq.heappush(heap, [lo[0] + hi[0]] + lo[1:] + hi[1:])
return {symbol: code for symbol, code in heap[0][1:]}
codes = build_huffman_codes(probs)
avg_bits = sum(probs[symbol] * len(code) for symbol, code in codes.items())
plt.imshow(photo_like, cmap='gray')
plt.title(f'Image with noise ({sp[1]}x{sp[0]}) = {sp[0]*sp[1]} bytes')
plt.axis('off')
plt.show()
print(f'naive fixed-width encoding: 8.00 bits/pixel')
print(f'Shannon entropy (theoretical floor): {entropy:.2f} bits/pixel')
print(f'Huffman coding achieves: {avg_bits:.2f} bits/pixel')
Huffman coding gets close to the entropy bound (it can only match it exactly when every probability happens to be a power of 2). Real lossless formats like PNG combine an idea like this (Huffman/arithmetic coding) with a predictive filter first — predicting each pixel from its neighbors and encoding only the (usually small) prediction error, which has much lower entropy than the raw pixel values.
Lossy compression: the DCT¶
JPEG's core idea is to take advantage of the human visual system's tendency to overlook fine details. JPEG divides the image into small 8x8 blocks, then transforms each block into the frequency domain so that it can discard the high-frequency components that are barely noticeable to the human eye. JPEG uses the discrete cosine transform (DCT), which is a close relative of the Fourier transform from Lesson 15. Like the DFT, the DCT re-expresses a block as a sum of frequency components — but it uses only cosines (no imaginary part) to avoid boundary artifacts.
block = photo_like[6:14, 6:14].astype(np.float64) # a block straddling the rectangle's sharp edge
dct_block = cv2.dct(block)
energy = dct_block**2
print(f'fraction of block energy in top-left 2x2 coefficients: {energy[:2, :2].sum() / energy.sum():.2%}')
print(f'fraction of block energy in top-left 4x4 coefficients: {energy[:4, :4].sum() / energy.sum():.2%}')
fig, axes = plt.subplots(1, 2, figsize=(6, 3))
axes[0].imshow(block, cmap='gray')
axes[0].set_title('8x8 pixel block')
axes[1].imshow(np.log1p(np.abs(dct_block)), cmap='gray')
axes[1].set_title('log|DCT coefficients|\n(top-left = low frequency)')
for ax in axes:
ax.axis('off')
plt.tight_layout()
plt.show()
Quantization: throwing away the coefficients that matter least¶
Compression happens by quantizing the DCT coefficients — dividing by some value and rounding. The values to be divided by are taken from a quantization table, which are designed based on human visual sensitivity. Many coefficients (especially high-frequency ones, which carry the least energy) become exactly zero and take almost no space to store.
def keep_top_left(dct_block, n):
"""Zero out every coefficient outside the top-left nxn corner (crude stand-in for JPEG-style quantization)."""
masked = np.zeros_like(dct_block)
masked[:n, :n] = dct_block[:n, :n]
return masked
fig, axes = plt.subplots(1, 4, figsize=(10, 3))
for ax, n in zip(axes, [8, 4, 2, 1]):
kept = keep_top_left(dct_block, n)
reconstructed = cv2.idct(kept)
nonzero = np.count_nonzero(kept)
ax.imshow(reconstructed, cmap='gray', vmin=block.min(), vmax=block.max())
ax.set_title(f'{nonzero}/64 coeffs kept', fontsize=9)
ax.axis('off')
plt.tight_layout()
plt.show()
Even keeping just the top-left 2x2 (4 out of 64 coefficients — a 16x reduction) preserves most of the block's overall structure; only the finest detail is lost.
Real JPEG: the rate-distortion tradeoff¶
Real JPEG encoding follows this same DCT-then-quantize recipe per block (with a carefully designed, frequency-dependent quantization table, plus run-length and Huffman coding of the resulting sparse coefficients — combining the lossy and lossless ideas from this lesson). We use OpenCV's built-in encoder directly, sweeping the quality setting, and measure both file size and reconstruction error.
def psnr(original, reconstructed):
mse = np.mean((original.astype(np.float64) - reconstructed.astype(np.float64))**2)
return float('inf') if mse == 0 else 10 * np.log10(255**2 / mse)
color_img = np.zeros((200, 200, 3), dtype=np.uint8)
cv2.rectangle(color_img, (20, 20), (180, 180), (255, 120, 30), -1)
cv2.circle(color_img, (100, 100), 60, (30, 200, 255), -1)
color_img = np.clip(color_img.astype(np.float64) + rng.normal(0, 8, color_img.shape), 0, 255).astype(np.uint8)
sp = color_img.shape
qualities = [10, 30, 50, 80, 95, 100]
sizes, psnrs = [], []
for q in qualities:
ok, encoded = cv2.imencode('.jpg', color_img, [cv2.IMWRITE_JPEG_QUALITY, q])
decoded = cv2.imdecode(encoded, cv2.IMREAD_COLOR)
sizes.append(len(encoded))
psnrs.append(psnr(color_img, decoded))
plt.imshow(color_img)
plt.title(f'Image with noise ({sp[1]}x{sp[0]}) = {sp[0]*sp[1]*sp[2]} bytes')
plt.axis('off')
plt.show()
raw_size = color_img.size
print(f'{"quality":>8} {"file size (bytes)":>18} {"compression ratio":>18} {"PSNR (dB)":>10}')
for q, s, p in zip(qualities, sizes, psnrs):
print(f'{q:>8} {s:>18} {raw_size / s:>17.1f}x {p:>10.2f}')
fig, axes = plt.subplots(1, 2, figsize=(9, 3.5))
axes[0].plot(qualities, sizes, marker='o')
axes[0].set_xlabel('JPEG quality')
axes[0].set_ylabel('file size (bytes)')
axes[0].set_title('Size vs. quality')
axes[1].plot(sizes, psnrs, marker='o')
axes[1].set_xlabel('file size (bytes)')
axes[1].set_ylabel('PSNR (dB)')
axes[1].set_title('Rate-distortion curve')
plt.tight_layout()
plt.show()
This is the fundamental tradeoff of lossy compression: every extra byte you're willing to spend buys diminishing returns in quality (the rate-distortion curve flattens out at high quality), and there is no free lunch — only a choice of where on the curve to sit. For a typical photograph, JPEG quality of 75 will not be noticeable unless you zoom in, and it will reduce the file size by 6--8x.
JPEG or PNG? Choosing the right format¶
The trade-offs above translate into a practical rule:
Use JPEG for photographs. Photographs are dominated by smooth gradients — exactly what the DCT concentrates into a few low-frequency coefficients, so discarding the rest is barely visible and JPEG buys a large size reduction for a small perceptual loss.
Use PNG for graphics (screenshots, diagrams, text, icons, anything with flat regions or sharp edges). Graphics are the opposite case: a sharp edge spreads energy across all DCT frequencies, so the same quantization that photos tolerate shows up as visible ringing and blocking around edges and text. Lossless compression handles flat regions and sharp edges almost for free (that's exactly what RLE and the predictive filter above are good at), so PNG is often both artifact-free and smaller than JPEG for this kind of content.
photo = cv2.cvtColor(cv2.imread('../img/house.png'), cv2.COLOR_BGR2RGB)
graphic = np.full((100, 200, 3), 255, dtype=np.uint8)
cv2.rectangle(graphic, (20, 20), (180, 100), (40, 40, 40), -1)
cv2.putText(graphic, 'REPORT', (30, 70), cv2.FONT_HERSHEY_SIMPLEX, 1.1, (255, 255, 255), 2)
cv2.line(graphic, (30, 80), (160, 80), (255, 255, 255), 2)
def encode_sizes(image, jpeg_quality=75):
ok, jpg = cv2.imencode('.jpg', image, [cv2.IMWRITE_JPEG_QUALITY, jpeg_quality])
ok, png = cv2.imencode('.png', image, [cv2.IMWRITE_PNG_COMPRESSION, 9])
return jpg, png
jpg_photo, png_photo = encode_sizes(photo)
jpg_graphic, png_graphic = encode_sizes(graphic)
print(f'{"":12} {"Original":>11} {"PNG":>11} {"JPEG (q75)":>11}')
print(f'{"photograph":12} {int(np.prod(photo.shape)):>9} B {len(png_photo):>9} B {len(jpg_photo):>9} B (PNG is bigger than JPEG by {len(png_photo)/len(jpg_photo):.1f}x)')
print(f'{"graphic":12} {int(np.prod(graphic.shape)):>9} B {len(png_graphic):>9} B {len(jpg_graphic):>9} B (PNG is smaller than JPEG by {len(png_graphic)/len(jpg_graphic):.1f}x)')
fig, axes = plt.subplots(2, 2, figsize=(7, 7))
axes[0, 0].imshow(photo); axes[0, 0].set_title(f'Photo: original ({int(np.prod(photo.shape))} B)')
axes[0, 1].imshow(cv2.imdecode(jpg_photo, cv2.IMREAD_COLOR)); axes[0, 1].set_title(f'Photo: JPEG q75 ({len(jpg_photo)} B)')
axes[1, 0].imshow(graphic); axes[1, 0].set_title(f'Graphic: original ({int(np.prod(graphic.shape))} B)')
axes[1, 1].imshow(cv2.imdecode(jpg_graphic, cv2.IMREAD_COLOR)); axes[1, 1].set_title(f'Graphic: JPEG q75 ({len(jpg_graphic)} B)')
for ax in axes.ravel():
ax.axis('off')
plt.tight_layout()
plt.show()
Photo by Peter Herrmann on Unsplash
A closer look at JPEG artifacts¶
Zooming into a high-contrast region (the roofline against the sky) at a much lower quality setting makes the cost of aggressive quantization visible directly: blocky 8x8 squares and ringing around the sharp edge, exactly the failure mode described above.
y0, y1, x0, x1 = 10, 60, 50, 150 # roofline against the sky
low_q = 10
jpg_low_q, _ = encode_sizes(photo, jpeg_quality=low_q)
decoded_low_q = cv2.imdecode(jpg_low_q, cv2.IMREAD_COLOR)
zoom = 4
original_crop = cv2.resize(photo[y0:y1, x0:x1], None, fx=zoom, fy=zoom, interpolation=cv2.INTER_NEAREST)
jpeg_crop = cv2.resize(decoded_low_q[y0:y1, x0:x1], None, fx=zoom, fy=zoom, interpolation=cv2.INTER_NEAREST)
fig, axes = plt.subplots(1, 2, figsize=(8, 4))
axes[0].imshow(original_crop)
axes[0].set_title('Original crop (zoomed 4x)')
axes[1].imshow(jpeg_crop)
axes[1].set_title(f'Same crop, JPEG q{low_q}\n(full image: {len(jpg_low_q)} B)')
for ax in axes:
ax.axis('off')
plt.tight_layout()
plt.show()
Exercise¶
- Run the RLE experiment on
photo_like(which has Gaussian noise) instead offlat_regions_img. How many runs does it produce, and what does that say about RLE's suitability for noisy natural images? keep_top_leftis a deliberately crude stand-in for real quantization (which shrinks all coefficients gradually rather than hard-zeroing a block of them). Replace it withnp.round(dct_block / step) * stepfor a few differentstepvalues, and compare the visual artifacts to the top-left-corner version.- At very low JPEG quality, you may notice a visible grid pattern aligned with the 8x8 block boundaries ("blocking artifacts"). Given what block boundaries have to do with the DCT, explain why quantizing each block independently would produce exactly this kind of artifact.