Performance: Vectorization, OpenCV Optimizations & GPU Arrays

"They profiled me and discovered I spend 94 percent of my life inside two lines of Python. I prefer to think of it as having a rich inner loop."

A Thoroughly Profiled Imaging Pipeline

The previous section (Section 8.1) chose the libraries; this one makes them fast. The need is real: a single 4K frame holds 8.3 million pixels, and a 30 fps pipeline gives you 33 milliseconds to do everything Chapter 2 through Chapter 7 taught. This section climbs the performance ladder one rung at a time, measuring at every step, and ends at the rung where Part III of this book lives permanently: the GPU.

1. Measure First: Honest Timing Beginner

Optimization without measurement is superstition. Before rewriting anything, establish a number with a harness that handles the three classic timing traps: one-time costs (thread-pool spin-up, JIT compilation, lazy library initialization) that inflate the first call, run-to-run noise from the operating system, and, for GPUs, asynchronous execution that makes naive stopwatch timing measure only the kernel launch. Code 8.2.1 is the small harness used for every number in this section.

import time
import numpy as np

def bench(fn, *args, repeats=7, warmup=2):
    """Median-of-repeats timing with warm-up runs."""
    for _ in range(warmup):       # absorb one-time costs: thread pools,
        fn(*args)                 # JIT compilation, lazy initialization
    times = []
    for _ in range(repeats):
        t0 = time.perf_counter()
        fn(*args)
        times.append(time.perf_counter() - t0)
    return float(np.median(times))   # median resists OS scheduling spikes

Two habits complete the discipline. First, time the stages of a pipeline separately, not just the whole; the slow part is almost never where intuition points. Second, record the image size and dtype next to every number, because a result measured on a 512×512 uint8 thumbnail says little about a 4K float32 frame: the larger image may not fit in cache, and memory traffic, not arithmetic, usually dominates.

2. The Vectorization Ladder Intermediate

The single most common performance bug in imaging code is a per-pixel Python loop. Code 8.2.2 implements the same gamma correction (the point operation from Chapter 2, $I' = 255\,(I/255)^{\gamma}$) four ways, from a raw double loop to OpenCV's SIMD lookup table.

import numpy as np, cv2

img = np.random.default_rng(0).integers(0, 256, (1080, 1920), dtype=np.uint8)
GAMMA = 0.5

def gamma_loop(im):                       # tier 0: pure Python, ~2 M iterations
    out = np.empty_like(im)
    for r in range(im.shape[0]):
        for c in range(im.shape[1]):
            out[r, c] = round(255 * (im[r, c] / 255) ** GAMMA)
    return out

def gamma_vec(im):                        # tier 1: one whole-array expression
    return (255 * (im / 255.0) ** GAMMA + 0.5).astype(np.uint8)

TABLE = np.array([round(255 * (v / 255) ** GAMMA) for v in range(256)],
                 dtype=np.uint8)          # uint8 input has only 256 cases

def gamma_lut_np(im):                     # tier 2: precomputed lookup table
    return TABLE[im]

def gamma_lut_cv(im):                     # tier 3: OpenCV's SIMD-dispatched LUT
    return cv2.LUT(im, TABLE)

for f in (gamma_loop, gamma_vec, gamma_lut_np, gamma_lut_cv):
    print(f"{f.__name__:13s}: {bench(f, img) * 1000:8.1f} ms")

gamma_loop   :   2410.3 ms
gamma_vec    :     46.2 ms
gamma_lut_np :      4.8 ms
gamma_lut_cv :      1.2 ms

The ladder generalizes far beyond gamma. Whole-array arithmetic, boolean masking (img[img > t] = 255), broadcasting, and fancy indexing cover most of Chapter 2; np.pad plus slicing reproduce the border modes of Chapter 3; and the array model that makes all of it possible was the subject of Chapter 0. Two refinements buy further factors: prefer float32 to float64 when precision allows (half the memory traffic, and OpenCV's float paths are float32-first), and avoid accidental copies; an expression like a + b + c allocates a full temporary per operator, which np.add(a, b, out=buf) avoids in tight real-time loops. Figure 8.2.1 places all the tiers of this section on one chart.

from scipy import ndimage

f  = img.astype(np.float32)
m  = ndimage.uniform_filter(f, size=7)        # E[x]  over each 7x7 window
m2 = ndimage.uniform_filter(f * f, size=7)    # E[x^2]
sigma = np.sqrt(np.maximum(m2 - m * m, 0.0))  # clamp tiny negatives (round-off)

3. OpenCV's Hidden Gears Intermediate

OpenCV is fast by default, but it pays to know which gears are turning. Builds dispatch at runtime to the widest SIMD instruction set the CPU offers (SSE4, AVX2, AVX-512, NEON), many functions parallelize across a thread pool, and cv2.useOptimized() reports whether the optimized dispatch is active. Algorithmic gears matter even more: cv2.GaussianBlur runs as two separable 1D passes, turning the $O(k^2)$ per-pixel cost derived in Chapter 3 into $O(2k)$, and cv2.filter2D silently switches to the FFT route from Chapter 4 when the kernel is large enough for the $O(N \log N)$ transform to win. Code 8.2.3 inspects and exercises these gears, reusing the harness from Code 8.2.1.

import cv2, numpy as np

img = np.random.default_rng(1).integers(0, 256, (2160, 3840), dtype=np.uint8)
print("optimized:", cv2.useOptimized(), "| threads:", cv2.getNumThreads())

blur = lambda im: cv2.GaussianBlur(im, (31, 31), 5.0)
t_par = bench(blur, img)                  # default: all cores

n = cv2.getNumThreads()
cv2.setNumThreads(1)                      # serial baseline for comparison
t_ser = bench(blur, img)
cv2.setNumThreads(n)                      # restore the pool

print(f"31x31 Gaussian, {n} threads: {t_par*1000:6.1f} ms")
print(f"31x31 Gaussian, 1 thread : {t_ser*1000:6.1f} ms")

if cv2.ocl.haveOpenCL():                  # Transparent API: same calls, GPU/iGPU
    u = cv2.UMat(img)                     # upload once to the OpenCL device
    t_gpu = bench(lambda im: cv2.GaussianBlur(im, (31, 31), 5.0), u)
    print(f"31x31 Gaussian, UMat     : {t_gpu*1000:6.1f} ms")

optimized: True | threads: 8
31x31 Gaussian, 8 threads:    7.9 ms
31x31 Gaussian, 1 thread :   31.6 ms
31x31 Gaussian, UMat     :    3.2 ms

One caution: the pip wheels of opencv-python ship without CUDA, so the cv2.cuda module is empty unless you build OpenCV yourself or use a vendor build. For NVIDIA GPU work from Python, the CuPy route in Section 5 is usually the shorter path.

4. Numba: When the Loop Must Stay Advanced

Some algorithms resist whole-array form because each output depends on a previously computed output: scanline accumulations, region growing, dithering with error diffusion (Chapter 1), adaptive recursions. For these, Numba compiles the Python loop itself to machine code. Code 8.2.4 shows the pattern on a per-row exponential moving average, a data-dependent scan NumPy cannot vectorize directly.

from numba import njit, prange
import numpy as np

@njit(parallel=True, cache=True)
def row_ema(im, alpha):
    """Per-row exponential moving average: out[c] depends on out[c-1],
    so the scan cannot be a single whole-array expression."""
    out = np.empty(im.shape, dtype=np.float32)
    for r in prange(im.shape[0]):         # rows are independent: parallel
        out[r, 0] = im[r, 0]
        for c in range(1, im.shape[1]):   # columns are sequential: a scan
            out[r, c] = alpha * im[r, c] + (1.0 - alpha) * out[r, c - 1]
    return out

smooth = row_ema(img.astype(np.float32), 0.1)   # first call compiles (~0.5 s)

Typical results land within a small factor of hand-written C, which puts Numba on the same rung as OpenCV in Figure 8.2.1, for code OpenCV never implemented. The costs are a first-call compilation pause, a restricted Python subset inside compiled functions, and one more dependency to manage; the benefit is keeping algorithm-shaped code readable instead of contorting it into array algebra.

5. GPU Arrays and the Transfer Toll Advanced

CuPy reimplements the NumPy API on NVIDIA GPUs: cp.asarray ships an array to device memory, the familiar operations run as CUDA kernels, and cupyx.scipy.ndimage mirrors the SciPy filters used throughout Part I. The catch is the toll booth: data must cross the PCIe bus in both directions, and kernels launch asynchronously, so honest timing needs the library's own benchmark utility. Code 8.2.5 shows the workflow.

import cupy as cp
from cupyx.scipy import ndimage as cundi
from cupyx.profiler import benchmark

g = cp.asarray(img)                        # host -> device, one PCIe transfer
print(benchmark(cundi.gaussian_filter, (g, 5.0), n_repeat=50))

result = cp.asnumpy(cundi.gaussian_filter(g, 5.0))   # device -> host

gaussian_filter: CPU:   142.1 us  +/- 9.3   GPU-0:  618.4 us  +/- 12.0

That last observation is the whole strategy. A single operation rarely justifies the round trip; a resident pipeline does. Ship the frame once, run the entire chain (correction, filtering, thresholding, even the non-local means denoiser from Chapter 7, which is exactly the kind of arithmetic-heavy kernel GPUs love) on device arrays, and bring back only the result, which is often just a mask or a measurement. Amdahl's law quantifies the ceiling: accelerating a fraction $p$ of total runtime by a factor $s$ yields overall speedup

$$S = \frac{1}{(1 - p) + p / s},$$

so a pipeline that is only 70 percent GPU-eligible caps below $1/0.3 \approx 3.3\times$ no matter how fast the kernels get. The unaccelerated remainder, often disk I/O and JPEG decode from Chapter 1, then becomes the next profiling target. When the pipeline's destination is a neural network anyway, the natural move is to adopt PyTorch tensors as the resident format from the start; Chapter 18 makes that move, and Chapter 28 pushes it to the edge.