scikit-image — Scientific Image Processing

Overview

scikit-image is a Python library for image processing in the SciPy ecosystem. It provides algorithms for reading/writing images, filtering (noise reduction, edge detection), geometric transforms, segmentation (thresholding, watershed, active contours), object measurement (area, intensity, shape descriptors), and feature detection. Images are represented as NumPy arrays, enabling seamless integration with NumPy, SciPy, matplotlib, and pandas. Widely used for fluorescence microscopy, histology, and general bioimage analysis.

When to Use

• Preprocessing fluorescence microscopy images: background subtraction, denoising, illumination correction
• Segmenting cells, nuclei, or organelles using thresholding or watershed
• Measuring object properties: area, perimeter, intensity statistics, shape descriptors
• Applying morphological operations: erosion, dilation, opening, closing, fill holes
• Detecting keypoints or local features in biological images
• Converting between image formats and color spaces
• Use OpenCV instead for real-time video processing or GPU-accelerated operations
• For deep-learning cell segmentation, use CellPose instead (better accuracy for touching cells)
• Use napari instead for interactive multi-dimensional image visualization and annotation
• For whole-slide image tiling, use PathML or histolab instead

Prerequisites

• Python packages: scikit-image, numpy, scipy, matplotlib
• Input requirements: Images as files (TIFF, PNG, JPEG) or NumPy arrays; fluorescence images as 2D/3D grayscale arrays
• Environment: Python 3.9+

pip install scikit-image numpy scipy matplotlib

# For reading proprietary microscopy formats
pip install tifffile aicsimageio

# Verify
python -c "import skimage; print(skimage.__version__)"

Quick Start

from skimage import io, filters, measure
import numpy as np

# Load → denoise → threshold → measure
img = io.imread("cells.tif")
img_smooth = filters.gaussian(img, sigma=1.5)
threshold = filters.threshold_otsu(img_smooth)
binary = img_smooth > threshold

regions = measure.regionprops(measure.label(binary))
print(f"Found {len(regions)} objects")
print(f"Mean area: {np.mean([r.area for r in regions]):.1f} px²")

Core API

Module 1: Image I/O and Data Types

from skimage import io, img_as_float, img_as_uint
import numpy as np

# Read single image
img = io.imread("nuclei.tif")
print(f"Shape: {img.shape}, dtype: {img.dtype}")  # (512, 512), uint16

# Read image collection from directory
from skimage import io as ski_io
images = ski_io.ImageCollection("data/*.tif")
print(f"Loaded {len(images)} images")

# Type conversions (critical for correct arithmetic)
img_f = img_as_float(img)      # uint16 → float64, range [0, 1]
img_u8 = (img_f * 255).astype(np.uint8)  # → 8-bit

# Save image
io.imsave("output.tif", img_u8)

# Multi-channel fluorescence (TIFF with CZYX or ZCYX dims)
import tifffile

stack = tifffile.imread("multichannel.tif")  # shape: (C, Z, Y, X)
dapi = stack[0]   # DAPI channel
gfp = stack[1]    # GFP channel
print(f"DAPI: {dapi.shape}, GFP: {gfp.shape}")

# Maximum intensity projection along Z
mip = dapi.max(axis=0)
io.imsave("dapi_mip.tif", mip)

Module 2: Filters and Preprocessing

from skimage import filters, restoration
import numpy as np

# Gaussian blur (denoising, smoothing)
from skimage.filters import gaussian
smoothed = gaussian(img, sigma=2.0)

# Median filter (salt-and-pepper noise removal)
from skimage.filters import median
from skimage.morphology import disk
denoised = median(img, footprint=disk(3))

# Top-hat transform (background subtraction for uneven illumination)
from skimage.morphology import white_tophat, disk
background_removed = white_tophat(img, footprint=disk(50))
print(f"Background removed: range [{background_removed.min()}, {background_removed.max()}]")

# Edge detection
from skimage.filters import sobel, laplace, prewitt

edges_sobel = sobel(img_as_float(img))
edges_laplace = laplace(img_as_float(img))

# Difference of Gaussians (blob-like structure detection)
from skimage.filters import difference_of_gaussians
blob_enhanced = difference_of_gaussians(img_as_float(img), low_sigma=1, high_sigma=3)

# Contrast enhancement (CLAHE: local histogram equalization)
from skimage.exposure import equalize_adapthist
enhanced = equalize_adapthist(img_as_float(img), clip_limit=0.03)

Module 3: Thresholding and Segmentation

from skimage import filters, morphology, segmentation
from skimage.color import label2rgb
import numpy as np

# Automatic thresholding methods
from skimage.filters import (threshold_otsu, threshold_li,
                              threshold_triangle, threshold_yen)

img_f = img_as_float(img)
print(f"Otsu: {threshold_otsu(img_f):.3f}")
print(f"Li: {threshold_li(img_f):.3f}")

# Apply threshold and clean binary mask
binary = img_f > threshold_otsu(img_f)
binary_clean = morphology.remove_small_objects(binary, min_size=50)
binary_filled = morphology.remove_small_holes(binary_clean, area_threshold=100)

# Watershed segmentation (separate touching objects)
from skimage.segmentation import watershed
from skimage.feature import peak_local_max
from scipy import ndimage as ndi

# Distance transform → local maxima → watershed
distance = ndi.distance_transform_edt(binary_filled)
coords = peak_local_max(distance, min_distance=20, labels=binary_filled)
mask = np.zeros(distance.shape, dtype=bool)
mask[tuple(coords.T)] = True
markers = ndi.label(mask)[0]
labels = watershed(-distance, markers, mask=binary_filled)

print(f"Segmented objects: {labels.max()}")
overlay = label2rgb(labels, image=img_f, bg_label=0)

Module 4: Morphological Operations

from skimage.morphology import (erosion, dilation, opening, closing,
                                 disk, ball, binary_erosion, binary_dilation)

# Erosion and dilation
eroded = erosion(binary, footprint=disk(3))
dilated = dilation(binary, footprint=disk(5))

# Opening: erosion then dilation (removes small objects, smooths edges)
opened = opening(binary, footprint=disk(3))

# Closing: dilation then erosion (fills small holes)
closed = closing(binary, footprint=disk(5))

# Skeletonization
from skimage.morphology import skeletonize
skeleton = skeletonize(binary)
print(f"Skeleton pixels: {skeleton.sum()}")

Module 5: Measurement and Region Properties

from skimage import measure
import pandas as pd

# Label connected components
labeled = measure.label(binary_filled)

# Extract region properties
props = measure.regionprops(labeled, intensity_image=img_as_float(img))

# Convert to DataFrame
data = []
for r in props:
    data.append({
        "label": r.label,
        "area": r.area,
        "perimeter": r.perimeter,
        "eccentricity": r.eccentricity,
        "mean_intensity": r.mean_intensity,
        "max_intensity": r.max_intensity,
        "centroid_y": r.centroid[0],
        "centroid_x": r.centroid[1],
        "bbox": r.bbox,
    })

df = pd.DataFrame(data)
print(f"Objects: {len(df)}")
print(df[["area", "mean_intensity", "eccentricity"]].describe().round(2))

# Filter by property thresholds
cells = df[(df["area"] > 100) & (df["area"] < 5000) & (df["eccentricity"] < 0.9)]
print(f"Valid cells: {len(cells)}")

# Measure co-localization: fraction of channel-1 signal in channel-2 positive mask
from skimage.measure import regionprops_table
import numpy as np

# For multi-channel images
table = regionprops_table(
    labeled, intensity_image=np.stack([dapi, gfp], axis=-1),
    properties=["label", "area", "mean_intensity"]
)

Module 6: Feature Detection and Transforms

from skimage.feature import blob_log, blob_dog, corner_harris, corner_peaks
from skimage import transform

# Laplacian of Gaussian blob detection (nuclei, puncta)
blobs = blob_log(img_as_float(img), min_sig