Image Processing with Python
Scientific image processing differs from photographic image editing in a critical way: the goal is quantification, not aesthetics. When a biologist processes microscopy images, the objective is to count cells, measure their areas, and quantify fluorescence intensity, not to make the image look better. This means every processing step must be scientifically justified (not just visually pleasing), reproducible (same parameters produce same results), and documented (another researcher can understand and replicate the analysis). Python's scriptable approach to image processing satisfies all three requirements, while manual editing in software like ImageJ or Photoshop satisfies none of them reliably.
Step 1: Load and Display Images
scikit-image loads images as NumPy arrays with consistent conventions. from skimage import io, then image = io.imread('sample.tif'). Grayscale images are 2D arrays with shape (height, width). Color images are 3D arrays with shape (height, width, 3) for RGB or (height, width, 4) for RGBA. Values range from 0 to 255 for 8-bit images, 0 to 65535 for 16-bit images (common in microscopy), or 0.0 to 1.0 for float images. Use skimage.img_as_float(image) to convert any format to float [0, 1] for consistent processing.
Display images with matplotlib: plt.imshow(image, cmap='gray') for grayscale (specify cmap='gray' explicitly, otherwise matplotlib applies a color map that misrepresents intensity). plt.colorbar() adds an intensity scale. For side-by-side comparison: fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5)), ax1.imshow(original, cmap='gray'), ax2.imshow(processed, cmap='gray'). Always include a scale bar or pixel size annotation when publishing microscopy images: ax.plot([10, 10 + scale_pixels], [height - 10, height - 10], 'w-', linewidth=2) draws a white line representing a known physical distance.
OpenCV uses BGR color order (not RGB), which is a constant source of bugs when mixing OpenCV with matplotlib. import cv2, then image = cv2.imread('photo.jpg') loads in BGR. Convert with cv2.cvtColor(image, cv2.COLOR_BGR2RGB) before displaying with matplotlib. For scientific work, scikit-image is generally preferred because it uses standard RGB ordering, represents images consistently as NumPy arrays, and provides algorithms specifically designed for quantitative analysis rather than computer vision applications.
Handle multi-channel and stack images common in microscopy. TIFF files from confocal microscopes often contain multiple channels (DAPI, GFP, mCherry) and multiple Z-slices. io.imread('stack.tif') loads the entire stack. Access individual channels: channel_0 = stack[:, :, 0] for multi-channel, or slice_5 = stack[5, :, :] for Z-stacks. tifffile.imread('microscopy.tif') (pip install tifffile) handles complex TIFF formats including OME-TIFF with metadata about pixel sizes, channel names, and acquisition parameters.
Step 2: Apply Filters and Enhancements
Gaussian blurring reduces noise while preserving edges. from skimage.filters import gaussian, then smoothed = gaussian(image, sigma=2). The sigma parameter controls the blur amount: sigma=1 for mild smoothing, sigma=3 to 5 for heavy smoothing. Choose sigma based on the noise level and the size of features you want to preserve: sigma should be smaller than the features of interest but larger than the noise pattern. For microscopy, sigma=1 to 2 usually removes shot noise without blurring cell boundaries.
Median filtering removes salt-and-pepper noise (random bright or dark pixels) while preserving edges better than Gaussian blur. from skimage.filters import median, from skimage.morphology import disk, then filtered = median(image, disk(3)). The disk(3) footprint defines a circular neighborhood of radius 3 pixels. Median filtering is preferred for images with impulse noise (hot pixels in CCD cameras, dead pixels in CMOS sensors) because it replaces each pixel with the median of its neighborhood, eliminating outlier pixels without smoothing edges.
Contrast enhancement makes features visible in low-contrast images. from skimage import exposure. exposure.rescale_intensity(image, in_range=(low, high)) stretches the intensity range to use the full display range. exposure.equalize_hist(image) redistributes intensities to achieve a flat histogram, maximizing visual contrast. exposure.equalize_adapthist(image, clip_limit=0.03) applies contrast-limited adaptive histogram equalization (CLAHE), which enhances local contrast without amplifying noise, ideal for images with uneven illumination such as large microscopy fields.
Background subtraction corrects uneven illumination. Create a background estimate by heavy Gaussian blurring (sigma=50 or larger): background = gaussian(image, sigma=50). Subtract: corrected = image - background (for additive background) or corrected = image / background (for multiplicative background, i.e., vignetting). The rolling-ball algorithm (from skimage.morphology import white_tophat, disk; corrected = white_tophat(image, disk(50))) provides a more robust background estimate that handles both bright and dark variations.
Step 3: Segment Objects of Interest
Thresholding separates foreground objects from background based on intensity. from skimage.filters import threshold_otsu, then thresh = threshold_otsu(image), binary = image > thresh. Otsu's method automatically finds the intensity value that best separates two classes (foreground and background) by maximizing between-class variance. The result is a binary image where True pixels are foreground. For images where a single global threshold fails (uneven illumination, varying background), use local thresholding: from skimage.filters import threshold_local, then binary = image > threshold_local(image, block_size=51).
Edge detection finds boundaries between regions. from skimage.feature import canny, then edges = canny(image, sigma=2). The sigma parameter controls edge sensitivity: lower sigma detects more edges (including noise), higher sigma detects only strong edges. The Canny edge detector applies Gaussian smoothing, computes gradients, suppresses non-maximum responses, and links edges using hysteresis thresholding, producing thin, connected edge contours. For gradient-based analysis without the full Canny pipeline, skimage.filters.sobel(image) computes edge magnitude directly.
Watershed segmentation separates touching objects. After thresholding produces a binary image where adjacent objects are merged, watershed splits them along intensity valleys. from skimage.segmentation import watershed, from skimage.feature import peak_local_max. Find markers at object centers: markers = peak_local_max(distance_transform, min_distance=20, labels=binary_image). Apply watershed: labels = watershed(-distance_transform, markers, mask=binary_image). Each connected object gets a unique label integer. This is the standard approach for separating touching cells, particles, or grains in microscopy images.
Region-based segmentation groups pixels by similarity rather than edges. from skimage.segmentation import slic, then segments = slic(image, n_segments=200, compactness=10). SLIC superpixels group pixels into perceptually uniform regions, useful as a preprocessing step for more complex analysis. from skimage.segmentation import felzenszwalb provides graph-based segmentation that adapts region size to local image structure. For interactive segmentation refinement, use the random walker algorithm: from skimage.segmentation import random_walker with user-provided seed points for foreground and background.
Step 4: Measure and Quantify Features
regionprops extracts measurements from labeled images. from skimage.measure import label, regionprops. labels = label(binary_image) assigns a unique integer to each connected component. props = regionprops(labels, intensity_image=original_image) computes properties for each labeled region. Access measurements: for region in props: print(region.area, region.perimeter, region.eccentricity, region.mean_intensity). Common measurements include area (pixel count), perimeter (boundary length), centroid (center position), bounding_box, major_axis_length, minor_axis_length, eccentricity (0 for circle, near 1 for elongated), and solidity (area / convex area).
Convert pixel measurements to physical units using the known pixel size. If the microscope metadata or a calibration image indicates that 1 pixel = 0.5 micrometers, then: area_um2 = region.area * pixel_size**2, perimeter_um = region.perimeter * pixel_size, diameter_um = region.equivalent_diameter_area * pixel_size. Always report measurements in physical units (micrometers, millimeters, square micrometers) rather than pixels, because pixel units depend on magnification and resolution settings that vary between acquisitions.
Compile measurements into a DataFrame for statistical analysis. data = [{'label': r.label, 'area': r.area * px**2, 'perimeter': r.perimeter * px, 'intensity': r.mean_intensity, 'eccentricity': r.eccentricity, 'centroid_y': r.centroid[0], 'centroid_x': r.centroid[1]} for r in props]. df = pd.DataFrame(data). Now you can filter (df[df['area'] > min_area] removes objects smaller than a threshold), compute statistics (df['area'].describe()), create plots (df['area'].hist()), and export results (df.to_csv('measurements.csv')).
Intensity measurements quantify fluorescence, staining, or absorption. region.mean_intensity is the average pixel value within the object. region.max_intensity and region.min_intensity give the range. For total fluorescence, multiply mean intensity by area. For background-corrected measurements, subtract the mean intensity of a background region from each object's measurement. For multi-channel images, compute regionprops on each channel separately using the same label image: props_ch1 = regionprops(labels, channel_1), props_ch2 = regionprops(labels, channel_2), then merge measurements by label.
Step 5: Batch Process Image Sets
Wrap your processing pipeline into a function that takes an image path and returns measurements. def process_image(path): image = io.imread(path), apply filters, segment, measure, return DataFrame of measurements. Apply to all images: results = pd.concat([process_image(f) for f in Path('images').glob('*.tif')], ignore_index=True). Add a source column: each sub-DataFrame includes the filename so you can trace measurements back to their source image. Save: results.to_csv('all_measurements.csv', index=False).
Parameter consistency is essential for batch processing. Use the same sigma, threshold method, minimum size filter, and measurement parameters for every image in a dataset. Define parameters once at the top of the script: SIGMA = 2, MIN_AREA = 100, BLOCK_SIZE = 51. Pass them to the processing function. Never adjust parameters per-image to "get better results" because this introduces experimenter bias: the measurements are no longer comparable across images. If a single set of parameters does not work for all images, the images need better standardization (consistent illumination, exposure, staining), not per-image parameter tuning.
Quality control verifies that batch processing succeeded. After processing, generate a summary: how many objects were detected per image? What is the distribution of areas? Are there images with zero detections (potential processing failures) or abnormally many detections (potential over-segmentation)? Create a visual QC sheet: for each image, save a side-by-side figure showing the original image and the segmentation overlay (labeled regions colored on top of the original). Review these QC figures to catch systematic failures before analyzing the measurements.
For large image sets (thousands of images, multi-gigabyte files), use memory-efficient processing. Load images one at a time rather than loading all into memory. Use dask-image for lazy loading and parallel processing of large image arrays. For 3D stacks (confocal Z-stacks, CT scans), process slice by slice when possible. Write results to disk incrementally rather than collecting all results in memory. Monitor memory usage with psutil: if memory exceeds a threshold, force garbage collection with gc.collect() and warn about potential issues.
Scientific image processing is about quantification, not aesthetics. Every processing step must use consistent parameters across all images in a dataset, and measurements should always be reported in physical units with documented processing parameters.