Integral Image

An integral image, also called a summed-area table, stores at every pixel the sum of all pixels above and to the left of it. From it, the sum over any axis-aligned rectangle takes four lookups and three additions, whatever the rectangle’s size, and the same construction in $d$ dimensions gives the sum over any box from $2^d$ lookups. Box sums of this kind are the inner loop of many classical methods: the rectangle features of the Viola–Jones detector, the box filters of SURF, the local means and variances of adaptive thresholding and normalized cross-correlation, and the patch distances of fast non-local means.

The name “summed-area table” comes from computer graphics: Lewis (1995) credits the technique to Crow (1984), who introduced it to low-pass filter texture images quickly, and the expanded version of Lewis’s paper notes that it had been used independently in several fields. As Haar-like features describes, Viola and Jones noted the relation to Crow’s tables and chose the name “integral image” for its use in image analysis. That article covers the integral image as Viola and Jones used it: the one-pass recurrences, the four-value rectangle sum, the second table of squared intensities for variance normalization and Lienhart and Maydt’s rotated table for 45° rectangles. This article treats summed-area tables in general: in any number of dimensions, as a filtering tool, and above all as a numerical object whose entries grow with the area they cover, which decides the integer width or floating-point precision a table needs.

Definition

For an image $i(x, y)$ , the integral image is

ii(x, y) = \sum_{x' \le x} \; \sum_{y' \le y} i(x', y'),

and the sum of $i$ over the rectangle with inclusive corners $(x_0, y_0)$ and $(x_1, y_1)$ is

\sum_{x = x_0}^{x_1} \sum_{y = y_0}^{y_1} i(x, y) = ii(x_1, y_1) - ii(x_0 - 1, y_1) - ii(x_1, y_0 - 1) + ii(x_0 - 1, y_0 - 1),

with $ii = 0$ at negative coordinates. The table can instead be stored with one extra row and column of zeros in front, so that entry $(Y, X)$ holds the sum over $y < Y$ and $x < X$ and the rectangle formula needs no special case at the border; OpenCV’s cv::integral returns such a $(W + 1) \times (H + 1)$ table for a $W \times H$ image. The examples below use this padded form, index arrays as [row, column], and build the table with one cumulative sum down the columns followed by one along the rows.

Intuition

A cumulative sum turns a sum over an interval into a difference of two numbers: the running total at the end of the interval minus the running total just before it. In one dimension that is all there is. In two dimensions the table entry at the bottom-right corner of a rectangle counts everything above and to the left of that corner; subtracting the entries to the left of the rectangle and above it removes the unwanted strips, but removes the region above and to the left twice, so the entry at the top-left corner is added back. The work of summing has moved from query time to construction time: one pass over the image builds the table, and afterwards every rectangle costs the same.

The price is that the table’s entries are large. The entry at the far corner holds the sum of the whole image, and a box sum is the small difference of four such large numbers. With exact integer arithmetic this does not matter. With floating-point numbers it decides how accurate the box sums are, and it decides it by where the box lies, not by how big it is.

Formal Definition

For a $d$ -dimensional array $f$ on the grid $\{0, \dots, n_1 - 1\} \times \dots \times \{0, \dots, n_d - 1\}$ , the padded summed-area table is

S(k_1, \dots, k_d) = \sum_{j_1 < k_1} \cdots \sum_{j_d < k_d} f(j_1, \dots, j_d), \qquad 0 \le k_a \le n_a .

Because the sum factorizes over the axes, $S$ is obtained by $d$ successive one-dimensional cumulative sums, one along each axis, in any order; scikit-image’s integral_image computes it this way for arrays of any dimension. The sum over the box $[a_1, b_1) \times \dots \times [a_d, b_d)$ follows by inclusion–exclusion:

\sum_{a_1 \le j_1 < b_1} \cdots \sum_{a_d \le j_d < b_d} f(\mathbf{j}) = \sum_{\boldsymbol{\epsilon} \in \{0, 1\}^d} (-1)^{\,d - |\boldsymbol{\epsilon}|} \, S\big(c_1(\epsilon_1), \dots, c_d(\epsilon_d)\big), \qquad c_a(1) = b_a, \; c_a(0) = a_a ,

where $|\boldsymbol{\epsilon}|$ is the number of far corners chosen. A corner enters with a plus sign when it uses an even number of near coordinates and with a minus sign otherwise, so the box sum takes $2^d$ lookups: 2 in one dimension, 4 for an image and 8 for a video volume or a 3D scan. scikit-image’s integrate enumerates the corners as the binary numbers from 0 to $2^d - 1$ in the same way.

In two dimensions Viola and Jones computed the table with a cumulative sum down each column (which they call the cumulative row sum) and a running sum along each row, as Haar-like features describes, and Lewis (1995) wrote it as the single recurrence

s(u, v) = f(u, v) + s(u - 1, v) + s(u, v - 1) - s(u - 1, v - 1),

which needs three additions and subtractions per pixel instead of two; template matching gives Lewis’s use of it. In the continuous limit, the table is the double integral of the image, which explains the name: a box sum is the integral of the image against a box function, and the second mixed derivative of a box is four signed impulses at its corners (Haar-like features gives this view from the journal version of Viola and Jones).

Properties

Cost

Building the table takes $d$ additions per element with separate cumulative sums. Each box sum then takes $2^d - 1$ additions and subtractions, independent of the box’s size, against $w^d - 1$ additions for summing a $w^d$ box directly. For a box filter, which needs the sum at every position, the comparison is with the two other standard ways of computing it: a separable filter applies a one-dimensional box of $w$ taps along each axis, $d(w - 1)$ additions per output, and running sums along each axis update the previous window by adding the sample that enters and subtracting the one that leaves, about $2d$ operations per output, also independent of $w$ (image filtering and convolution mention this running-sum box filter). For a filter of one fixed size, running sums are as cheap and need no large table. The summed-area table pays off when boxes of many different sizes and positions are needed from the same image, which is the situation in a detector that scales its features rather than the image, in SURF’s box-filter scale space, or in a filter whose window size varies from pixel to pixel.

The table also needs memory: one entry per pixel, with more bits per entry than the image, plus a second table when squared values are needed.

Integer range

For an image of $N$ pixels with values at most $m$ , the largest entry is at most $mN$ , so a table of a $b$ -bit unsigned integer type holds every entry exactly when $mN \le 2^b - 1$ . A table of squared values needs $m^2 N \le 2^b - 1$ . For 8-bit images these limits are small: a 32-bit signed table of a white image overflows at 2,902 × 2,902 pixels, and a 32-bit signed table of its squared values at 182 × 182 pixels (the example prints these sizes for each type). Thresholding reports that Shafait, Keysers and Breuel (2008) point out the overflow of the squared table in 32-bit integers.

Overflow of the table is not by itself fatal. In NumPy, whose cumsum documentation states that arithmetic is modular for integer types and raises no error on overflow, a table that wraps around stores each entry modulo $2^b$ . Inclusion–exclusion uses only additions and subtractions, which are compatible with arithmetic modulo $2^b$ , so the four-corner formula returns the true box sum modulo $2^b$ , and that is the true box sum whenever the box sum itself is below $2^b$ . The width a table needs is therefore set by the largest box that will be queried, not by the image. This is our derivation, which the sources read here do not state; the example confirms it with a 16-bit table of an image whose entries wrapped tens of thousands of times. It depends on the arithmetic wrapping: a type that saturates, or a language in which overflow of signed integers is not defined to wrap, breaks it.

Floating-point exactness and error

A floating-point format with a $p$ -bit significand represents every integer up to $2^p$ exactly: $2^{24} = 16{,}777{,}216$ for 32-bit floats and $2^{53}$ for 64-bit floats. A table of integer pixel values in floating point is therefore exact as long as every partial sum that produced an entry stays below that bound. For a 64-bit float table of 8-bit images the bound is far away (beyond 5 million pixels on a side); for a 32-bit float table of a white image it is passed at 257 × 257 pixels. Beyond the bound each addition rounds to the spacing of floating-point numbers near the running total, which is proportional to the total, and the rounding errors of each entry accumulate along its row and column. A box sum is a difference of four entries; the large common part cancels, but the rounding errors of the four entries do not, so the error of a box sum grows with the magnitude of its corner entries, that is, with the area above and to the left of the box, and not with the size of the box. Relative to the box sum, the error is largest for small boxes far from the origin. NumPy’s documentation notes that the last element of cumsum may differ from sum for floating-point values, because sum may use pairwise summation, which reduces round-off; the mixed-precision training article compares sequential and pairwise accumulation of an image’s pixels in 16- and 32-bit floats.

Lewis gives the reason in signal-processing terms in the expanded version of his paper: the running-sum recursion has its pole at $z = 1$ , on the unit circle, so it is only marginally stable (in our reading, an error made in one entry is carried unchanged into every later entry); he recommends large-integer rather than floating-point arithmetic for the tables. scikit-image’s integral_image promotes floating-point input to at least 64-bit floats “for better accuracy and to avoid potential overflow”, and for integer input follows NumPy’s cumsum, which accumulates integer types smaller than the platform integer in the platform integer (64 bits on the machine of the example, as its printed dtype shows). OpenCV’s cv::integral offers 32-bit integer, 32-bit float and 64-bit float sums; its documentation describes the table of squared values as a 64-bit float array, although its sqdepth argument also accepts 32-bit floats.

Variance from two tables

The mean and variance of a window come from two tables, of $f$ and of $f^2$ :

\mu = \frac{S_1}{n}, \qquad \sigma^2 = \frac{S_2}{n} - \mu^2 ,

where $S_1$ and $S_2$ are the window sums of $f$ and $f^2$ over $n$ pixels. Viola and Jones normalized their detection windows this way, as Haar-like features describes, and Lewis used the same two tables for the denominator of normalized cross-correlation. The formula subtracts two nearly equal numbers when the window is bright and has little texture: for a window of a white page at grey level 235 with a noise standard deviation of 2, $S_2/n$ and $\mu^2$ are both about 55,000 and their difference is about 4. Any error in $S_2$ , of which a 32-bit float table of squares has plenty, appears magnified in $\sigma^2$ , and the computed variance can come out negative. scikit-image’s local thresholds keep both tables in 64-bit floats “for accuracy” and clip $\sigma^2$ at zero before the square root, because, as its code comments, the window mean of $f^2$ is not guaranteed to exceed the squared mean once rounding is considered.

The variance does not change when a constant $c$ is subtracted from every pixel, since $\sigma^2 = \overline{(f - c)^2} - (\overline{f - c})^2$ , so building the tables from $f - c$ , for instance with $c$ the image mean, reduces both the cancellation in the formula and the size of the table entries. In our reasoning it cannot make the table of squares small when the image has strong contrast, because $\sum (f - c)^2$ still grows in proportion to the area; the example measures how much it helps.

Other forms

Several extensions keep the same idea. A table per histogram bin, an integral histogram, gives the histogram of any rectangle in a fixed number of lookups per bin; histograms of oriented gradients describes its use by Zhu, Avidan, Yeh and Cheng (2006). Tables of other per-pixel quantities give window sums of products, as in the box-weighted patch distances of non-local means, or of whole descriptor vectors, as in VLAD’s MultiVLAD. Rotated tables for 45° rectangles are described in Haar-like features. Repeating a box filter approximates a Gaussian, and the Gaussian filter article describes Wells’s (1986) cascades of uniform filters; each pass can come from a summed-area table.

Examples

The example studies the numerics on a synthetic scanned page of 4,096 × 4,096 8-bit pixels: paper at grey level 235 with Gaussian noise of standard deviation 2, and lines of ring-shaped glyphs with 2-pixel strokes at grey level 30, about 9% of the pixels. Its pixels are integers, so 64-bit integer tables and direct sums give exact references. The page was designed once, before the example was run, to resemble the bright, low-variance background on which document thresholding computes local statistics; its parameters were not adjusted afterwards.

The script has five parts:

A. It computes the sums of all 65,536 tiles of 16 × 16 pixels from tables in 64-bit integers, 64-bit floats and 32-bit floats, and groups the errors by the exact value $S$ of the table entry at each tile’s far (bottom-right) corner, in powers of two above $2^{24}$ . It also prints the 32-bit error of single tiles along the diagonal and along the top row.
B. It builds integer tables that overflow: a 16-bit table of the page, read with tiles of 16 × 16 and 32 × 32 pixels, and a 32-bit unsigned table of the squared page, read with 15 × 15 windows.
C. It computes the standard deviation of 65,536 windows of 15 × 15 pixels from tables of $f$ and $f^2$ in 64-bit floats, in 32-bit floats, and in 32-bit floats after subtracting the page’s mean, and compares each with the exact value: the 64-bit result over all windows, the 32-bit results on the background windows (exact $\sigma \le 3$ ). A negative computed variance is counted and then clipped to zero, as scikit-image does.
D. For an image of maximal 8- or 16-bit values, it prints the first square size at which the largest table entry exceeds each type’s range: for the integer types the largest value they hold, for the float types the bound $2^{24}$ or $2^{53}$ up to which every integer is exact (beyond it exactness is no longer guaranteed). It checks five of the 8-bit cases by building the tables in NumPy just below and at that size.
E. It checks a general $d$ -dimensional table with $2^d$ corners against direct sums for the rows of scikit-image’s bundled page image ( $d = 1$ ), the image itself ( $d = 2$ ) and a synthetic 32-frame video of a moving bright square ( $d = 3$ ), compares it with scikit-image’s integral_image and integrate, and prints the operation counts of the cost section.

Before the run we expected exact results in the 64-bit tables and in the part of the 32-bit float table whose entries stay below $2^{24}$ (we derived this from the bound above), an error growing with $S$ elsewhere (in the share of tiles with a nonzero error, in the root-mean-square error from bin to bin, and along the diagonal), a root-mean-square relative tile error above 1% in the top bin of $S$ , exact box sums from the wrapped integer tables, standard-deviation errors within $10^{-6}$ grey levels from the 64-bit tables, median errors above one grey level and some negative variances from the 32-bit tables beyond $2^{24}$ , a smaller but nonzero standard-deviation error after subtracting the mean, and the operation counts of the cost section. The column of root-mean-square relative errors in part A was added after a pilot run, before the run shown, with the editor’s approval, because that expectation was stated on it and the pilot did not print it; nothing else changed. The script runs in about 4 seconds with NumPy 2.3.5 and scikit-image 0.25.2.

import os
os.environ["OMP_NUM_THREADS"] = "1"                 # one BLAS thread, set before NumPy loads
os.environ["OPENBLAS_NUM_THREADS"] = "1"
os.environ["MKL_NUM_THREADS"] = "1"
from itertools import product
from math import isqrt
import numpy as np
from numpy.lib.stride_tricks import sliding_window_view
from skimage import data
from skimage.transform import integral_image, integrate

H = W = 4096


def synthetic_page(seed=0):
    """4096 x 4096 uint8 'scanned page': paper at 235 with noise sigma 2, lines of
    ring-shaped glyphs (2-px strokes) at grey level 30."""
    rng = np.random.default_rng(seed)
    page = 235.0 + rng.normal(0.0, 2.0, (H, W))
    for y0 in range(128, H - 128, 64):            # one text line every 64 rows
        x = 128
        while x < W - 128 - 18:
            w = int(rng.integers(10, 19))          # glyph width 10..18, height 24
            g = np.zeros((24, w), bool)
            g[:2, :] = g[-2:, :] = True
            g[:, :2] = g[:, -2:] = True
            sl = page[y0:y0 + 24, x:x + w]
            sl[g] = 30.0 + rng.normal(0.0, 2.0, int(g.sum()))
            x += w + int(rng.integers(3, 7))
            if rng.random() < 0.15:                # word gap
                x += 12
    return np.clip(np.rint(page), 0, 255).astype(np.uint8)


def video_volume(seed=1, F=32, R=96):
    """F frames of R x R: dark noisy background (40, sigma 3) and a bright 20x20 square (200)
    moving 2 px right and 1 px down per frame."""
    rng = np.random.default_rng(seed)
    v = 40.0 + rng.normal(0.0, 3.0, (F, R, R))
    for f in range(F):
        y, x = 10 + f, 10 + 2 * f
        v[f, y:y + 20, x:x + 20] = 200.0 + rng.normal(0.0, 3.0, (20, 20))
    return np.clip(np.rint(v), 0, 255).astype(np.uint8)


def table(img, dtype):
    """Summed-area table with a leading zero row/column (and plane): entry [Y, X] is the
    sum over y < Y, x < X. One cumulative sum per axis, accumulated in `dtype`."""
    S = np.zeros(tuple(n + 1 for n in img.shape), dtype)
    acc = img.astype(dtype)
    for ax in range(img.ndim):
        acc = acc.cumsum(axis=ax, dtype=dtype)
    S[(slice(1, None),) * img.ndim] = acc
    return S


def box_sums(S, size, step=None):
    """Sums of all boxes of the given size (top-left corners every `step`), from 2^d
    signed corners of the table, in the table's own arithmetic. Returns (sums, n_corners)."""
    d = S.ndim
    step = step or (1,) * d
    out, corners = 0, 0                             # Python 0 takes the dtype of the first term
    for bits in product((0, 1), repeat=d):          # bit 1: corner at the box's far side
        idx = tuple(slice(b * k, S.shape[a] - (1 - b) * k if (1 - b) * k else None, st)
                    for a, (b, k, st) in enumerate(zip(bits, size, step)))
        plus = (d - sum(bits)) % 2 == 0             # even number of near corners: add
        out = out + S[idx] if plus else out - S[idx]  # add/subtract keep the table's dtype
        corners += 1
    return out, corners


page = synthetic_page()
I = page.astype(np.int64)
print(f"page {page.shape} {page.dtype}  min {page.min()} max {page.max()}  mean {I.mean():.4f}  "
      f"ink {np.mean(page < 128):.4f}  background std {I[:128].std():.4f}")

# ---- A. 16 x 16 tile sums from tables of I in three dtypes ---------------------------------
T = 16; n = H // T
direct = I.reshape(n, T, n, T).sum(axis=(1, 3))
S64i = table(page, np.int64)
far = S64i[T::T, T::T]                              # exact far-corner entry of each tile
err = {}
for dt in (np.int64, np.float64, np.float32):
    s, _ = box_sums(table(page, dt), (T, T), (T, T))
    s = s[:n, :n]                                   # drop the partial tiles at the border
    err[np.dtype(dt).name] = s.astype(np.float64) - direct
print("\nA. 16x16 tile sums: error against direct sums, by the far-corner entry S")
print(f"int64 mismatches {np.count_nonzero(err['int64'])}, float64 mismatches "
      f"{np.count_nonzero(err['float64'])}, of {direct.size} tiles")
edges = [0] + [2 ** k for k in range(24, 33)]
lab = np.digitize(far, edges, right=True)
e32 = err["float32"]
print(f"{'S in':<18}{'tiles':>7}{'f32 nonzero':>13}{'f32 RMS':>10}{'f32 max':>10}{'max rel':>10}{'RMS rel':>10}{'f64 nonzero':>13}")
for b in range(1, len(edges)):
    m = lab == b
    e = e32[m]
    lo = "0" if b == 1 else f"2^{22 + b}"
    print(f"({lo}, 2^{23 + b}]".ljust(18) + f"{m.sum():>7}{np.count_nonzero(e):>13}"
          f"{np.sqrt(np.mean(e ** 2)):>10.1f}{np.abs(e).max():>10.1f}"
          f"{np.max(np.abs(e) / direct[m]):>10.3g}{np.sqrt(np.mean((e / direct[m]) ** 2)):>10.3g}{np.count_nonzero(err['float64'][m]):>13}")
pos = [0, 32, 64, 96, 128, 160, 192, 224, 255]
print("float32 error, diagonal tiles (i = j):", [f"{e32[k, k]:.0f}" for k in pos])
print("float32 error, top-row tiles (i = 0): ", [f"{e32[0, k]:.0f}" for k in pos])

# ---- B. modular integer tables --------------------------------------------------------------
print("\nB. integer tables that overflow (NumPy integer arithmetic wraps modulo 2^bits)")
S16 = table(page, np.uint16)
print(f"uint16 table of I: {np.mean(S64i >= 2 ** 16):.4f} of entries wrapped")
for k in (16, 32):
    m = H // k
    ref = I.reshape(m, k, m, k).sum(axis=(1, 3))
    s, _ = box_sums(S16, (k, k), (k, k))
    print(f"  {k}x{k} tile sums: mismatches {np.count_nonzero(s[:m, :m] != ref)} of {ref.size}; "
          f"tiles with true sum >= 2^16: {np.count_nonzero(ref >= 2 ** 16)}")
WIN = 15
I2 = I * I
w1 = I.reshape(n, T, n, T)[:, :WIN, :, :WIN].sum(axis=(1, 3))      # windows at (16i, 16j)
w2 = I2.reshape(n, T, n, T)[:, :WIN, :, :WIN].sum(axis=(1, 3))
S2_64i = table(I2, np.int64)
s, _ = box_sums(table(I2.astype(np.uint32), np.uint32), (WIN, WIN), (T, T))
print(f"uint32 table of I^2: {np.mean(S2_64i >= 2 ** 32):.4f} of entries wrapped; 15x15 window sums "
      f"mismatches {np.count_nonzero(s[:n, :n] != w2)} of {w2.size}")

# ---- C. window standard deviation from tables of I and I^2 ----------------------------------
print("\nC. 15x15 window sigma = sqrt(S2/n - (S1/n)^2), against the exact value")
N = WIN * WIN
sig = np.sqrt(N * w2 - w1 * w1) / N                 # exact integers inside the root
mu = I.mean()


def window_sigma(x, dt):
    s1, _ = box_sums(table(x, dt), (WIN, WIN), (T, T))
    s2, _ = box_sums(table(x.astype(dt) ** 2, dt), (WIN, WIN), (T, T))
    m1, m2 = s1[:n, :n] / dt(N), s2[:n, :n] / dt(N)
    var = m2 - m1 * m1                              # in the table's own precision
    return np.sqrt(np.clip(var, 0, None)).astype(np.float64), np.count_nonzero(var < 0, axis=None), var < 0


s64, neg64, _ = window_sigma(page, np.float64)
s32, _, neg32 = window_sigma(page, np.float32)
s32c, _, neg32c = window_sigma((I - mu).astype(np.float32), np.float32)
print(f"float64 tables: max |d sigma| over all {sig.size} windows {np.abs(s64 - sig).max():.2e}, "
      f"negative variances {neg64}")
bgw = sig <= 3                                      # background windows
labw = np.digitize(S64i[WIN::T, WIN::T][:n, :n], edges, right=True)
print(f"background windows (exact sigma <= 3): {bgw.sum()}; exact sigma median {np.median(sig[bgw]):.3f}")
print(f"{'S in':<18}{'windows':>8}{'f32 med|d|':>12}{'f32 max|d|':>12}{'f32 neg':>9}"
      f"{'f32-mu med|d|':>15}{'f32-mu max|d|':>15}{'f32-mu neg':>12}")
for b in range(1, len(edges)):
    m = (labw == b) & bgw
    d32, d32c = np.abs(s32 - sig)[m], np.abs(s32c - sig)[m]
    lo = "0" if b == 1 else f"2^{22 + b}"
    print(f"({lo}, 2^{23 + b}]".ljust(18) + f"{m.sum():>8}{np.median(d32):>12.3g}{d32.max():>12.3g}"
          f"{np.count_nonzero(neg32[m]):>9}{np.median(d32c):>15.3g}{d32c.max():>15.3g}"
          f"{np.count_nonzero(neg32c[m]):>12}")

# ---- D. sizes at which tables overflow or stop being exact ----------------------------------
print("\nD. first square side n (all pixels at the maximum) that breaks each table")
limits = [("int32", 2 ** 31 - 1), ("uint32", 2 ** 32 - 1), ("float32 exact", 2 ** 24),
          ("float64 exact", 2 ** 53), ("int64", 2 ** 63 - 1)]
print(f"{'':<12}" + "".join(f"{t:>15}" for t, _ in limits))
for bits, top in ((8, 255), (16, 65535)):
    for what, v in (("I", top), ("I^2", top * top)):
        print(f"{bits}-bit {what:<6}" + "".join(f"{isqrt(lim // v) + 1:>15,}" for _, lim in limits))
print("NumPy check, 8-bit all-255 images: does the table equal the int64 table?")
for what, dt, side in (("I", np.uint32, 4105), ("I", np.int32, 2902), ("I", np.float32, 257),
                       ("I^2", np.int32, 182), ("I^2", np.uint32, 258)):
    res = []
    for k in (side - 1, side):
        x = np.full((k, k), 255, np.int64)
        x = x * x if what == "I^2" else x
        ref = x.cumsum(0).cumsum(1)
        got = x.astype(dt).cumsum(0, dtype=dt).cumsum(1, dtype=dt)
        res.append(bool(np.array_equal(got.astype(np.int64), ref)))
    print(f"  {what:<4}{np.dtype(dt).name:>8}: side {side - 1} {res[0]}, side {side} {res[1]}")

# ---- E. n dimensions, the library, and operation counts -------------------------------------
print("\nE. n-D tables: 5^d box sums from 2^d corners against direct sums")
pg, vol = data.page(), video_volume()
print(f"bundled page {pg.shape} {pg.dtype}; volume {vol.shape} {vol.dtype} min {vol.min()} max {vol.max()}")
for d, x in ((1, pg.astype(np.int64)), (2, pg.astype(np.int64)), (3, vol.astype(np.int64))):
    if d == 1:                                      # each row of the page as a 1-D signal
        worst = max(int(np.abs(box_sums(table(r, np.int64), (5,))[0]
                               - sliding_window_view(r, 5).sum(-1)).max()) for r in x)
        corners = box_sums(table(x[0], np.int64), (5,))[1]
    else:
        s, corners = box_sums(table(x, np.int64), (5,) * d)
        ref = sliding_window_view(x, (5,) * d).sum(axis=tuple(range(-d, 0)))
        worst = int(np.abs(s - ref).max())
    print(f"  d={d}: corners {corners}, max |table - direct| {worst}")
ii = integral_image(pg)
print(f"skimage integral_image(page): dtype {ii.dtype}, equals our table "
      f"{np.array_equal(ii, table(pg, np.int64)[1:, 1:])}")
rng = np.random.default_rng(2)
ok = []
for x in (pg, vol):
    iix = integral_image(x)
    for _ in range(5):
        a = [int(rng.integers(0, s - 1)) for s in x.shape]
        b = [int(rng.integers(lo, s)) for lo, s in zip(a, x.shape)]
        box = tuple(slice(lo, hi + 1) for lo, hi in zip(a, b))
        ok.append(integrate(iix, [tuple(a)], [tuple(b)])[0] == x[box].astype(np.int64).sum())
print(f"skimage integrate equals the direct sum for {sum(ok)} of {len(ok)} random boxes")
print("additions per output pixel of a w^d box filter:")
print(f"{'d':>3}{'w':>5}{'direct':>9}{'separable':>11}{'running':>9}{'table':>7}")
for d in (1, 2, 3):
    for w in (3, 15, 63):
        print(f"{d:>3}{w:>5}{w ** d - 1:>9}{d * (w - 1):>11}{2 * d:>9}{d + 2 ** d - 1:>7}")

Output:

page (4096, 4096) uint8  min 21 max 246  mean 216.1918  ink 0.0917  background std 2.0234

A. 16x16 tile sums: error against direct sums, by the far-corner entry S
int64 mismatches 0, float64 mismatches 0, of 65536 tiles
S in                tiles  f32 nonzero   f32 RMS   f32 max   max rel   RMS rel  f64 nonzero
(0, 2^24]            1575            0       0.0       0.0         0         0            0
(2^24, 2^25]         1430         1268       3.7      14.0  0.000237  6.52e-05            0
(2^25, 2^26]         2535         2394       6.7      23.0  0.000422   0.00012            0
(2^26, 2^27]         4389         4237      12.8      48.0   0.00101  0.000236            0
(2^27, 2^28]         7133         7028      25.6     111.0   0.00221  0.000476            0
(2^28, 2^29]        10873        10790      55.6     227.0   0.00394   0.00103            0
(2^29, 2^30]        14937        14878     108.0     481.0   0.00897   0.00202            0
(2^30, 2^31]        16259        16221     232.4     897.0    0.0169   0.00424            0
(2^31, 2^32]         6405         6386     439.6    1577.0    0.0286   0.00778            0
float32 error, diagonal tiles (i = j): ['0', '0', '-4', '-35', '19', '372', '-77', '373', '233']
float32 error, top-row tiles (i = 0):  ['0', '0', '0', '0', '0', '0', '0', '0', '0']

B. integer tables that overflow (NumPy integer arithmetic wraps modulo 2^bits)
uint16 table of I: 0.9994 of entries wrapped
  16x16 tile sums: mismatches 0 of 65536; tiles with true sum >= 2^16: 0
  32x32 tile sums: mismatches 16384 of 16384; tiles with true sum >= 2^16: 16384
uint32 table of I^2: 0.9700 of entries wrapped; 15x15 window sums mismatches 0 of 65536

C. 15x15 window sigma = sqrt(S2/n - (S1/n)^2), against the exact value
float64 tables: max |d sigma| over all 65536 windows 3.58e-12, negative variances 0
background windows (exact sigma <= 3): 37302; exact sigma median 2.013
S in               windows  f32 med|d|  f32 max|d|  f32 neg  f32-mu med|d|  f32-mu max|d|  f32-mu neg
(0, 2^24]             1502        0.15        2.02        5        0.00198         0.0456           0
(2^24, 2^25]          1163        1.43        3.42      357        0.00804          0.164           0
(2^25, 2^26]          1787        1.96        5.52      747         0.0712          0.315           0
(2^26, 2^27]          2461        2.05        8.36     1099         0.0921          0.802           0
(2^27, 2^28]          3714         2.1        12.3     1815          0.227           2.25          32
(2^28, 2^29]          5685        2.17        17.7     2872          0.539            2.3         670
(2^29, 2^30]          7938        2.15        26.5     3985           1.23           5.32        2294
(2^30, 2^31]          8952        2.19        43.1     4507           1.98           8.34        3606
(2^31, 2^32]          4100        2.19        62.8     2089           2.05           8.68        1925

D. first square side n (all pixels at the maximum) that breaks each table
                      int32         uint32  float32 exact  float64 exact          int64
8-bit I               2,902          4,105            257      5,943,261    190,184,349
8-bit I^2               182            258             17        372,182     11,909,806
16-bit I                 182            257             17        370,731     11,863,374
16-bit I^2                 1              2              1          1,449         46,342
NumPy check, 8-bit all-255 images: does the table equal the int64 table?
  I     uint32: side 4104 True, side 4105 False
  I      int32: side 2901 True, side 2902 False
  I    float32: side 256 True, side 257 False
  I^2    int32: side 181 True, side 182 False
  I^2   uint32: side 257 True, side 258 False

E. n-D tables: 5^d box sums from 2^d corners against direct sums
bundled page (191, 384) uint8; volume (32, 96, 96) uint8 min 27 max 212
  d=1: corners 2, max |table - direct| 0
  d=2: corners 4, max |table - direct| 0
  d=3: corners 8, max |table - direct| 0
skimage integral_image(page): dtype uint64, equals our table True
skimage integrate equals the direct sum for 10 of 10 random boxes
additions per output pixel of a w^d box filter:
  d    w   direct  separable  running  table
  1    3        2          2        2      2
  1   15       14         14        2      2
  1   63       62         62        2      2
  2    3        8          4        4      5
  2   15      224         28        4      5
  2   63     3968        124        4      5
  3    3       26          6        6     10
  3   15     3374         42        6     10
  3   63   250046        186        6     10

The integer and 64-bit tables are exact throughout: every tile sum agrees with the direct sum, so a 64-bit float table of an 8-bit image of this size is not merely accurate but exact, because no entry approaches $2^{53}$ . The expectation of an exact 32-bit float table below $2^{24}$ held: all 1,575 tiles whose far corner stays at or below $2^{24}$ have zero error, and so do all tiles of the top row, which fall in that bin because their far corners cannot exceed $255 \times 16 \times 4{,}096 = 16{,}711{,}680 < 2^{24}$ although they reach across the whole width; on the diagonal the error is zero at the origin and nonzero from tile 64 on, as expected. Error is therefore a matter of how much image lies above and to the left of a box, not of the distance from the origin along a row. Above $2^{24}$ the root-mean-square error roughly doubles with each doubling of $S$ , from 3.7 to 439.6 grey-level units, as the spacing of 32-bit floats near the corner entries doubles. The share of tiles with a nonzero error also grows, from 89% in the first bin above $2^{24}$ to 99.8%, but not at the last bin, where it falls slightly (16,221 of 16,259 tiles against 6,386 of 6,405). The expectation of a growing share of nonzero errors therefore failed at the last bin, while the expectations of a root-mean-square error rising from each bin to the next and of the diagonal profile held. The expectation of a root-mean-square relative error above 1% in the top bin failed as well: the printed value is 0.78%. The worst tile in that bin is off by 2.9%, and the worst tile exceeds 0.1% of its sum in every bin from $(2^{26}, 2^{27}]$ on, so a 32-bit table of a large 8-bit image can give box sums that are wrong in the third or second significant digit, a large error for an operation that is exact in principle, but on average not by more than 1% at this image size.

The wrapped integer tables confirm the modular argument, and the expectation of exact box sums from them held. More than 99.9% of the 16-bit table’s entries overflowed, yet every 16 × 16 tile sum is exact, because no such sum reaches $2^{16}$ ; with 32 × 32 tiles every sum exceeds $2^{16}$ and every one is wrong, exactly the tiles the condition predicts. The 32-bit table of squares wrapped in 97% of its entries and still gives every 15 × 15 window sum exactly.

The variance shows how much worse cancellation makes the precision problem. From 64-bit tables the window standard deviation is within $4 \times 10^{-12}$ grey levels of the exact value. From 32-bit tables it is already off by a median of 0.15 grey levels in the first bin, where the table of $f$ is still exact, because the table of $f^2$ passes $2^{24}$ after a few hundred pixels; in every bin beyond, the median error exceeds 1.4 grey levels for a background whose true standard deviation is about 2, and from 31% to 51% of the background windows’ variances come out negative and are clipped to zero. Subtracting the mean reduces the median error by one to two orders of magnitude where the tables are small, with median errors of 0.008 to 0.09 grey levels in the three bins above $2^{24}$ , and the shifted tables give no negative variance up to $2^{27}$ . It does not rescue the far part of the page: in the last bins the median error is close to that of the unshifted tables (2.05 against 2.19 in the top bin); in our reading, because the dark text makes $(f - \mu)^2$ large, so its table still grows with the area. Every expectation about the variance held: the 64-bit error is far below $10^{-6}$ , the raw 32-bit medians exceed one grey level with negative variances in every bin beyond $2^{24}$ , and the shifted medians are smaller than the raw ones in each of those bins but still above 0.1 grey levels in the top bin. For window statistics on large images the remedies that work are a 64-bit float table or exact integer tables wide enough for the window sums.

Part D gives the sizes behind the integer-range section, and the NumPy checks show each table type exact at the last size below the limit and wrong at the first size above it. For 16-bit images the limits are drastic: a single pixel of maximal value squared already overflows a 32-bit signed integer, and even a 64-bit float table of squares passes $2^{53}$ at 1,449 × 1,449 pixels. Part E confirms that the general $d$ -dimensional table with $2^d$ signed corners reproduces direct box sums exactly in one, two and three dimensions and agrees with scikit-image: its integral_image equals the unpadded table and returned it as 64-bit unsigned integers for the 8-bit image, as NumPy’s cumsum promotes small integer types, and its integrate matches the direct sum for all 10 random boxes. The operation counts show the table’s cost per output fixed at 2, 5 and 10 additions for $d = 1, 2, 3$ , against direct sums that grow as $w^d$ ; the expected counts held, as they must for counts derived from the formulas.

Common Misconceptions

“A 32-bit float table is precise enough for 8-bit images.” It is exact only while its entries stay below $2^{24}$ , which for a white image means about 257 × 257 pixels. Beyond that, box sums from a large image can carry errors in the third or second significant digit, and variances from a table of squares can be meaningless, as the example shows. Use 64-bit floats or integers.
“The error of a box sum depends on the size of the box.” It depends on the magnitude of the four corner entries, that is, on the area above and to the left of the box. A small box far from the origin has the largest relative error; a box along the top row of the example is exact anywhere.
“If the table overflows, the box sums are wrong.” With wrapping integer arithmetic the box sums are right as long as each box sum fits in the type; the table needs enough bits for the largest box queried, not for the whole image. Saturating arithmetic, or signed overflow in a language that does not define it, does not have this property.
“Subtracting the mean solves the precision problem.” It removes most of the cancellation in the variance formula and shrinks the table of the image itself, but the table of squared deviations still grows with area wherever the image has contrast.
“An integral image is always the fastest box filter.” For a single, fixed window size, running sums along each axis cost about as much per pixel and need no large table. The table wins when many different box sizes and positions are needed from one image.

Where It Is Used

Rectangle features and cascades. As Haar-like features and Viola–Jones describe, the Viola–Jones detector evaluates every feature from an integral image and normalizes windows with a second table of squares; pedestrian detection describes channel features summed the same way, and multi-block local binary patterns compare rectangle sums.
Box-filter scale spaces. SURF approximates Gaussian derivatives with box filters read from an integral image, so that its filters can grow instead of the image shrinking; scikit-image’s blob_doh, described in blob detection, does the same, and binary descriptors such as BRIEF and ORB compare box-smoothed values.
Normalized cross-correlation. Template matching computes the correlation by the Fourier transform and the window means and energies from tables of $f$ and $f^2$ , following Lewis (1995).
Local thresholding. Niblack’s and Sauvola’s methods in thresholding need the local mean and standard deviation at every pixel; with tables of $f$ and $f^2$ the cost no longer depends on the window size. scikit-image’s implementation cites Shafait, Keysers and Breuel (2008).
Denoising. Fast non-local means computes box-weighted patch distances for every pixel and offset from a table of squared differences.
Histograms and descriptors. Integral histograms give the HOG histogram of any rectangle in constant time, and VLAD’s MultiVLAD aggregates descriptors over grid tiles from one table.
Reference values. Because a table gives exact means over arbitrary rectangles, RoI pooling uses one to compute the ideal feature of each pooling bin in its example.
Graphics. Texture filtering in computer graphics, where Lewis places the technique’s origin in Crow’s (1984) paper.