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：画像の入出力とデータ型

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 しきい値処理でクラスが 1 つしか見つからない	二峰性の分布がない	比較のため、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 スタックが検出されていない	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 et al.（2014）「scikit-image: image processing in Python」— PeerJ 2:e453
• scikit-image のサンプルギャラリー — カテゴリ別の注釈付きコード例