Skip to content

Convolution and Filtering with scipy.ndimage

Convolution is one of the most fundamental operations in image processing and scientific computing. It allows us to compute weighted sums over neighborhoods in arrays, enabling edge detection, smoothing, and feature extraction.

What is Convolution?

Convolution combines two functions to produce a third function that expresses how one function modifies the other. In discrete form, for a 1D signal, convolution is defined as:

\[s'(t) = \sum_{j=t-\tau}^{t} s(j) f(t-j)\]

where \(s(t)\) is the input signal, \(f\) is the filter kernel, and \(s'(t)\) is the output. The kernel is typically small and is "flipped" as it slides across the signal.

For 2D images, this generalizes naturally:

\[I'(x,y) = \sum_{i,j} I(x+i, y+j) \cdot K(i, j)\]

where \(K\) is a 2D kernel and \(I\) is the image.

Why Convolution Matters

Convolution is the mathematical foundation for: - Blurring and smoothing - Averaging neighboring values - Edge detection - Finding boundaries where intensity changes - Feature extraction - Highlighting patterns in images - Denoising - Removing noise while preserving structure

1D Convolution Example

Let's start with a simple 1D example to understand how ndi.convolve() works:

import numpy as np
from scipy import ndimage as ndi
import matplotlib.pyplot as plt

# Create a simple 1D signal with a spike
signal = np.array([0, 0, 0, 1, 0, 0, 0])

# Define a simple averaging kernel
kernel = np.array([0.25, 0.5, 0.25])

# Apply convolution
result = ndi.convolve(signal, kernel)
print("Signal:", signal)
print("Kernel:", kernel)
print("Result:", result)

The convolve() function slides the kernel across the signal, computing weighted sums at each position. At boundaries, it uses a default boundary handling mode (we'll discuss these next).

2D Convolution and Boundary Modes

When working with 2D images, boundary handling becomes important—what happens at the edges where the kernel extends beyond the image? scipy.ndimage provides several boundary modes:

import numpy as np
from scipy import ndimage as ndi

# Create a small test image
image = np.array([[1, 2, 3],
                  [4, 5, 6],
                  [7, 8, 9]], dtype=float)

# Define a simple kernel
kernel = np.array([[1, 0, -1],
                   [2, 0, -2],
                   [1, 0, -1]]) / 8.0

# Different boundary modes
modes = ['nearest', 'constant', 'wrap', 'reflect']

for mode in modes:
    result = ndi.convolve(image, kernel, mode=mode)
    print(f"Mode: {mode}")
    print(result)
    print()

Boundary Modes Explained

  • nearest: Repeat the edge pixel values (default)
  • constant: Pad with a constant value (default 0)
  • wrap: Wrap around as if the array is periodic
  • reflect: Mirror values at the boundary

Choosing the Right Mode

  • Use reflect for natural images where boundaries shouldn't repeat
  • Use wrap for periodic data (e.g., spherical coordinates)
  • Use nearest when edge pixels matter
  • Use constant (with cval parameter) for synthetic backgrounds

Difference Filters for Edge Detection

Difference filters compute local variations in intensity, revealing edges and boundaries. The simplest is the first-order difference:

import numpy as np
from scipy import ndimage as ndi
from PIL import Image

# Create a simple synthetic image with a step edge
image = np.zeros((5, 10))
image[:, 5:] = 1  # Right half is brighter

# Horizontal difference filter (detects vertical edges)
h_kernel = np.array([[-1, 1]])
h_edges = ndi.convolve(image, h_kernel)

# Vertical difference filter (detects horizontal edges)
v_kernel = np.array([[-1], [1]])
v_edges = ndi.convolve(image, v_kernel)

print("Horizontal edges:\n", h_edges)
print("Vertical edges:\n", v_edges)

Difference filters highlight sharp transitions in intensity, making them ideal for boundary detection in images.

Gaussian Smoothing Kernels

Gaussian smoothing is one of the most important preprocessing steps in image analysis. It reduces noise while preserving edges (better than uniform averaging). The 2D Gaussian kernel is:

\[G(x, y) = \frac{1}{2\pi\sigma^2} e^{-\frac{x^2+y^2}{2\sigma^2}}\]
import numpy as np
from scipy import ndimage as ndi
import matplotlib.pyplot as plt

# Create an image with noise
np.random.seed(42)
image = np.random.rand(50, 50)
# Add a bright spot in the center
image[20:30, 20:30] += 2

# Apply Gaussian filter with sigma=2
smoothed = ndi.gaussian_filter(image, sigma=2)

# Compare
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].imshow(image, cmap='gray')
axes[0].set_title('Noisy Image')
axes[1].imshow(smoothed, cmap='gray')
axes[1].set_title('After Gaussian Smoothing (σ=2)')
plt.tight_layout()
plt.show()

The gaussian_filter() function is optimized using separable convolution—it applies 1D Gaussian filters horizontally and vertically, which is much faster than 2D convolution.

Why Gaussian Smoothing?

  • Natural weighting: nearby pixels have more influence
  • Frequency response: removes high-frequency noise gradually
  • Separable: can be computed efficiently in 1D passes
  • Parameter σ controls the amount of smoothing

Sobel Edge Detection

The Sobel operator combines horizontal and vertical difference filters to detect edges while reducing noise:

\[S_x = \begin{bmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{bmatrix}, \quad S_y = \begin{bmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ 1 & 2 & 1 \end{bmatrix}\]

The edge magnitude at each point is \(|S| = \sqrt{S_x^2 + S_y^2}\), and the edge direction is \(\theta = \arctan(S_y / S_x)\).

import numpy as np
from scipy import ndimage as ndi
import matplotlib.pyplot as plt

# Create a synthetic image with features
image = np.zeros((100, 100))
image[25:75, 25:75] = 1  # White square
image[40:60, 40:60] = 0.5  # Gray inner square

# Compute Sobel edges
sx = ndi.sobel(image, axis=0)  # Horizontal edges
sy = ndi.sobel(image, axis=1)  # Vertical edges
magnitude = np.sqrt(sx**2 + sy**2)
direction = np.arctan2(sy, sx)

# Visualize
fig, axes = plt.subplots(2, 2, figsize=(10, 10))

axes[0, 0].imshow(image, cmap='gray')
axes[0, 0].set_title('Original Image')

axes[0, 1].imshow(sx, cmap='RdBu_r')
axes[0, 1].set_title('Sobel-X (Horizontal Edges)')

axes[1, 0].imshow(sy, cmap='RdBu_r')
axes[1, 0].set_title('Sobel-Y (Vertical Edges)')

axes[1, 1].imshow(magnitude, cmap='gray')
axes[1, 1].set_title('Edge Magnitude')

for ax in axes.flat:
    ax.axis('off')

plt.tight_layout()
plt.show()

# Find strong edges
threshold = 0.1
strong_edges = magnitude > threshold
print(f"Pixels with strong edges: {strong_edges.sum()}")

Understanding the Sobel Kernel

The Sobel kernel weights the central row/column more heavily (coefficient 2 instead of 1) because the center pixel is closer and more reliable. The operator combines:

  1. Smoothing in one direction (reduces noise)
  2. Differencing in the perpendicular direction (detects edges)

This combination makes Sobel more robust to noise than simple difference filters.

Key Takeaway

Sobel edge detection is a practical tool that: - Combines smoothing and differencing for noise robustness - Provides both magnitude (edge strength) and direction - Forms the basis for more advanced edge detection (Canny, etc.) - Runs efficiently on modern hardware

Advanced: Custom Kernels

You can create custom kernels for specialized applications:

import numpy as np
from scipy import ndimage as ndi

# Laplacian kernel (detects both edges and blobs)
laplacian_kernel = np.array([[0, -1, 0],
                             [-1, 4, -1],
                             [0, -1, 0]])

# Mexican Hat (Ricker wavelet) - useful for blob detection
size = 7
y, x = np.ogrid[-size//2:size//2+1, -size//2:size//2+1]
sigma = 1.5
mexican_hat = -(1 - (x**2 + y**2) / (2*sigma**2)) * np.exp(-(x**2 + y**2) / (2*sigma**2))
mexican_hat /= mexican_hat.sum()  # Normalize

# Test on image
image = np.random.rand(50, 50)
image[20:30, 20:30] += 1  # Add a blob

result_laplacian = ndi.convolve(image, laplacian_kernel)
result_hat = ndi.convolve(image, mexican_hat)

print("Laplacian output range:", result_laplacian.min(), "-", result_laplacian.max())
print("Mexican Hat output range:", result_hat.min(), "-", result_hat.max())

Performance Considerations

For large images, consider these optimizations:

import numpy as np
from scipy import ndimage as ndi
import time

image = np.random.rand(1000, 1000)

# Gaussian filter (highly optimized)
start = time.time()
result1 = ndi.gaussian_filter(image, sigma=3)
t1 = time.time() - start

# Manual 2D convolution with same kernel
gaussian_kernel = np.exp(-(np.arange(-3, 4)**2) / 6.0)
gaussian_kernel_2d = gaussian_kernel[:, None] * gaussian_kernel[None, :]
start = time.time()
result2 = ndi.convolve(image, gaussian_kernel_2d)
t2 = time.time() - start

print(f"Gaussian filter: {t1:.4f}s")
print(f"Manual convolution: {t2:.4f}s")
print(f"Speedup: {t2/t1:.1f}x")

Performance Tips

  • Use gaussian_filter() instead of convolve() with Gaussian kernels
  • For separable kernels, apply 1D convolutions sequentially
  • Use order parameter in gaussian_filter() for derivatives
  • Consider frequency-domain methods (FFT) for very large kernels

Practical Example: Image Enhancement

Here's a complete example combining multiple techniques:

import numpy as np
from scipy import ndimage as ndi
import matplotlib.pyplot as plt

# Create synthetic image
np.random.seed(42)
image = np.zeros((100, 100))
image[20:80, 20:80] = 1
image += 0.2 * np.random.randn(100, 100)  # Add noise

# Step 1: Denoise with Gaussian filter
denoised = ndi.gaussian_filter(image, sigma=1)

# Step 2: Detect edges with Sobel
sx = ndi.sobel(denoised, axis=0)
sy = ndi.sobel(denoised, axis=1)
edges = np.sqrt(sx**2 + sy**2)

# Step 3: Sharpen using difference of Gaussians
blurred1 = ndi.gaussian_filter(image, sigma=1)
blurred2 = ndi.gaussian_filter(image, sigma=3)
sharpened = blurred1 - 0.5 * blurred2

# Visualize
fig, axes = plt.subplots(2, 2, figsize=(10, 10))

axes[0, 0].imshow(image, cmap='gray')
axes[0, 0].set_title('Noisy Original')

axes[0, 1].imshow(denoised, cmap='gray')
axes[0, 1].set_title('Denoised')

axes[1, 0].imshow(edges, cmap='gray')
axes[1, 0].set_title('Edge Detection')

axes[1, 1].imshow(sharpened, cmap='gray')
axes[1, 1].set_title('Sharpened')

for ax in axes.flat:
    ax.axis('off')

plt.tight_layout()
plt.show()

Summary

Convolution is the fundamental building block for image processing and filtering. Key points:

  • Convolution applies a weighted sum over neighborhoods using a kernel
  • Different boundary modes handle edges differently
  • Gaussian smoothing is efficient and noise-preserving
  • Sobel edge detection combines smoothing and differencing
  • scipy.ndimage provides optimized implementations

See also: - Generic Filter Operations - Arbitrary neighborhood computations - Binary Structures and Morphology - Connected components and morphological operations