scikit-image-processing

scikit-image — 科学图像处理

概述

scikit-image 是 SciPy 生态系统中的 Python 图像处理库。它提供用于图像读写、滤波（降噪、边缘检测）、几何变换、分割（阈值处理、分水岭、主动轮廓）、对象测量（面积、强度、形状描述符）和特征检测的算法。图像以 NumPy 数组表示，可与 NumPy、SciPy、matplotlib 和 pandas 无缝集成。广泛用于荧光显微镜、组织学和一般生物图像分析。

适用场景

• 预处理荧光显微镜图像：背景扣除、降噪、照明校正
• 使用阈值处理或分水岭分割细胞、细胞核或细胞器
• 测量对象属性：面积、周长、强度统计数据、形状描述符
• 应用形态学操作：腐蚀、膨胀、开运算、闭运算、填充孔洞
• 检测生物图像中的关键点或局部特征
• 在图像格式和色彩空间之间转换
• 实时视频处理或 GPU 加速操作请改用 OpenCV
• 深度学习细胞分割请改用 CellPose（对相互接触的细胞准确度更高）
• 交互式多维图像可视化和标注请改用 napari
• 全切片图像切块请改用 PathML 或 histolab

前置条件

• Python 包：scikit-image、numpy、scipy、matplotlib
• 输入要求：图像文件（TIFF、PNG、JPEG）或 NumPy 数组；荧光图像应为 2D/3D 灰度数组
• 环境：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__)"

快速入门

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²")

核心 API

模块 1：图像 I/O 和数据类型

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)

模块 2：滤波器和预处理

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)

模块 3：阈值处理和分割

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)

模块 4：形态学操作

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()}")

模块 5：测量和区域属性

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"]
)

模块 6：特征检测和变换

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_sigma=5, max_sigma=20,
                  num_sigma=5, threshold=0.05)
print(f"Blobs detected: {len(blobs)}")
# blobs columns: [y, x, sigma] where radius = sqrt(2) * sigma

# Difference of Gaussians (faster alternative)
blobs_dog = blob_dog(img_as_float(img), min_sigma=5, max_sigma=20, threshold=0.02)

# Geometric transforms
from skimage import transform

# Rescale
img_small = transform.rescale(img_as_float(img), 0.5)

# Rotate
img_rotated = transform.rotate(img_as_float(img), angle=15, resize=True)

# Affine registration (align two images)
from skimage.registration import phase_cross_correlation
shift, error, _ = phase_cross_correlation(ref_img, moving_img)
print(f"Alignment shift: {shift} px, error: {error:.4f}")

核心概念

图像数组和约定

scikit-image 使用 NumPy 数组表示图像。形状约定如下：

图像类型	形状	dtype
2D 灰度图像	(H, W)	uint8、uint16、float64
RGB 彩色图像	(H, W, 3)	uint8
多通道图像	(H, W, C)	任意
Z 轴堆栈	(Z, H, W)	任意

dtype 很重要：大多数算法要求使用范围为 [0, 1] 的 float64。处理前使用 img_as_float(img)；保存时使用 img_as_uint(img) 转换回来。

常见工作流

工作流 1：荧光细胞分割和测量

目标：分割 DAPI 染色的细胞核，并测量每个细胞核的 GFP 荧光强度。

from skimage import io, filters, morphology, measure, img_as_float
from skimage.segmentation import watershed
from skimage.feature import peak_local_max
from scipy import ndimage as ndi
import pandas as pd
import numpy as np
import tifffile

# Load 2-channel image (DAPI=ch0, GFP=ch1)
img = tifffile.imread("cells.tif")
dapi = img_as_float(img[0])
gfp = img_as_float(img[1])

# Segment nuclei from DAPI channel
dapi_smooth = filters.gaussian(dapi, sigma=2)
threshold = filters.threshold_otsu(dapi_smooth)
binary = dapi_smooth > threshold
binary = morphology.remove_small_objects(binary, min_size=200)
binary = morphology.remove_small_holes(binary, area_threshold=500)

# Watershed to separate touching nuclei
distance = ndi.distance_transform_edt(binary)
coords = peak_local_max(distance, min_distance=30, labels=binary)
mask = np.zeros_like(distance, dtype=bool)
mask[tuple(coords.T)] = True
markers = ndi.label(mask)[0]
labels = watershed(-distance, markers, mask=binary)

# Measure GFP per nucleus
props = measure.regionprops(labels, intensity_image=gfp)
df = pd.DataFrame([{
    "nucleus_id": p.label,
    "area_px2": p.area,
    "gfp_mean": p.mean_intensity,
    "gfp_max": p.max_intensity,
} for p in props])

df.to_csv("nucleus_measurements.csv", index=False)
print(f"Nuclei: {len(df)}, mean GFP: {df['gfp_mean'].mean():.3f}")

工作流 2：批量图像处理

目标：对文件夹中的图像应用相同的预处理和测量流程。

from pathlib import Path
from skimage import io, filters, measure, img_as_float, morphology
import pandas as pd

results = []
for img_path in sorted(Path("data/").glob("*.tif")):
    img = img_as_float(io.imread(img_path))
    if img.ndim == 3:
        img = img.mean(axis=-1)  # convert RGB to grayscale

    # Preprocess
    smooth = filters.gaussian(img, sigma=1.5)
    thresh = filters.threshold_otsu(smooth)
    binary = morphology.remove_small_objects(smooth > thresh, min_size=50)

    # Measure
    labeled = measure.label(binary)
    props = measure.regionprops(labeled, intensity_image=img)
    for p in props:
        results.append({
            "image": img_path.stem,
            "object_id": p.label,
            "area": p.area,
            "mean_intensity": p.mean_intensity,
        })

df = pd.DataFrame(results)
df.to_csv("batch_results.csv", index=False)
print(f"Processed {df['image'].nunique()} images, {len(df)} objects total")

关键参数

函数	参数	默认值	范围/选项	效果
gaussian	sigma	1.0	0.5–10+	平滑核大小；数值越大，模糊程度越高
median	footprint	disk(1)	disk(1–10)	中值滤波器的邻域
threshold_otsu	—	—	—	返回自动阈值（无需调参）
remove_small_objects	min_size	64	任意整数	移除小于 N 个像素的对象
watershed	—	—	—	使用 peak_local_max 提供的标记位置
peak_local_max	min_distance	1	5–50	检测到的峰之间的最小间距（px）
blob_log	min_sigma / max_sigma	1/50	取决于情况	预期斑点半径范围
regionprops	intensity_image	None	数组	用于强度测量的图像

最佳实践

1. 处理前始终检查 dtype：对 uint8 执行减法等操作会静默地将结果截断为 0。先转换为浮点数：img = img_as_float(img)。

2. 可视化中间结果：对于每个分割流程，在测量前都要将二值掩膜叠加到原图上进行绘制。静默的分割错误是最常见的失败原因。

3. 使用具有代表性的样本调整阈值：对于双峰直方图，Otsu 效果良好。对于较难处理的图像，可使用 skimage 中的 try_all_threshold(img) 比较 Otsu/Li/Triangle。

4. 使用分水岭分割相接的对象：简单阈值法无法分开相接的细胞核。对于密集排列的细胞，始终使用距离变换标记应用分水岭算法。

5. 使用物理单位进行测量：使用显微镜元数据中的像素大小，将像素测量值转换为微米：area_um2 = area_px * pixel_size_um**2。

6. 使用已知样本进行验证：批量处理前，先用 3–5 张经过人工计数的图像验证流程。抽查对象计数是否符合预期。

常见用法

用法：测量荧光图像中的斑点强度

from skimage import io, filters, measure, img_as_float
from skimage.feature import blob_log
import numpy as np

img = img_as_float(io.imread("spots.tif"))
blobs = blob_log(img, min_sigma=2, max_sigma=8, num_sigma=5, threshold=0.05)

intensities = []
for y, x, sigma in blobs:
    r = int(np.ceil(np.sqrt(2) * sigma))
    region = img[max(0,int(y-r)):int(y+r), max(0,int(x-r)):int(x+r)]
    intensities.append({"y": y, "x": x, "radius": r, "mean_intensity": region.mean()})

import pandas as pd
df = pd.DataFrame(intensities)
print(f"Spots: {len(df)}, mean intensity: {df['mean_intensity'].mean():.4f}")

用法：从多个通道生成二值掩膜

from skimage import filters, morphology, img_as_float
import numpy as np

# Segment objects positive in BOTH channels
ch1 = img_as_float(dapi)
ch2 = img_as_float(gfp)

mask1 = ch1 > filters.threshold_otsu(ch1)
mask2 = ch2 > filters.threshold_otsu(ch2)
combined_mask = mask1 & mask2  # objects in both channels

combined_mask = morphology.binary_closing(combined_mask)
print(f"Co-positive pixels: {combined_mask.sum()}")

用法：将带标签的叠加图保存为出版级图片

from skimage.color import label2rgb
from skimage import io, img_as_ubyte
import matplotlib.pyplot as plt

overlay = label2rgb(labels, image=img_as_float(dapi), bg_label=0, alpha=0.5)

fig, axes = plt.subplots(1, 2, figsize=(12, 6))
axes[0].imshow(dapi, cmap="gray")
axes[0].set_title(f"DAPI (raw)")
axes[1].imshow(overlay)
axes[1].set_title(f"Segmented ({labels.max()} nuclei)")
for ax in axes:
    ax.axis("off")
plt.tight_layout()
plt.savefig("segmentation_overlay.png", dpi=300, bbox_inches="tight")
print("Saved segmentation_overlay.png")

故障排除

问题	原因	解决方案
算术运算中出现 OverflowError	uint8/uint16 运算发生整数溢出	先转换为浮点数：img = img_as_float(img)
Otsu 阈值法只找到一个类别	不存在双峰分布	尝试 threshold_li() 或 try_all_threshold(img) 进行比较
分水岭过度分割	标记点距离过近	增大 peak_local_max 中的 min_distance；平滑距离图
二值图像中的所有对象都合并在一起	阈值过低	使用 plt.hist(img.ravel()) 检查直方图；手动调整
regionprops 缺少强度统计数据	未提供 intensity_image	传入：regionprops(labels, intensity_image=img)
3D 图像被作为 2D 图像处理	未检测到 Z-stack	检查形状：img.shape；按切片处理或使用 3D 函数
二值图像中有微小噪声对象	阈值过于激进	应用 morphology.remove_small_objects(binary, min_size=50)

相关技能

• pathml — 使用底层 scikit-image 处理全切片图像
• matplotlib-scientific-plotting — 可视化分割结果和测量分布
• histolab-wsi-processing — 与 scikit-image 兼容的 H&E WSI 图块预处理

参考资料

• scikit-image 文档 — API 参考和教程
• GitHub：scikit-image/scikit-image — 源代码和示例
• van der Walt 等人（2014）“scikit-image：Python 中的图像处理” — PeerJ 2:e453
• scikit-image 示例图库 — 按类别整理并附有注释的代码示例