Tutorials

Image Processing with scipy.ndimage

Image processing is the manipulation of 2D (or higher-dimensional) arrays of values — from photographs to microscopy slides to satellite data. `scipy.ndimage` treats any NumPy array as an n-dimensional image and provides efficient implementations of the most important operations: blurring, edge detection, morphological transformations, and connected component labeling. No image file loading is needed — you work directly with arrays, which makes it easy to combine with the rest of the scientific Python stack.

### Gaussian Blur and Sharpening

`gaussian_filter` blurs an image by convolving it with a Gaussian kernel — it's the standard way to smooth noise before further processing. The inverse operation, sharpening, adds back the difference between the original and the blurred version (unsharp masking).


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

rng = np.random.default_rng(0)
# Synthetic image: white square on black background with noise
img = np.zeros((80, 80))
img[20:60, 20:60] = 1.0
img += 0.15 * rng.standard_normal((80, 80))

blurred  = gaussian_filter(img, sigma=3)
sharpened = img + 1.5 * (img - blurred)  # unsharp mask

fig, axes = plt.subplots(1, 3, figsize=(13, 4))
for ax, data, title in zip(axes, [img, blurred, sharpened],
                            ["Noisy original", "Gaussian blur (σ=3)", "Sharpened"]):
    ax.imshow(data, cmap="gray", origin="lower", vmin=-0.3, vmax=1.3)
    ax.set_title(title)
plt.tight_layout()
plt.show()
- `gaussian_filter(img, sigma=3)` applies a 2D Gaussian with standard deviation 3 pixels — larger sigma = more blur.
- `img + 1.5 * (img - blurred)` is unsharp masking: the high-frequency detail `(img - blurred)` is amplified and added back, making edges crisper.
- `vmin`/`vmax` are set to a fixed range so all three panels use the same color scale for comparison.

### Sobel Edge Detection

Edges are places where pixel values change sharply. `sobel` computes the derivative of an image along one axis; combining the x and y derivatives gives the gradient magnitude, which is large at edges regardless of direction.


from scipy.ndimage import gaussian_filter, sobel
import numpy as np
import matplotlib.pyplot as plt

img = np.zeros((80, 80))
img[20:60, 20:60] = 1.0

blurred = gaussian_filter(img, sigma=1)
sx = sobel(blurred, axis=1)   # horizontal derivative
sy = sobel(blurred, axis=0)   # vertical derivative
edges = np.hypot(sx, sy)       # gradient magnitude

fig, axes = plt.subplots(1, 3, figsize=(13, 4))
for ax, data, title in zip(axes, [img, blurred, edges],
                            ["Original", "Gaussian blur (σ=1)", "Sobel edges"]):
    ax.imshow(data, cmap="gray", origin="lower")
    ax.set_title(title)
plt.tight_layout()
plt.show()
- Applying `gaussian_filter` before `sobel` suppresses noise before differentiating — edges are cleaner with less spurious response.
- `sobel(img, axis=1)` differentiates along columns (left-right); `axis=0` differentiates along rows (up-down).
- `np.hypot(sx, sy)` computes `sqrt(sx² + sy²)` — the Euclidean gradient magnitude — highlighting edges regardless of orientation.

### Morphological Operations

Morphological operations alter the shape of objects in a binary image. Erosion shrinks objects by removing boundary pixels; dilation expands them by adding boundary pixels. Together they form "opening" (erosion then dilation) and "closing" (dilation then erosion) — standard tools for cleaning up binary masks.


from scipy.ndimage import binary_erosion, binary_dilation, binary_fill_holes
import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(1)
img = np.zeros((60, 60), dtype=bool)
img[15:45, 15:45] = True
img[27:33, 27:33] = False    # hole in the middle
img.ravel()[rng.choice(3600, 30, replace=False)] = True  # noise speckles

cleaned  = binary_erosion(img, iterations=2)
expanded = binary_dilation(img, iterations=2)
filled   = binary_fill_holes(img)

fig, axes = plt.subplots(1, 4, figsize=(14, 3.5))
for ax, data, title in zip(axes, [img, cleaned, expanded, filled],
                            ["Original", "Eroded (×2)", "Dilated (×2)", "Holes filled"]):
    ax.imshow(data, cmap="gray", origin="lower")
    ax.set_title(title)
plt.tight_layout()
plt.show()
- `binary_erosion(img, iterations=2)` shrinks the main square and removes the small noise speckles (which erode away entirely after 2 iterations).
- `binary_dilation(img, iterations=2)` expands every `True` region outward, merging close objects and enlarging the main square.
- `binary_fill_holes` fills the interior hole without changing the outer boundary — useful for cleaning up segmentation masks.

### Labeling Connected Components

`label` assigns a unique integer to each distinct connected region in a binary image — a foundational step in object counting, cell detection, and particle analysis.


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

rng = np.random.default_rng(3)
img = np.zeros((80, 80), dtype=bool)

# Three blobs at different positions
for cx, cy, r in [(20, 20, 10), (55, 30, 8), (35, 60, 12)]:
    y, x = np.ogrid[:80, :80]
    img[(x - cx)**2 + (y - cy)**2 < r**2] = True

labeled, n = label(img)
print(f"Detected {n} objects")

plt.figure(figsize=(5, 5))
plt.imshow(labeled, cmap="tab10", origin="lower", vmin=0, vmax=9)
plt.colorbar(label="Object label")
plt.title(f"Connected Components ({n} objects)")
plt.tight_layout()
plt.show()

# Measure each object
from scipy.ndimage import sum as ndi_sum
areas = [ndi_sum(img, labeled, index=i+1) for i in range(n)]
print("Object areas (pixels):", [int(a) for a in areas])
Detected 3 objects

Object areas (pixels): [305, 193, 437]

- `label(img)` returns `(labeled_array, n_features)` — the labeled array has the same shape as the input, with 0 for background and 1, 2, … for each connected region.
- `np.ogrid` creates coordinate grids for building circular masks without explicit loops — `(x - cx)**2 + (y - cy)**2 < r**2` is a disk of radius r centered at (cx, cy).
- `scipy.ndimage.sum(img, labeled, index=i+1)` counts the pixels belonging to object i+1 — a fast alternative to looping over the labeled array.

`scipy.ndimage` is well-suited for scientific image analysis directly on NumPy arrays. For detecting peaks in 1D signals derived from image profiles, see [finding peaks with SciPy](/tutorials/find-peaks-with-scipy), and for computing image gradients numerically, see [NumPy gradient](/tutorials/numpy-gradient).