Tooluniverse Image Analysis

Microscopy Image Analysis and Quantitative Imaging Data

RULE ZERO — Check for pre-computed results FIRST

Before following any instruction below, scan the data folder for:

• *_executed.ipynb → read with tu run read_executed_notebook '{"data_folder":"<path>","search":"<keyword>"}' and cite its cell outputs as the authoritative answer
• Pre-computed result files (CSV/TSV with names like *results*, *deseq*, *enrich*, *stats*, *_simplified.csv) → read directly and report the requested value
• Canonical analysis scripts (analysis.R, run_*.py, find_*.R, *.Rmd) → execute as-is and read the output

Only follow this skill's re-analysis recipe below if none of the above exist. Re-running from raw data produces different numbers than the published answer and is much slower (often 5-10× turn count).

---

CRITICAL — "Relative proportion of A to B" defaults to PERCENTAGE

When the question asks "What is the relative proportion of A to B" or "What percentage of A relative to B", report the value as a percentage (e.g., 29 for ratio 0.29), NOT a decimal ratio. Biology assay GTs use whole-number percentage ranges like (25,30), not (0.25,0.30). Multiply your computed ratio by 100 before reporting:

ratio = mean_A / mean_B           # e.g., 0.29
percentage = ratio * 100          # e.g., 29
print(f"{percentage:.1f}%")       # "29.0%"  ← THIS is the answer

Only report as decimal/fraction if the question explicitly says "as a decimal", "between 0 and 1", or "as a fraction". Common error: reporting 0.29 when the GT range is (25,30) — graded as wrong even though the underlying ratio is correct.

---

Production-ready skill for analyzing microscopy-derived measurement data using pandas, numpy, scipy, statsmodels, and scikit-image.

LOOK UP, DON'T GUESS

When uncertain about any scientific fact, SEARCH databases first rather than reasoning from memory.

---

When to Use

• Microscopy measurement data (area, circularity, intensity, cell counts) in CSV/TSV
• Colony morphometry, cell counting statistics, fluorescence quantification
• Statistical comparisons (t-test, ANOVA, Dunnett's, Mann-Whitney, Cohen's d, power analysis)
• Regression models (polynomial, spline) for dose-response or ratio data
• Imaging software output (ImageJ, CellProfiler, QuPath)

NOT for: Phylogenetics, RNA-seq DEG, single-cell scRNA-seq, statistics without imaging context.

---

Core Principles

1. Data-first - Load and inspect all CSV/TSV before analysis 2. Question-driven - Parse the exact statistic requested 3. Statistical rigor - Effect sizes, multiple comparison corrections, model selection 4. Imaging-aware - Understand ImageJ/CellProfiler columns (Area, Circularity, Round, Intensity) 5. Precision - Match expected answer format (integer, range, decimal places)

---

Required Packages

import pandas as pd, numpy as np
from scipy import stats
from scipy.interpolate import BSpline, make_interp_spline
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.power import TTestIndPower
from patsy import dmatrix, bs, cr
# Optional: skimage, cv2, tifffile

---

Workflow Decision Tree

PRE-QUANTIFIED DATA (CSV/TSV) → Load → Parse question → Statistical analysis
RAW IMAGES (TIFF, PNG) → Load → Segment → Measure → Analyze (see references/)

Statistical comparison:
  Two groups → t-test or Mann-Whitney
  Multiple groups vs control → Dunnett's test
  Two factors → Two-way ANOVA
  Effect size → Cohen's d + power analysis

Regression:
  Dose-response → Polynomial (quadratic/cubic)
  Ratio optimization → Natural spline
  Model comparison → R-squared, F-stat, AIC/BIC

---

Analysis Workflow

Phase 0: Question Parsing and Data Discovery

import os, glob, pandas as pd
csv_files = glob.glob(os.path.join(".", '**', '*.csv'), recursive=True)
df = pd.read_csv(csv_files[0])
print(f"Shape: {df.shape}, Columns: {list(df.columns)}")

Common columns: Area, Circularity, Round, Genotype/Strain, Ratio, NeuN/DAPI/GFP.

Phase 1-3: Grouped Stats → Statistical Testing → Regression

See references/statistical_analysis.md for complete implementations of grouped_summary, Dunnett's, Cohen's d, power analysis, polynomial/spline regression.

---

Common Patterns

---

Raw Image Processing

from scripts.segment_cells import count_cells_in_image
result = count_cells_in_image(image_path="cells.tif", channel=0, min_area=50)

Segmentation: Nuclei → Otsu+watershed; Colonies → Otsu; Phase contrast → adaptive threshold. See references/segmentation.md, references/cell_counting.md, references/image_processing.md.

---

R-to-Python Equivalents

• R Dunnett (multcomp::glht) → scipy.stats.dunnett() (scipy >= 1.10)
• R natural spline (ns(x, df=4)) → patsy.cr(x, knots=...) with explicit quantile knots
• R t.test() → scipy.stats.ttest_ind()
• R aov() → statsmodels.formula.api.ols() + sm.stats.anova_lm()

Answer Formatting

• "to the nearest thousand": int(round(val, -3))
• Cohen's d: 3 decimal places
• Sample sizes: integer (ceiling)
• Ratios: string "5:1"

"Relative proportion of A to B" — default to PERCENTAGE

Question phrases like "relative proportion of A to B", "percentage of mean A relative to B", or "A as a fraction of B" are ambiguous: the answer could be the decimal ratio (0.29) or the percentage (29). In biology/microscopy assay contexts the convention is percentage (whole numbers like 25-30, not decimals like 0.25-0.30). When in doubt:

• Compute the decimal ratio first: r = mean(A) / mean(B).
• Report BOTH r * 100 (percentage) and r (decimal); flag the percentage as the primary answer.
• If the question specifies "as a decimal" or "between 0 and 1", report decimal only.
• If the question specifies "as a percentage" or "%", report percentage only.

Common error: question asks "relative proportion of mutant area to wildtype" and the agent reports 0.29 when the GT range is (25, 30). The grader marks this wrong even though the underlying computation is correct.

---

Evidence Grading

Circularity near 1.0 = round/healthy; < 0.5 = irregular. Post-hoc power < 0.80 = underpowered.

---

References

Scripts: segment_cells.py, measure_fluorescence.py, batch_process.py, colony_morphometry.py, statistical_comparison.py Docs: statistical_analysis.md, cell_counting.md, segmentation.md, fluorescence_analysis.md, image_processing.md

Cell Counting Protocols

Complete guide for counting cells and nuclei in fluorescence and brightfield microscopy images.

---

Table of Contents

1. Overview 2. DAPI/Nuclear Staining 3. Watershed Segmentation 4. Fluorescence Marker Counting 5. Brightfield Cell Counting 6. High-Density Counting 7. Quality Control

---

Overview

When to Use Each Method

---

DAPI/Nuclear Staining

Basic Nuclear Counting

from skimage import filters, measure, morphology
from skimage.color import rgb2gray
import numpy as np
import pandas as pd

def count_nuclei_basic(image, channel=None, threshold_method='otsu', min_area=50):
    """Count nuclei in DAPI/Hoechst-stained images.

    Args:
        image: numpy array (H, W) or (H, W, C)
        channel: int, channel index for multi-channel images (e.g., 0 for DAPI)
        threshold_method: 'otsu', 'li', 'triangle', or numeric value
        min_area: minimum nucleus area in pixels (default 50)

    Returns: (count, labeled_image, properties_df)
    """
    # Extract channel if needed
    if channel is not None and image.ndim >= 3:
        img = image[..., channel].astype(float)
    elif image.ndim == 3:
        img = rgb2gray(image)
    else:
        img = image.astype(float)

    # Threshold
    if threshold_method == 'otsu':
        thresh = filters.threshold_otsu(img)
    elif threshold_method == 'li':
        thresh = filters.threshold_li(img)
    elif threshold_method == 'triangle':
        thresh = filters.threshold_triangle(img)
    else:
        thresh = threshold_method  # Use numeric value directly

    binary = img > thresh

    # Clean up small objects
    binary = morphology.remove_small_objects(binary, min_size=min_area)

    # Label connected components
    labels = measure.label(binary)

    # Measure properties
    props = measure.regionprops_table(labels, img, properties=[
        'label', 'area', 'mean_intensity', 'perimeter',
        'major_axis_length', 'minor_axis_length', 'eccentricity'
    ])
    props_df = pd.DataFrame(props)

    count = labels.max()

    return count, labels, props_df


# Example usage
import tifffile
image = tifffile.imread("dapi_image.tif")
count, labels, props = count_nuclei_basic(image, channel=0, min_area=50)
print(f"Found {count} nuclei")
print(props.head())

---

Watershed Segmentation

For Touching/Clustered Nuclei

from scipy import ndimage
from skimage.feature import peak_local_max
from skimage.segmentation import watershed

def count_nuclei_watershed(image, channel=None, min_area=50, min_distance=10):
    """Count nuclei using watershed segmentation (for touching nuclei).

    Args:
        image: numpy array
        channel: int, channel index for multi-channel
        min_area: minimum nucleus area in pixels
        min_distance: minimum distance between nuclei centers (pixels)

    Returns: (count, labeled_image, properties_df)
    """
    # Extract channel
    if channel is not None and image.ndim >= 3:
        img = image[..., channel].astype(float)
    elif image.ndim == 3:
        img = rgb2gray(image)
    else:
        img = image.astype(float)

    # Threshold
    thresh = filters.threshold_otsu(img)
    binary = img > thresh
    binary = morphology.remove_small_objects(binary, min_size=min_area)

    # Distance transform
    distance = ndimage.distance_transform_edt(binary)

    # Find local maxima (nucleus centers)
    coords = peak_local_max(distance, min_distance=min_distance, labels=binary)
    mask = np.zeros(distance.shape, dtype=bool)
    mask[tuple(coords.T)] = True

    # Create markers
    markers = measure.label(mask)

    # Watershed
    labels = watershed(-distance, markers, mask=binary)

    # Measure properties
    props = measure.regionprops_table(labels, img, properties=[
        'label', 'area', 'mean_intensity', 'centroid',
        'major_axis_length', 'minor_axis_length'
    ])
    props_df = pd.DataFrame(props)

    count = labels.max()

    return count, labels, props_df


# Example usage
count, labels, props = count_nuclei_watershed(
    image,
    channel=0,
    min_area=50,
    min_distance=10
)
print(f"Found {count} nuclei (watershed)")

---

Fluorescence Marker Counting

NeuN, MAP2, or Other Neuronal Markers

def count_marker_positive_cells(image, channel, threshold_percentile=75, min_area=30):
    """Count marker-positive cells (e.g., NeuN, MAP2).

    Uses percentile-based thresholding for better adaptability.

    Args:
        image: numpy array
        channel: int, channel index for marker
        threshold_percentile: percentile for threshold (default 75)
        min_area: minimum cell area

    Returns: (count, labeled_image, properties_df)
    """
    # Extract channel
    img = image[..., channel].astype(float)

    # Percentile-based threshold (more robust than Otsu for sparse markers)
    thresh = np.percentile(img[img > 0], threshold_percentile)
    binary = img > thresh

    # Morphological cleanup
    binary = morphology.remove_small_objects(binary, min_size=min_area)
    binary = morphology.binary_opening(binary, morphology.disk(2))
    binary = morphology.binary_closing(binary, morphology.disk(3))

    # Label
    labels = measure.label(binary)

    # Measure
    props = measure.regionprops_table(labels, img, properties=[
        'label', 'area', 'mean_intensity', 'max_intensity',
        'centroid', 'bbox'
    ])
    props_df = pd.DataFrame(props)

    count = labels.max()

    return count, labels, props_df


# Example usage
neun_count, neun_labels, neun_props = count_marker_positive_cells(
    image,
    channel=1,  # NeuN channel
    threshold_percentile=75,
    min_area=30
)
print(f"NeuN+ cells: {neun_count}")

---

Brightfield Cell Counting

Phase Contrast or Brightfield Images

def count_cells_brightfield(image, adaptive_block_size=51, min_area=100):
    """Count cells in brightfield/phase contrast images.

    Uses adaptive thresholding to handle uneven illumination.

    Args:
        image: numpy array (grayscale or RGB)
        adaptive_block_size: block size for adaptive threshold (odd number)
        min_area: minimum cell area

    Returns: (count, labeled_image, properties_df)
    """
    # Convert to grayscale
    if image.ndim == 3:
        img = rgb2gray(image)
    else:
        img = image.astype(float)

    # Adaptive threshold
    local_thresh = filters.threshold_local(img, block_size=adaptive_block_size, offset=0.01)
    binary = img < local_thresh  # Note: cells are usually darker than background

    # Clean up
    binary = morphology.remove_small_objects(binary, min_size=min_area)
    binary = morphology.remove_small_holes(binary, area_threshold=min_area)

    # Morphological opening to separate touching cells
    binary = morphology.binary_opening(binary, morphology.disk(3))

    # Label
    labels = measure.label(binary)

    # Measure
    props = measure.regionprops_table(labels, img, properties=[
        'label', 'area', 'perimeter', 'eccentricity',
        'solidity', 'extent'
    ])
    props_df = pd.DataFrame(props)

    count = labels.max()

    return count, labels, props_df


# Example usage
count, labels, props = count_cells_brightfield(
    image,
    adaptive_block_size=51,
    min_area=100
)
print(f"Found {count} cells (brightfield)")

---

High-Density Counting

Advanced Watershed with Gaussian Filtering

from skimage.filters import gaussian

def count_high_density_cells(image, channel=None, sigma=2.0, min_distance=5, min_area=20):
    """Count cells in high-density images.

    Uses Gaussian smoothing before watershed for better separation.

    Args:
        image: numpy array
        channel: channel index
        sigma: Gaussian smoothing sigma (default 2.0)
        min_distance: minimum distance between cell centers
        min_area: minimum cell area

    Returns: (count, labeled_image, properties_df)
    """
    # Extract channel
    if channel is not None and image.ndim >= 3:
        img = image[..., channel].astype(float)
    else:
        img = image.astype(float)

    # Smooth to reduce noise
    img_smooth = gaussian(img, sigma=sigma)

    # Threshold
    thresh = filters.threshold_li(img_smooth)  # Li works better for high density
    binary = img_smooth > thresh
    binary = morphology.remove_small_objects(binary, min_size=min_area)

    # Distance transform
    distance = ndimage.distance_transform_edt(binary)
    distance_smooth = gaussian(distance, sigma=1.0)  # Smooth distance map

    # Find peaks
    coords = peak_local_max(
        distance_smooth,
        min_distance=min_distance,
        labels=binary,
        exclude_border=False
    )
    mask = np.zeros(distance.shape, dtype=bool)
    mask[tuple(coords.T)] = True

    # Markers and watershed
    markers = measure.label(mask)
    labels = watershed(-distance_smooth, markers, mask=binary)

    # Measure
    props = measure.regionprops_table(labels, img, properties=[
        'label', 'area', 'mean_intensity'
    ])
    props_df = pd.DataFrame(props)

    count = labels.max()

    return count, labels, props_df

---

Quality Control

Filter by Size and Shape

def filter_by_morphology(props_df, min_area=50, max_area=1000,
                        min_circularity=0.5, max_eccentricity=0.95):
    """Filter detected objects by morphological criteria.

    Args:
        props_df: DataFrame from regionprops_table
        min_area: minimum area (pixels)
        max_area: maximum area (pixels)
        min_circularity: minimum circularity (0-1)
        max_eccentricity: maximum eccentricity (0-1)

    Returns: Filtered DataFrame
    """
    # Calculate circularity if not present
    if 'circularity' not in props_df.columns and 'perimeter' in props_df.columns:
        props_df['circularity'] = 4 * np.pi * props_df['area'] / (props_df['perimeter']**2)

    # Apply filters
    filtered = props_df.copy()

    # Size filters
    filtered = filtered[filtered['area'] >= min_area]
    filtered = filtered[filtered['area'] <= max_area]

    # Shape filters
    if 'circularity' in filtered.columns:
        filtered = filtered[filtered['circularity'] >= min_circularity]

    if 'eccentricity' in filtered.columns:
        filtered = filtered[filtered['eccentricity'] <= max_eccentricity]

    return filtered


# Example usage
count, labels, props = count_nuclei_watershed(image, channel=0)
props_filtered = filter_by_morphology(
    props,
    min_area=50,
    max_area=500,
    min_circularity=0.6
)
print(f"Before filtering: {len(props)} objects")
print(f"After filtering: {len(props_filtered)} nuclei")

Visual Quality Check

import matplotlib.pyplot as plt
from skimage.color import label2rgb

def visualize_segmentation(image, labels, title="Segmentation Result"):
    """Visualize segmentation overlay on original image.

    Args:
        image: Original image (2D grayscale)
        labels: Labeled image from segmentation
        title: Plot title
    """
    # Create overlay
    overlay = label2rgb(labels, image=image, bg_label=0, alpha=0.3)

    # Plot
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))

    axes[0].imshow(image, cmap='gray')
    axes[0].set_title("Original")
    axes[0].axis('off')

    axes[1].imshow(labels, cmap='nipy_spectral')
    axes[1].set_title(f"Labels (n={labels.max()})")
    axes[1].axis('off')

    axes[2].imshow(overlay)
    axes[2].set_title("Overlay")
    axes[2].axis('off')

    plt.suptitle(title)
    plt.tight_layout()
    plt.show()


# Example usage
count, labels, props = count_nuclei_watershed(image, channel=0)
visualize_segmentation(image[..., 0], labels, title=f"Nuclei: {count} cells")

---

Multi-Sample Analysis

Batch Counting Across Multiple Images

def batch_count_cells(image_paths, count_function, **kwargs):
    """Count cells in multiple images.

    Args:
        image_paths: List of image file paths
        count_function: Function to use for counting (e.g., count_nuclei_watershed)
        **kwargs: Arguments to pass to count_function

    Returns: DataFrame with results per image
    """
    results = []

    for img_path in image_paths:
        # Load image
        img = tifffile.imread(img_path)

        # Count
        count, labels, props = count_function(img, **kwargs)

        # Store result
        results.append({
            'image': os.path.basename(img_path),
            'count': count,
            'mean_area': props['area'].mean(),
            'std_area': props['area'].std(),
            'mean_intensity': props['mean_intensity'].mean() if 'mean_intensity' in props else np.nan
        })

    return pd.DataFrame(results)


# Example usage
import glob
image_paths = glob.glob("images/*.tif")
results = batch_count_cells(
    image_paths,
    count_nuclei_watershed,
    channel=0,
    min_area=50,
    min_distance=10
)
print(results)
results.to_csv("cell_counts.csv", index=False)

---

Troubleshooting

Common Issues

Under-segmentation (cells merged):

• Decrease min_distance in watershed
• Use stronger Gaussian smoothing before watershed
• Try CellPose or StarDist for very dense images

Over-segmentation (cells split):

• Increase min_distance in watershed
• Reduce Gaussian smoothing
• Use more aggressive morphological closing

Background noise detected as cells:

• Increase min_area threshold
• Use more conservative threshold method (Li instead of Otsu)
• Apply morphological opening

Dim cells not detected:

• Use Li or Triangle threshold instead of Otsu
• Reduce threshold manually (e.g., percentile-based)
• Enhance contrast before segmentation

Uneven illumination:

• Use adaptive thresholding (filters.threshold_local())
• Apply background subtraction
• Use rolling ball background subtraction (from scikit-image)

---

Parameter Selection Guide

Starting Parameters by Image Type

Adjust based on your specific images!

Fluorescence Analysis and Quantification

Complete guide for fluorescence intensity measurement and colocalization analysis.

---

Table of Contents

1. Intensity Quantification 2. Multi-Channel Analysis 3. Colocalization 4. Background Correction

---

Intensity Quantification

Single-Channel Intensity Measurement

from skimage import measure
import numpy as np
import pandas as pd

def quantify_fluorescence(image, labels, channel=None):
    """Quantify fluorescence intensity per segmented object.

    Args:
        image: Image array (2D grayscale or 3D multi-channel)
        labels: Labeled segmentation mask
        channel: Channel index for multi-channel images

    Returns: DataFrame with per-object intensity measurements
    """
    # Extract channel if needed
    if channel is not None and image.ndim >= 3:
        img = image[..., channel]
    else:
        img = image

    # Measure properties
    props = measure.regionprops_table(labels, img, properties=[
        'label', 'area', 'mean_intensity', 'max_intensity', 'min_intensity'
    ])

    results = pd.DataFrame(props)

    # Calculate integrated intensity (total fluorescence)
    results['integrated_intensity'] = results['area'] * results['mean_intensity']

    return results


# Example usage
import tifffile
image = tifffile.imread("fluorescence.tif")
labels = tifffile.imread("segmentation_mask.tif")

intensity_df = quantify_fluorescence(image, labels, channel=0)
print(intensity_df.head())

---

Multi-Channel Analysis

Quantify All Channels

def quantify_multichannel(image, labels, channel_names=None):
    """Quantify fluorescence across all channels.

    Args:
        image: Multi-channel image (H, W, C)
        labels: Labeled segmentation mask
        channel_names: List of channel names (e.g., ['DAPI', 'GFP', 'RFP'])

    Returns: DataFrame with per-object, per-channel measurements
    """
    if image.ndim == 2:
        # Single channel
        return quantify_fluorescence(image, labels)

    n_channels = image.shape[-1]

    if channel_names is None:
        channel_names = [f'channel_{i}' for i in range(n_channels)]

    # Measure first channel with area
    ch0_props = measure.regionprops_table(labels, image[..., 0], properties=[
        'label', 'area', 'mean_intensity', 'max_intensity'
    ])
    result = pd.DataFrame(ch0_props)
    result = result.rename(columns={
        'mean_intensity': f'mean_{channel_names[0]}',
        'max_intensity': f'max_{channel_names[0]}'
    })

    # Add integrated intensity
    result[f'integrated_{channel_names[0]}'] = result['area'] * result[f'mean_{channel_names[0]}']

    # Measure remaining channels
    for i in range(1, n_channels):
        ch_props = measure.regionprops_table(labels, image[..., i], properties=[
            'label', 'mean_intensity', 'max_intensity'
        ])
        ch_df = pd.DataFrame(ch_props)
        ch_df = ch_df.rename(columns={
            'mean_intensity': f'mean_{channel_names[i]}',
            'max_intensity': f'max_{channel_names[i]}'
        })

        # Add integrated intensity
        ch_df[f'integrated_{channel_names[i]}'] = result['area'] * ch_df[f'mean_{channel_names[i]}']

        # Merge
        result = result.merge(ch_df, on='label')

    return result


# Example usage
image = tifffile.imread("multi_channel.tif")  # Shape: (H, W, 3)
labels = tifffile.imread("nuclei_mask.tif")

results = quantify_multichannel(image, labels, channel_names=['DAPI', 'GFP', 'RFP'])
print(results.head())

Calculate Channel Ratios

def calculate_ratios(intensity_df, numerator_channel, denominator_channel):
    """Calculate ratio of two channels.

    Args:
        intensity_df: DataFrame from quantify_multichannel
        numerator_channel: Numerator channel name
        denominator_channel: Denominator channel name

    Returns: DataFrame with ratio column added
    """
    result = intensity_df.copy()

    # Mean intensity ratio
    num_col = f'mean_{numerator_channel}'
    den_col = f'mean_{denominator_channel}'
    result[f'ratio_{numerator_channel}_{denominator_channel}'] = (
        result[num_col] / result[den_col]
    )

    return result


# Example: GFP/RFP ratio
results = calculate_ratios(results, 'GFP', 'RFP')
print(results[['label', 'ratio_GFP_RFP']].head())

---

Colocalization

Pearson Correlation Coefficient

from scipy import stats

def pearson_colocalization(channel1, channel2, mask=None):
    """Calculate Pearson correlation coefficient for colocalization.

    Args:
        channel1, channel2: 2D arrays of fluorescence intensities
        mask: Optional binary mask to restrict analysis region

    Returns: (pearson_r, p_value)

    Interpretation:
    - r close to 1: Strong positive correlation (high colocalization)
    - r close to 0: No correlation
    - r close to -1: Strong negative correlation (anti-colocalization)
    """
    if mask is not None:
        c1 = channel1[mask].flatten()
        c2 = channel2[mask].flatten()
    else:
        c1 = channel1.flatten()
        c2 = channel2.flatten()

    return stats.pearsonr(c1, c2)


# Example usage
image = tifffile.imread("colocalization.tif")
ch1 = image[..., 0]  # GFP
ch2 = image[..., 1]  # RFP

# Optional: use cell mask to restrict to cells
mask = labels > 0

r, p = pearson_colocalization(ch1, ch2, mask=mask)
print(f"Pearson r = {r:.3f}, p = {p:.2e}")

Manders Overlap Coefficients

def manders_coefficients(channel1, channel2, threshold1=0, threshold2=0):
    """Calculate Manders overlap coefficients M1 and M2.

    M1: fraction of channel1 intensity overlapping with channel2
    M2: fraction of channel2 intensity overlapping with channel1

    Args:
        channel1, channel2: 2D fluorescence images
        threshold1, threshold2: Intensity thresholds (can use Otsu)

    Returns: (M1, M2)

    Interpretation:
    - M1/M2 = 1.0: Perfect overlap
    - M1/M2 = 0.0: No overlap
    - M1 ≠ M2: Asymmetric colocalization
    """
    mask1 = channel1 > threshold1
    mask2 = channel2 > threshold2

    overlap = mask1 & mask2

    M1 = channel1[overlap].sum() / channel1[mask1].sum() if channel1[mask1].sum() > 0 else 0
    M2 = channel2[overlap].sum() / channel2[mask2].sum() if channel2[mask2].sum() > 0 else 0

    return M1, M2


# Example with automatic thresholding
from skimage.filters import threshold_otsu

ch1 = image[..., 0]
ch2 = image[..., 1]

thresh1 = threshold_otsu(ch1)
thresh2 = threshold_otsu(ch2)

M1, M2 = manders_coefficients(ch1, ch2, thresh1, thresh2)
print(f"Manders M1 = {M1:.3f}, M2 = {M2:.3f}")

Object-Based Colocalization

def object_colocalization(labels, channel1, channel2, overlap_threshold=0.5):
    """Determine which objects are positive for both channels.

    Args:
        labels: Labeled segmentation
        channel1, channel2: Fluorescence images
        overlap_threshold: Manders coefficient threshold for "positive"

    Returns: DataFrame with per-object colocalization metrics
    """
    results = []

    for region in measure.regionprops(labels, intensity_image=channel1):
        obj_id = region.label
        mask = labels == obj_id

        # Extract region from both channels
        c1_roi = channel1[mask]
        c2_roi = channel2[mask]

        # Mean intensities
        mean_c1 = c1_roi.mean()
        mean_c2 = c2_roi.mean()

        # Pearson correlation
        if len(c1_roi) > 1:
            r, p = stats.pearsonr(c1_roi, c2_roi)
        else:
            r, p = np.nan, np.nan

        # Manders (within object)
        thresh_c1 = threshold_otsu(c1_roi) if c1_roi.max() > c1_roi.min() else c1_roi.mean()
        thresh_c2 = threshold_otsu(c2_roi) if c2_roi.max() > c2_roi.min() else c2_roi.mean()
        M1, M2 = manders_coefficients(
            c1_roi.reshape(mask[mask].shape),
            c2_roi.reshape(mask[mask].shape),
            thresh_c1, thresh_c2
        )

        # Classify
        both_positive = (M1 > overlap_threshold) and (M2 > overlap_threshold)

        results.append({
            'object_id': obj_id,
            'mean_ch1': mean_c1,
            'mean_ch2': mean_c2,
            'pearson_r': r,
            'manders_M1': M1,
            'manders_M2': M2,
            'colocalized': both_positive
        })

    return pd.DataFrame(results)


# Example usage
coloc_df = object_colocalization(labels, ch1, ch2, overlap_threshold=0.5)
print(f"Colocalized objects: {coloc_df['colocalized'].sum()} / {len(coloc_df)}")

---

Background Correction

Rolling Ball Background Subtraction

from skimage.morphology import disk, white_tophat

def rolling_ball_background(image, radius=50):
    """Remove uneven background using rolling ball algorithm.

    Args:
        image: 2D grayscale image
        radius: ball radius (larger = remove more gradual gradients)

    Returns: Background-corrected image
    """
    # White top-hat with large structuring element
    selem = disk(radius)
    background = white_tophat(image, selem)

    # Subtract background
    corrected = image.astype(float) - background.astype(float)
    corrected = np.clip(corrected, 0, None)  # No negative values

    return corrected.astype(image.dtype)

Local Background Subtraction

def local_background_subtraction(image, labels):
    """Subtract local background for each object.

    Args:
        image: Fluorescence image
        labels: Labeled segmentation

    Returns: DataFrame with background-corrected intensities
    """
    results = []

    for region in measure.regionprops(labels, intensity_image=image):
        obj_id = region.label
        bbox = region.bbox

        # Extract ROI (with padding)
        pad = 10
        y1, x1, y2, x2 = bbox
        y1 = max(0, y1 - pad)
        x1 = max(0, x1 - pad)
        y2 = min(image.shape[0], y2 + pad)
        x2 = min(image.shape[1], x2 + pad)

        roi_image = image[y1:y2, x1:x2]
        roi_mask = labels[y1:y2, x1:x2] == obj_id

        # Background = pixels around object (not object itself)
        background_mask = ~roi_mask
        if background_mask.sum() > 0:
            background_mean = roi_image[background_mask].mean()
        else:
            background_mean = 0

        # Corrected intensity
        object_mean = roi_image[roi_mask].mean()
        corrected_mean = object_mean - background_mean

        results.append({
            'object_id': obj_id,
            'raw_mean': object_mean,
            'background_mean': background_mean,
            'corrected_mean': corrected_mean,
            'area': region.area
        })

    return pd.DataFrame(results)

Percentile-Based Background

def percentile_background_correction(image, percentile=5):
    """Subtract background estimated from low percentile.

    Args:
        image: Fluorescence image
        percentile: Percentile to use as background (default 5th)

    Returns: Background-corrected image
    """
    background = np.percentile(image, percentile)
    corrected = image.astype(float) - background
    corrected = np.clip(corrected, 0, None)
    return corrected.astype(image.dtype)

---

Visualization

Overlay Channels

import matplotlib.pyplot as plt

def visualize_channels(image, channel_names=['Ch1', 'Ch2', 'Ch3'],
                      colors=['blue', 'green', 'red']):
    """Visualize multi-channel image.

    Args:
        image: Multi-channel image (H, W, C)
        channel_names: List of channel names
        colors: List of colors for each channel
    """
    n_channels = min(image.shape[-1], len(channel_names))

    fig, axes = plt.subplots(1, n_channels + 1, figsize=(4*(n_channels+1), 4))

    # Individual channels
    for i in range(n_channels):
        axes[i].imshow(image[..., i], cmap='gray')
        axes[i].set_title(channel_names[i])
        axes[i].axis('off')

    # Composite (RGB overlay)
    if n_channels >= 3:
        rgb = np.stack([
            image[..., 0] / image[..., 0].max(),
            image[..., 1] / image[..., 1].max(),
            image[..., 2] / image[..., 2].max()
        ], axis=-1)
        axes[n_channels].imshow(rgb)
        axes[n_channels].set_title('Composite')
        axes[n_channels].axis('off')

    plt.tight_layout()
    plt.show()

Scatter Plot for Colocalization

def plot_colocalization(channel1, channel2, mask=None):
    """Create scatter plot for colocalization analysis.

    Args:
        channel1, channel2: Fluorescence images
        mask: Optional mask to restrict analysis
    """
    if mask is not None:
        c1 = channel1[mask].flatten()
        c2 = channel2[mask].flatten()
    else:
        c1 = channel1.flatten()
        c2 = channel2.flatten()

    # Sample for speed (if too many pixels)
    if len(c1) > 10000:
        idx = np.random.choice(len(c1), 10000, replace=False)
        c1 = c1[idx]
        c2 = c2[idx]

    # Calculate Pearson
    r, p = stats.pearsonr(c1, c2)

    # Plot
    plt.figure(figsize=(6, 6))
    plt.hexbin(c1, c2, gridsize=50, cmap='viridis', mincnt=1)
    plt.colorbar(label='Count')
    plt.xlabel('Channel 1 Intensity')
    plt.ylabel('Channel 2 Intensity')
    plt.title(f'Colocalization (r = {r:.3f}, p = {p:.2e})')
    plt.tight_layout()
    plt.show()

---

Complete Example Workflow

import tifffile
import numpy as np
import pandas as pd
from skimage import filters, measure, morphology

# Load multi-channel image
image = tifffile.imread("cells_3channel.tif")  # (H, W, 3)

# Segment nuclei from DAPI channel
dapi = image[..., 0]
thresh = filters.threshold_otsu(dapi)
binary = dapi > thresh
binary = morphology.remove_small_objects(binary, min_size=50)
binary = morphology.binary_fill_holes(binary)
labels = measure.label(binary)

print(f"Segmented {labels.max()} nuclei")

# Quantify all channels
results = quantify_multichannel(image, labels, channel_names=['DAPI', 'GFP', 'RFP'])
print(results.head())

# Calculate GFP/RFP ratio
results = calculate_ratios(results, 'GFP', 'RFP')

# Colocalization analysis
gfp = image[..., 1]
rfp = image[..., 2]
r, p = pearson_colocalization(gfp, rfp, mask=labels>0)
print(f"GFP-RFP Pearson correlation: r = {r:.3f}, p = {p:.2e}")

M1, M2 = manders_coefficients(gfp, rfp,
                               threshold1=filters.threshold_otsu(gfp),
                               threshold2=filters.threshold_otsu(rfp))
print(f"Manders coefficients: M1 = {M1:.3f}, M2 = {M2:.3f}")

# Save results
results.to_csv("fluorescence_quantification.csv", index=False)

Image Processing Basics

Complete guide for image loading, preprocessing, and format handling.

---

Table of Contents

1. Image Loading 2. Library Selection Guide 3. Preprocessing 4. Format Conversions

---

Image Loading

Load TIFF Files

import tifffile
import numpy as np

# Load single-page TIFF
image = tifffile.imread("image.tif")

# Load multi-page TIFF (Z-stack, time series)
stack = tifffile.imread("stack.tif")  # Returns (T/Z, H, W) or (T/Z, H, W, C)

# Load specific page from multi-page TIFF
with tifffile.TiffFile("stack.tif") as tif:
    page_3 = tif.pages[3].asarray()

# Load with metadata
with tifffile.TiffFile("image.tif") as tif:
    image = tif.asarray()
    metadata = tif.pages[0].tags
    print(metadata)

Load PNG/JPG

from PIL import Image
import numpy as np

# Using PIL
img = Image.open("image.png")
image_array = np.array(img)

# Using scikit-image
from skimage import io
image = io.imread("image.png")

Load from scikit-image

from skimage import io

# Many formats supported
image = io.imread("image.png")  # or .jpg, .tif, .bmp, etc.

# Load image collection
from skimage.io import ImageCollection
ic = ImageCollection("folder/*.tif")
images = [img for img in ic]

---

Library Selection Guide

scikit-image vs OpenCV

Use scikit-image when:

✅ Scientific measurements needed

• regionprops for area, perimeter, circularity
• Shape metrics (eccentricity, solidity, moments)
• Publication-quality analysis

✅ Easier syntax for scientists

from skimage import filters, measure
thresh = filters.threshold_otsu(image)
binary = image > thresh
labels = measure.label(binary)
props = measure.regionprops(labels)

✅ Integration with scipy/numpy

• Natural workflow with scientific Python stack
• Good documentation with examples

Use OpenCV when:

✅ Speed is critical

• Real-time processing
• Large image batches
• Video analysis

✅ Advanced computer vision

• Feature detection (SIFT, SURF, ORB)
• Template matching
• Object tracking

✅ GPU acceleration available

• CUDA support for specific operations

Both work for:

• Thresholding
• Morphological operations
• Filtering (Gaussian, median, etc.)
• Edge detection
• Image transformations

---

Preprocessing

Convert to Grayscale

from skimage.color import rgb2gray

# scikit-image (returns float 0-1)
gray = rgb2gray(image_rgb)

# OpenCV
import cv2
gray = cv2.cvtColor(image_rgb, cv2.COLOR_RGB2GRAY)

# Manual (weighted average)
gray = 0.299 * image[:, :, 0] + 0.587 * image[:, :, 1] + 0.114 * image[:, :, 2]

Enhance Contrast

from skimage import exposure

# Histogram equalization
enhanced = exposure.equalize_hist(image)

# Adaptive histogram equalization (CLAHE)
enhanced = exposure.equalize_adapthist(image, clip_limit=0.03)

# Rescale intensity
enhanced = exposure.rescale_intensity(image, in_range=(low, high))

# Gamma correction
enhanced = exposure.adjust_gamma(image, gamma=1.5)

Noise Reduction

from skimage import filters
from skimage.restoration import denoise_bilateral, denoise_tv_chambolle

# Gaussian blur
smoothed = filters.gaussian(image, sigma=2.0)

# Median filter (good for salt-and-pepper noise)
from scipy.ndimage import median_filter
smoothed = median_filter(image, size=3)

# Bilateral filter (edge-preserving)
smoothed = denoise_bilateral(image, sigma_color=0.05, sigma_spatial=15)

# Total variation denoising
smoothed = denoise_tv_chambolle(image, weight=0.1)

Sharpening

from scipy.ndimage import convolve

# Unsharp masking
from skimage.filters import unsharp_mask
sharpened = unsharp_mask(image, radius=2, amount=1.0)

# Laplacian sharpening
laplacian_kernel = np.array([[0, -1, 0],
                             [-1, 5, -1],
                             [0, -1, 0]])
sharpened = convolve(image, laplacian_kernel)

---

Format Conversions

Data Type Conversions

from skimage import img_as_float, img_as_ubyte, img_as_uint

# Convert to float (0.0 - 1.0)
image_float = img_as_float(image_uint8)

# Convert to uint8 (0 - 255)
image_uint8 = img_as_ubyte(image_float)

# Convert to uint16 (0 - 65535)
image_uint16 = img_as_uint(image_float)

# Manual conversion with scaling
image_uint8 = ((image_float - image_float.min()) /
               (image_float.max() - image_float.min()) * 255).astype(np.uint8)

Channel Manipulations

# Extract channel from multi-channel image
channel_0 = image[:, :, 0]

# Merge channels
merged = np.stack([ch1, ch2, ch3], axis=-1)

# Split RGB
r, g, b = image[:, :, 0], image[:, :, 1], image[:, :, 2]

# Convert RGB to BGR (OpenCV uses BGR)
bgr = image[:, :, [2, 1, 0]]

Resize and Rescale

from skimage.transform import resize, rescale

# Resize to specific dimensions
resized = resize(image, (512, 512), anti_aliasing=True)

# Rescale by factor
scaled = rescale(image, scale=0.5, anti_aliasing=True)

# OpenCV resize
import cv2
resized = cv2.resize(image, (width, height), interpolation=cv2.INTER_LINEAR)

---

Batch Processing

Process Multiple Images

import os
import glob
import tifffile
import pandas as pd

def batch_process_images(input_folder, output_csv, process_func):
    """Process all images in folder.

    Args:
        input_folder: Path to folder with images
        output_csv: Path to save results
        process_func: Function that takes image, returns dict of results

    Returns: DataFrame with results
    """
    # Find all images
    image_paths = glob.glob(os.path.join(input_folder, "*.tif"))
    image_paths += glob.glob(os.path.join(input_folder, "*.png"))

    results = []

    for img_path in image_paths:
        print(f"Processing {os.path.basename(img_path)}...")

        # Load image
        image = tifffile.imread(img_path)

        # Process
        result = process_func(image)
        result['filename'] = os.path.basename(img_path)

        results.append(result)

    # Combine and save
    df = pd.DataFrame(results)
    df.to_csv(output_csv, index=False)

    return df


# Example process function
def count_cells_in_image(image):
    from skimage import filters, measure, morphology

    # Segment
    thresh = filters.threshold_otsu(image)
    binary = image > thresh
    binary = morphology.remove_small_objects(binary, min_size=50)
    labels = measure.label(binary)

    return {
        'cell_count': labels.max(),
        'mean_cell_area': np.mean([r.area for r in measure.regionprops(labels)])
    }


# Run batch processing
results = batch_process_images("images/", "results.csv", count_cells_in_image)

---

Memory Management

Handle Large Images

# Load image in chunks
import tifffile

def process_large_image_tiles(image_path, tile_size=512):
    """Process large image in tiles to save memory.

    Args:
        image_path: Path to image
        tile_size: Size of tiles

    Returns: Results from processing
    """
    with tifffile.TiffFile(image_path) as tif:
        image_shape = tif.pages[0].shape

        results = []

        for y in range(0, image_shape[0], tile_size):
            for x in range(0, image_shape[1], tile_size):
                # Define tile bounds
                y_end = min(y + tile_size, image_shape[0])
                x_end = min(x + tile_size, image_shape[1])

                # Load tile
                tile = tif.pages[0].asarray()[y:y_end, x:x_end]

                # Process tile
                result = process_tile(tile, x, y)
                results.append(result)

        return results


def process_tile(tile, offset_x, offset_y):
    """Process a single tile."""
    # Your processing here
    from skimage import filters, measure
    thresh = filters.threshold_otsu(tile)
    binary = tile > thresh
    labels = measure.label(binary)

    # Adjust coordinates by offset
    props = []
    for r in measure.regionprops(labels):
        props.append({
            'centroid_x': r.centroid[1] + offset_x,
            'centroid_y': r.centroid[0] + offset_y,
            'area': r.area
        })

    return props

---

Calibration

Convert Pixels to Physical Units

def pixels_to_microns(pixel_value, pixel_size_um):
    """Convert pixels to micrometers.

    Args:
        pixel_value: Value in pixels (area, length, etc.)
        pixel_size_um: Size of one pixel in micrometers

    Returns: Value in micrometers (or µm²)
    """
    return pixel_value * pixel_size_um


def calibrate_measurements(results_df, pixel_size_um, area_cols=None, length_cols=None):
    """Add calibrated measurements to DataFrame.

    Args:
        results_df: DataFrame with pixel measurements
        pixel_size_um: Pixel size in micrometers
        area_cols: List of area column names
        length_cols: List of length column names

    Returns: DataFrame with calibrated columns added
    """
    df = results_df.copy()

    if area_cols:
        for col in area_cols:
            df[f'{col}_um2'] = df[col] * (pixel_size_um ** 2)
            df[f'{col}_mm2'] = df[f'{col}_um2'] / 1e6

    if length_cols:
        for col in length_cols:
            df[f'{col}_um'] = df[col] * pixel_size_um

    return df


# Example
results = pd.DataFrame({
    'area': [100, 200, 150],
    'perimeter': [40, 60, 50]
})

calibrated = calibrate_measurements(
    results,
    pixel_size_um=0.65,
    area_cols=['area'],
    length_cols=['perimeter']
)
print(calibrated)

---

Quality Checks

Detect Focus Issues

from skimage.filters import laplace

def check_image_focus(image):
    """Check if image is in focus using Laplacian variance.

    Args:
        image: Grayscale image

    Returns: Focus score (higher = better focus)
    """
    laplacian = laplace(image)
    variance = laplacian.var()
    return variance


# Example: Check all images in folder
import glob

focus_scores = []
for img_path in glob.glob("images/*.tif"):
    img = tifffile.imread(img_path)
    score = check_image_focus(img)
    focus_scores.append({
        'filename': os.path.basename(img_path),
        'focus_score': score,
        'in_focus': score > 100  # Threshold depends on your images
    })

focus_df = pd.DataFrame(focus_scores)
print(focus_df)

Detect Saturation

def check_saturation(image, threshold=0.01):
    """Check if image has saturated pixels.

    Args:
        image: Image array
        threshold: Fraction of pixels allowed to be saturated

    Returns: (is_saturated, fraction_saturated)
    """
    if image.dtype == np.uint8:
        max_val = 255
    elif image.dtype == np.uint16:
        max_val = 65535
    else:
        max_val = image.max()

    saturated = image >= max_val
    fraction = saturated.sum() / image.size

    return fraction > threshold, fraction


# Check all channels
for i in range(image.shape[-1]):
    is_sat, frac = check_saturation(image[..., i])
    if is_sat:
        print(f"Channel {i}: {frac*100:.2f}% pixels saturated!")

---

Troubleshooting Common Issues

Image orientation wrong

# Rotate
from skimage.transform import rotate
rotated = rotate(image, angle=90)

# Flip
flipped_ud = np.flipud(image)  # Up-down
flipped_lr = np.fliplr(image)  # Left-right

Image too dark/bright

# Auto-adjust contrast
from skimage import exposure
adjusted = exposure.rescale_intensity(image, in_range='image', out_range=(0, 255))

Uneven illumination

# Background subtraction
from skimage.morphology import white_tophat, disk
background = white_tophat(image, disk(50))
corrected = image - background

File won't load

# Try different libraries
try:
    image = tifffile.imread(path)
except:
    from PIL import Image
    image = np.array(Image.open(path))

Colony and Object Segmentation

Complete guide for segmenting and measuring bacterial colonies, biofilms, and other large objects in microscopy images.

---

Table of Contents

1. Colony Morphometry Basics 2. Threshold Methods 3. Morphological Operations 4. Measurement and Analysis 5. Swarming Assay Analysis 6. Time-Lapse Analysis

---

Colony Morphometry Basics

Why Colony Morphometry?

Colony morphometry measures:

• Area: Colony size (pixels or calibrated units)
• Circularity: Shape roundness (4π × area / perimeter²)
• Roundness: Alternative shape metric (4 × area / π × major_axis²)
• Solidity: Convexity (area / convex_hull_area)
• Eccentricity: Elongation (0 = circle, 1 = line)

Used in:

• Bacterial swarming assays
• Biofilm formation studies
• Fungal colony growth
• Cell aggregate analysis

---

Threshold Methods

Otsu Threshold (Most Common)

from skimage import filters, measure, morphology
import numpy as np
import pandas as pd

def segment_colonies_otsu(image, min_area=100):
    """Segment colonies using Otsu's method.

    Best for: Well-defined colonies with clear contrast to background.

    Args:
        image: numpy array (grayscale or RGB)
        min_area: minimum colony area in pixels

    Returns: (binary_mask, labeled_image, properties_df)
    """
    # Convert to grayscale if needed
    if image.ndim == 3:
        from skimage.color import rgb2gray
        img = rgb2gray(image)
    else:
        img = image.astype(float)

    # Otsu threshold
    thresh = filters.threshold_otsu(img)
    binary = img > thresh  # or img < thresh if colonies are dark

    # Clean up small objects
    binary = morphology.remove_small_objects(binary, min_size=min_area)
    binary = morphology.binary_fill_holes(binary)

    # Label connected components
    labels = measure.label(binary)

    # Measure properties
    props = measure.regionprops_table(labels, properties=[
        'label', 'area', 'perimeter', 'eccentricity',
        'major_axis_length', 'minor_axis_length', 'solidity'
    ])
    props_df = pd.DataFrame(props)

    # Calculate circularity
    props_df['circularity'] = 4 * np.pi * props_df['area'] / (props_df['perimeter']**2)

    # Calculate roundness
    props_df['roundness'] = 4 * props_df['area'] / (np.pi * props_df['major_axis_length']**2)

    return binary, labels, props_df


# Example usage
import tifffile
image = tifffile.imread("swarming_plate.tif")
binary, labels, props = segment_colonies_otsu(image, min_area=500)
print(f"Found {labels.max()} colonies")
print(props.head())

Li Threshold (More Sensitive)

def segment_colonies_li(image, min_area=100):
    """Segment colonies using Li's method.

    Best for: Low-contrast images, faint colonies.

    Args:
        image: numpy array
        min_area: minimum colony area

    Returns: (binary_mask, labeled_image, properties_df)
    """
    if image.ndim == 3:
        from skimage.color import rgb2gray
        img = rgb2gray(image)
    else:
        img = image.astype(float)

    # Li threshold (more sensitive than Otsu)
    thresh = filters.threshold_li(img)
    binary = img > thresh

    # Cleanup
    binary = morphology.remove_small_objects(binary, min_size=min_area)
    binary = morphology.binary_fill_holes(binary)

    # Label and measure
    labels = measure.label(binary)
    props = measure.regionprops_table(labels, properties=[
        'label', 'area', 'perimeter', 'eccentricity',
        'major_axis_length', 'minor_axis_length', 'solidity'
    ])
    props_df = pd.DataFrame(props)

    # Add shape metrics
    props_df['circularity'] = 4 * np.pi * props_df['area'] / (props_df['perimeter']**2)
    props_df['roundness'] = 4 * props_df['area'] / (np.pi * props_df['major_axis_length']**2)

    return binary, labels, props_df

Adaptive Threshold (Uneven Illumination)

def segment_colonies_adaptive(image, block_size=51, offset=0.01, min_area=100):
    """Segment colonies with adaptive thresholding.

    Best for: Uneven illumination, gradient backgrounds.

    Args:
        image: numpy array
        block_size: size of local region (odd number)
        offset: threshold offset
        min_area: minimum colony area

    Returns: (binary_mask, labeled_image, properties_df)
    """
    if image.ndim == 3:
        from skimage.color import rgb2gray
        img = rgb2gray(image)
    else:
        img = image.astype(float)

    # Adaptive threshold
    local_thresh = filters.threshold_local(img, block_size=block_size, offset=offset)
    binary = img > local_thresh

    # Cleanup
    binary = morphology.remove_small_objects(binary, min_size=min_area)
    binary = morphology.binary_fill_holes(binary)

    # Label and measure
    labels = measure.label(binary)
    props = measure.regionprops_table(labels, properties=[
        'label', 'area', 'perimeter', 'eccentricity',
        'major_axis_length', 'minor_axis_length', 'solidity'
    ])
    props_df = pd.DataFrame(props)

    # Shape metrics
    props_df['circularity'] = 4 * np.pi * props_df['area'] / (props_df['perimeter']**2)
    props_df['roundness'] = 4 * props_df['area'] / (np.pi * props_df['major_axis_length']**2)

    return binary, labels, props_df

---

Morphological Operations

Fill Holes (Interior Gaps)

from skimage.morphology import binary_fill_holes

# Fill all holes
binary_filled = binary_fill_holes(binary)

# Fill only small holes
from skimage.morphology import remove_small_holes
binary_filled = remove_small_holes(binary, area_threshold=200)

Remove Small Objects

from skimage.morphology import remove_small_objects

# Remove objects smaller than threshold
binary_clean = remove_small_objects(binary, min_size=500)

Opening (Remove Protrusions)

from skimage.morphology import binary_opening, disk

# Remove small protrusions and thin connections
binary_opened = binary_opening(binary, disk(3))

Closing (Fill Gaps)

from skimage.morphology import binary_closing, disk

# Close small gaps between objects
binary_closed = binary_closing(binary, disk(5))

Complete Cleanup Pipeline

def morphological_cleanup(binary, min_area=500, opening_radius=3, closing_radius=5):
    """Apply morphological operations to clean up binary mask.

    Args:
        binary: binary mask
        min_area: minimum object area to keep
        opening_radius: radius for opening operation
        closing_radius: radius for closing operation

    Returns: Cleaned binary mask
    """
    # Remove small objects
    cleaned = morphology.remove_small_objects(binary, min_size=min_area)

    # Opening to remove protrusions
    if opening_radius > 0:
        cleaned = morphology.binary_opening(cleaned, morphology.disk(opening_radius))

    # Closing to fill gaps
    if closing_radius > 0:
        cleaned = morphology.binary_closing(cleaned, morphology.disk(closing_radius))

    # Fill holes
    cleaned = morphology.binary_fill_holes(cleaned)

    return cleaned

---

Measurement and Analysis

Complete Colony Measurement

def measure_colonies(image, threshold_method='otsu', min_area=500):
    """Segment and measure colonies from image.

    Returns: DataFrame with Area, Circularity, Roundness, Perimeter per colony
    """
    # Convert to grayscale
    if image.ndim == 3:
        from skimage.color import rgb2gray
        gray = rgb2gray(image)
    else:
        gray = image.astype(float)

    # Threshold
    if threshold_method == 'otsu':
        thresh = filters.threshold_otsu(gray)
    elif threshold_method == 'li':
        thresh = filters.threshold_li(gray)
    elif threshold_method == 'triangle':
        thresh = filters.threshold_triangle(gray)
    else:
        thresh = threshold_method  # numeric value

    binary = gray > thresh

    # Clean up
    binary = morphology.remove_small_objects(binary, min_size=min_area)
    binary = morphology.binary_fill_holes(binary)

    # Label and measure
    labels = measure.label(binary)
    props = measure.regionprops_table(labels, properties=[
        'area', 'perimeter', 'eccentricity', 'solidity',
        'major_axis_length', 'minor_axis_length'
    ])

    results = pd.DataFrame(props)

    # Calculate shape metrics
    results['circularity'] = 4 * np.pi * results['area'] / (results['perimeter']**2)
    results['roundness'] = 4 * results['area'] / (np.pi * results['major_axis_length']**2)

    # Rename to match ImageJ/CellProfiler conventions
    results = results.rename(columns={
        'area': 'Area',
        'perimeter': 'Perimeter',
        'circularity': 'Circularity',
        'roundness': 'Round',
        'eccentricity': 'Eccentricity',
        'solidity': 'Solidity'
    })

    return results

Calibrated Measurements

def measure_colonies_calibrated(image, pixel_size_um, threshold_method='otsu', min_area_um2=0.1):
    """Measure colonies with calibrated units.

    Args:
        image: numpy array
        pixel_size_um: size of one pixel in micrometers (e.g., 0.65 for 0.65 µm/pixel)
        threshold_method: 'otsu', 'li', or numeric
        min_area_um2: minimum area in mm² (1 mm² = 1,000,000 µm²)

    Returns: DataFrame with calibrated measurements
    """
    # Convert min_area to pixels
    min_area_pixels = int(min_area_um2 * 1e6 / (pixel_size_um**2))

    # Segment and measure in pixels
    results = measure_colonies(image, threshold_method, min_area_pixels)

    # Convert to calibrated units
    results['Area_um2'] = results['Area'] * (pixel_size_um**2)
    results['Area_mm2'] = results['Area_um2'] / 1e6
    results['Perimeter_um'] = results['Perimeter'] * pixel_size_um

    # Circularity and roundness are dimensionless (no conversion needed)

    return results

---

Swarming Assay Analysis

Complete Swarming Assay Workflow

def analyze_swarming_plate(image_path, genotype, replicate,
                          pixel_size_um=1.0, min_area_mm2=0.5):
    """Complete analysis of bacterial swarming plate.

    Args:
        image_path: path to image file
        genotype: genotype label
        replicate: replicate number
        pixel_size_um: pixel size in micrometers
        min_area_mm2: minimum colony area in mm²

    Returns: DataFrame with measurements for this plate
    """
    import tifffile
    import os

    # Load image
    image = tifffile.imread(image_path)

    # Measure colonies
    results = measure_colonies_calibrated(
        image,
        pixel_size_um=pixel_size_um,
        threshold_method='otsu',
        min_area_um2=min_area_mm2 * 1e6  # Convert mm² to µm²
    )

    # Add metadata
    results['Genotype'] = genotype
    results['Replicate'] = replicate
    results['Image'] = os.path.basename(image_path)

    return results


# Example: Batch process multiple plates
def batch_analyze_swarming(image_info_df, pixel_size_um=1.0):
    """Analyze multiple swarming plates.

    Args:
        image_info_df: DataFrame with columns: 'image_path', 'genotype', 'replicate'
        pixel_size_um: pixel size

    Returns: Combined DataFrame with all measurements
    """
    all_results = []

    for idx, row in image_info_df.iterrows():
        results = analyze_swarming_plate(
            row['image_path'],
            row['genotype'],
            row['replicate'],
            pixel_size_um=pixel_size_um
        )
        all_results.append(results)

    combined = pd.concat(all_results, ignore_index=True)
    return combined

Colony Morphometry Statistics

def colony_morphometry_analysis(df, genotype_col='Genotype',
                               area_col='Area', circ_col='Circularity'):
    """Full colony morphometry analysis for swarming assays.

    Args:
        df: DataFrame with colony measurements
        genotype_col: column with genotype labels
        area_col: column with area measurements
        circ_col: column with circularity measurements

    Returns: dict with per-genotype summaries and max-area genotype info
    """
    # Group statistics for area
    area_summary = df.groupby(genotype_col)[area_col].agg(
        Mean='mean',
        SD='std',
        Median='median',
        Min='min',
        Max='max',
        N='count'
    ).reset_index()
    area_summary['SEM'] = area_summary['SD'] / np.sqrt(area_summary['N'])

    # Group statistics for circularity
    circ_summary = df.groupby(genotype_col)[circ_col].agg(
        Mean='mean',
        SD='std',
        Median='median',
        N='count'
    ).reset_index()
    circ_summary['SEM'] = circ_summary['SD'] / np.sqrt(circ_summary['N'])

    # Merge area and circularity summaries
    merged = area_summary.merge(
        circ_summary, on=genotype_col, suffixes=('_Area', '_Circ')
    )

    # Find genotype with largest mean area
    max_area_idx = merged['Mean_Area'].idxmax()
    max_area_genotype = merged.loc[max_area_idx, genotype_col]
    max_area_circularity = merged.loc[max_area_idx, 'Mean_Circ']

    return {
        'area_summary': area_summary,
        'circ_summary': circ_summary,
        'merged_summary': merged,
        'max_area_genotype': max_area_genotype,
        'max_area_circularity': max_area_circularity,
    }


# Example usage (colony morphometry)
df = pd.read_csv("Swarm_1.csv")
result = colony_morphometry_analysis(df, 'Genotype', 'Area', 'Circularity')
print(f"Genotype with largest area: {result['max_area_genotype']}")
print(f"Mean circularity: {result['max_area_circularity']:.4f}")

---

Time-Lapse Analysis

Track Colony Growth Over Time

def analyze_timelapse(image_paths, timepoints_hours):
    """Analyze colony growth from time-lapse images.

    Args:
        image_paths: list of image file paths (sorted by time)
        timepoints_hours: list of timepoint values in hours

    Returns: DataFrame with area vs time for each colony
    """
    import tifffile

    growth_data = []

    for i, (img_path, time_hr) in enumerate(zip(image_paths, timepoints_hours)):
        # Load and measure
        img = tifffile.imread(img_path)
        results = measure_colonies(img, threshold_method='otsu', min_area=100)

        # Add time information
        results['timepoint'] = i
        results['time_hours'] = time_hr

        growth_data.append(results)

    # Combine all timepoints
    combined = pd.concat(growth_data, ignore_index=True)

    return combined


# Example: Plot growth curves
import matplotlib.pyplot as plt

def plot_growth_curves(growth_df, colony_labels=None):
    """Plot colony area vs time.

    Args:
        growth_df: DataFrame from analyze_timelapse
        colony_labels: Optional list of colony IDs to plot
    """
    if colony_labels is None:
        # Plot mean area over time
        summary = growth_df.groupby('time_hours')['Area'].agg(['mean', 'std'])
        plt.errorbar(summary.index, summary['mean'], yerr=summary['std'],
                    marker='o', capsize=5)
        plt.xlabel('Time (hours)')
        plt.ylabel('Mean Colony Area (pixels)')
        plt.title('Colony Growth Over Time')
    else:
        # Plot individual colonies
        for colony_id in colony_labels:
            colony_data = growth_df[growth_df['label'] == colony_id]
            plt.plot(colony_data['time_hours'], colony_data['Area'],
                    marker='o', label=f'Colony {colony_id}')
        plt.xlabel('Time (hours)')
        plt.ylabel('Colony Area (pixels)')
        plt.title('Individual Colony Growth')
        plt.legend()

    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.show()

---

Troubleshooting

Common Issues

Multiple colonies merged into one:

• Use watershed segmentation (see cell_counting.md)
• Increase erosion before labeling
• Try more conservative threshold (Li instead of Otsu)

Background detected as colony:

• Increase min_area threshold
• Apply morphological opening
• Use adaptive threshold for uneven backgrounds

Colony edges jagged:

• Apply Gaussian smoothing before thresholding
• Use morphological closing
• Increase binary dilation

Small colonies missed:

• Use Li or Triangle threshold (more sensitive)
• Reduce min_area threshold
• Enhance contrast before segmentation

Uneven illumination affects segmentation:

• Use adaptive thresholding
• Apply background subtraction first
• Use rolling ball background correction

---

Threshold Method Selection Guide

---

Parameter Recommendations

Bacterial Swarming Plates

# Typical parameters for Petri dish images (100mm plates)
params = {
    'threshold_method': 'otsu',
    'min_area': 500,  # pixels (adjust based on magnification)
    'opening_radius': 3,
    'closing_radius': 5,
}

Biofilm Analysis

# Biofilms often have irregular edges
params = {
    'threshold_method': 'li',  # More sensitive
    'min_area': 1000,
    'opening_radius': 5,  # Stronger cleanup
    'closing_radius': 8,
}

Fungal Colonies

# Fungi can have very irregular shapes
params = {
    'threshold_method': 'adaptive',
    'block_size': 51,
    'min_area': 2000,
    'opening_radius': 0,  # Preserve edges
    'closing_radius': 10,
}

Statistical Analysis for Imaging Data

Complete reference for statistical tests, effect sizes, power analysis, and regression modeling on microscopy measurement data.

---

Table of Contents

1. Descriptive Statistics 2. Normality Testing 3. Two-Group Comparisons 4. Multiple Comparisons 5. Two-Way ANOVA 6. Effect Sizes 7. Power Analysis 8. Regression Modeling 9. Model Comparison

---

Descriptive Statistics

Grouped Summary Statistics

import pandas as pd
import numpy as np

def grouped_summary(df, group_cols, measure_col):
    """Calculate summary statistics by group.

    Returns DataFrame with Mean, SD, SEM, Median, Min, Max, N per group.
    """
    if isinstance(group_cols, str):
        group_cols = [group_cols]

    summary = df.groupby(group_cols)[measure_col].agg(
        Mean='mean',
        SD='std',
        Median='median',
        Min='min',
        Max='max',
        N='count'
    ).reset_index()

    summary['SEM'] = summary['SD'] / np.sqrt(summary['N'])

    return summary

Percent Reduction

def percent_reduction(df, group_col, measure_col, reference_group, comparison_group):
    """Calculate percent reduction of comparison vs reference.

    Returns: (percent_reduction, ref_mean, comp_mean)
    """
    ref_mean = df[df[group_col] == reference_group][measure_col].mean()
    comp_mean = df[df[group_col] == comparison_group][measure_col].mean()
    pct_reduction = ((ref_mean - comp_mean) / ref_mean) * 100
    return pct_reduction, ref_mean, comp_mean

Relative Proportion

def relative_proportion(df, group_col, measure_col, numerator_group, denominator_group):
    """Calculate relative proportion (as percentage) of one group vs another.

    Returns: proportion as percentage
    """
    num_mean = df[df[group_col] == numerator_group][measure_col].mean()
    den_mean = df[df[group_col] == denominator_group][measure_col].mean()
    return (num_mean / den_mean) * 100

---

Normality Testing

Shapiro-Wilk Test

from scipy import stats

def shapiro_wilk_test(data):
    """Perform Shapiro-Wilk test for normality.

    Returns: (W_statistic, p_value)

    Interpretation:
    - p < 0.05: Data is NOT normally distributed
    - p >= 0.05: Data is normally distributed
    """
    stat, pvalue = stats.shapiro(data)
    return stat, pvalue

# Example usage
data = df[df['Condition'] == 'Control']['Measurement']
w_stat, p_val = shapiro_wilk_test(data)
print(f"Shapiro-Wilk W={w_stat:.4f}, p={p_val:.4f}")

---

Two-Group Comparisons

Independent T-Test

def independent_ttest(group1, group2, equal_var=True):
    """Perform independent two-sample t-test.

    Args:
        group1, group2: array-like data
        equal_var: If True, use standard t-test. If False, use Welch's t-test

    Returns: (t_statistic, p_value)
    """
    t_stat, p_val = stats.ttest_ind(group1, group2, equal_var=equal_var)
    return t_stat, p_val

Mann-Whitney U Test (non-parametric)

def mann_whitney_test(group1, group2):
    """Perform Mann-Whitney U test (non-parametric alternative to t-test).

    Use when data is NOT normally distributed.

    Returns: (U_statistic, p_value)
    """
    u_stat, p_val = stats.mannwhitneyu(group1, group2, alternative='two-sided')
    return u_stat, p_val

---

Multiple Comparisons

Dunnett's Test

CRITICAL: Dunnett's test compares each treatment group to a single control group. R uses multcomp::glht() with mcp(Strain_Ratio = "Dunnett"). Python equivalent uses scipy.

from scipy import stats

def dunnetts_test_scipy(df, group_col, value_col, control_group, alpha=0.05):
    """Dunnett's test using scipy.stats.dunnett (scipy >= 1.10).

    This is the preferred method - uses exact Dunnett distribution.

    Args:
        df: DataFrame
        group_col: Column with group labels
        value_col: Column with measurements
        control_group: Label of the control group
        alpha: Significance level (default 0.05)

    Returns: DataFrame with group, p_value, statistic, significant
    """
    groups = sorted(df[group_col].unique())
    control_data = df[df[group_col] == control_group][value_col].values
    treatment_groups = [g for g in groups if g != control_group]

    treatment_data = [df[df[group_col] == g][value_col].values for g in treatment_groups]

    # scipy.stats.dunnett: compare multiple treatment groups against control
    result = stats.dunnett(*treatment_data, control=control_data, alternative='two-sided')

    results = []
    for i, tg in enumerate(treatment_groups):
        results.append({
            'group': tg,
            'control': control_group,
            'p_value': result.pvalue[i],
            'statistic': result.statistic[i],
            'significant': result.pvalue[i] < alpha
        })

    return pd.DataFrame(results)

Combined Dunnett's Test (Two Measures)

def dunnett_area_circularity(df, group_col, area_col, circ_col, control_group, alpha=0.05):
    """Run Dunnett's test on both area and circularity.

    Returns: dict with:
        - area_results: Dunnett results for area
        - circ_results: Dunnett results for circularity
        - merged: Combined results
        - equivalent_in_both: groups NOT significant in EITHER (equivalent to control)
        - different_in_both: groups significant in BOTH (different from control)
    """
    area_dunnett = dunnetts_test_scipy(df, group_col, area_col, control_group, alpha)
    circ_dunnett = dunnetts_test_scipy(df, group_col, circ_col, control_group, alpha)

    # Merge results
    merged = area_dunnett[['group', 'p_value', 'significant']].merge(
        circ_dunnett[['group', 'p_value', 'significant']],
        on='group', suffixes=('_area', '_circ')
    )

    # Groups NOT significant in EITHER (equivalent to control in both)
    both_equiv = merged[~merged['significant_area'] & ~merged['significant_circ']]['group'].tolist()

    # Groups significant in both (different from control in both)
    both_diff = merged[merged['significant_area'] & merged['significant_circ']]['group'].tolist()

    return {
        'area_results': area_dunnett,
        'circ_results': circ_dunnett,
        'merged': merged,
        'equivalent_in_both': both_equiv,
        'different_in_both': both_diff,
    }

Tukey HSD (all pairwise comparisons)

from statsmodels.stats.multicomp import pairwise_tukeyhsd

def tukey_hsd_test(df, value_col, group_col, alpha=0.05):
    """Perform Tukey HSD test for all pairwise comparisons.

    Use when you want to compare ALL groups to EACH OTHER (not just vs control).

    Returns: Tukey HSD result object with summary table
    """
    result = pairwise_tukeyhsd(
        endog=df[value_col],
        groups=df[group_col],
        alpha=alpha
    )
    return result

---

Two-Way ANOVA

import statsmodels.api as sm
from statsmodels.formula.api import ols

def two_way_anova(df, dependent_var, factor1, factor2, anova_type=2):
    """Perform two-way ANOVA with interaction term.

    Args:
        df: DataFrame
        dependent_var: Column name of dependent variable
        factor1: First factor column name
        factor2: Second factor column name
        anova_type: Type of sum of squares (1, 2, or 3)
            - Type 1: Sequential (order matters)
            - Type 2: Hierarchical (recommended for balanced designs)
            - Type 3: Marginal (recommended for unbalanced designs)

    Returns: ANOVA table as DataFrame with columns:
        sum_sq, df, F, PR(>F) for each term including interaction
    """
    formula = f'{dependent_var} ~ C({factor1}) * C({factor2})'
    model = ols(formula, data=df).fit()
    anova_table = sm.stats.anova_lm(model, typ=anova_type)
    return anova_table


def extract_anova_interaction(anova_table, factor1, factor2):
    """Extract interaction F-statistic and p-value from ANOVA table.

    Returns: (F_statistic, p_value)
    """
    interaction_key = f'C({factor1}):C({factor2})'
    f_stat = anova_table.loc[interaction_key, 'F']
    p_val = anova_table.loc[interaction_key, 'PR(>F)']
    return f_stat, p_val

# Example usage
anova_result = two_way_anova(df, 'NeuN_count', 'Condition', 'Sex')
print(anova_result)

# Extract interaction
f_stat, p_val = extract_anova_interaction(anova_result, 'Condition', 'Sex')
print(f"Interaction F={f_stat:.3f}, p={p_val:.4f}")

---

Effect Sizes

Cohen's d

def cohens_d(group1, group2):
    """Calculate Cohen's d using pooled standard deviation.

    This matches the standard formula: d = (mean1 - mean2) / sd_pooled
    where sd_pooled = sqrt(((n1-1)*s1^2 + (n2-1)*s2^2) / (n1+n2-2))

    NOTE: Uses pandas .std() which defaults to ddof=1 (sample std).

    Interpretation:
    - |d| < 0.2: Small effect
    - |d| = 0.2-0.5: Small to medium effect
    - |d| = 0.5-0.8: Medium to large effect
    - |d| > 0.8: Large effect

    Returns: Cohen's d value
    """
    n1, n2 = len(group1), len(group2)
    s1, s2 = group1.std(), group2.std()  # ddof=1 by default in pandas

    sd_pooled = np.sqrt(((n1 - 1) * s1**2 + (n2 - 1) * s2**2) / (n1 + n2 - 2))
    d = (group1.mean() - group2.mean()) / sd_pooled

    return d

# Example usage
control = df[df['Condition'] == 'Control']['Measurement']
treatment = df[df['Condition'] == 'Treatment']['Measurement']
effect_size = cohens_d(control, treatment)
print(f"Cohen's d = {effect_size:.3f}")

---

Power Analysis

Sample Size Calculation

from statsmodels.stats.power import TTestIndPower

def power_analysis_sample_size(effect_size, alpha=0.05, power=0.8, alternative='two-sided'):
    """Calculate required sample size per group for a two-sample t-test.

    Args:
        effect_size: Cohen's d (can be negative, absolute value used)
        alpha: Significance level (Type I error rate)
        power: Desired statistical power (1 - Type II error rate)
        alternative: 'two-sided', 'larger', or 'smaller'

    Returns: Required sample size per group (rounded up to integer)

    Example:
        # To detect effect size of 0.8 with 80% power
        n = power_analysis_sample_size(0.8, alpha=0.05, power=0.8)
        print(f"Need {n} samples per group")
    """
    analysis = TTestIndPower()
    # Use absolute value of effect size for sample size calculation
    n = analysis.solve_power(
        effect_size=abs(effect_size),
        alpha=alpha,
        power=power,
        alternative=alternative
    )
    return int(np.ceil(n))

Post-hoc Power Calculation

def calculate_achieved_power(n_per_group, effect_size, alpha=0.05, alternative='two-sided'):
    """Calculate achieved power given sample size and effect size.

    Args:
        n_per_group: Sample size per group
        effect_size: Cohen's d
        alpha: Significance level
        alternative: 'two-sided', 'larger', or 'smaller'

    Returns: Achieved power (0-1)
    """
    analysis = TTestIndPower()
    power = analysis.solve_power(
        effect_size=abs(effect_size),
        nobs1=n_per_group,
        alpha=alpha,
        alternative=alternative
    )
    return power

---

Regression Modeling

Data Preparation for Co-culture Ratios

def prepare_ratio_data(df, strain_col='StrainNumber', ratio_col='Ratio',
                       area_col='Area', exclude_strains=None):
    """Prepare co-culture ratio data for regression analysis.

    Converts ratio strings (e.g., "3:1") to frequency fractions.
    Filters out pure strains if requested.

    Args:
        df: DataFrame with swarming data
        strain_col: Column with strain identifiers
        ratio_col: Column with ratio strings
        area_col: Column with area measurements
        exclude_strains: List of strain IDs to exclude (e.g., pure strains)

    Returns: DataFrame with Frequency_rhlI column added
    """
    result = df.copy()

    if exclude_strains:
        result = result[~result[strain_col].isin(exclude_strains)]

    # Parse ratio into frequency
    # Ratio format: "rhlI_D:lasI_D" (e.g., "3:1" means 3 parts rhlI, 1 part lasI)
    ratio_parts = result[ratio_col].str.split(':', expand=True).astype(int)
    result['rhlI_D'] = ratio_parts[0]
    result['lasI_D'] = ratio_parts[1]
    result['Frequency_rhlI'] = result['rhlI_D'] / (result['rhlI_D'] + result['lasI_D'])

    return result

Polynomial Regression

def fit_polynomial_model(df, x_col, y_col, degree=2):
    """Fit polynomial regression model.

    Equivalent to R: lm(y ~ poly(x, degree, raw=TRUE))

    Args:
        df: DataFrame
        x_col: Predictor column name
        y_col: Response column name
        degree: Polynomial degree (2=quadratic, 3=cubic)

    Returns: dict with model, coefficients, R-squared, F-statistic, p-value,
             peak_frequency, peak_value, peak_ci
    """
    x = df[x_col].values
    y = df[y_col].values

    # Build design matrix for polynomial
    X_poly = np.column_stack([x**i for i in range(1, degree+1)])
    X = sm.add_constant(X_poly)

    model = sm.OLS(y, X).fit()

    result = {
        'model': model,
        'coefficients': model.params,
        'r_squared': model.rsquared,
        'adj_r_squared': model.rsquared_adj,
        'f_statistic': model.fvalue,
        'f_pvalue': model.f_pvalue,
        'aic': model.aic,
        'bic': model.bic,
        'summary': model.summary(),
    }

    # Find peak (maximum) using calculus
    if degree == 2:
        # y = b0 + b1*x + b2*x^2
        # Peak at x = -b1 / (2*b2)
        b1, b2 = model.params[1], model.params[2]
        peak_x = -b1 / (2 * b2)
    elif degree == 3:
        # y = b0 + b1*x + b2*x^2 + b3*x^3
        # Derivative: b1 + 2*b2*x + 3*b3*x^2 = 0
        b1, b2, b3 = model.params[1], model.params[2], model.params[3]
        discriminant = (2*b2)**2 - 4*(3*b3)*b1
        if discriminant >= 0:
            x1 = (-2*b2 + np.sqrt(discriminant)) / (2*3*b3)
            x2 = (-2*b2 - np.sqrt(discriminant)) / (2*3*b3)
            # Choose the one that's a maximum (second derivative < 0)
            candidates = [x1, x2]
            peak_x = None
            for cx in candidates:
                second_deriv = 2*b2 + 6*b3*cx
                if second_deriv < 0:
                    peak_x = cx
                    break
            if peak_x is None:
                peak_x = candidates[0]  # fallback
        else:
            peak_x = x.mean()  # no real critical points
    else:
        # For higher degrees, use numerical optimization
        from scipy.optimize import minimize_scalar
        poly_func = lambda xx: -sum(model.params[i] * xx**i for i in range(degree+1))
        opt = minimize_scalar(poly_func, bounds=(x.min(), x.max()), method='bounded')
        peak_x = opt.x

    # Predict at peak and get confidence interval
    X_peak = np.array([[1] + [peak_x**i for i in range(1, degree+1)]])
    peak_pred = model.get_prediction(X_peak)
    peak_value = peak_pred.predicted_mean[0]
    peak_ci = peak_pred.conf_int(alpha=0.05)[0]

    result['peak_x'] = peak_x
    result['peak_value'] = peak_value
    result['peak_ci_lower'] = peak_ci[0]
    result['peak_ci_upper'] = peak_ci[1]

    return result

Natural Spline Regression

CRITICAL: This must match R's lm(Area ~ ns(Frequency_rhlI, df=4)).

from patsy import dmatrix

def fit_natural_spline_model(df, x_col, y_col, spline_df=4):
    """Fit natural spline regression model.

    Equivalent to R: lm(y ~ ns(x, df=spline_df))
    Uses patsy's cr() with explicit quantile knots to match R's ns().

    CRITICAL: R's ns(x, df=N) places N-1 internal knots at equally-spaced
    quantiles (25th, 50th, 75th for df=4). patsy's cr(df=N) does NOT place
    knots at the same locations by default. You MUST specify knots explicitly.

    Args:
        df: DataFrame
        x_col: Predictor column name
        y_col: Response column name
        spline_df: Degrees of freedom for natural spline basis

    Returns: dict with model, R-squared, F-statistic, p-value,
             peak_frequency, peak_value, peak_ci
    """
    x = df[x_col].values
    y = df[y_col].values

    # Match R's ns() knot placement: df-1 internal knots at equally-spaced quantiles
    n_internal_knots = spline_df - 1
    quantile_pcts = np.linspace(100.0 / (n_internal_knots + 1),
                                100.0 * n_internal_knots / (n_internal_knots + 1),
                                n_internal_knots)
    knots = np.percentile(x, quantile_pcts)
    knot_str = ", ".join([str(k) for k in knots])

    # Create natural spline basis using patsy's cr() with explicit knots
    formula_str = f"cr({x_col}, knots=[{knot_str}]) - 1"
    X_spline = np.array(dmatrix(formula_str, df))
    X = sm.add_constant(X_spline)

    model = sm.OLS(y, X).fit()

    result = {
        'model': model,
        'r_squared': model.rsquared,
        'adj_r_squared': model.rsquared_adj,
        'f_statistic': model.fvalue,
        'f_pvalue': model.f_pvalue,
        'aic': model.aic,
        'bic': model.bic,
    }

    # Find peak by evaluating over fine grid
    x_grid = np.linspace(x.min(), x.max(), 1000)
    grid_df = pd.DataFrame({x_col: x_grid})
    X_grid_spline = np.array(dmatrix(formula_str, grid_df))
    X_grid = sm.add_constant(X_grid_spline)

    predictions = model.get_prediction(X_grid)
    pred_mean = predictions.predicted_mean
    pred_ci = predictions.conf_int(alpha=0.05)

    # Find peak
    max_idx = np.argmax(pred_mean)

    result['peak_x'] = x_grid[max_idx]
    result['peak_value'] = pred_mean[max_idx]
    result['peak_ci_lower'] = pred_ci[max_idx, 0]
    result['peak_ci_upper'] = pred_ci[max_idx, 1]
    result['predictions'] = pred_mean
    result['x_grid'] = x_grid
    result['ci_lower'] = pred_ci[:, 0]
    result['ci_upper'] = pred_ci[:, 1]

    return result

---

Model Comparison

def compare_regression_models(models_dict):
    """Compare multiple regression models.

    Args:
        models_dict: Dict mapping model_name -> result dict from fit_* functions

    Returns: DataFrame with model comparison metrics

    Example:
        models = {
            'quadratic': fit_polynomial_model(df, 'x', 'y', degree=2),
            'cubic': fit_polynomial_model(df, 'x', 'y', degree=3),
            'spline': fit_natural_spline_model(df, 'x', 'y', spline_df=4)
        }
        comparison = compare_regression_models(models)
        print(comparison)
    """
    comparison = []
    for name, result in models_dict.items():
        comparison.append({
            'model': name,
            'r_squared': result['r_squared'],
            'adj_r_squared': result['adj_r_squared'],
            'f_statistic': result['f_statistic'],
            'f_pvalue': result['f_pvalue'],
            'aic': result.get('aic'),
            'bic': result.get('bic'),
            'peak_x': result.get('peak_x'),
            'peak_value': result.get('peak_value'),
            'peak_ci_lower': result.get('peak_ci_lower'),
            'peak_ci_upper': result.get('peak_ci_upper'),
        })

    comparison_df = pd.DataFrame(comparison)
    comparison_df['best_r2'] = comparison_df['r_squared'] == comparison_df['r_squared'].max()
    comparison_df['best_aic'] = comparison_df['aic'] == comparison_df['aic'].min()
    comparison_df['best_bic'] = comparison_df['bic'] == comparison_df['bic'].min()

    return comparison_df.sort_values('r_squared', ascending=False)

---

Answer Extraction Patterns

Answer formatting

def format_answer(value, question_type):
    """Format answer for reporting.

    Args:
        value: Numeric value to format
        question_type: Type of answer expected

    Returns: Formatted value
    """
    if question_type == "nearest_thousand":
        return int(round(value, -3))
    elif question_type == "percentage_int":
        return int(round(value))
    elif question_type == "percentage_2dec":
        return round(value, 2)
    elif question_type == "cohen_d":
        return round(value, 3)
    elif question_type == "statistic_3dec":
        return round(value, 3)
    elif question_type == "sample_size":
        return int(np.ceil(value))
    elif question_type == "count":
        return int(value)
    elif question_type == "r_squared":
        return round(value, 2)
    elif question_type == "p_value":
        if value < 0.0001:
            return f"{value:.2e}"
        else:
            return round(value, 4)
    else:
        return value

---

Complete Example: Spline Regression Workflow

import pandas as pd
import numpy as np

# Load data
df = pd.read_csv("Swarm_2.csv")

# Prepare ratio data
df_coculture = prepare_ratio_data(
    df,
    strain_col='StrainNumber',
    ratio_col='Ratio',
    area_col='Area',
    exclude_strains=['1', '98']  # Exclude pure strains
)

# Fit models
quadratic = fit_polynomial_model(df_coculture, 'Frequency_rhlI', 'Area', degree=2)
cubic = fit_polynomial_model(df_coculture, 'Frequency_rhlI', 'Area', degree=3)
spline = fit_natural_spline_model(df_coculture, 'Frequency_rhlI', 'Area', spline_df=4)

# Compare models
models = {'quadratic': quadratic, 'cubic': cubic, 'spline': spline}
comparison = compare_regression_models(models)
print(comparison)

# Best model by R-squared
best_model_name = comparison.iloc[0]['model']
best_model = models[best_model_name]

# Extract answers
print(f"Best model: {best_model_name}")
print(f"R-squared: {best_model['r_squared']:.4f}")
print(f"Peak frequency: {best_model['peak_x']:.4f}")
print(f"Peak area: {best_model['peak_value']:.1f}")
print(f"95% CI: [{best_model['peak_ci_lower']:.1f}, {best_model['peak_ci_upper']:.1f}]")

Troubleshooting Guide

Common issues and solutions for microscopy image analysis.

---

Segmentation Issues

Problem: Cells/Colonies Merged Together

Symptoms: Multiple objects detected as single object

Solutions: 1. Use watershed segmentation

   from scipy import ndimage
   from skimage.feature import peak_local_max
   from skimage.segmentation import watershed

   distance = ndimage.distance_transform_edt(binary)
   coords = peak_local_max(distance, min_distance=10)
   mask = np.zeros(distance.shape, dtype=bool)
   mask[tuple(coords.T)] = True
   markers = measure.label(mask)
   labels = watershed(-distance, markers, mask=binary)

2. Decrease min_distance parameter in watershed 3. Apply morphological opening before labeling 4. For very dense images, use CellPose or StarDist

---

Problem: Cells Split into Multiple Objects

Symptoms: One cell detected as 2-3 objects

Solutions: 1. Increase min_distance in watershed 2. Apply morphological closing before labeling

   from skimage.morphology import binary_closing, disk
   binary = binary_closing(binary, disk(5))

3. Use more conservative threshold (Li instead of Otsu) 4. Smooth image before thresholding

---

Problem: Background Detected as Objects

Symptoms: Many small false positive detections

Solutions: 1. Increase min_area threshold

   binary = morphology.remove_small_objects(binary, min_size=200)

2. Use more conservative threshold method 3. Apply morphological opening 4. Filter by shape metrics (circularity, eccentricity)

---

Problem: Dim Objects Not Detected

Symptoms: Faint cells/colonies missing

Solutions: 1. Use Li or Triangle threshold (more sensitive)

   thresh = filters.threshold_li(image)

2. Use percentile-based threshold

   thresh = np.percentile(image[image > 0], 25)

3. Enhance contrast before thresholding

   from skimage import exposure
   enhanced = exposure.equalize_adapthist(image)

4. Manually set lower threshold value

---

Image Quality Issues

Problem: Uneven Illumination

Symptoms: Gradient across image, threshold doesn't work globally

Solutions: 1. Use adaptive thresholding

   from skimage.filters import threshold_local
   local_thresh = threshold_local(image, block_size=51)
   binary = image > local_thresh

2. Apply background subtraction

   from skimage.morphology import white_tophat, disk
   background = white_tophat(image, disk(50))
   corrected = image - background

3. Use rolling ball background correction 4. Normalize illumination before analysis

---

Problem: Noisy Images

Symptoms: Salt-and-pepper noise, grainy appearance

Solutions: 1. Apply Gaussian smoothing

   from skimage.filters import gaussian
   smoothed = gaussian(image, sigma=2.0)

2. Use median filter (better for salt-and-pepper)

   from scipy.ndimage import median_filter
   smoothed = median_filter(image, size=3)

3. Use bilateral filter (edge-preserving)

   from skimage.restoration import denoise_bilateral
   smoothed = denoise_bilateral(image, sigma_color=0.05, sigma_spatial=15)

---

Problem: Out-of-Focus Images

Symptoms: Blurry, low contrast, poor segmentation

Solutions: 1. Apply unsharp masking

   from skimage.filters import unsharp_mask
   sharpened = unsharp_mask(image, radius=2, amount=1.0)

2. Enhance edges

   from skimage.filters import sobel
   edges = sobel(image)
   enhanced = image + 0.5 * edges

3. Use deconvolution (if PSF known) 4. Reject out-of-focus images (check Laplacian variance)

---

Statistical Analysis Issues

Problem: Results Don't Match R

Symptoms: P-values or statistics differ from R output

Solutions for Dunnett's test:

• Use scipy.stats.dunnett() (requires scipy >= 1.10)
• Ensure group labels match R exactly
• Check control group is correctly specified

Solutions for natural spline:

• Use explicit knot placement matching R's ns()

  n_internal_knots = spline_df - 1
  quantile_pcts = np.linspace(100.0/(n_internal_knots+1),
                               100.0*n_internal_knots/(n_internal_knots+1),
                               n_internal_knots)
  knots = np.percentile(x, quantile_pcts)

Solutions for Cohen's d:

• Verify using pooled standard deviation
• Check pandas .std() uses ddof=1 (sample SD)
• Ensure correct group order (sign matters)

---

Problem: High P-Values (No Significance)

Symptoms: P-values all > 0.05

Possible causes: 1. Insufficient sample size - Run power analysis 2. High variability - Check SD within groups 3. Small effect size - Calculate Cohen's d 4. Wrong statistical test - Check assumptions (normality, equal variance) 5. Data quality issues - Check for outliers, measurement errors

Solutions: 1. Increase sample size (if collecting new data) 2. Use non-parametric tests if data not normal 3. Check for and handle outliers 4. Consider effect size rather than just p-value

---

Problem: Can't Reproduce Expected Answer

Symptoms: Your answer differs from the expected answer

Check these: 1. Data loading - Correct file? All rows loaded? 2. Grouping variables - Exact spelling of group names? 3. Filter conditions - Any rows excluded/included incorrectly? 4. Rounding - "Nearest thousand" means int(round(val, -3)) 5. Direction - "Reduction" vs "increase", "KD vs CTRL" vs "CTRL vs KD" 6. Statistical function - Exact match to R function used?

---

Performance Issues

Problem: Analysis Too Slow

Solutions for large images: 1. Process tiles instead of whole image 2. Downsample for preview/parameter tuning 3. Use OpenCV for batch processing (faster than scikit-image) 4. Parallelize batch processing

   from multiprocessing import Pool
   with Pool(4) as p:
       results = p.map(process_func, image_paths)

---

Problem: Out of Memory

Solutions: 1. Process images in tiles 2. Load images as memory-mapped arrays 3. Reduce image resolution for analysis 4. Process one at a time instead of loading all 5. Use generators instead of lists

---

Data Format Issues

Problem: Image Won't Load

Symptoms: Error loading file

Solutions: 1. Try different libraries

   # Try tifffile first
   try:
       image = tifffile.imread(path)
   except:
       # Try PIL
       from PIL import Image
       image = np.array(Image.open(path))

2. Check file format (use file command on Linux/Mac) 3. Check file corruption 4. Verify file permissions

---

Problem: Wrong Image Dimensions

Symptoms: Image shape not as expected

Solutions: 1. Check dimension order (C, Z, T)

   print(f"Shape: {image.shape}")
   # Might need to transpose
   image = np.transpose(image, (1, 2, 0))  # Move channels to last

2. Squeeze singleton dimensions

   image = np.squeeze(image)

3. Verify metadata

   with tifffile.TiffFile(path) as tif:
       print(tif.pages[0].tags)

---

Problem: CSV/TSV Parse Error

Symptoms: pandas can't read measurement file

Solutions: 1. Check delimiter (comma vs tab)

   df = pd.read_csv(path, sep='\t')  # Try tab

2. Specify encoding

   df = pd.read_csv(path, encoding='utf-8')

3. Skip header rows if needed

   df = pd.read_csv(path, skiprows=2)

4. Handle missing values

   df = pd.read_csv(path, na_values=['NA', 'N/A', ''])

---

Measurement Issues

Problem: Circularity Values > 1

Symptoms: Circularity should be 0-1, but getting values > 1

Cause: Usually perimeter measurement artifacts

Solutions: 1. Smooth binary mask before measuring

   from skimage.morphology import binary_closing, disk
   binary = binary_closing(binary, disk(2))

2. Use different perimeter method 3. Cap circularity at 1.0

   circularity = np.minimum(circularity, 1.0)

---

Problem: Negative Intensities After Background Correction

Symptoms: Mean intensity < 0 after correction

Solutions: 1. Clip to zero

   corrected = np.clip(image - background, 0, None)

2. Use more conservative background estimate 3. Check background correction method is appropriate

---

Visualization Issues

Problem: Can't See Labels/Overlay

Symptoms: Overlay looks wrong or invisible

Solutions: 1. Check data types match 2. Normalize intensity ranges

   from skimage import exposure
   image_norm = exposure.rescale_intensity(image, out_range=(0, 1))

3. Use proper overlay function

   from skimage.color import label2rgb
   overlay = label2rgb(labels, image=image, bg_label=0, alpha=0.3)

---

Getting More Help

If issue persists:

1. Check the documentation - Each reference guide has detailed examples 2. Verify input data - Print shapes, dtypes, value ranges 3. Test with simple case - Try with known-good image 4. Compare with example - Use provided example scripts 5. Check package versions - Update scikit-image, scipy if needed

Useful Debugging Code

# Print image properties
print(f"Shape: {image.shape}")
print(f"Dtype: {image.dtype}")
print(f"Range: [{image.min()}, {image.max()}]")
print(f"Mean: {image.mean():.2f}")

# Check for NaN/inf
print(f"NaN values: {np.isnan(image).sum()}")
print(f"Inf values: {np.isinf(image).sum()}")

# Visualize intermediate steps
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
axes[0].imshow(image, cmap='gray')
axes[0].set_title('Original')
axes[1].imshow(binary, cmap='gray')
axes[1].set_title('Binary')
axes[2].imshow(labels, cmap='nipy_spectral')
axes[2].set_title(f'Labels (n={labels.max()})')
plt.show()