Scikit Image Processing

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

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)

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)

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)

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

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

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

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

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