Project 2: Fun with
Filters and Frequencies

E Harrison  —  CS180, Fall 2026

CS180 Projects  /  Home  /  Assignment

Overview

This project explores what we can do with 2D convolution and the frequency content of images. In Part 1 I implement convolution from scratch and use finite-difference and derivative-of-Gaussian filters to detect edges. In Part 2 I use low- and high-pass filters to sharpen images, build hybrid images whose interpretation changes with viewing distance, and blend images seamlessly with Gaussian and Laplacian stacks.

Part 1: Fun with Filters

1.1  Convolutions from Scratch

I wrote convolution two ways using only NumPy. The four-loop version iterates over every output pixel and, within that, over every kernel entry, accumulating the weighted sum by hand. The two-loop version keeps the loops over output pixels but replaces the inner two loops with a single dot product between the flattened kernel and the image patch beneath it.

Both versions flip the kernel along both axes first (that is what makes this a convolution rather than a cross-correlation) and support zero padding in "same" mode (output matches the input size, padding by (k−1)/2) and "full" mode (padding by k−1, output grows by k−1). With no padding they return the smaller “valid” output.

import numpy as np

def zero_padding(image, pad_x, pad_y):
    w, h = image.shape
    padded_image = np.zeros((w + 2 * pad_x, h + 2 * pad_y))
    padded_image[pad_x:pad_x + w, pad_y:pad_y + h] = image
    return padded_image

def add_padding(image, kernel, padding):
    kw, kh = kernel.shape
    if padding == "same":
        image = zero_padding(image, (kw - 1) // 2, (kh - 1) // 2)
    elif padding == "full":
        image = zero_padding(image, kw - 1, kh - 1)
    return image

def output_dims(image, kernel, padding):
    w, h = image.shape
    kw, kh = kernel.shape
    if padding == "same":
        return w, h
    elif padding == "full":
        return w + kw - 1, h + kh - 1
    else:
        return w - kw + 1, h - kh + 1

def convolve_4loop(image, kernel, padding=None):
    ow, oh = output_dims(image, kernel, padding)
    if padding is not None:
        image = add_padding(image, kernel, padding)
    kw, kh = kernel.shape
    kernel = np.flip(kernel, (0, 1))
    convolved = np.zeros((ow, oh))
    for i in range(ow):
        for j in range(oh):
            output = 0
            for x in range(kw):
                for y in range(kh):
                    output += kernel[x, y] * image[i + x, j + y]
            convolved[i, j] = output
    return convolved

def convolve_2loop(image, kernel, padding=None):
    ow, oh = output_dims(image, kernel, padding)
    if padding is not None:
        image = add_padding(image, kernel, padding)
    kw, kh = kernel.shape
    kernel = np.flip(kernel, (0, 1))
    convolved = np.zeros((ow, oh))
    for i in range(ow):
        for j in range(oh):
            kernel_flat = kernel.flatten()
            image_flat = image[i:i + kw, j:j + kh].flatten()
            convolved[i, j] = np.dot(kernel_flat, image_flat)
    return convolved

I tested all three implementations on a grayscale photo with a 9×9 box filter (every entry 1/81) and the finite-difference operators Dx = [1, 0, −1] and Dy = [1, 0, −1]T. The four-loop, two-loop, and scipy.signal.convolve2d outputs are identical, so only one set is shown below.

Runtime. I timed each implementation on the 542×540 cameraman image with the 9×9 box filter. The full-resolution photo above takes minutes with four loops.

Boundary handling. My implementation zero-pads the image, which is exactly what convolve2d does by default (boundary="fill", fillvalue=0). That is why the outputs agree to floating-point precision. Zero padding pulls black into the border pixels, so the box-filtered image has a slightly darkened frame.

1.2  Finite Difference Operator

Convolving the cameraman image with Dx and Dy gives the partial derivatives ∂I/∂x and ∂I/∂y, which respond to vertical and horizontal edges respectively. Combining them gives the gradient magnitude √((∂I/∂x)² + (∂I/∂y)²), which is large wherever intensity changes quickly in any direction. Thresholding the magnitude produces a binary edge map.

def edge_detection(image, threshold, dx_kernel=None, dy_kernel=None):
    if dx_kernel is None:
        dx_kernel = np.array([[1, 0, -1]])
    if dy_kernel is None:
        dy_kernel = np.array([[1], [0], [-1]])
    dx_image = signal.convolve2d(image, dx_kernel, mode="same")
    dy_image = signal.convolve2d(image, dy_kernel, mode="same")
    grad_mag_image = np.sqrt(np.square(dx_image) + np.square(dy_image))
    edge_image = np.where(grad_mag_image > threshold, 255.0, 0.0)
    return grad_mag_image, edge_image

grad_mag_image, edge_image = edge_detection(cameraman, threshold=90)

Partial derivatives are shown with 0 mapped to mid-gray: bright means a positive response and dark means negative.

Choosing the threshold. The threshold trades noise for missed edges. The grass in the lower half is full of small, high-frequency texture that produces many weak gradient responses. A low threshold keeps those responses and speckles the ground with false edges. A high threshold erases the thinner tripod legs and the faint buildings in the background. I chose 90 by eye. At that value the cameraman, camera, and tripod keep clear outlines, and most of the grass noise is suppressed.

1.3  Derivative of Gaussian (DoG) Filter

The finite-difference edges are noisy because differentiation amplifies high frequencies. The fix is to smooth first. I built a 2D Gaussian (ksize=9, σ=3) by taking the outer product of a 1D cv2.getGaussianKernel with itself, blurred the cameraman, and then repeated the gradient pipeline from 1.2.

def gaussian_kernel(ksize, sigma):
    gaussian_1d = cv2.getGaussianKernel(ksize=ksize, sigma=sigma)
    return gaussian_1d @ gaussian_1d.T

gaussian = gaussian_kernel(ksize=9, sigma=3)
dx_kernel, dy_kernel = np.array([[1, 0, -1]]), np.array([[1], [0], [-1]])

# (a) two convolutions: blur, then take finite differences
blurred_image = signal.convolve2d(image, gaussian, mode="same")
grad_mag, edges = edge_detection(blurred_image, threshold=20)

# (b) one convolution: fold the derivative into the Gaussian (DoG filters)
gaussian_dx_kernel = signal.convolve2d(gaussian, dx_kernel)
gaussian_dy_kernel = signal.convolve2d(gaussian, dy_kernel)
grad_mag_1conv, edges_1conv = edge_detection(
    image, threshold=20, dx_kernel=gaussian_dx_kernel, dy_kernel=gaussian_dy_kernel
)

(a) Blur, then differentiate

What differences do you see? Compared with the plain finite differences in 1.2, the edges are far cleaner. The grass noise is almost entirely gone, and the outlines are thicker, smoother, and more continuous because the blur spreads each edge over several pixels. The price is lost fine detail: small structures such as the thin tripod features and the distant buildings become weaker or merge together. Smoothing also lowers the overall gradient values, so the threshold had to drop from 90 to 20.

(b) A single convolution with DoG filters

Convolution is associative, so Dx * (G * I) = (Dx * G) * I. Convolving the Gaussian with Dx and Dy once produces two derivative of Gaussian filters that blur and differentiate in a single pass.

The DoG results match the two-step results. Away from the image border, the gradient magnitudes from (a) and (b) differ by at most about 10⁻¹³, which is floating-point round-off. Small differences can appear only within a few pixels of the frame edge, where zero padding is applied once in (b) but twice in (a).

Part 2: Fun with Frequencies

2.1  Image “Sharpening”

A Gaussian is a low-pass filter, so subtracting a blurred image from the original leaves only its high frequencies: the edges and fine texture. Adding a scaled copy of those high frequencies back makes the image look sharper:

sharpened = I + α(I − G * I) = ((1 + α)e − αG) * I

Here e is the unit impulse (identity kernel). The right-hand side shows that the whole operation collapses into a single convolution with the unsharp mask filter (1 + α)e − αG.

def unsharp_mask_filter(ksize, sigma, alpha):
    blur_kernel = gaussian_kernel(ksize=ksize, sigma=sigma)
    identity_kernel = np.zeros((ksize, ksize))
    identity_kernel[(ksize - 1) // 2, (ksize - 1) // 2] = 1.0
    return (1 + alpha) * identity_kernel - alpha * blur_kernel

unsharp_kernel = unsharp_mask_filter(ksize=15, sigma=5, alpha=3)
sharpened_image = signal.convolve2d(image, unsharp_kernel, mode="same")

Gaussian: ksize=15, σ=5. The high-frequency image is amplified 4× and offset to mid-gray for visibility.

Varying the sharpening amount. Increasing α boosts the high frequencies further. Small values make the carvings on the dome crisper. Large values over-sharpen: bright/dark halos appear around every edge (see the dome against the sky and the dark cypress trees), and noise in the sky starts to show.

Another image. Next I applied the unsharp mask to the photo from Part 1.1. This image is much larger than the Taj Mahal (1200 px wide as processed), so the ksize=15, σ=5 Gaussian only picks up very fine detail, and the sharpening is hard to see at normal viewing size. I increased the Gaussian to ksize=31, σ=10. That moves more mid-frequency structure into the “high-frequency” band, so faces, jersey lettering, and the folds of clothing get visibly more contrast.

Gaussian: ksize=31, σ=10. The high-frequency image is amplified 2× and offset to mid-gray for visibility.

As with the Taj Mahal, small α gives a subtle boost in clarity. By α = 4 the image is clearly over-processed: the white shoes and towels clip to pure white, the dark clothing clips to black, and bright halos surround the players.

Blur, then re-sharpen. To test whether sharpening actually restores lost detail, I took the sharp Taj Mahal, blurred it with the same Gaussian, and then applied the unsharp mask (α = 3) again.

The re-sharpened image is not the original. Re-sharpening restores contrast at large-scale edges such as the dome silhouette, the minarets, and the tree line, but the fine detail (the carvings and window tracery) does not come back. It also adds strong ringing halos around every edge. The σ = 5 blur drove those high frequencies to nearly zero, and sharpening can only amplify frequencies that are still present. So “sharpening” is a perceptual trick that raises edge contrast; it does not recover lost information.

2.2  Hybrid Images

A hybrid image combines the low frequencies of one image with the high frequencies of another. Up close, the high frequencies dominate what we see. From far away (or when the image is small), only the low frequencies survive, so the interpretation flips.

After aligning the two images on two hand-picked points, I low-pass the first image with a Gaussian of standard deviation σlow and high-pass the second with an impulse-minus-Gaussian filter of σhigh, then average the two.

def highpass_kernel(ksize, sigma):
    blur_kernel = gaussian_kernel(ksize=ksize, sigma=sigma)
    identity_kernel = np.zeros((ksize, ksize))
    identity_kernel[(ksize - 1) // 2, (ksize - 1) // 2] = 1.0
    return identity_kernel - blur_kernel

def hybrid_image(im1, im2, low_sig, high_sig, ksize=9):
    low_freq_kernel = gaussian_kernel(ksize=ksize, sigma=low_sig)
    high_freq_kernel = highpass_kernel(ksize=ksize, sigma=high_sig)
    low_freq_img = signal.convolve2d(im1, low_freq_kernel, mode="same")
    high_freq_img = signal.convolve2d(im2, high_freq_kernel, mode="same")

    w = min(low_freq_img.shape[0], high_freq_img.shape[0])
    h = min(low_freq_img.shape[1], high_freq_img.shape[1])
    low_freq_img = crop_image(low_freq_img, w, h)
    high_freq_img = crop_image(high_freq_img, w, h)
    return (low_freq_img + high_freq_img) / 2

def fourier(img):
    return np.log(np.abs(np.fft.fftshift(np.fft.fft2(img))))

Full process: Oski + Stanford Tree. Up close you see the Stanford Tree; step back and it turns into Oski. Oski supplies the low frequencies and the Tree supplies the high frequencies, both with σ = 4 (ksize = 9).

The log-magnitude spectra show the filtering directly. Low-passing Oski dims everything away from the center of the spectrum, keeping only the coarse shape. High-passing the Tree removes the bright central peak while keeping the outer frequencies. The hybrid's spectrum is the sum of the two: a strong center from Oski surrounded by the Tree's high-frequency energy.

Choosing the cutoffs. A Gaussian with spatial standard deviation σ has a frequency-domain standard deviation of roughly 1/(2πσ) cycles per pixel, so a larger σ means a lower cutoff. I tuned σ_low and σ_high for each pair by eye. σ_low has to be large enough that the low-frequency image disappears up close, but not so large that it is unrecognizable from a distance. σ_high has to be large enough that the high-pass image keeps enough structure to be readable, while still removing its low-frequency shading, which would otherwise compete with the other image. Larger, higher-resolution images such as Derek + Nutmeg needed larger σ values (and a larger kernel) than the small Oski and Tree images.

More hybrids. In each row, the left image provides the low frequencies and the middle image provides the high frequencies.

2.3  Gaussian and Laplacian Stacks

A Gaussian stack is like a Gaussian pyramid without the downsampling. Each level is the previous level blurred again, so every level keeps the original resolution. A Laplacian stack takes the difference of consecutive Gaussian levels, which isolates one frequency band per level. The last Gaussian level is kept as a low-frequency residual, so summing all Laplacian levels plus the residual reconstructs the original image exactly. Both are implemented from scratch with convolve2d on each color channel.

def gaussian_stack(img, sigma, stack_height):
    w, h, c = img.shape
    stack = np.zeros((stack_height, w, h, c))
    stack[0] = img
    for i in range(1, stack_height):
        for j in range(c):
            layer_sigma = 2 ** (i - 1) * sigma   # sigma doubles at each level
            kernel = gaussian_kernel(ksize=layer_sigma * 3 + 1, sigma=layer_sigma)
            stack[i, ..., j] = signal.convolve2d(stack[i-1, ..., j], kernel, mode="same")
    return stack

def laplacian_stack(img, sigma, stack_height):
    gauss_stack = gaussian_stack(img, sigma, stack_height + 1)
    stack = gauss_stack[:-1] - gauss_stack[1:]
    return stack, gauss_stack[-1]   # band-pass levels + low-frequency residual

Below is my recreation of Figure 3.42 from Szeliski. I built 5-level Laplacian stacks of the apple and orange and a Gaussian stack of a vertical step mask, all starting at σ=3. The blur doubles at each level (σi = 2i−1 · 3, with ksize = 3σi + 1), so each level captures a band roughly an octave lower than the one before. Each row shows level i of the apple weighted by mask level i, the orange weighted by 1 − mask level i, and their sum. The bottom row shows the collapsed stacks. The mask becomes blurrier at each level, so low frequencies are blended over a wide seam and high frequencies over a narrow one.

Laplacian levels are min-max normalized to [0, 255] so the zero-centered band-pass detail is visible. The mask stack is built at twice the image size and center-cropped so zero padding at the borders doesn't bleed into the mask.

2.4  Multiresolution Blending

Following Burt and Adelson (1983), I blend two images by combining their Laplacian stacks band by band, using a Gaussian stack of the mask as the per-band weights:

LSi = GMi · LAi + (1 − GMi) · LBi

Summing the blended levels (plus the blended low-frequency residuals) gives the final image. The mask can be any shape: a vertical step for the oraple, or a hand-drawn irregular region.

def blended_image(im1, im2, mask, sigma, stack_height):
    lstack1, res1 = laplacian_stack(im1, sigma, stack_height)
    lstack2, res2 = laplacian_stack(im2, sigma, stack_height)
    mstack = gaussian_stack(mask, sigma, stack_height + 1)

    # start with the blended low-frequency residuals
    m = mstack[-1]
    blended = m * res1 + (1 - m) * res2

    # add each blended band-pass level
    for i in range(stack_height):
        m = mstack[i]
        blended += m * lstack1[i] + (1 - m) * lstack2[i]
    return blended

The Oraple. With a vertical step mask, simply pasting the halves together leaves a hard, obvious seam. The multiresolution blend (5 levels, σ=3) hides it: the colors fade smoothly while the textures still change over a narrow band.

XStation 5. Here I merge an Xbox controller into a PS5 controller with an irregular mask. The mask selects the Xbox's face buttons, thumbsticks, and grips and places them on the PS5 body. The blend smooths the differences in shading and white balance between the two controllers along the seam.

Favorite: the LEGO Skyline. My favorite result replaces the front of a real Nissan Skyline with the front of a LEGO Skyline, using a hand-drawn, irregular mask around the hood, bumper, and windshield. With a hard paste, the LEGO section looks cut out: its lighting and color don't match the real car. The 5-level blend softens that transition into the real car's paint.

The figure below mirrors Figure 10 of the Burt and Adelson paper. Each row is one Laplacian level: the masked LEGO band, the masked real-car band, and their sum. The bottom row collapses each column into the masked input images and the final blend.

Website template from Jon Barron.