scikit-image - Scientific Image Processing

scikit-image treats images as NumPy arrays. It provides a comprehensive suite of algorithms for filtering, feature detection, and object measurement, making it the standard for research-grade image analysis.

When to Use

• Preprocessing scientific images (noise reduction, contrast enhancement).
• Image segmentation (separating cells, particles, or regions of interest).
• Feature extraction (detecting edges, corners, blobs, or textures).
• Geometric transformations (rescaling, rotating, warping).
• Morphological operations (thinning, skeletonization, hole filling).
• Measuring object properties (area, perimeter, eccentricity).
• Restoring degraded images (deconvolution, inpainting).

Reference Documentation

Official docs: https://scikit-image.org/
User Guide: https://scikit-image.org/docs/stable/user_guide.html
Search patterns: skimage.filters, skimage.segmentation, skimage.feature, skimage.morphology

Core Principles

Images are NumPy Arrays

A grayscale image is a 2D array (M, N). A color image is a 3D array (M, N, 3). A multichannel 3D volume is (P, M, N, C).

Coordinate System

The origin (0, 0) is at the top-left corner. Coordinates are always represented as (row, column).

Data Types and Ranges

scikit-image handles various dtypes with specific ranges:

• uint8: 0 to 255
• uint16: 0 to 65535
• float: -1 to 1 or 0 to 1

Quick Reference

Installation

pip install scikit-image

Standard Imports

import numpy as np
import matplotlib.pyplot as plt
from skimage import io, filters, segmentation, feature, measure, morphology, color, util

Basic Pattern - Load and Filter

from skimage import io, filters, color

# Load image
image = io.imread('data.png')

# Convert to grayscale if needed
gray_image = color.rgb2gray(image)

# Apply a filter (e.g., Gaussian blur)
blurred = filters.gaussian(gray_image, sigma=2.0)

# Display
io.imshow(blurred)
io.show()

Critical Rules

✅ DO

• Use utility functions for conversion - Use util.img_as_float, util.img_as_ubyte to handle rescaling automatically when changing dtypes.
• Grayscale for analysis - Most feature extraction and segmentation algorithms expect 2D grayscale arrays.
• Check image shape - Always verify image.shape and image.dtype before processing.
• Vectorize - Use NumPy operations instead of loops over pixels.
• Apply filters before segmentation - Denoising (Gaussian, Median) significantly improves segmentation results.
• Use label for object counting - measure.label is the standard way to identify connected components.

❌ DON'T

• Manually rescale dtypes - Avoid image / 255.0; use util.img_as_float.
• Ignore the "Coordinate Warning" - Be careful with (x, y) vs (row, col). scikit-image uses (row, col).
• Modify the input image - Most functions return a new array; work with copies if you need to mutate.
• Apply color-sensitive filters to grayscale - Some filters behave differently on 3D vs 2D arrays.

Anti-Patterns (NEVER)

# ❌ BAD: Manual dtype conversion (doesn't handle scaling)
image_float = image.astype(float) # Values stay 0-255 but as floats

# ✅ GOOD: Safe utility conversion
from skimage import util
image_float = util.img_as_float(image) # Values scaled to 0.0-1.0

# ❌ BAD: Looping over pixels
for r in range(rows):
    for c in range(cols):
        if image[r, c] > 128:
            image[r, c] = 255

# ✅ GOOD: Vectorized thresholding
image[image > 0.5] = 1.0

# ❌ BAD: Plotting multiple images without a shared axis
plt.imshow(img1); plt.show()
plt.imshow(img2); plt.show()

# ✅ GOOD: Using subplots for comparison
fig, ax = plt.subplots(1, 2)
ax[0].imshow(img1, cmap='gray')
ax[1].imshow(img2, cmap='gray')

Filtering and Restoration (skimage.filters)

Denoising and Edge Detection

from skimage import filters

# Edge detection
edges_sobel = filters.sobel(image)
edges_canny = feature.canny(image, sigma=3)

# Denoising
denoised = filters.median(image, morphology.disk(3))

# Thresholding (Automatic)
val = filters.threshold_otsu(image)
binary = image > val

Morphology (skimage.morphology)

Shaping and Structural Analysis

from skimage import morphology

# Binary morphology
struct_element = morphology.disk(5)
eroded = morphology.erosion(binary, struct_element)
dilated = morphology.dilation(binary, struct_element)
opened = morphology.opening(binary, struct_element) # Erosion then Dilation

# Skeletonization (Reducing objects to 1-pixel width)
skeleton = morphology.skeletonize(binary)

# Removing small objects
clean_binary = morphology.remove_small_objects(binary, min_size=64)

Segmentation (skimage.segmentation)

Separating Objects

from skimage import segmentation, color

# Watershed Segmentation (Marker-based)
from scipy import ndimage as ndi
distance = ndi.distance_transform_edt(binary)
coords = feature.peak_local_max(distance, footprint=np.ones((3, 3)), labels=binary)
mask = np.zeros(distance.shape, dtype=bool)
mask[tuple(coords.T)] = True
markers, _ = ndi.label(mask)
labels = segmentation.watershed(-distance, markers, mask=binary)

# SLIC (Superpixels)
segments = segmentation.slic(image, n_segments=100, compactness=10)
out = color.label2rgb(segments, image, kind='avg')

Feature Detection (skimage.feature)

Keypoints and Textures

from skimage import feature

# Local Binary Patterns (Texture analysis)
lbp = feature.local_binary_pattern(image, P=8, R=1)

# Corner detection (Harris)
coords = feature.corner_peaks(feature.corner_harris(image), min_distance=5)

# Blob detection (Difference of Gaussian)
blobs = feature.blob_dog(image, max_sigma=30, threshold=.1)
# blobs: array of [row, col, sigma]

Measurements (skimage.measure)

Quantifying Results

from skimage import measure

# Label connected components
labels = measure.label(binary)

# Calculate properties
props = measure.regionprops(labels)

for prop in props:
    print(f"Label: {prop.label}")
    print(f"Area: {prop.area}")
    print(f"Centroid: {prop.centroid}")
    print(f"Eccentricity: {prop.eccentricity}")

# Find contours
contours = measure.find_contours(binary, 0.8)

Practical Workflows

1. Particle Counting Pipeline

def count_particles(image_path):
    # 1. Load and grayscale
    img = color.rgb2gray(io.imread(image_path))
    
    # 2. Denoise
    img_denoised = filters.gaussian(img, sigma=1)
    
    # 3. Threshold
    thresh = filters.threshold_otsu(img_denoised)
    binary = img_denoised < thresh # Assuming dark particles
    
    # 4. Morphological cleaning
    binary = morphology.remove_small_objects(binary, 50)
    binary = morphology.closing(binary, morphology.disk(3))
    
    # 5. Label and measure
    labels = measure.label(binary)
    return measure.regionprops(labels), labels

# properties, label_img = count_particles('samples.tif')

2. Micrograph Analysis (Nuclei Segmentation)

def segment_nuclei(dna_image):
    """Identify nuclei in a DAPI/Hoechst stained image."""
    # Local thresholding for uneven illumination
    local_thresh = filters.threshold_local(dna_image, block_size=35)
    binary = dna_image > local_thresh
    
    # Fill holes
    filled = ndi.binary_fill_holes(binary)
    
    # Separate touching nuclei
    distance = ndi.distance_transform_edt(filled)
    local_maxi = feature.peak_local_max(distance, indices=False, footprint=np.ones((15, 15)), labels=filled)
    markers = measure.label(local_maxi)
    labels = segmentation.watershed(-distance, markers, mask=filled)
    
    return labels

3. Change Detection (Image Subtraction)

def detect_change(img_before, img_after, threshold=0.1):
    # Ensure float
    im1 = util.img_as_float(color.rgb2gray(img_before))
    im2 = util.img_as_float(color.rgb2gray(img_after))
    
    # Difference
    diff = np.abs(im1 - im2)
    
    # Filter noise in difference
    diff_clean = filters.median(diff, morphology.disk(2))
    
    return diff_clean > threshold

Performance Optimization

Using skimage.util.view_as_windows

from skimage.util import view_as_windows

# Create overlapping patches without copying data (using strides)
# Fast for local neighborhood calculations
patches = view_as_windows(image, (64, 64), step=32)
# Shape: (n_patches_r, n_patches_c, 64, 64)

Parallel Processing

# Many skimage functions are already Cython-optimized.
# For batch processing, use standard joblib or multiprocessing.
from joblib import Parallel, delayed

results = Parallel(n_jobs=-1)(delayed(filters.gaussian)(img) for img in image_list)

Common Pitfalls and Solutions

Handling Multichannel Images

# ❌ Problem: Filters failing on RGB images
# result = filters.median(rgb_image) # Raises error

# ✅ Solution: Filter per channel or convert
# filters.median(color.rgb2gray(rgb_image))
# OR apply to each channel
# result = np.stack([filters.median(rgb_image[..., i]) for i in range(3)], axis=-1)

Memory issues with label

# ❌ Problem: Large noisy images produce millions of 1-pixel labels
# ✅ Solution: Clear small noise before labeling
binary = morphology.remove_small_objects(binary, min_size=10)
labels = measure.label(binary)

Coordinates: (X, Y) vs (Row, Col) confusion

# ❌ Problem: Points appear swapped on plot
# Matplotlib: plt.plot(x, y)
# scikit-image: image[row, col]

# ✅ Solution: Explicitly swap for visualization
coords = feature.peak_local_max(img) # Returns [row, col]
plt.plot(coords[:, 1], coords[:, 0], 'r.') # Plot as [x, y]

scikit-image provides the mathematical rigor needed for scientific discovery. By building on the NumPy ecosystem, it allows for a seamless workflow from raw sensor data to quantifiable insights.