Implementing 2D Cross-Correlation From Scratch in NumPy

Implementing 2D Cross-Correlation From Scratch in NumPy#

Twitter Handle LinkedIn Profile GitHub Profile Tag Tag

The cleanest way to confirm you understand an operator is to write it in twenty lines. This page does exactly that: it re-implements the author’s two helper functions — cross_correlation_2d and calculate_feature_map_shape — in idiomatic, fully type-annotated NumPy, then runs them on a small numeric case and on a real edge-detection kernel.

It is for you if you have read the concept chapter and want the sliding-window operator under your fingers. By the end you will have a cross_correlation_2d you can step through line by line, a calculate_feature_map_shape that specialises the general output-shape formula (162) to stride \(1\) and zero padding, and a printed feature map showing a Sobel-style vertical-edge detector firing exactly on the edge — the whole story of the intro in twenty executable lines.

Prerequisites

This page implements, as code, the operator defined in the concept chapter; you will get more out of the implementation if you have read the definition and the output-shape formula there first. The notation follows the deep-learning notation page.

The output shape, straight from the formula#

The concept chapter gives the general shape rule as \(H_{out} = \lfloor (H_{in} + 2P - K_h) / S \rfloor + 1\). The author’s helper covers the common case this chapter starts from — no padding, unit stride — so we set \(P = 0\) and \(S = 1\) and the floor disappears:

(165)#\[ H_{out} = H_{in} - K_h + 1, \qquad W_{out} = W_{in} - K_w + 1. \]

That is the whole of calculate_feature_map_shape. It is deliberately tiny: it exists to keep the sliding window in cross_correlation_2d from recomputing the output size by hand.

 1def calculate_feature_map_shape(
 2    x: NDArray[np.floating],
 3    kernel: NDArray[np.floating],
 4) -> tuple[int, int]:
 5    """Feature-map shape for stride 1 and zero padding.
 6
 7    Specialises the general formula
 8    ``(H_in + 2 * P - K_h) // S + 1`` to ``P = 0`` and ``S = 1`` — the
 9    convention this chapter uses until stride and padding are introduced.
10
11    Args:
12        x: Input grid of shape ``(H_in, W_in)``.
13        kernel: Kernel grid of shape ``(K_h, K_w)``.
14
15    Returns:
16        The output shape ``(H_in - K_h + 1, W_in - K_w + 1)``.
17    """
18    h_in, w_in = x.shape
19    k_h, k_w = kernel.shape
20    return h_in - k_h + 1, w_in - k_w + 1

The sliding-window operator#

With the shape in hand, the operator writes itself. Allocate an output grid of that shape, then for every anchor position slice the corresponding patch out of the input and write its dot product with the kernel into one cell. This is Definition 79 translated line for line from mathematics into NumPy — no flipping, no padding, unit stride, exactly as deep learning means when it says “convolution”.

 1def cross_correlation_2d(
 2    x: NDArray[np.floating],
 3    kernel: NDArray[np.floating],
 4) -> NDArray[np.floating]:
 5    """Sliding-window 2D cross-correlation (the ``Conv2d`` operator).
 6
 7    Slide ``kernel`` across ``x`` in unit strides with no padding. At each
 8    anchor position, return the dot product of the kernel with the
 9    ``K_h x K_w`` patch of ``x`` it covers.
10
11    Args:
12        x: Input grid of shape ``(H_in, W_in)``.
13        kernel: Kernel grid of shape ``(K_h, K_w)``.
14
15    Returns:
16        Feature map of shape ``(H_in - K_h + 1, W_in - K_w + 1)``.
17    """
18    out_h, out_w = calculate_feature_map_shape(x, kernel)
19    k_h, k_w = kernel.shape
20    out = np.zeros((out_h, out_w), dtype=x.dtype)
21    for i in range(out_h):
22        for j in range(out_w):
23            patch = x[i : i + k_h, j : j + k_w]   # K_h x K_w window
24            out[i, j] = np.sum(patch * kernel)    # scalar dot product
25    return out

A worked numeric call#

Before trusting the function on anything interesting, reproduce Example 31 from the concept chapter: the \(5 \times 5\) ramp under a \(3 \times 3\) all-ones kernel should return the \(3 \times 3\) feature map whose top-left entry is \(63\) and whose centre entry is \(117\). If those two match, the sliding window is anchored correctly.

 1x = np.array(
 2    [
 3        [1.0, 2.0, 3.0, 4.0, 5.0],
 4        [6.0, 7.0, 8.0, 9.0, 10.0],
 5        [11.0, 12.0, 13.0, 14.0, 15.0],
 6        [16.0, 17.0, 18.0, 19.0, 20.0],
 7        [21.0, 22.0, 23.0, 24.0, 25.0],
 8    ]
 9)
10box = np.ones((3, 3))
11
12feature_map = cross_correlation_2d(x, box)
13print("shape:  ", feature_map.shape)      # (3, 3)
14print("Y[0,0]: ", feature_map[0, 0])      # 63.0
15print("Y[1,1]: ", feature_map[1, 1])      # 117.0
16print(feature_map)

shape:   (3, 3)
Y[0,0]:  63.0
Y[1,1]:  117.0
[[ 63.  72.  81.]
 [108. 117. 126.]
 [153. 162. 171.]]

An edge-detection demo#

Counting sums is fine for checking the anchor logic, but it hides what kernels are for. A more honest demo is an edge detector: a kernel whose dot product is large exactly when the patch straddles a boundary and near zero over flat regions. Build a \(6 \times 6\) image that is bright (\(10\)) on the left half and dark (\(0\)) on the right — a single vertical edge down the middle — and slide a classic vertical Sobel-style filter over it ([Zhang et al., 2023], Ch. 7.1, discusses exactly this family of hand-crafted detectors before the network learns its own).

 1# A 6x6 image: bright block on the left, dark block on the right.
 2image = np.zeros((6, 6), dtype=np.float64)
 3image[:, :3] = 10.0
 4
 5# Vertical-edge detector: +1 on the left column, -1 on the right.
 6vertical_edge = np.array(
 7    [
 8        [1.0, 0.0, -1.0],
 9        [1.0, 0.0, -1.0],
10        [1.0, 0.0, -1.0],
11    ]
12)
13
14edge_map = cross_correlation_2d(image, vertical_edge)
15print("feature-map shape:", edge_map.shape)   # (4, 4)
16print(edge_map)

feature-map shape: (4, 4)
[[ 0. 30. 30.  0.]
 [ 0. 30. 30.  0.]
 [ 0. 30. 30.  0.]
 [ 0. 30. 30.  0.]]

The feature map is \(4 \times 4\), exactly as (165) predicts (\((6 - 3 + 1) \times (6 - 3 + 1)\)). Its two middle columns read \(30\) and the outer columns read \(0\): the detector fires only where its \(3 \times 3\) window straddles the boundary between the bright and dark halves, and stays silent over the flat regions on either side. That is the whole job of a convolutional layer in one picture — light up where the pattern is present, stay dark where it is not — except that in a real network the kernel entries are learned from data rather than hand-set as they are here.

What this page does not do is multi-channel input, padding, stride, or backpropagation. Each is a small extension of the loop above: multi-channel sums one cross-correlation per input channel (see (163)), padding surrounds x with zeros before the loop, stride skips indices in the range, and backpropagation is the same dot product run in reverse. For all of those at production speed, use torch.nn.Conv2d — but now you know what it is doing underneath.

Summary

If this page had to be one sentence: 2D cross-correlation is twenty lines of NumPy — one shape helper and a double loop that takes a patch–kernel dot product at every anchor — and a hand-set Sobel filter on a synthetic edge shows exactly the “fire on the pattern, stay dark elsewhere” behaviour that a learned kernel generalises.

  • calculate_feature_map_shape specialises the general output-shape formula to stride \(1\) and zero padding: \((H - K + 1) \times (W - K + 1)\).

  • cross_correlation_2d is Definition 79 translated line for line: slice a patch, multiply, sum, advance.

  • The Sobel demo fires (\(30\)) only where the window straddles the edge; learned kernels do the same thing, but discover their weights from data.

For the formal operator, see the concept chapter; for the three priors that make this worth doing, see the intro.

Further reading