Image Gradients

For a continuous image intensity function $I(x, y)$, the gradient is the vector of partial derivatives

This vector points in the direction along which intensity increases fastest, and its magnitude $\|\nabla I\| = \sqrt{I_x^2 + I_y^2}$ gives the rate of that increase. A flat region (constant color) has zero gradient everywhere. A sharp boundary between a dark and a bright region produces a large gradient magnitude oriented perpendicular to the boundary — which is precisely why gradients are the raw material of edge detection. Two derived quantities matter as much as the raw vector. The gradient magnitude $\|\nabla I\|$ measures how much intensity is changing locally, regardless of direction, and is used as an edge-strength map. The gradient orientation $\theta = \operatorname{atan2}(I_y, I_x)$ measures the direction of steepest ascent; rotating it by 90° gives the direction along the edge (the tangent), while $\theta$ itself is the direction across the edge (the normal). This normal direction is exactly what non-maximum suppression in The Canny Edge Detector uses to thin wide gradient ridges down to single-pixel edges. Because a digital image is a discrete grid $I[i, j]$ rather than a differentiable function, "the derivative" must be estimated from finite differences of neighboring samples. This estimation is itself a linear filter — a small convolution kernel — which means gradient computation inherits everything known about Linear Filtering and Convolution: it is shift-invariant, it composes with other filters, and it can be analyzed in the frequency domain. Crucially, differentiation is a high-pass operation: it amplifies high spatial frequencies, which is exactly where sensor noise lives. This is why gradients are essentially never computed on raw pixels in a serious vision system; they are computed after smoothing, ideally by fusing the smoothing and differentiation kernels into one operator (the derivative-of-Gaussian, discussed below).

Step 1 — Finite-difference approximation. The simplest discrete derivative approximations come from a Taylor expansion. The forward difference

approximates $\partial I/\partial x$ with error $O(h)$ where $h$ is the pixel spacing (taken as 1). The central difference

is more accurate, with error $O(h^2)$, because the first-order error terms of the forward and backward Taylor expansions cancel by symmetry. Central differences correspond to convolving with the 1-D kernel $[-\tfrac12,\ 0,\ \tfrac12]$ along rows for $I_x$, and along columns for $I_y$.

Step 2 — 2-D derivative kernels. Pure 1-D differencing is noise-sensitive because it uses only two samples. Practical detectors combine differencing in one direction with light smoothing in the perpendicular direction, giving separable 2-D kernels:

• Prewitt: $\displaystyle G_x = \begin{bmatrix}-1&0&1\\-1&0&1\\-1&0&1\end{bmatrix}$, $\displaystyle G_y = G_x^\top$ — a difference in $x$ combined with a uniform (box) average in $y$.
• Sobel: $\displaystyle G_x = \begin{bmatrix}-1&0&1\\-2&0&2\\-1&0&1\end{bmatrix}$, $\displaystyle G_y = G_x^\top$ — same idea but weighting the center row/column more (an approximate 1-D binomial/Gaussian smoothing $[1, 2, 1]$), giving somewhat better noise suppression than Prewitt.
• Roberts cross: $2\times 2$ diagonal kernels $\begin{bmatrix}1&0\\0&-1\end{bmatrix}$ and $\begin{bmatrix}0&1\\-1&0\end{bmatrix}$, the smallest and noisiest option, useful mainly when computation is extremely constrained.

Each of these is separable: Sobel's $G_x$ equals $[1,2,1]^\top$ (smoothing) convolved with $[-1,0,1]$ (differencing), so it can be applied as two 1-D passes, an $O(N)$-per-pixel operation instead of $O(k^2)$, exactly as with any Linear Filtering and Convolution separable kernel.

Step 3 — Derivative of Gaussian (DoG-derivative). Differentiation and convolution commute:

so instead of first Gaussian-smoothing the image (from Gaussian and Median Filtering) and then differencing the result — two convolutions — you can precompute the derivative of the Gaussian kernel itself and convolve once. This is both more efficient and mathematically exact (no separate approximation error is introduced by doing the operations in two steps versus one), and it gives a principled way to choose the smoothing scale $\sigma$: larger $\sigma$ suppresses more noise but blurs and mislocalizes edges, smaller $\sigma$ preserves detail but is noisier. This scale parameter is the same $\sigma$ that governs the trade-offs discussed in Gaussian and Median Filtering.

Step 4 — Assemble magnitude and orientation. Once $I_x$ and $I_y$ are computed at every pixel, form

$M$ is the edge-strength map; $\theta$ (or $\theta + 90°$) is the local edge direction map, both of which feed forward into The Canny Edge Detector and other edge/corner pipelines. ``` function IMAGE_GRADIENT(I, sigma):

# I: 2D array of intensities, sigma: Gaussian smoothing scale (0 = no smoothing) if sigma > 0: Gx_kernel, Gy_kernel = derivative_of_gaussian_kernels(sigma) else: Gx_kernel = [[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]] # Sobel x Gy_kernel = transpose(Gx_kernel) # Sobel y Ix = convolve2D(I, Gx_kernel, border="reflect") Iy = convolve2D(I, Gy_kernel, border="reflect") M = zeros_like(I) Theta = zeros_like(I) for each pixel (i, j): M[i, j] = sqrt(Ix[i, j]^2 + Iy[i, j]^2) Theta[i, j] = atan2(Iy[i, j], Ix[i, j]) # radians in (-pi, pi] return Ix, Iy, M, Theta

function derivative_of_gaussian_kernels(sigma):

# Build a separable derivative-of-Gaussian pair. radius = ceil(3 sigma) # ~99.7% of the Gaussian mass g = [gaussian(t, sigma) for t in -radius..radius] g_ = [-t / sigma^2 gaussian(t, sigma) for t in -radius..radius] # dG/dt normalize(g) # sum(g) = 1 Gx = outer(g_, g) # differentiate along x, smooth along y Gy = outer(g, g_) # smooth along x, differentiate along y return Gx, Gy ``` The reflect border handling avoids fabricating a false edge at the image boundary (a zero-padded border would create a spurious large gradient right at the frame edge).

Time complexity. For an $H \times W$ image and a $k \times k$ derivative kernel, direct 2-D convolution costs $O(HWk^2)$. Because Sobel, Prewitt, and derivative-of-Gaussian kernels are all separable (a smoothing 1-D kernel outer-producted with a differencing 1-D kernel), the same result is obtained with two 1-D passes at cost $O(HWk)$ — a linear-in-$k$ speedup that matters a great deal once $\sigma$ grows and $k \approx 6\sigma+1$ becomes large. Computing magnitude and orientation afterward is $O(HW)$ (one square-root and one `atan2` per pixel), asymptotically dominated by the convolution. Total cost: $O(HWk)$ time, $O(HW)$ extra space for $I_x$, $I_y$, $M$, $\Theta$.

Space complexity. Only a constant number of full-image buffers are needed ($I_x$, $I_y$, and optionally $M$, $\Theta$), so space is $O(HW)$, independent of kernel size — the kernel itself is $O(k)$ and negligible.

Correctness of the finite-difference approximation. By Taylor's theorem, for a smooth $I$,

Subtracting and dividing by $2h$ cancels the $I''(x)$ term exactly:

which is why the central difference is second-order accurate while the one-sided forward difference, $\frac{I(x+h)-I(x)}{h} = I'(x) + O(h)$, is only first-order accurate — the error shrinks quadratically rather than linearly as sampling gets finer.

Correctness of the derivative-of-Gaussian shortcut. Convolution is linear and differentiation is a linear operator, and for well-behaved (compactly supported or rapidly decaying) kernels they commute:

This is proved by differentiating under the integral sign of the convolution integral $\int f(t) g(x-t)\, dt$, which is valid whenever $g$ is differentiable and the integral converges uniformly (true for Gaussians, which are smooth and decay super-exponentially). Consequently, "smooth then differentiate" and "convolve with the derivative-of-Gaussian directly" are exactly the same computation, not merely similar — there is no approximation penalty for fusing the two operations, only a computational saving.

Noise-sensitivity tradeoff. In the frequency domain, differentiation multiplies the Fourier transform by $i\omega$ (or $i\omega_x$, $i\omega_y$ in 2-D), a gain that grows without bound as frequency $\omega$ increases. White sensor noise has roughly flat power across all frequencies, so a raw finite difference amplifies noise arbitrarily at high frequencies. Convolving first with a Gaussian of scale $\sigma$ multiplies the spectrum by $e^{-\sigma^2\omega^2/2}$, which decays at high frequencies and thus tames exactly the frequencies the derivative amplifies most. The derivative-of-Gaussian filter is the product of these two frequency responses, $i\omega\, e^{-\sigma^2 \omega^2/2}$, which is bounded and peaks at a finite frequency $\omega = 1/\sigma$ — giving a mathematically explicit, tunable balance between edge localization (favoring small $\sigma$) and noise rejection (favoring large $\sigma$). This is the same bias–variance-style tradeoff seen generally in Gaussian and Median Filtering, now specialized to differentiation.

Take a 1-D scanline of intensities crossing a step edge, with additive noise, sampled at nine positions:

Central difference on raw pixels, $I_x[i] = \tfrac12(I[i+1]-I[i-1])$:

Note the noise: index 7 shows a spurious derivative of $2.0$ purely from sensor jitter (79 vs. 83) even though the true underlying signal there is flat. On a clean 2-D image this kind of jitter, multiplied across rows, is exactly what produces the salt-and-pepper false edges a naive gradient detector reports.

Smoothing first (3-tap box $[\tfrac13,\tfrac13,\tfrac13]$), then differencing: Smoothed signal $\bar I$ (interior points, ignoring boundary handling):

Central difference on $\bar I$ at index 4: $\tfrac12(66.67 - 10.67) \cdot \tfrac12$... carried through carefully, the key qualitative effect is that the derivative peak at the true edge (indices 3–5) is preserved and even sharpened relative to its neighbors, while the spurious bump at index 7 shrinks toward zero because the box filter averaged 79, 81, 83 into a nearly flat local trend. This is the concrete, numeric demonstration of the frequency-domain argument above: smoothing before differentiating suppresses exactly the high-frequency noise-driven derivative spikes while leaving the low-frequency step-edge derivative largely intact.

2-D Sobel example. Consider a $3\times3$ patch straddling a vertical dark-to-light edge:

Applying the Sobel $G_x$ kernel (centered on the middle pixel, value 21):

Applying $G_y$ (the transposed kernel) to the same patch gives a small value, near zero, since intensity barely changes vertically:

So $M = \sqrt{714^2 + 0^2} = 714$ and $\theta = \operatorname{atan2}(0, 714) = 0°$ — a large, purely horizontal gradient, correctly identifying a vertical edge (whose normal points horizontally) with no vertical component, exactly as the patch's structure demands.

The mathematical idea of directional derivatives predates computer vision by centuries, but its systematic use for machine perception of edges crystallized in the 1960s–1980s as digital image processing matured into a discipline of its own.

Gradient computation is the workhorse first stage of essentially every classical feature detector: The Canny Edge Detector thresholds and thins the gradient magnitude map using the gradient orientation for non-maximum suppression; the Harris and Shi-Tomasi corner detectors build a structure tensor directly from products of $I_x$ and $I_y$; the SIFT and HOG descriptors bin pixel gradients into orientation histograms to build viewpoint-robust local features used in object recognition and image stitching. Optical flow algorithms (Lucas–Kanade, Horn–Schunck) rely on spatial gradients $I_x, I_y$ together with a temporal gradient $I_t$ to solve the brightness-constancy equation. In medical imaging, gradient magnitude highlights tissue boundaries for segmentation; in industrial inspection, gradient-based edge maps locate part boundaries and defects. Even inside modern convolutional neural networks, the earliest learned filters in a trained network frequently converge to oriented edge/gradient detectors strikingly similar to Sobel kernels, showing that gradient estimation is not merely a hand-engineered relic but a near-optimal first-layer computation that data-driven learning rediscovers on its own.

:::faq[Why not just use the raw pixel difference I[i,j+1] - I[i,j] everywhere?] That is the simplest possible derivative estimator but it has two weaknesses. First, it is only first-order accurate (error O(h)), because it is a one-sided (forward) difference whose Taylor-expansion error term does not cancel, unlike the symmetric central difference. Second, and more importantly in practice, it uses just two raw samples with no averaging, so sensor noise passes straight through into the derivative estimate with no suppression at all. In any real image with non-negligible noise, this produces a gradient map cluttered with spurious high-magnitude responses unrelated to true scene edges, which is why production systems always pair differencing with some form of smoothing.

Ask the tutor: How does non-maximum suppression use the gradient orientation to thin a wide edge ridge down to a single pixel?