Laplacian Operator

The Laplacian of an image is the sum of its second derivatives in the two coordinate directions. It is zero on flat regions and linear ramps and large at thin lines, spots, and on either side of an edge. Because it treats all directions alike, it was central to early theories of edge detection; smoothed with a Gaussian, it detects blobs across scales; and it appears in sharpening, the Laplacian pyramid, and variational methods for optical flow and image editing.

Definition

The Laplacian is a second-order differential operator. Where the gradient is a vector giving the direction and rate of intensity change, the Laplacian is one number per pixel: the curvature of the intensity surface, summed over two perpendicular directions.

Intuition

Picture the image as a landscape with brightness as height. The Laplacian is negative at the top of a hill, positive at the bottom of a valley, and zero on flat ground or a uniform slope, so a bright spot gives a negative response ringed by a positive one.

Across a blurred edge, the intensity profile is an S-shaped curve. Its slope peaks at the middle of the edge, and its second derivative changes sign there: positive on the dark side, where the profile curves upward, negative on the bright side. The Laplacian therefore marks an edge not by a maximum but by a zero crossing between a positive and a negative band. In discrete form, the Laplacian at a pixel is proportional to the difference between the mean of its neighbors and the pixel itself.

Formal Definition

For an image $I(x, y)$,

$$\nabla^2 I = \frac{\partial^2 I}{\partial x^2} + \frac{\partial^2 I}{\partial y^2}.$$

It is also written $\Delta I$, and it equals the divergence of the gradient, $\nabla \cdot \nabla I$. In the Fourier domain it multiplies each spatial frequency $(\omega_x, \omega_y)$ by $-(\omega_x^2 + \omega_y^2)$: it removes the mean and amplifies fine detail in proportion to the square of its frequency.

Discrete kernels

A Taylor expansion along one axis gives the second difference

$$I(x + 1) - 2 I(x) + I(x - 1) = \frac{\partial^2 I}{\partial x^2} + \frac{1}{12} \frac{\partial^4 I}{\partial x^4} + \dots,$$

the 1D kernel $[1, -2, 1]$. Adding the same kernel along $y$ gives the 4-neighbor (five-point) Laplacian. Using the two diagonals instead gives a kernel with spacing $\sqrt{2}$, which approximates $2 \nabla^2 I$; adding it to the 4-neighbor kernel gives the 8-neighbor kernel, which approximates $3 \nabla^2 I$:

$$L_4 = \begin{bmatrix} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0 \end{bmatrix}, \qquad L_8 = \begin{bmatrix} 1 & 1 & 1 \\ 1 & -8 & 1 \\ 1 & 1 & 1 \end{bmatrix} \approx 3 \nabla^2, \qquad L_9 = \frac{1}{6} \begin{bmatrix} 1 & 4 & 1 \\ 4 & -20 & 4 \\ 1 & 4 & 1 \end{bmatrix}.$$

Extending the Taylor expansion to fourth order, $L_4$ has the error $\tfrac{1}{12}(I_{xxxx} + I_{yyyy})$, which is not rotationally symmetric, and $L_8/3$ has an anisotropic error of the same size and opposite sign. $L_9$, two-thirds axis kernel and one-third diagonal kernel, cancels it, leaving the isotropic error $\tfrac{1}{12}\nabla^4 I$; it is the two-dimensional isotropic stencil discussed by Thampi et al. (2013). All three kernels are symmetric under a 180-degree rotation, so convolution and correlation give the same result and, unlike for the Sobel operator, no flip convention is needed. Libraries differ in sign: many sharpening kernels use the negated form with a positive center.

Laplacian of Gaussian and difference of Gaussians

Because convolution commutes with differentiation, the Laplacian of a Gaussian-smoothed image equals the image convolved with the Laplacian of the Gaussian (LoG):

$$\nabla^2 (G_\sigma * I) = (\nabla^2 G_\sigma) * I, \qquad \nabla^2 G_\sigma(x, y) = \frac{x^2 + y^2 - 2\sigma^2}{2\pi\sigma^6}\, e^{-(x^2 + y^2)/2\sigma^2}.$$

The LoG kernel is a negative center surrounded by a positive ring, often called the Mexican-hat filter. The Gaussian satisfies the diffusion equation $\partial G_\sigma / \partial \sigma = \sigma \nabla^2 G_\sigma$, and replacing the derivative with a difference between two nearby scales gives, as Lowe (2004) derives,

$$G_{k\sigma} - G_{\sigma} \approx (k - 1)\, \sigma^2 \nabla^2 G_\sigma .$$

A difference of Gaussians (DoG) thus approximates the LoG with the factor $\sigma^2$ needed for comparisons across scales already included, and is cheap when a stack of blurred images is computed anyway.

Properties

• Rotational invariance. Rotating the image and then applying $\nabla^2$ gives the same result as applying $\nabla^2$ and then rotating. Every rotation-invariant linear differential operator with constant coefficients is a polynomial in $\nabla^2$, so the Laplacian is the lowest-order one, and a single filter built on it responds to edges of every orientation. The discrete kernels are only approximately isotropic.
• Null space. It gives zero on constants and on linear ramps $a x + b y + c$.
• Second order, hence noisy. The quadratic frequency response amplifies noise much more than a first derivative does. On unit-variance white noise, $L_4$ produces output with standard deviation $\sqrt{20} \approx 4.5$, against about 0.43 for the normalized Sobel derivative.
• Zero crossings and curvature. In coordinates aligned with the gradient ($n$) and the edge ($t$), $\nabla^2 I = I_{nn} + I_{tt}$. On a straight edge $I_{tt} = 0$ and the zero crossing of $\nabla^2 I$ lies exactly where the gradient magnitude peaks. On a curved edge $I_{tt}$ is proportional to the curvature of the edge times the gradient magnitude, so the zero crossings move away from the edge’s center of curvature. For the disk of radius 8 used in the example below, filtered with a LoG of $\sigma = 4$, the zero crossing lies at a radius of about 9.1 pixels.
• Closed contours. Zero crossings are level curves of a smooth function, so they form closed contours or end at the border.

Examples

Edges: Marr–Hildreth

Marr and Hildreth (1980) proposed detecting intensity changes separately at several scales. They argued that the right filter at each scale is a second derivative of a Gaussian which, under simple conditions, need not depend on orientation, so edges are the zero crossings of $\nabla^2 G * I$. Zero-crossing segments from channels of different sizes are then combined into what they called the raw primal sketch, and they related the center-surround $\nabla^2 G$ filter to cells in the visual system.

In practice, LoG zero crossings must be filtered by the slope of the response across them, since weak crossings arise from noise and from gradient minima. Gradient-based detectors, which look for maxima of the first derivative along the gradient direction, such as the Canny edge detector, are the main alternative.

Blobs and scale selection

The raw LoG response falls with scale, because larger Gaussians flatten everything. Lindeberg (1998) studied the scale-normalized Laplacian $t \nabla^2 L$, where $L$ is the image smoothed with a Gaussian of variance $t = \sigma^2$, and detected blobs as points that are local extrema over both space and scale. For a Gaussian blob of variance $t_0$ he showed that the normalized response at the center peaks exactly at $t = t_0$, so the selected scale measures the size of the blob.

A uniform disk of radius $r$ can be worked out directly. Its smoothed center value is $L(0, 0; t) = 1 - e^{-r^2/2t}$. The diffusion equation $\partial_t L = \tfrac{1}{2} \nabla^2 L$ gives $t \nabla^2 L = -(r^2/t)\, e^{-r^2/2t}$ at the center, whose magnitude peaks at $t = r^2/2$, that is at $\sigma = r/\sqrt{2}$, with value $-2/e \approx -0.736$.

Lowe’s SIFT detector (2004) builds on this: it finds keypoints as extrema of difference-of-Gaussian images, which approximate the scale-normalized LoG, in a scale space organized into octaves.

Code: blob scale, edge zero crossings, kernel isotropy

The example measures the scale-normalized LoG at the center of a disk of radius 8 over a range of $\sigma$, counts zero crossings at a blurred step edge with and without noise, and compares how the three discrete kernels respond to a plane wave in different directions.

import numpy as np
from scipy import ndimage
from scipy.special import erf

# 1. Scale selection: a bright disk of radius r on a dark background.
r = 8.0
y, x = np.mgrid[:129, :129] - 64.0
# 4x4 supersampling gives the disk soft, unaliased boundary pixels
sub = (np.arange(4) + 0.5) / 4 - 0.5
disk = np.mean([np.hypot(x + dx, y + dy) < r for dx in sub for dy in sub], axis=0)

print(f"disk radius r = {r}, predicted best scale r/sqrt(2) = {r / np.sqrt(2):.2f}")
print(" sigma   LoG at center   sigma^2 * LoG")
sigmas = np.arange(2.0, 12.01, 0.25)
norm = []
for s in sigmas:
    log = ndimage.gaussian_laplace(disk, s)[64, 64]
    norm.append(s**2 * log)
    if s in (2, 4, 5, 5.5, 5.75, 6, 7, 8, 10, 12):
        print(f"{s:6.2f}   {log:13.4f}   {s**2 * log:13.4f}")
best = sigmas[np.argmin(norm)]          # bright blob -> negative response
print(f"strongest normalized response at sigma = {best:.2f}")

# 2. Zero crossings at a step edge, with and without noise.
rng = np.random.default_rng(1)
n = 200
cols = np.arange(n) - 99.5               # the edge lies between columns 99 and 100
step = np.tile(0.5 * (1 + erf(cols / (np.sqrt(2) * 1.0))), (64, 1))  # blur 1 px
noisy = step + 0.02 * rng.standard_normal(step.shape)

def zero_crossings(row):
    """Positions (in columns) where the row changes sign, by linear interpolation."""
    i = np.where(np.sign(row[:-1]) * np.sign(row[1:]) < 0)[0]
    return i + row[i] / (row[i] - row[i + 1])

four = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]], dtype=float)
print("\nzero crossings along row 32 (true edge at column 99.5)")
print("  filter            image  crossings  strong  strongest at")
filters = [("4-neighbor", lambda im: ndimage.correlate(im, four))]
filters += [(f"LoG sigma={s:g}", lambda im, s=s: ndimage.gaussian_laplace(im, s))
            for s in (1, 2, 4)]
for name, f in filters:
    for label, img in [("clean", step), ("noisy", noisy)]:
        resp = f(img)[32, 15:-15]          # skip the borders
        zc = zero_crossings(resp)
        i = zc.astype(int)
        jump = np.abs(resp[i + 1] - resp[i])   # contrast across each crossing
        strong = jump > 0.25 * jump.max()
        print(f"  {name:16s}  {label}  {len(zc):9d}  {strong.sum():6d}"
              f"  {zc[jump.argmax()] + 15:9.2f}")

# 3. Isotropy of discrete Laplacian kernels: response to a plane wave
#    cos(w (x cos t + y sin t)), whose true Laplacian is -w^2 times the wave.
kernels = {
    "4-neighbor": four,
    "8-neighbor/3": np.array([[1, 1, 1], [1, -8, 1], [1, 1, 1]]) / 3.0,
    "9-point (1/6)": np.array([[1, 4, 1], [4, -20, 4], [1, 4, 1]]) / 6.0,
}
w = np.pi / 4                            # a wavelength of 8 pixels
print(f"\nresponse / exact for a plane wave with wavelength {2 * np.pi / w:.0f} px")
print("angle   " + "".join(f"{k:>15s}" for k in kernels))
for t_deg in (0, 22.5, 45):
    t = np.deg2rad(t_deg)
    wx, wy = w * np.cos(t), w * np.sin(t)
    row = f"{t_deg:5.1f}   "
    for k in kernels.values():
        # frequency response of a 3x3 kernel at (wx, wy)
        H = sum(k[i + 1, j + 1] * np.cos(wx * j + wy * i)
                for i in (-1, 0, 1) for j in (-1, 0, 1))
        row += f"{H / -w**2:15.4f}"
    print(row)

Output:

disk radius r = 8.0, predicted best scale r/sqrt(2) = 5.66
 sigma   LoG at center   sigma^2 * LoG
  2.00         -0.0015         -0.0060
  4.00         -0.0337         -0.5396
  5.00         -0.0284         -0.7092
  5.50         -0.0242         -0.7323
  5.75         -0.0222         -0.7333
  6.00         -0.0203         -0.7290
  7.00         -0.0139         -0.6788
  8.00         -0.0095         -0.6062
 10.00         -0.0047         -0.4651
 12.00         -0.0025         -0.3564
strongest normalized response at sigma = 5.75

zero crossings along row 32 (true edge at column 99.5)
  filter            image  crossings  strong  strongest at
  4-neighbor        clean          1       1      99.50
  4-neighbor        noisy         97      52     123.56
  LoG sigma=1       clean          1       1      99.50
  LoG sigma=1       noisy         66       1      99.53
  LoG sigma=2       clean          1       1      99.50
  LoG sigma=2       noisy         28       1      99.48
  LoG sigma=4       clean          1       1      99.49
  LoG sigma=4       noisy         11       1      99.47

response / exact for a plane wave with wavelength 8 px
angle        4-neighbor   8-neighbor/3  9-point (1/6)
  0.0            0.9496         0.9496         0.9496
 22.5            0.9621         0.9377         0.9499
 45.0            0.9746         0.9257         0.9501

The normalized response peaks at $\sigma = 5.75$, the grid value closest to the predicted $r/\sqrt{2} = 5.66$, and its value there, $-0.733$, is within 0.5% of the analytic $-2/e$. The raw LoG, by contrast, peaks at $\sigma = 4$, which is $r/2$ (the maximum of $(r^2/t^2)\, e^{-r^2/2t}$ lies at $t = r^2/4$), so comparing raw responses across scales would select a scale that is too small.

At the clean edge every filter has a single zero crossing at the true position. With noise of 2% of the edge contrast, the unsmoothed 4-neighbor Laplacian produces 97 crossings, 52 of them strong, the strongest far from the edge. Gaussian smoothing reduces the count to 66, 28, and 11 for $\sigma = 1, 2, 4$, and the slope test leaves one strong crossing within 0.05 pixels of the edge.

The plane-wave test shows the anisotropy of the discrete kernels at a wavelength of 8 pixels. As the Taylor expansion predicts, the 4-neighbor kernel responds 2.6% more strongly along the diagonal than along an axis, the 8-neighbor kernel 2.5% less, and $L_9$ varies by less than 0.1%.

Sharpening

Subtracting a multiple of the Laplacian boosts whatever differs from its surroundings:

$$I_{\text{sharp}} = I - c\, \nabla^2 I .$$

With $L_4$ and $c = 1$ this is a single convolution with a kernel that has 5 at the center and $-1$ at the four neighbors, the common sharpening kernel that convolution compares with unsharp masking. More generally, $G_\sigma * I \approx I + \tfrac{\sigma^2}{2} \nabla^2 I$ for small $\sigma$, so unsharp masking, $I + \alpha (I - G_\sigma * I)$, is approximately $I - \alpha \tfrac{\sigma^2}{2} \nabla^2 I$, and it amplifies noise too.

The Laplacian pyramid

Each level of a Laplacian image pyramid is the difference between a Gaussian pyramid level and the expanded next-coarser level, a band-pass image much like a DoG. Burt and Adelson noted the resemblance to a Laplacian operator whose scale doubles from level to level, hence the name; the levels are not literally discrete Laplacians.

Variational methods and the Poisson equation

Energies that penalize the squared gradient of an unknown field lead, through the Euler–Lagrange equations, to Laplacians. In the Horn–Schunck method for optical flow, the smoothness term produces $\nabla^2 u$ and $\nabla^2 v$, each approximated by the difference between a local average and the value at the pixel. Where the image gradient vanishes, the flow satisfies Laplace’s equation $\nabla^2 u = 0$ and fills in smoothly from the surroundings.

In Poisson image editing, Pérez, Gangnet and Blake (2003) paste a region by solving a Poisson equation inside it: the Laplacian of the unknown pixels must equal the divergence of a guidance field, such as the source image’s gradient, while the values on the region’s boundary come from the destination image. The pasted content keeps its detail but blends in without a visible seam.

Common Misconceptions

• “The 8-neighbor kernel is the isotropic one.” Its direction-dependent error is as large as that of the 4-neighbor kernel, with the opposite sign; $L_9$ cancels it.
• “Zero crossings of the Laplacian are edges.” They include crossings caused by noise and by gradient minima, and on curved edges they are displaced. They mark edge candidates that need a strength test.
• “LoG and DoG are the same filter.” A DoG with scale ratio $k$ only approximates $(k - 1)\sigma^2 \nabla^2 G_\sigma$, more closely as $k$ approaches 1.
• “Raw LoG responses can be compared across scales.” Without the factor $\sigma^2$ the selected scale is too small, as the example shows.

Where It Is Used

• Edge detection by Marr–Hildreth zero crossings, alongside gradient-based detectors built on operators such as Sobel.
• Blob and keypoint detection with scale-normalized LoG extrema and SIFT’s difference of Gaussians.
• Sharpening and the band-pass levels of the Laplacian image pyramid.
• Variational and PDE methods: smoothness terms in optical flow such as Horn–Schunck, Poisson image editing, and diffusion: smoothing with a Gaussian filter of width $\sigma$ solves the heat equation $\partial_t I = \nabla^2 I$ up to time $t = \sigma^2/2$.