Block Matching

Block matching estimates motion between two frames by cutting the current frame into rectangular blocks and finding, for each block, the displacement at which it best matches a block of a reference frame. Each block gets one motion vector. It is the simplest matching form of motion estimation, and it is how video encoders find the vectors for motion-compensated prediction: the encoder sends each block’s vector and the difference between the block and its displaced match, instead of the block itself. The same search, with the displacement restricted to a row, is the local core of stereo matching. Block matching is a form of template matching in which every block of one frame serves as a template searched for in a small window of another.

The algorithm has three ingredients: a criterion that scores a candidate match, a set of candidate displacements, and a strategy for searching that set. This article covers the criteria and block size, exhaustive search and the fast search patterns that video coding produced (three-step, new three-step, four-step, diamond and hexagon), the assumption those patterns share, sub-pixel refinement, and the role of block matching in the ITU-T video coding standards from H.261 to H.265 (HEVC).

Problem

Given a current frame $I_t$ and a reference frame $I_r$, usually the previous frame, find for each $B \times B$ block of $I_t$ the displacement $\mathbf{v} = (v_x, v_y)$ such that the block of $I_r$ at that offset looks most like it. The answer is a field of vectors, one per block, sampled on the block grid, rather than the per-pixel field of optical flow estimation.

What “best” means depends on the use. A video encoder wants the vector that makes the prediction residual cheapest to code, whether or not it follows the true motion; a vision system wants the vector of the physical motion. Motion Estimation notes the difference; this article returns to it under Properties and Behavior.

Inputs and Outputs

• Input: a current frame, one or more reference frames, a block size $B$ (or a set of block shapes), a search range $p$, a matching criterion and a search strategy.
• Output: one motion vector per block, in whole pixels or in fractions of a pixel after refinement, and usually the matching cost at that vector, which an encoder uses to choose between coding modes.

Intuition

Slide the block over the reference frame within a window of $\pm p$ pixels, compute a difference score at every offset, and keep the offset with the lowest score. If the image content inside the block moved rigidly by $\mathbf{v}$ between the frames, the difference at $\mathbf{v}$ is only noise, and everywhere else it is larger.

The cost of this exhaustive search grows with the square of the search range, and motion estimation is the most time-consuming part of a video encoder (Chen et al., 2002). Fast searches test a few positions, move toward the best one, and repeat with a finer step. They work when the score decreases steadily toward the correct offset, like walking downhill in a single valley. Repeated texture and flat regions break that picture: the score surface then has several valleys or none.

Mathematical Formulation

Let $\mathcal{B}$ be the set of pixel positions $\mathbf{x}$ in a block of $I_t$. A candidate vector $\mathbf{v}$ predicts each pixel $I_t(\mathbf{x})$ from $I_r(\mathbf{x} + \mathbf{v})$; the difference is the displaced frame difference. The usual criteria compare the block with its candidate match (the similarity and distance measures article treats the underlying measures as distances between vectors):

$$\mathrm{SAD}(\mathbf{v}) = \sum_{\mathbf{x} \in \mathcal{B}} \left| I_t(\mathbf{x}) - I_r(\mathbf{x} + \mathbf{v}) \right|, \qquad \mathrm{SSD}(\mathbf{v}) = \sum_{\mathbf{x} \in \mathcal{B}} \left( I_t(\mathbf{x}) - I_r(\mathbf{x} + \mathbf{v}) \right)^2,$$ $$\mathrm{NCC}(\mathbf{v}) = \frac{\sum_{\mathbf{x}} \left( I_t(\mathbf{x}) - \bar{I}_t \right) \left( I_r(\mathbf{x} + \mathbf{v}) - \bar{I}_r(\mathbf{v}) \right)}{\sqrt{\sum_{\mathbf{x}} \left( I_t(\mathbf{x}) - \bar{I}_t \right)^2} \sqrt{\sum_{\mathbf{x}} \left( I_r(\mathbf{x} + \mathbf{v}) - \bar{I}_r(\mathbf{v}) \right)^2}},$$

where $\bar{I}_t$ and $\bar{I}_r(\mathbf{v})$ are the means of the two blocks. The estimate is

$$\hat{\mathbf{v}} = \arg\min_{\mathbf{v} \in \mathcal{W}} \mathrm{SAD}(\mathbf{v}), \qquad \mathcal{W} = \{-p, \dots, p\}^2,$$

or the minimum of SSD, or the maximum of NCC.

The criteria differ in cost and in what they tolerate:

• SAD, the $\ell_1$ distance between the blocks, needs only subtractions, absolute values and additions. Szeliski (2022, §9.1) notes that video coding uses it for its speed, and that it is unsuited to gradient-based refinement because it is not differentiable at zero.
• SSD, the squared $\ell_2$ distance, is the least-squares criterion, optimal for Gaussian noise, and differentiable, which is what Lucas–Kanade refinement minimizes. It weights large differences more heavily than SAD, so a few outlier pixels, such as an occluding edge crossing the block, count for more.
• NCC subtracts each block’s mean and divides by the product of their norms, so it is unchanged by a gain and an offset in brightness. Szeliski (2022, §9.1) notes that it lies in $[-1, 1]$, which makes scores comparable across blocks, that it is undefined for a block of zero variance, and that it degrades in noisy, low-contrast regions.

Video encoders add two refinements. The H.264 and HEVC reference encoders can measure distortion at fractional positions with the sum of absolute transformed differences (SATD), a SAD computed after a Hadamard transform (the default in JM; in HM when its Hadamard flag is on), because the residual is transform-coded and SATD estimates its coding cost better (Kim et al., 2013, §4.1.3; Tourapis et al., 2009). And they minimize a Lagrangian cost rather than the distortion alone. The HEVC test model, HM10, chooses motion parameters by

$$J = D(\mathbf{v}) + \lambda \, R(\mathbf{v}),$$

with $D$ the SAD for the integer search and SA(T)D for the fractional search, $R$ the bits needed to code the vector difference and the related indices, and $\lambda$ a multiplier derived from the quantization parameter (Kim et al., 2013, §4.1.4, §4.6.2.2). The rate term favours vectors close to the predicted vector, which is cheap to code.

The unimodal assumption

Exhaustive search finds the minimum of the criterion over $\mathcal{W}$ by construction. Fast searches find it only if the criterion has one minimum in the window and grows with distance from it. Chen et al. (2002) state the assumption directly: fast integer-pel searches assume that the matching error surface in the search window is unimodal, decreasing monotonically as the tested point approaches the global optimum. They call this assumption far from reality for complex video and show an integer-pel error surface on which a search is easily trapped in a local minimum, while finding the assumption usually true in the much smaller fractional-pel window around the best integer vector.

Algorithm

Exhaustive search

Exhaustive search, also called full search, evaluates every candidate:

1. For each block of the current frame, at position $(x_0, y_0)$:
2. For every $\mathbf{v} \in \{-p, \dots, p\}^2$ whose block lies inside the reference frame, compute the criterion between the block and the reference block at $(x_0 + v_x, y_0 + v_y)$.
3. Keep the $\mathbf{v}$ with the best score. Ties need a rule: FFmpeg’s exhaustive search keeps the zero vector if nothing beats it strictly, and otherwise the first best candidate in raster order.
4. Optionally refine the result to sub-pixel precision.

With $p = 7$ there are $15^2 = 225$ candidates per block; with $p = 16$, the range Szeliski (2022, §9.1.1) gives as an example for video compression, there are 1,089.

Fast search patterns

Fast searches replace the full grid with a sequence of small patterns. The original papers of most of these methods were published by IEEE and could not be read for this article. The patterns are described here as implemented in FFmpeg’s motion estimation code (libavfilter/motion_estimation.c, release 7.1), which its documentation lists under the methods’ names, and the attributions follow the reference list of Chen et al. (2002), except for the new three-step search, which that list does not include. The counts of tested positions follow from the geometry of each pattern (our reasoning).

• Two-dimensional logarithmic search. Chen et al. (2002) cite Jain and Jain (1981) for it. The FFmpeg version tests the four points of a cross at distance $s$ around the current centre, moves to the best, and halves $s$ only when the centre wins, starting from $s = \lceil p/2 \rceil$.
• Three-step search (TSS). Credited by Chen et al. (2002) to Koga and colleagues at the 1981 National Telecommunications Conference. It tests the eight neighbours at distance $s$ on a $3 \times 3$ grid around the centre, moves to the best of the nine, halves $s$, and repeats. For $p = 7$ the steps are 4, 2 and 1, and it always tests $1 + 3 \times 8 = 25$ positions.
• New three-step search (NTSS). The name matches the title of a 1994 paper by Li, Zeng and Liou. In the FFmpeg version the first step tests the eight neighbours at distance 1 as well as the eight at distance $s$, 17 positions with the centre. The search stops there if the centre wins, finishes with one $3 \times 3$ step if an inner neighbour wins, and otherwise continues as TSS. The extra inner points make small motions cheap to find.
• Four-step search (4SS) of Po and Ma (1996). FFmpeg repeats a $3 \times 3$ pattern with step 2 until its centre wins, then repeats it with step 1 until the centre wins again.
• Diamond search (DS). Chen et al. (2002) cite Zhu and Ma’s 1997 paper at the International Conference on Information, Communications and Signal Processing; the same authors published a journal article under the same title in 2000 (Zhu and Ma, 2000). A large diamond of nine points (the centre, four points at distance 2 along the axes and four diagonal neighbours) is moved to its best point until the centre wins; a small diamond of five points (the centre and its four axis neighbours) then gives the final vector. The first step tests 9 positions; each move adds 5 new ones after a move along an axis and 3 after a diagonal move, because the new diamond overlaps the old; the final step adds 4. A search that stops at once costs 13.
• Hexagon-based search (HEXBS) of Zhu, Lin and Chau (2002). A hexagon of seven points, the centre and six points at distance 2 horizontally or offset $(\pm 1, \pm 2)$, moves until its centre wins, then the four axis neighbours are checked. Each move adds only 3 new positions, and a search that stops at once costs 11. Chen et al. (2002) note that it was shown to find the same vector as diamond search with fewer points.

All of these start at the zero vector and move greedily. They inherit the unimodal assumption: once a pattern’s centre is better than all its neighbours, the search stops, even if a lower minimum lies elsewhere in the window.

Predictive searches

Later encoders added a second idea: start where the answer probably is. Neighbouring blocks usually move alike, so the vectors of already-coded blocks and of the co-located block in the previous frame predict the current one. Chen et al. (2002) describe how searches that start at the zero vector were trapped in local minima on sequences with large motion, such as Bus and Stefan, and how starting from a predicted vector avoids this unless the neighbours belong to another moving object. Their Unsymmetrical-cross Multi-Hexagon-grid Search (UMHexagonS) combines a predicted start, a cross search longer horizontally than vertically (they argue that horizontal motion dominates in natural video), a full search over a $5 \times 5$ area, a set of 16-point hexagons at growing scales to escape local minima, and a final hexagon search. Integrated, together with their fractional-pel search and an early termination test, into version 5.0 of the H.264 reference software and compared with its fast full search (search range 32), it saved, on the common test sequences with one reference frame and CAVLC entropy coding, 91.23% of the integer-pel search time and 46.81% of the total encoding time, at an average bit-rate increase of 0.8% (equivalently, a loss of 0.04 dB), according to the authors’ Table 1.

The H.264 reference software, JM, kept several such options. Its manual lists full search, fast full search (the default, which the JVT common test conditions specify), Uneven Multi-Hexagon Search (UMHex), a simplified hexagon search, and Enhanced Predictive Zonal Search (EPZS) with diamond, square and other refinement patterns (Tourapis et al., 2009, SearchMode). Its default search range is 16, and with rate-distortion optimization enabled the window is centred on the median predicted vector rather than on zero.

Sub-Pixel Refinement

Objects do not move by whole pixels, and a whole-pixel vector leaves a prediction error of up to half a pixel in each direction. Two families of refinement are used.

Interpolated search. The reference frame is interpolated at fractional positions and the search continues on the finer grid. The HM10 encoder first finds the integer vector, then tests the 8 half-pel positions around it, then the 8 quarter-pel positions around the new best, using SAD for the integer search and SA(T)D for the fractional ones (Kim et al., 2013, §4.6.2.2). The JM manual describes the same three layers, with SAD as the default metric at full-pel and SATD at half- and quarter-pel (Tourapis et al., 2009). The interpolation filter is fixed by the standard, because the decoder must form the same prediction (see the section on video coding standards).

Curve fitting. Fit a function to the costs at the integer minimum and its neighbours and take the fitted function’s minimum. A separable parabola through the costs $c_{-1}, c_0, c_{+1}$ at offsets $-1, 0, +1$ along one axis gives the fractional correction

$$\delta = \frac{c_{-1} - c_{+1}}{2 \left( c_{-1} - 2 c_0 + c_{+1} \right)},$$

added to the integer vector. Szeliski (2022, §9.1.3) lists interpolating the matching score around the best value among the ways to get sub-pixel estimates, beside gradient-based refinement in the Lucas–Kanade manner; in stereo (§12.4.1) he notes that such fitting works only where intensities vary smoothly and the window lies on one surface, and that the advisability of fitting curves to integer-sampled costs has been questioned. It needs no interpolated image, only the costs already computed. A parabola suits SSD, which is locally quadratic near a smooth minimum, better than SAD, whose surface near the minimum is closer to a V (our reasoning; the example below compares both).

Parameters

• Block size $B$. Larger blocks average more pixels, so the minimum stands out more clearly from noise, and they need fewer vectors; smaller blocks follow motion boundaries and non-rigid motion more closely but are more often ambiguous. The H.264 recommendation frames the choice as a trade-off between the coding gain of motion compensation with smaller blocks and the amount of data needed to represent the motion (ITU-T H.264, 2003, §0.5.3). H.261 used one vector per $16 \times 16$ macroblock; H.264 lets the encoder split a macroblock down to $4 \times 4$.
• Search range $p$. It must exceed the largest expected motion. Exhaustive cost grows as $(2p+1)^2$, and a larger window also offers more chances for a spurious match in repeated or flat texture. HM10’s encoder description sets the integer search range to 64 (Kim et al., 2013, §4.6.2.2).
• Criterion. SAD by default; SSD where refinement is gradient-based; NCC where brightness changes between frames; SATD with a rate term in encoders.
• Start point. Zero for the classical fast patterns; a predicted vector in JM’s rate-distortion mode, in HM and in predictive searches.
• Search strategy. Exhaustive when accuracy or hardware regularity matters; fast patterns when computation is limited and motion is small; predictive searches for large motion at low cost.

Complexity

For a frame of $N$ pixels with blocks of $B \times B$, exhaustive search evaluates $(2p+1)^2$ candidates per block at $B^2$ absolute differences each, so about $N (2p+1)^2$ operations per frame, independent of the block size. For CIF video (352 × 288 pixels) and $p = 15$, that is about $9.7 \times 10^7$ absolute differences per frame. Memory is small: the block and its search window. Exhaustive search can be computed for all candidates at once as a cost volume, and the regular access pattern suits hardware.

The fast patterns reduce the candidate count from $(2p+1)^2$ to roughly $8 \log_2 (p+1) + 1$ for TSS and to a number that grows with the length of the vector for diamond and hexagon searches. Chen et al. (2002) report that motion estimation took 60% (one reference frame) to 80% (five reference frames) of the encoding time of the H.264 reference codec in their experiments, which is why the search strategy mattered.

Implementation

The example compares exhaustive search (ES), three-step search (TSS), diamond search (DS) and hexagon-based search (HEXBS) on frame pairs with known motion. Each current frame samples a reference image at $\mathbf{x} + \mathbf{v}$, so every block’s true vector is $\mathbf{v}$, and both frames get independent Gaussian noise ($\sigma = 3$ grey levels). Three images are used: camera, a natural photograph with textured areas and a flat sky; gravel, a dense irregular texture; and a synthetic repeated texture, a sum of two sinusoidal gratings with a 6-pixel period plus a weak copy of the gravel texture, so that the best match is unique but has near copies every 6 pixels. Blocks are $16 \times 16$ over the central $256 \times 256$ pixels (256 blocks), the search range is $\pm 7$, and four motions are tested. For each method the code prints the mean endpoint error, the fraction of blocks whose vector is exactly right, and the mean number of distinct positions evaluated per block; a row with $8 \times 8$ blocks shows the effect of block size. The fast patterns follow the FFmpeg implementations described above and skip candidates outside the window. The run used Python 3.13, NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2, took about 3 seconds, and printed the same output twice. The three code blocks run in sequence in one session.

import numpy as np
from numpy.lib.stride_tricks import sliding_window_view
from scipy import ndimage
from skimage import data

P = 7                                    # search range: vectors in [-P, P]^2
rng = np.random.default_rng(0)

def textures():
    g = data.gravel().astype(float)
    y, x = np.mgrid[0:512, 0:512]
    stripes = 128 + 50 * (np.sin(2*np.pi*x/6) + np.sin(2*np.pi*y/6)) + 0.3 * (g - g.mean())
    return {"camera": data.camera().astype(float), "gravel": g, "stripes": stripes}

def frame_pair(img, v, noise=3.0):
    """Current frame samples the reference at x + v, so every block's true vector is v = (dx, dy)."""
    cur = ndimage.shift(img, (-v[1], -v[0]), order=3, mode="reflect")
    return img + rng.normal(0, noise, img.shape), cur + rng.normal(0, noise, img.shape)

def sad(a, b): return np.abs(a - b).sum()

def search(ref, block, by, bx, method):
    """Integer motion search for one block; returns (vector, number of distinct positions tested)."""
    B = block.shape[0]
    cache = {}
    def cost(dx, dy):
        if (dx, dy) not in cache:
            cache[(dx, dy)] = sad(block, ref[by+dy:by+dy+B, bx+dx:bx+dx+B])
        return cache[(dx, dy)]
    def best_of(c, pts):                  # move to the best in-window candidate (centre kept on ties)
        cands = [c] + [(c[0]+a, c[1]+b) for a, b in pts if abs(c[0]+a) <= P and abs(c[1]+b) <= P]
        return min(cands, key=lambda q: cost(*q))
    sq = [(a, b) for a in (-1, 0, 1) for b in (-1, 0, 1) if (a, b) != (0, 0)]
    small_dia = [(-1, 0), (1, 0), (0, -1), (0, 1)]
    if method == "ES":
        win = sliding_window_view(ref[by-P:by+P+B, bx-P:bx+P+B], (B, B))
        s = np.abs(win - block).sum(axis=(2, 3))
        dy, dx = np.unravel_index(np.argmin(s), s.shape)
        return (dx - P, dy - P), s.size
    c = (0, 0)
    if method == "TSS":                   # steps 4, 2, 1 on a 3x3 grid
        for step in (4, 2, 1):
            c = best_of(c, [(a*step, b*step) for a, b in sq])
    elif method == "DS":                  # large diamond until the centre wins, then small diamond
        large_dia = [(-2, 0), (2, 0), (0, -2), (0, 2), (-1, -1), (-1, 1), (1, -1), (1, 1)]
        while (n := best_of(c, large_dia)) != c: c = n
        c = best_of(c, small_dia)
    elif method == "HEXBS":               # hexagon until the centre wins, then small diamond
        hexagon = [(-2, 0), (2, 0), (-1, -2), (1, -2), (-1, 2), (1, 2)]
        while (n := best_of(c, hexagon)) != c: c = n
        c = best_of(c, small_dia)
    return c, len(cache)

def evaluate(ref, cur, v, method, B=16):
    errs, tests = [], []
    for by in range(128, 384, B):
        for bx in range(128, 384, B):
            (dx, dy), n = search(ref, cur[by:by+B, bx:bx+B], by, bx, method)
            errs.append(np.hypot(dx - v[0], dy - v[1])); tests.append(n)
    errs = np.array(errs)
    return errs.mean(), (errs == 0).mean(), np.mean(tests)

motions = [(1, -1), (3, 2), (-5, 4), (6, -6)]
imgs = textures()
print("Integer search, 16x16 blocks, range +-7: mean EPE | fraction exact | mean distinct positions tested")
for name, img in imgs.items():
    pairs = [frame_pair(img, v) for v in motions]
    for method, B in [("ES", 16), ("ES", 8), ("TSS", 16), ("DS", 16), ("HEXBS", 16)]:
        res = [evaluate(r, c, v, method, B) for (r, c), v in zip(pairs, motions)]
        cells = "  ".join(f"{e:4.2f} {f:4.2f} {n:5.1f}" for e, f, n in res)
        print(f"{name:8s} {method:5s} B={B:2d}  {cells}")

Integer search, 16x16 blocks, range +-7: mean EPE | fraction exact | mean distinct positions tested
camera   ES    B=16  0.20 0.94 225.0  0.17 0.95 225.0  0.25 0.95 225.0  0.39 0.92 225.0
camera   ES    B= 8  0.91 0.78 225.0  0.82 0.78 225.0  0.93 0.78 225.0  1.10 0.78 225.0
camera   TSS   B=16  1.21 0.70  25.0  1.66 0.64  25.0  0.84 0.71  25.0  3.33 0.57  25.0
camera   DS    B=16  0.12 0.95  16.0  0.73 0.79  22.7  3.34 0.46  27.0  5.09 0.40  27.8
camera   HEXBS B=16  0.64 0.68  13.8  0.65 0.83  16.6  3.35 0.41  18.7  5.07 0.36  18.4
gravel   ES    B=16  0.00 1.00 225.0  0.00 1.00 225.0  0.00 1.00 225.0  0.00 1.00 225.0
gravel   ES    B= 8  0.00 1.00 225.0  0.00 1.00 225.0  0.00 1.00 225.0  0.00 1.00 225.0
gravel   TSS   B=16  0.00 1.00  25.0  0.03 1.00  25.0  0.00 1.00  25.0  1.65 0.86  25.0
gravel   DS    B=16  0.00 1.00  16.0  0.02 1.00  22.7  1.55 0.82  29.8  4.93 0.55  30.5
gravel   HEXBS B=16  0.19 0.82  13.8  0.00 1.00  17.0  1.23 0.84  20.5  5.30 0.53  19.4
stripes  ES    B=16  0.00 1.00 225.0  0.00 1.00 225.0  0.00 1.00 225.0  0.00 1.00 225.0
stripes  ES    B= 8  0.00 1.00 225.0  0.00 1.00 225.0  0.00 1.00 225.0  0.00 1.00 225.0
stripes  TSS   B=16  1.41 0.79  25.0  7.31 0.00  25.0  1.22 0.80  25.0  8.49 0.00  25.0
stripes  DS    B=16  0.00 1.00  16.0  1.45 0.76  22.0  8.49 0.00  18.9  8.49 0.00  13.0
stripes  HEXBS B=16  0.00 1.00  14.0  2.67 0.55  17.0  8.49 0.00  14.0  8.49 0.00  11.0

Exhaustive search tests 225 positions per block and recovers every vector on the two textured images. On camera it fails on 5–8% of the blocks, most of them presumably in the sky, where the block’s content is flat and the noise decides where the minimum falls (our reasoning; the code does not print which blocks fail). With $8 \times 8$ blocks the failures rise to 22%, since a smaller block contains fewer pixels to outweigh the noise and more blocks fall wholly in flat regions.

The fast searches test 11 to 31 positions, 5–14% of the exhaustive cost, and their accuracy depends on where the error surface is unimodal. For the small motion $(1, -1)$ the diamond search is right on every gravel and repeated-texture block, and the hexagon search on 82% of the gravel blocks. Three-step search has the largest mean error of the fast searches on camera at the two small motions (64–70% of blocks right): its first step jumps 4 pixels from the start, onto a part of the surface that need not slope toward the answer. On camera the diamond search even beats exhaustive search for that motion (mean error 0.12 against 0.20 pixels): starting at zero and stopping as soon as the centre wins keeps its vectors short in the flat sky, where the noisy global minimum is often far away (our reasoning). That bias toward small vectors helps when motion is small and hurts when it is not. For the larger motions $(-5, 4)$ and $(6, -6)$, the diamond and hexagon searches are right on fewer than half of the camera blocks and on 53–84% of the gravel blocks. Gravel grains are a few pixels across, so a few pixels away from the true vector the surface becomes a noisy plateau that gives the search no direction (our reasoning).

The repeated texture breaks the assumption by construction. The criterion has a deep local minimum every 6 pixels, and for $(6, -6)$ one of them lies exactly at the zero vector, where every fast search starts: all of them end at the zero vector on every block, the diamond and hexagon searches after their first pattern, with an error of $6\sqrt{2} \approx 8.49$ pixels. For $(-5, 4)$ the diamond and hexagon searches walk into the replica at $(1, -2)$, also $6\sqrt{2}$ from the truth. Three-step search behaves differently: it is wrong on every block for $(3, 2)$ but right on 80% for $(-5, 4)$. Its first step samples a grid with spacing 4, and whether that grid lands nearer the true minimum or a replica depends on how the grid falls against the 6-pixel period, not on the size of the motion (our reasoning). Exhaustive search is unaffected, because it compares the replicas directly and the weak irregular texture makes the true minimum the lowest.

The second block refines the exhaustive integer vector for two motions that are not whole pixels, $(2.5, -1.25)$ and $(-3.25, 1.75)$. It compares no refinement, the separable parabola fitted to SAD and to SSD costs (the sources leave the fitting function open, so both are printed), and a half-pel then quarter-pel search with 8 neighbours at each level, as HM10 does, on a bilinearly interpolated reference. Bilinear interpolation is a simplification: H.264 and HEVC specify longer filters.

# ---- sub-pixel refinement and matching criteria ----
def bilinear(img, y, x):
    return ndimage.map_coordinates(img, [y, x], order=1)

def refine(ref, block, by, bx, v0, mode):
    B = block.shape[0]
    yy, xx = np.mgrid[by:by+B, bx:bx+B].astype(float)
    if mode == "integer":
        return np.array(v0, float)
    if mode in ("parabola-SAD", "parabola-SSD"):
        f = (lambda d: np.abs(d).sum()) if mode == "parabola-SAD" else (lambda d: (d**2).sum())
        out = []
        for ax in (0, 1):                 # x then y, separately
            c = [f(block - ref[by+v0[1]+(k if ax else 0):by+v0[1]+(k if ax else 0)+B,
                                  bx+v0[0]+(0 if ax else k):bx+v0[0]+(0 if ax else k)+B]) for k in (-1, 0, 1)]
            den = c[0] - 2*c[1] + c[2]
            out.append(v0[ax] + (0.5 * (c[0] - c[2]) / den if den > 0 else 0.0))
        return np.array(out)
    v = np.array(v0, float)               # half-pel then quarter-pel: 8 neighbours each
    for step in (0.5, 0.25):
        cands = [v + step*np.array([a, b]) for a in (-1, 0, 1) for b in (-1, 0, 1)]
        costs = [sad(block, bilinear(ref, yy + c[1], xx + c[0])) for c in cands]
        v = cands[int(np.argmin(costs))]
    return v

print("\nSub-pixel motion, ES integer search then refinement: mean EPE")
for name in ("camera", "gravel"):
    for v in [(2.5, -1.25), (-3.25, 1.75)]:
        ref, cur = frame_pair(imgs[name], v)
        row = []
        for mode in ("integer", "parabola-SAD", "parabola-SSD", "quarter-pel search"):
            e = []
            for by in range(128, 384, 16):
                for bx in range(128, 384, 16):
                    blk = cur[by:by+16, bx:bx+16]
                    v0, _ = search(ref, blk, by, bx, "ES")
                    e.append(np.hypot(*(refine(ref, blk, by, bx, v0, mode) - v)))
            row.append(f"{mode} {np.mean(e):.3f}")
        print(f"{name:7s} v={v}:  " + "  ".join(row))

Sub-pixel motion, ES integer search then refinement: mean EPE
camera  v=(2.5, -1.25):  integer 0.800  parabola-SAD 0.523  parabola-SSD 0.511  quarter-pel search 0.292
camera  v=(-3.25, 1.75):  integer 0.607  parabola-SAD 0.457  parabola-SSD 0.433  quarter-pel search 0.346
gravel  v=(2.5, -1.25):  integer 0.560  parabola-SAD 0.186  parabola-SSD 0.151  quarter-pel search 0.000
gravel  v=(-3.25, 1.75):  integer 0.354  parabola-SAD 0.156  parabola-SSD 0.107  quarter-pel search 0.000

The integer vectors leave at least the rounding error, about 0.56 and 0.35 pixels for the two motions, plus the failed sky blocks on camera. Both parabolas reduce the error, and the SSD parabola does slightly better than the SAD parabola on every case, as expected from the shapes of the two cost surfaces. The quarter-pel search recovers both motions exactly on gravel. That result is favourable by design: both motions lie on the quarter-pel grid, so the search can land exactly on them, while a motion between grid points would leave up to an eighth of a pixel per component (our reasoning). On camera the quarter-pel search leaves 0.29–0.35 pixels, well above its result on gravel and against the hypothesis written before the run, which expected the search to cut the error “to <0.2”: it searches only around the integer vector, so blocks whose integer vector was already wrong, such as flat sky blocks, keep their error under every method (our reasoning; the code does not print errors by block).

The third block repeats the exhaustive search for motion $(3, 2)$ after changing the brightness of the current frame by a gain of 1.3 and an offset of $-20$ grey levels, clipped to $[0, 255]$, and compares SAD, SSD and zero-mean NCC.

def es_criterion(ref, block, by, bx, crit):
    B = block.shape[0]
    win = sliding_window_view(ref[by-P:by+P+B, bx-P:bx+P+B], (B, B))
    if crit == "SAD":
        s = np.abs(win - block).sum(axis=(2, 3))
    elif crit == "SSD":
        s = ((win - block)**2).sum(axis=(2, 3))
    else:                                 # zero-mean NCC, maximised -> minimise its negative
        w = win - win.mean(axis=(2, 3), keepdims=True); b = block - block.mean()
        s = -(w * b).sum(axis=(2, 3)) / (np.sqrt((w**2).sum(axis=(2, 3)) * (b**2).sum()) + 1e-9)
    dy, dx = np.unravel_index(np.argmin(s), s.shape)
    return dx - P, dy - P

print("\nGain 1.3 and offset -20 on the current frame, v=(3, 2): fraction of blocks exact")
for name in ("camera", "gravel"):
    ref, cur = frame_pair(imgs[name], (3, 2))
    cur = np.clip(1.3 * cur - 20, 0, 255)
    row = []
    for crit in ("SAD", "SSD", "NCC"):
        ok = [es_criterion(ref, cur[by:by+16, bx:bx+16], by, bx, crit) == (3, 2)
              for by in range(128, 384, 16) for bx in range(128, 384, 16)]
        row.append(f"{crit} {np.mean(ok):.2f}")
    print(f"{name:7s} " + "  ".join(row))

Gain 1.3 and offset -20 on the current frame, v=(3, 2): fraction of blocks exact
camera  SAD 0.50  SSD 0.69  NCC 0.91
gravel  SAD 0.98  SSD 1.00  NCC 1.00

NCC is unaffected by the gain and offset, apart from clipped pixels, and keeps 91% of the camera blocks right, close to the 95% that exhaustive SAD search found for the same motion without the change in the first block. SAD and SSD are not invariant: on low-contrast camera blocks the brightness change outweighs the differences between candidate positions, and SAD loses half the blocks. On gravel the texture contrast is large enough that all three criteria still find the true vector on at least 98% of the blocks.

Properties and Behavior

• Exactness within its model. Exhaustive search returns the global minimum of the criterion over the window, which is the true vector when the block moves rigidly by an integer amount, stays inside the frame and is distinctive enough to stand out from the noise. Fast searches return a local minimum.
• The vector minimizes a cost, not an error in motion. In flat regions, on edges and in repeated texture, many vectors fit nearly equally well; on a straight edge only the component across it is determined (the aperture problem). An encoder is indifferent to which of the near-equal vectors it gets, since any of them predicts the block well, and with a rate term it prefers the cheapest. A vision system that needs physical motion must detect or regularize these cases.
• Centre bias. Fast searches that start at zero, and the rate term of encoder costs, favour short vectors, which suits video with slow motion and fails with fast motion, as the example shows.
• Brightness constancy. SAD and SSD assume that a point keeps its intensity between frames (brightness constancy) without linearizing it; NCC relaxes this to an affine change per block.
• Regularity. Exhaustive search has a fixed number of operations and a regular memory access pattern, while the fast patterns have data-dependent control flow; the JM manual notes that the JVT common test conditions specified fast full search (Tourapis et al., 2009).

Block Matching in Video Coding Standards

The ITU-T recommendations specify the decoder, that is, how a vector and a residual reconstruct a block. They do not specify how the encoder finds the vector: H.264 says explicitly that the motion estimation process in the encoder is not specified (ITU-T H.264, 2003, §0.5.3), and H.265 leaves the encoding algorithms unspecified as well (ITU-T H.265, 2013, §0.6). Block matching is therefore the encoder’s choice, and the standards constrain only what a vector may be.

• H.261 (1993 edition). Motion compensation is optional in the encoder. The decoder accepts one vector per $16 \times 16$ luminance macroblock, used for all four of its $8 \times 8$ luminance blocks, with integer components of at most $\pm 15$; the chrominance vector is the luminance vector halved and truncated. A positive component means the prediction comes from pixels to the right of or below the predicted ones, the same convention as $I_r(\mathbf{x} + \mathbf{v})$ above. An optional loop filter smooths the prediction (§3.2.2, §3.2.3).
• H.263 (1996 edition). Uses half-pixel precision, “as opposed to” H.261’s full-pixel vectors and loop filter (§3.4); half-pixel values come from bilinear interpolation, and in the default mode vector components lie in $[-16, 15.5]$.
• H.264 (2003 edition). A macroblock can be split into $16 \times 16$, $16 \times 8$, $8 \times 16$ or $8 \times 8$ partitions, and each $8 \times 8$ into $8 \times 4$, $4 \times 8$ or $4 \times 4$, each with its own vector (Figure 6-5). Vectors have quarter-sample accuracy for luma; half-sample values come from a 6-tap filter with taps $(1, -5, 20, 20, -5, 1)$ and quarter-sample values from averaging two neighbouring samples at full- or half-sample positions (§8.4.2.2.1). The encoder can choose the reference picture from several stored decoded pictures, and vectors are coded as differences from a prediction formed from neighbouring vectors (§0.5.3).
• H.265, HEVC (2013 edition). Quarter-sample vectors with longer filters: an 8-tap filter $(-1, 4, -11, 40, 40, -11, 4, -1)$ for half-sample and 7-tap filters such as $(-1, 4, -10, 58, 17, -5, 1)$ for quarter-sample luma positions (§8.5.3.3.3.2). Prediction blocks are carved from coding tree blocks of $16 \times 16$ to $64 \times 64$ luma samples in the Main profile, including asymmetric partitions.

The search range, the criterion and the pattern are thus left to the encoder, and the reference encoders show the choices made in practice. JM offers full search, fast full search, UMHexagonS, a simplified hexagon search and EPZS (Tourapis et al., 2009). HM10 calls its integer search a “3 step motion search strategy”, a different procedure from Koga’s three-step search despite the name: it starts at the predicted vector, runs a diamond (the default) or square pattern search, adds a raster search when the result lies far from the start, refines around the best point with diamond or square patterns until two rounds pass without improvement, and then performs the half- and quarter-pel searches described above (Kim et al., 2013, §4.6.2.2).

Block matching serves outside coding as well. FFmpeg’s minterpolate filter converts frame rates by motion-compensated interpolation, with a choice of exhaustive, three-step, two-dimensional logarithmic, new three-step, four-step, diamond, hexagon-based, EPZS and uneven multi-hexagon searches over $16 \times 16$ blocks by default (FFmpeg 7.1 documentation). The example in the two-stream networks article computes the flow for the temporal stream by exhaustive block matching.

Limitations

• One vector per block. A block that straddles a motion boundary, or contains rotation, zoom or deformation, has no single correct vector; its pixels on one side are predicted badly. Variable block sizes in H.264 and HEVC reduce the problem at motion boundaries.
• Ambiguity. Flat blocks, blocks on straight edges and blocks of repeated texture have several near-equal minima. Exhaustive search picks one of them, which may be wrong in motion terms; fast searches may stop at a wrong one even when the true minimum is distinct.
• Search range. Motion beyond $\pm p$ cannot be found at all, while a larger window costs quadratically more and offers more spurious matches. Hierarchical search over an image pyramid addresses large motion (below).
• Integer grid. Without refinement the vector is off by up to half a pixel in each component; refinement depends on interpolation and on smooth image content.
• Appearance change. Lighting change, shadows, blur and occlusion violate the constancy that SAD and SSD assume. NCC tolerates an affine brightness change but not occlusion or deformation.
• Not true motion. Vectors chosen for compression can be erratic in flat regions and need smoothing or a different cost before they can serve as optical flow.

Variants

• Hierarchical block matching. The search runs first on a downsampled level of an image pyramid, over a correspondingly smaller range, and the result, scaled up, seeds a small search at the next finer level. Szeliski (2022, §9.1.1) describes this scheme and notes that it is not guaranteed to give the full-search result but usually comes close and is much faster. Coarse-to-fine estimation treats it in general.
• Predictive and zonal searches. Starting points from neighbouring and previous vectors, as in EPZS and UMHexagonS, followed by small patterns, described above.
• Gradient-based refinement. After the integer search, Lucas–Kanade iterations on the SSD give continuous sub-pixel estimates (Szeliski, 2022, §9.1.3).
• Frequency-domain matching. Phase correlation, described in the Fourier transform article, finds a translation from the cross-power spectrum instead of a spatial search; it suits global or large translations.
• Stereo. In a rectified stereo pair the correspondence of a pixel lies on the same row, so the search is one-dimensional over disparity. Szeliski (2022, §9.1.1) notes that stereo matching almost always searches all disparities explicitly, since the one-dimensional search has far fewer hypotheses. Window-based matching with SAD or SSD, a winner-take-all choice and parabolic sub-pixel interpolation are the local stereo pipeline that Stereo Matching and Cost Volume describe.
• Dense matching. Computing a block match centred on every pixel, rather than on a block grid, gives a dense but noisy flow field; optical flow estimation adds smoothness and occlusion reasoning to such data terms.