"""Antigen Puncta Quantification for fluorescence microscopy. Gradio app deployable on Hugging Face Spaces. Supports four antigens with antigen-specific analysis goals: CDA - identify and quantify puncta with similar fluorescence intensity while minimizing background signal. PTAU - detect and quantify puncta with similar fluorescence characteristics while reducing background noise. MBP - identify MBP fluorescence within cells while excluding background signal (area / intensity based, not puncta counting). SYN - identify and quantify smaller, easily missed puncta with similar fluorescence intensity. The "Calibrate on gold standard" tab tunes the detector for each antigen from the user's gold-standard ROI images with known manual counts (a supervised parameter search). Calibrated settings are used automatically by the Analyze and Batch tabs and can be exported / imported as JSON. """ import json import os import re import tempfile import gradio as gr import numpy as np import pandas as pd import tifffile from PIL import Image from scipy import ndimage as ndi from skimage import draw, feature, filters, segmentation # The LoG detection threshold never drops below this many robust noise # sigmas (in peak-amplitude terms), so very sensitive threshold settings on # dim images cannot flood the detector with noise peaks. Kept constant (not # tied to min_snr) so calibration and analysis see identical candidates. NOISE_FLOOR_SNR = 2.0 # The intensity scale used to normalize the image for detection is capped at # background + this many robust noise sigmas. Without the cap, a single # large saturated artifact (dust, a hair, a dead-pixel patch) inflates the # 99.9th-percentile scale and silently suppresses every real punctum. INTENSITY_CAP_SNR = 50.0 # Bimodality test for the intensity-similarity filter: the log-intensity # modes must differ by at least log(2) (a 2x brightness ratio) AND by this # many robust within-mode sigmas before the filter treats the candidate # population as two distinct populations (signal vs background/artifacts). SIM_BIMODAL_SEP = 4.0 SIM_MIN_LOG_RATIO = float(np.log(2.0)) # A bright mode smaller than this fraction of candidates is treated as # artifacts (dust); otherwise the brighter mode is taken as the puncta. SIM_MIN_BRIGHT_FRAC = 0.15 # A bright mode more than this far above the dim mode (log units; 30x in # linear terms) is never plausible same-stain puncta - treat it as # saturated artifacts regardless of its size (e.g. rim detections around a # clipped dust blob can outnumber 15 percent of candidates). SIM_MAX_LOG_RATIO = float(np.log(30.0)) # MBP cell masks must be at least this many robust background sigmas # brighter than the unmasked area, otherwise the field is treated as having # no cells (prevents fabricated signal on negative controls). Applied when # at least MBP_MIN_BG_FRAC of the field is outside the mask. MBP_MIN_CELL_CONTRAST = 4.0 MBP_MIN_BG_FRAC = 0.10 # When the cell mask covers nearly the whole field (no background sample to # validate against - either a confluent culture or a blank field with a low # threshold), require the within-mask positive/negative intensity split to # be genuinely bimodal: the class medians must differ by at least this many # robust within-class sigmas. Real staining separates by 8-10+; noise and # smooth illumination gradients stay below about 3.5. MBP_MIN_POS_SEPARATION = 4.0 APP_TITLE = "Antigen Puncta Quantification (CDA, PTAU, MBP, SYN)" ALL_ANTIGENS = ["CDA", "PTAU", "MBP", "SYN"] PUNCTA_ANTIGENS = ["CDA", "PTAU", "SYN"] ANTIGEN_GOALS = { "CDA": "Identify and quantify puncta with similar fluorescence intensity " "while minimizing background signal.", "PTAU": "Detect and quantify puncta with similar fluorescence " "characteristics while reducing background noise.", "MBP": "Identify MBP fluorescence within cells while excluding " "background signal (area and intensity based).", "SYN": "Identify and quantify smaller, easily missed puncta with " "similar fluorescence intensity.", } # Unified parameter schema. Puncta keys are ignored in MBP (area) mode and # the MBP-only keys (cell_sigma, cell_sensitivity, min_area) are ignored in # puncta mode, so every antigen carries a full dict. PARAM_KEYS = [ "bg_sigma", "denoise_sigma", "min_sigma", "max_sigma", "log_threshold", "min_snr", "sim_filter", "sim_tol", "cell_sigma", "cell_sensitivity", "min_area", ] DEFAULT_PARAMS = { "CDA": { "bg_sigma": 20.0, "denoise_sigma": 0.6, "min_sigma": 1.6, "max_sigma": 5.0, "log_threshold": 0.08, "min_snr": 5.0, "sim_filter": True, "sim_tol": 3.5, "cell_sigma": 8.0, "cell_sensitivity": 1.0, "min_area": 25, }, "PTAU": { "bg_sigma": 20.0, "denoise_sigma": 1.0, "min_sigma": 1.6, "max_sigma": 6.0, "log_threshold": 0.06, "min_snr": 5.0, "sim_filter": True, "sim_tol": 3.5, "cell_sigma": 8.0, "cell_sensitivity": 1.0, "min_area": 25, }, "MBP": { "bg_sigma": 40.0, "denoise_sigma": 1.5, "min_sigma": 1.6, "max_sigma": 5.0, "log_threshold": 0.08, "min_snr": 5.0, "sim_filter": False, "sim_tol": 3.5, "cell_sigma": 8.0, "cell_sensitivity": 1.0, "min_area": 25, }, "SYN": { "bg_sigma": 15.0, "denoise_sigma": 0.5, "min_sigma": 0.9, "max_sigma": 4.0, "log_threshold": 0.04, "min_snr": 4.0, "sim_filter": True, "sim_tol": 3.5, "cell_sigma": 8.0, "cell_sensitivity": 1.0, "min_area": 25, }, } CHANNEL_CHOICES = [ "Auto (strongest channel)", "Grayscale / first channel", "Red", "Green", "Blue", "Max projection of channels", ] # --------------------------------------------------------------------------- # Image loading # --------------------------------------------------------------------------- def load_image(path): """Load TIFF (via tifffile, keeps 16-bit depth) or PNG/JPG (via PIL).""" ext = os.path.splitext(str(path))[1].lower() if ext in (".tif", ".tiff"): arr = tifffile.imread(path) else: with Image.open(path) as im: arr = np.array(im) return np.asarray(arr) def _reduce_to_2d(arr, channel_choice): a = np.squeeze(np.asarray(arr)) # Higher-dimensional stacks (e.g. CZYX, TCZYX): max-project the LARGEST # leading axis first (Z or T), preserving small axes (channels) so the # channel selection below still applies. Projecting axis 0 blindly would # blend the channels of common CZYX exports. while a.ndim > 3: lead_sizes = a.shape[:-2] a = np.squeeze(a.max(axis=int(np.argmax(lead_sizes)))) if a.ndim < 2: raise ValueError("Image must be at least 2-dimensional.") if a.ndim == 2: return a.astype(np.float64) # ndim == 3: RGB(A) channel-last first, so color names keep their # usual meaning and Grayscale means the RGB mean. if a.shape[-1] in (3, 4) and a.shape[-1] <= min(a.shape[:2]): rgb = a[..., :3].astype(np.float64) named = {"Red": 0, "Green": 1, "Blue": 2} if channel_choice in named: return rgb[..., named[channel_choice]] if channel_choice == "Grayscale / first channel": return rgb.mean(axis=-1) if channel_choice == "Max projection of channels": return rgb.max(axis=-1) sums = rgb.reshape(-1, 3).sum(axis=0) return rgb[..., int(np.argmax(sums))] # Other channel-last layouts (2 or 5-6 channels): move channels first # so they are selected, not max-projected over image rows. if a.shape[-1] <= 6 and a.shape[-1] < min(a.shape[:2]): a = np.moveaxis(a, -1, 0) if a.shape[0] <= 6: # channel-first stack stack = a.astype(np.float64) named = {"Red": 0, "Green": 1, "Blue": 2} if channel_choice in named: return stack[min(named[channel_choice], stack.shape[0] - 1)] if channel_choice == "Grayscale / first channel": return stack[0] if channel_choice == "Max projection of channels": return stack.max(axis=0) sums = stack.reshape(stack.shape[0], -1).sum(axis=1) return stack[int(np.argmax(sums))] return a.max(axis=0).astype(np.float64) # z-stack MIP def to_2d(arr, channel_choice=CHANNEL_CHOICES[0]): """Reduce any loaded array to a single, finite, 2D float64 channel. Handles grayscale, RGB(A), channel-last and channel-first multichannel layouts, z-stacks (maximum intensity projection) and 4D+ stacks. Non-finite pixels (NaN/Inf padding from stitching or registration software) are replaced with the darkest finite value so they cannot silently poison percentiles and thresholds downstream. """ a = _reduce_to_2d(arr, channel_choice).astype(np.float64) if not np.isfinite(a).all(): finite = a[np.isfinite(a)] if finite.size == 0: raise ValueError("Image contains no finite pixel values.") lo, hi = float(finite.min()), float(finite.max()) a = np.nan_to_num(a, nan=lo, posinf=hi, neginf=lo) return a # --------------------------------------------------------------------------- # Core analysis # --------------------------------------------------------------------------- def _remove_small_regions(mask, min_px): """Drop connected regions smaller than min_px (version-stable helper).""" labels, n = ndi.label(mask) if n == 0: return mask sizes = np.bincount(labels.ravel()) keep = sizes >= int(min_px) keep[0] = False return keep[labels] def _robust_sigma(x): """Robust noise estimate: 1.4826 * median absolute deviation.""" x = np.asarray(x, dtype=np.float64).ravel() if x.size == 0: return 0.0 med = np.median(x) return 1.4826 * float(np.median(np.abs(x - med))) def preprocess(img2d, params): """Denoise and subtract smooth background. Returns (smoothed, background_subtracted, noise_sigma). The noise estimate is computed on the unclipped residual so it is not biased by the non-negativity clip. """ img = np.asarray(img2d, dtype=np.float64) if img.ndim != 2: raise ValueError("preprocess expects a 2D image.") d = float(params.get("denoise_sigma", 0.0)) sm = filters.gaussian(img, sigma=d, preserve_range=True) if d > 0 else img bg = filters.gaussian(sm, sigma=float(params["bg_sigma"]), preserve_range=True) resid = sm - bg noise = _robust_sigma(resid) bgsub = np.clip(resid, 0.0, None) return sm, bgsub, noise def _similarity_population(peaks): """Pick the punctum population for the similarity filter. If the candidate peak intensities are clearly bimodal (two modes at least 2x apart and well separated relative to their spread), the brighter mode is taken as the puncta and the dimmer mode as background confusers - unless the bright mode is a small minority, in which case it is treated as artifacts (dust) and the dimmer mode is kept, matching the intent 'quantify puncta with similar intensity while minimizing background'. Returns (population boolean mask, note string). """ peaks = np.asarray(peaks, dtype=np.float64) n = len(peaks) all_mask = np.ones(n, dtype=bool) if n < 8: return all_mask, "all" logp = np.log(np.maximum(peaks, 1e-9)) try: t = filters.threshold_otsu(logp) except ValueError: return all_mask, "all" bright = logp > t nb = int(bright.sum()) if nb < 2 or n - nb < 2: return all_mask, "all" mb = float(np.median(logp[bright])) md = float(np.median(logp[~bright])) spread = max(_robust_sigma(logp[bright]), _robust_sigma(logp[~bright]), 0.05) if (mb - md) < SIM_MIN_LOG_RATIO or (mb - md) / spread < SIM_BIMODAL_SEP: return all_mask, "all" if (mb - md) > SIM_MAX_LOG_RATIO: return ~bright, "dimmer mode" if nb >= max(4, SIM_MIN_BRIGHT_FRAC * n): return bright, "brighter mode" return ~bright, "dimmer mode" def _apply_similarity(df, tol): """Robust intensity-similarity filter on peak_bgsub. Returns (kept_df, rejected_df, population_note). Used by both detect_puncta and the calibration sweep so their counts always agree. """ if len(df) < 4: return df, df.iloc[0:0], "all" peaks = df["peak_bgsub"].to_numpy() pop_mask, note = _similarity_population(peaks) pop = peaks[pop_mask] med = float(np.median(pop)) sig = _robust_sigma(pop) if sig > 0: keep = np.abs(peaks - med) <= float(tol) * sig if note != "all": keep &= pop_mask else: keep = pop_mask if note != "all" else np.ones(len(df), dtype=bool) return (df[keep].reset_index(drop=True), df[~keep].reset_index(drop=True), note) def _blob_measurements(raw, bgsub, y, x, r): """Per-punctum measurements in a local window (fast for many blobs).""" h, w = raw.shape pad = int(np.ceil(r + 5)) y0, y1 = max(0, int(y) - pad), min(h, int(y) + pad + 1) x0, x1 = max(0, int(x) - pad), min(w, int(x) + pad + 1) yy, xx = np.mgrid[y0:y1, x0:x1] dist = np.hypot(yy - y, xx - x) inside = dist <= r annulus = (dist > r + 1) & (dist <= r + 4) win_raw = raw[y0:y1, x0:x1] win_bg = bgsub[y0:y1, x0:x1] if not inside.any(): return None local_bg = float(np.median(win_raw[annulus])) if annulus.any() else float( np.median(win_raw)) return { "peak_bgsub": float(win_bg[inside].max()), "mean_raw": float(win_raw[inside].mean()), "integrated_raw": float(win_raw[inside].sum()), "area_px": int(inside.sum()), "local_bg_raw": local_bg, } def detect_puncta(img2d, params): """Detect and quantify puncta. Returns (DataFrame, info dict). Pipeline: Gaussian denoise -> Gaussian background subtraction -> Laplacian-of-Gaussian blob detection -> background/noise rejection by signal-to-noise ratio -> robust intensity-similarity filter (keeps puncta whose background-subtracted peak intensity lies within sim_tol robust standard deviations of the population median). """ img = np.asarray(img2d, dtype=np.float64) sm, bgsub, noise = preprocess(img, params) info = {"noise_sigma": noise, "background_median": float(np.median(img)), "n_candidates": 0, "n_rejected_snr": 0, "n_rejected_sim": 0, "sim_population": "all", "sim_warning": ""} # Detection scale: the 99.9th percentile of the background-subtracted # image, capped at a robust ceiling so one large saturated artifact # cannot inflate the scale and suppress every real punctum. p999 = float(np.percentile(bgsub, 99.9)) med_bgsub = float(np.median(bgsub)) hi = min(p999, med_bgsub + INTENSITY_CAP_SNR * noise) if noise > 0 \ else p999 empty = pd.DataFrame(columns=[ "punctum_id", "y_px", "x_px", "radius_px", "area_px", "peak_bgsub", "mean_raw", "integrated_raw", "local_bg_raw", "snr"]) if not np.isfinite(hi) or hi <= 0: return empty, info norm = np.clip(bgsub / hi, 0.0, 1.0) # A LoG response of a matched Gaussian blob is roughly half its peak # amplitude, so NOISE_FLOOR_SNR noise sigmas of peak amplitude map to # about 0.5 * NOISE_FLOOR_SNR * noise in LoG-response terms. noise_floor = 0.5 * NOISE_FLOOR_SNR * noise / hi eff_threshold = max(float(params["log_threshold"]), noise_floor) info["effective_threshold"] = eff_threshold blobs = feature.blob_log( norm, min_sigma=float(params["min_sigma"]), max_sigma=float(params["max_sigma"]), num_sigma=6, threshold=eff_threshold, overlap=0.5, ) info["n_candidates"] = len(blobs) if len(blobs) == 0: return empty, info rows = [] for (y, x, s) in blobs: r = max(1.0, float(s) * np.sqrt(2.0)) m = _blob_measurements(img, bgsub, float(y), float(x), r) if m is None: continue snr = m["peak_bgsub"] / noise if noise > 0 else np.inf rows.append({"y_px": float(y), "x_px": float(x), "radius_px": r, "snr": float(snr), **m}) df = pd.DataFrame(rows) if df.empty: return empty, info # Background / noise rejection ("minimizing background signal"). keep = df["snr"] >= float(params["min_snr"]) info["n_rejected_snr"] = int((~keep).sum()) df = df[keep].reset_index(drop=True) # Intensity-similarity filter ("puncta with similar fluorescence # intensity"): robust median band, bimodal-population aware. if params.get("sim_filter", False): df, rejected, info["sim_population"] = \ _apply_similarity(df, params["sim_tol"]) info["n_rejected_sim"] = len(rejected) if len(rejected) > 0 and len(df) > 0: if float(rejected["peak_bgsub"].median()) > \ float(df["peak_bgsub"].median()): info["sim_warning"] = ( "The similarity filter rejected candidates brighter " "than the kept puncta; verify on the overlay that " "these are artifacts, not true puncta.") df.insert(0, "punctum_id", np.arange(1, len(df) + 1)) cols = ["punctum_id", "y_px", "x_px", "radius_px", "area_px", "peak_bgsub", "mean_raw", "integrated_raw", "local_bg_raw", "snr"] return df[cols], info def analyze_mbp(img2d, params): """MBP fluorescence within cells, excluding background. A coarse cell mask is built from a strongly smoothed copy (triangle threshold scaled by cell_sensitivity); MBP-positive pixels are then thresholded inside the cell mask only. A robust fallback threshold guards against Otsu failure when foreground dominates. Returns (metrics dict, positive_mask, cell_mask). """ img = np.asarray(img2d, dtype=np.float64) d = float(params.get("denoise_sigma", 1.0)) sm = filters.gaussian(img, sigma=d, preserve_range=True) if d > 0 else img coarse = filters.gaussian(img, sigma=float(params["cell_sigma"]), preserve_range=True) t_cell = filters.threshold_triangle(coarse) cell_mask = coarse > (t_cell * float(params["cell_sensitivity"])) cell_mask = _remove_small_regions(cell_mask, 200) cell_mask = ndi.binary_fill_holes(cell_mask) if not cell_mask.any(): raise ValueError( "No cell region found. Lower the cell sensitivity parameter.") # Validate that the mask is genuinely brighter than the rest of the # field. Triangle thresholding always returns a threshold, so without # this check a cell-free negative control gets a fabricated 'cell mask' # carved out of background structure, with roughly half of it then # called MBP-positive. no_cells_msg = ( "No clear cell region or MBP signal detected. This looks like a " "negative control or an empty field, so no MBP measurement is " "reported (a fabricated measurement would be worse). To force a " "measurement anyway, switch to Manual parameters and lower 'MBP: " "cell mask sensitivity'.") bg_frac = 1.0 - float(cell_mask.mean()) if bg_frac >= MBP_MIN_BG_FRAC: bg_coarse = coarse[~cell_mask] contrast = float(np.median(coarse[cell_mask]) - np.median(bg_coarse)) bg_spread = max(_robust_sigma(bg_coarse), 1e-9) if contrast < MBP_MIN_CELL_CONTRAST * bg_spread: raise ValueError( f"{no_cells_msg} (Cell mask contrast {contrast:.1f} vs " f"background variation {bg_spread:.1f}.)") vals = sm[cell_mask] try: t_pos = float(filters.threshold_otsu(vals)) except ValueError: t_pos = float(np.median(vals)) pos = (sm > t_pos) & cell_mask # Guard against threshold collapse (nearly all cell pixels "positive"). frac = pos.sum() / cell_mask.sum() if frac > 0.85: t_pos = float(np.median(vals) + 2.0 * _robust_sigma(vals)) pos = (sm > t_pos) & cell_mask # When the mask swallowed nearly the whole field there was no # background sample to validate against (either a confluent culture or # a blank with a low threshold). Real staining splits into two well # separated intensity classes inside the mask; noise and illumination # gradients do not. if bg_frac < MBP_MIN_BG_FRAC: below = vals[vals <= t_pos] above = vals[vals > t_pos] if below.size < 10 or above.size < 10: raise ValueError(no_cells_msg) separation = (float(np.median(above)) - float(np.median(below))) / \ max(_robust_sigma(below), _robust_sigma(above), 1e-9) if separation < MBP_MIN_POS_SEPARATION: raise ValueError( f"{no_cells_msg} (Within-field intensity separation " f"{separation:.1f}, below the required " f"{MBP_MIN_POS_SEPARATION:.1f}.)") pos = _remove_small_regions(pos, int(params["min_area"])) cell_area = int(cell_mask.sum()) pos_area = int(pos.sum()) bg_pixels = img[~cell_mask] metrics = { "cell_area_px": cell_area, "cell_mask_coverage_of_field": cell_area / img.size, "mbp_positive_area_px": pos_area, "mbp_area_fraction_of_cells": pos_area / cell_area if cell_area else 0.0, "mbp_mean_intensity_raw": float(img[pos].mean()) if pos_area else 0.0, "mbp_integrated_intensity_raw": float(img[pos].sum()), "background_median_raw": float(np.median(bg_pixels)) if bg_pixels.size else float(np.median(img)), "positive_threshold": t_pos, } return metrics, pos, cell_mask # --------------------------------------------------------------------------- # Overlays # --------------------------------------------------------------------------- def to_display_rgb(img2d): lo, hi = np.percentile(img2d, (1.0, 99.7)) if hi <= lo: hi = lo + 1.0 g = np.clip((np.asarray(img2d, dtype=np.float64) - lo) / (hi - lo), 0, 1) g8 = (g * 255).astype(np.uint8) return np.stack([g8, g8, g8], axis=-1) def overlay_puncta(img2d, df): rgb = to_display_rgb(img2d) h, w = img2d.shape for _, row in df.iterrows(): cy, cx = int(round(row["y_px"])), int(round(row["x_px"])) rad = int(round(row["radius_px"])) + 2 for rr_off in (0, 1): rr, cc = draw.circle_perimeter(cy, cx, rad + rr_off, shape=(h, w)) rgb[rr, cc] = (255, 64, 64) return rgb def overlay_mbp(img2d, pos_mask, cell_mask): rgb = to_display_rgb(img2d) cell_edge = segmentation.find_boundaries(cell_mask, mode="outer") pos_edge = segmentation.find_boundaries(pos_mask, mode="outer") rgb[cell_edge] = (90, 170, 255) rgb[pos_edge] = (255, 64, 64) return rgb # --------------------------------------------------------------------------- # Parameter resolution and shared analysis entry point # --------------------------------------------------------------------------- def sliders_to_params(values): p = dict(zip(PARAM_KEYS, values)) p["sim_filter"] = bool(p["sim_filter"]) p["min_area"] = int(p["min_area"]) return p def resolve_params(antigen, source, calibrated_state, slider_values): if source.startswith("Manual"): return sliders_to_params(slider_values), "manual sliders" if isinstance(calibrated_state, dict) and antigen in calibrated_state: return dict(calibrated_state[antigen]), "calibrated" return dict(DEFAULT_PARAMS[antigen]), "antigen defaults" def analyze_one(img2d, antigen, params, pixel_size_um): """Run the antigen-appropriate analysis. Returns (summary_rows, df, overlay_rgb, status_text).""" h, w = img2d.shape px = float(pixel_size_um) if pixel_size_um else 0.0 if antigen == "MBP": metrics, pos, cells = analyze_mbp(img2d, params) overlay = overlay_mbp(img2d, pos, cells) summary = [ ("Cell area (px)", metrics["cell_area_px"]), ("Cell mask coverage of field", round(metrics["cell_mask_coverage_of_field"], 4)), ("MBP-positive area (px)", metrics["mbp_positive_area_px"]), ("MBP area fraction of cells", round(metrics["mbp_area_fraction_of_cells"], 4)), ("MBP mean intensity (raw)", round(metrics["mbp_mean_intensity_raw"], 2)), ("MBP integrated intensity (raw)", round(metrics["mbp_integrated_intensity_raw"], 1)), ("Background median (raw)", round(metrics["background_median_raw"], 2)), ("Positive threshold used", round(metrics["positive_threshold"], 2)), ] if px > 0: summary.insert(3, ("MBP-positive area (um^2)", round(metrics["mbp_positive_area_px"] * px * px, 2))) df = pd.DataFrame(summary, columns=["metric", "value"]) status = ( "MBP mode: fluorescence measured inside the detected cell mask " "only; pixels outside cells are treated as background and " "excluded. Cell outline shown in blue, MBP-positive regions in " "red.") if metrics["mbp_area_fraction_of_cells"] > 0.5: status += ( " Warning: most of the detected cell area is MBP-positive. " "If this field is confluent (cells everywhere), the cell " "mask may have collapsed onto the bright signal - check the " "blue outline in the overlay and, if needed, lower 'MBP: " "cell mask sensitivity' in Manual mode.") return summary, df, overlay, status df, info = detect_puncta(img2d, params) overlay = overlay_puncta(img2d, df) n = len(df) area_px = h * w summary = [ ("Puncta count", n), ("Density (per megapixel)", round(n / area_px * 1e6, 2)), ("Mean peak intensity (bg-subtracted)", round(float(df["peak_bgsub"].mean()), 2) if n else 0), ("Median peak intensity (bg-subtracted)", round(float(df["peak_bgsub"].median()), 2) if n else 0), ("Mean radius (px)", round(float(df["radius_px"].mean()), 2) if n else 0), ("Total integrated intensity (raw)", round(float(df["integrated_raw"].sum()), 1) if n else 0), ("Candidates before filtering", info["n_candidates"]), ("Rejected as background/noise (SNR)", info["n_rejected_snr"]), ("Rejected by intensity-similarity filter", info["n_rejected_sim"]), ("Image noise sigma (robust)", round(info["noise_sigma"], 3)), ("Background median (raw)", round(info["background_median"], 2)), ] if px > 0: area_mm2 = area_px * (px ** 2) / 1e6 summary.insert(2, ("Density (per mm^2)", round(n / area_mm2, 2) if area_mm2 else 0)) if n: summary.insert(6, ("Mean radius (um)", round(float(df["radius_px"].mean()) * px, 3))) summary_df = pd.DataFrame(summary, columns=["metric", "value"]) status = ( f"{antigen}: {n} puncta kept. " f"{info['n_rejected_snr']} candidates rejected as background/noise, " f"{info['n_rejected_sim']} rejected by the intensity-similarity " f"filter.") if info.get("sim_population") == "brighter mode": status += (" The candidate intensities were bimodal; the brighter " "population was quantified as puncta and the dimmer " "population was treated as background.") elif info.get("sim_population") == "dimmer mode": status += (" A small number of much brighter outliers were treated " "as artifacts and excluded.") if info.get("sim_warning"): status += " " + info["sim_warning"] return summary, summary_df, overlay, status, df # --------------------------------------------------------------------------- # Gradio callbacks # --------------------------------------------------------------------------- def _write_csv(df, prefix): out_dir = tempfile.mkdtemp(prefix="puncta_") path = os.path.join(out_dir, f"{prefix}.csv") df.to_csv(path, index=False) return path def cb_analyze(file_path, antigen, channel_choice, source, pixel_size, calibrated_state, *slider_values): if not file_path: raise gr.Error("Upload an image first.") try: img = to_2d(load_image(file_path), channel_choice) except Exception as exc: raise gr.Error(f"Could not read the image: {exc}") if min(img.shape) < 32: raise gr.Error( f"Image is too small after channel reduction (got shape " f"{img.shape}; minimum 32 x 32 pixels). If the original image " f"is larger, its channel layout may be unusual - try a " f"different Channel setting.") params, origin = resolve_params(antigen, source, calibrated_state, slider_values) try: result = analyze_one(img, antigen, params, pixel_size) except ValueError as exc: raise gr.Error(str(exc)) if antigen == "MBP": _, summary_df, overlay, status = result csv_path = _write_csv(summary_df, f"{antigen}_metrics") puncta_df = pd.DataFrame() else: _, summary_df, overlay, status, puncta_df = result csv_path = _write_csv(puncta_df, f"{antigen}_puncta") status = f"{status} Parameters: {origin}." shown = puncta_df.head(200).round(3) if not puncta_df.empty else puncta_df return overlay, summary_df, shown, csv_path, status def cb_batch(file_paths, antigen, channel_choice, source, pixel_size, calibrated_state, *slider_values): if not file_paths: raise gr.Error("Upload one or more images first.") params, origin = resolve_params(antigen, source, calibrated_state, slider_values) rows = [] for fp in file_paths: name = os.path.basename(str(fp)) try: img = to_2d(load_image(fp), channel_choice) if min(img.shape) < 32: raise ValueError( f"Degenerate image after channel reduction " f"(shape {img.shape}); check the channel setting and " f"the file's layout.") result = analyze_one(img, antigen, params, pixel_size) summary = result[0] row = {"file": name, "antigen": antigen, "error": ""} row.update({metric: value for metric, value in summary}) except Exception as exc: row = {"file": name, "antigen": antigen, "error": str(exc)} rows.append(row) df = pd.DataFrame(rows) csv_path = _write_csv(df, f"{antigen}_batch_results") status = (f"Processed {len(rows)} file(s) with {origin} parameters " f"for {antigen}.") return df, csv_path, status CAL_MIN_SIGMAS = [0.8, 1.2, 1.6] CAL_THRESHOLDS = [0.02, 0.035, 0.055, 0.08, 0.12, 0.18] CAL_SNRS = [2.0, 3.0, 4.0, 5.0, 7.0] CAL_SIM = [(False, 3.5), (True, 2.5), (True, 3.5), (True, 5.0)] def _count_with_filters(df_all, noise, min_snr, sim_filter, sim_tol): """Apply the SNR and similarity filters to a precomputed candidate table, using the same _apply_similarity as detect_puncta so calibrated settings reproduce exactly at analysis time.""" if df_all.empty: return 0 snr = df_all["peak_bgsub"] / noise if noise > 0 else np.full( len(df_all), np.inf) df = df_all[snr >= min_snr].reset_index(drop=True) if sim_filter: df, _, _ = _apply_similarity(df, sim_tol) return len(df) def _strictness(p): """Ordering key so calibration ties resolve to the LEAST permissive settings. Without this, a clean gold-standard ROI (where many combinations reproduce the count exactly) would select the most permissive combination and flood routine images with false positives.""" return (p["log_threshold"], p["min_snr"], 1 if p["sim_filter"] else 0, p["min_sigma"]) def cb_calibrate(file_paths, counts_text, antigen, channel_choice, calibrated_state, progress=gr.Progress()): if antigen not in PUNCTA_ANTIGENS: raise gr.Error("Calibration applies to the puncta antigens " "(CDA, PTAU, SYN). MBP uses area-based analysis.") if not file_paths: raise gr.Error("Upload gold-standard ROI images first.") tokens = [t for t in re.split(r"[,\s]+", (counts_text or "").strip()) if t] try: counts = [int(float(t)) for t in tokens] except ValueError: raise gr.Error("Expected counts must be numbers separated by commas " "or spaces.") if len(counts) != len(file_paths): raise gr.Error( f"Got {len(file_paths)} image(s) but {len(counts)} count(s). " "Enter one expected count per image, in upload order.") base = dict(DEFAULT_PARAMS[antigen]) images = [] for fp in file_paths: try: images.append(to_2d(load_image(fp), channel_choice)) except Exception as exc: raise gr.Error(f"Could not read {os.path.basename(str(fp))}: {exc}") # Precompute candidate tables for each (image, min_sigma, threshold) # so the cheap SNR / similarity filters can be swept without re-running # blob detection. total = len(images) * len(CAL_MIN_SIGMAS) * len(CAL_THRESHOLDS) step = 0 candidates = {} noises = {} for i, img in enumerate(images): for ms in CAL_MIN_SIGMAS: for th in CAL_THRESHOLDS: step += 1 progress(step / total, desc=f"Scanning detector settings ({step}/{total})") p = dict(base, min_sigma=ms, log_threshold=th, min_snr=0.0, sim_filter=False) df, info = detect_puncta(images[i], p) candidates[(i, ms, th)] = df noises[i] = info["noise_sigma"] def mape(pred): return float(np.mean([abs(p - c) / max(c, 1) for p, c in zip(pred, counts)])) best = None for ms in CAL_MIN_SIGMAS: for th in CAL_THRESHOLDS: for snr in CAL_SNRS: for sim_f, sim_t in CAL_SIM: pred = [_count_with_filters(candidates[(i, ms, th)], noises[i], snr, sim_f, sim_t) for i in range(len(images))] err = mape(pred) cand = dict(base, min_sigma=ms, log_threshold=th, min_snr=snr, sim_filter=sim_f, sim_tol=sim_t) if (best is None or err < best["err"] - 1e-9 or (abs(err - best["err"]) <= 1e-9 and _strictness(cand) > _strictness(best["params"]))): best = {"err": err, "pred": pred, "params": cand} default_pred = [] for i in range(len(images)): df, info = detect_puncta(images[i], base) default_pred.append(len(df)) state = dict(calibrated_state) if isinstance(calibrated_state, dict) else {} state[antigen] = best["params"] report = pd.DataFrame({ "file": [os.path.basename(str(fp)) for fp in file_paths], "expected_count": counts, "default_params_count": default_pred, "calibrated_count": best["pred"], }) p = best["params"] msg = ( f"Calibration complete for {antigen}. " f"Mean absolute count error: {mape(default_pred) * 100:.1f} percent " f"with defaults, {best['err'] * 100:.1f} percent after calibration. " f"Selected: min_sigma={p['min_sigma']}, " f"log_threshold={p['log_threshold']}, min_snr={p['min_snr']}, " f"similarity filter=" f"{'on, tol ' + str(p['sim_tol']) if p['sim_filter'] else 'off'}. " f"The Analyze and Batch tabs will now use these settings for " f"{antigen} automatically.") return state, report, msg def cb_export_params(calibrated_state): merged = {} for a in ALL_ANTIGENS: if isinstance(calibrated_state, dict) and a in calibrated_state: merged[a] = calibrated_state[a] else: merged[a] = DEFAULT_PARAMS[a] out_dir = tempfile.mkdtemp(prefix="puncta_params_") path = os.path.join(out_dir, "antigen_parameters.json") with open(path, "w", encoding="utf-8") as f: json.dump(merged, f, indent=2) return path def cb_import_params(file_path, calibrated_state): if not file_path: raise gr.Error("Upload a parameters JSON file first.") try: with open(file_path, "r", encoding="utf-8") as f: data = json.load(f) except Exception as exc: raise gr.Error(f"Could not read JSON: {exc}") state = dict(calibrated_state) if isinstance(calibrated_state, dict) else {} loaded = [] for a in ALL_ANTIGENS: if a in data and isinstance(data[a], dict): p = dict(DEFAULT_PARAMS[a]) for k in PARAM_KEYS: if k in data[a]: p[k] = data[a][k] p["sim_filter"] = bool(p["sim_filter"]) p["min_area"] = int(p["min_area"]) state[a] = p loaded.append(a) if not loaded: raise gr.Error("No antigen parameter blocks (CDA, PTAU, MBP, SYN) " "found in the JSON file.") return state, f"Loaded parameters for: {', '.join(loaded)}." def cb_antigen_change(antigen, calibrated_state): if isinstance(calibrated_state, dict) and antigen in calibrated_state: p = calibrated_state[antigen] else: p = DEFAULT_PARAMS[antigen] updates = [gr.update(value=p[k]) for k in PARAM_KEYS] return updates + [gr.update(value=f"Goal: {ANTIGEN_GOALS[antigen]}")] # --------------------------------------------------------------------------- # UI # --------------------------------------------------------------------------- def _make_sliders(): """Advanced parameter controls, in PARAM_KEYS order.""" p = DEFAULT_PARAMS["CDA"] return [ gr.Slider(5, 100, value=p["bg_sigma"], step=1, label="Background sigma (px)", info="Scale of the smooth background that is subtracted."), gr.Slider(0, 3, value=p["denoise_sigma"], step=0.1, label="Denoise sigma (px)", info="Gaussian pre-smoothing to reduce noise."), gr.Slider(0.5, 5, value=p["min_sigma"], step=0.1, label="Min punctum sigma (px)", info="Lower this to catch smaller puncta (as for SYN)."), gr.Slider(2, 12, value=p["max_sigma"], step=0.5, label="Max punctum sigma (px)"), gr.Slider(0.005, 0.4, value=p["log_threshold"], step=0.005, label="Detection threshold", info="Lower is more sensitive; higher suppresses " "background."), gr.Slider(0, 10, value=p["min_snr"], step=0.5, label="Minimum signal-to-noise ratio", info="Rejects candidates close to the background noise " "level."), gr.Checkbox(value=p["sim_filter"], label="Apply intensity-similarity filter", info="Keep only puncta with similar fluorescence " "intensity (robust median band)."), gr.Slider(1, 10, value=p["sim_tol"], step=0.5, label="Similarity tolerance (robust sigmas)"), gr.Slider(2, 30, value=p["cell_sigma"], step=1, label="MBP: cell mask smoothing sigma (px)"), gr.Slider(0.3, 3, value=p["cell_sensitivity"], step=0.05, label="MBP: cell mask sensitivity", info="Lower includes dimmer cell regions."), gr.Slider(5, 500, value=p["min_area"], step=5, label="MBP: minimum region area (px)"), ] with gr.Blocks(title=APP_TITLE) as demo: gr.Markdown(f"# {APP_TITLE}") gr.Markdown( "Upload fluorescence images (8/16-bit TIFF, PNG or JPG), pick the " "antigen, and quantify. Calibrate the detector on your gold-standard " "ROIs in the Calibrate tab; calibrated settings are then used " "automatically.") calibrated_state = gr.State({}) with gr.Row(): antigen_dd = gr.Dropdown(ALL_ANTIGENS, value="CDA", label="Antigen") channel_dd = gr.Dropdown(CHANNEL_CHOICES, value=CHANNEL_CHOICES[0], label="Channel") source_radio = gr.Radio( ["Auto (calibrated when available, else antigen defaults)", "Manual (advanced sliders below)"], value="Auto (calibrated when available, else antigen defaults)", label="Parameter source") pixel_size_num = gr.Number(value=0, label="Pixel size (um/px)", info="0 = report in pixels only.") goal_md = gr.Markdown(f"Goal: {ANTIGEN_GOALS['CDA']}") with gr.Accordion("Advanced parameters (used in Manual mode)", open=False): sliders = _make_sliders() with gr.Tabs(): with gr.Tab("Analyze"): with gr.Row(): with gr.Column(): an_file = gr.File(label="Image", type="filepath", file_types=["image", ".tif", ".tiff"]) an_btn = gr.Button("Run analysis", variant="primary") an_status = gr.Markdown("") with gr.Column(): an_overlay = gr.Image(label="Detection overlay", type="numpy") an_summary = gr.Dataframe(label="Summary metrics", interactive=False) an_table = gr.Dataframe(label="Per-punctum measurements " "(first 200 shown)", interactive=False) an_csv = gr.File(label="Download full results (CSV)") with gr.Tab("Calibrate on gold standard"): gr.Markdown( "Train the detector on your gold-standard ROIs: upload ROI " "images for one puncta antigen (CDA, PTAU or SYN), enter the " "manually verified punctum count for each image (same order, " "separated by commas), and run calibration. The app searches " "detector settings to match your counts and stores them for " "that antigen.") cal_files = gr.File(label="Gold-standard ROI images", type="filepath", file_count="multiple", file_types=["image", ".tif", ".tiff"]) cal_counts = gr.Textbox( label="Expected counts (one per image, comma separated)", placeholder="e.g. 42, 37, 51") cal_btn = gr.Button("Run calibration", variant="primary") cal_report = gr.Dataframe(label="Calibration report", interactive=False) cal_status = gr.Markdown("") with gr.Row(): exp_btn = gr.Button("Export parameters (JSON)") exp_file = gr.File(label="Exported parameters") with gr.Row(): imp_file = gr.File(label="Import parameters (JSON)", type="filepath", file_types=[".json"]) imp_btn = gr.Button("Load imported parameters") imp_status = gr.Markdown("") with gr.Tab("Batch"): gr.Markdown("Process many images with the current settings for " "the selected antigen.") b_files = gr.File(label="Images", type="filepath", file_count="multiple", file_types=["image", ".tif", ".tiff"]) b_btn = gr.Button("Run batch", variant="primary") b_table = gr.Dataframe(label="Batch results", interactive=False) b_csv = gr.File(label="Download batch results (CSV)") b_status = gr.Markdown("") with gr.Tab("Help"): gr.Markdown( "\n".join( ["### Antigen-specific behavior"] + [f"- **{a}**: {ANTIGEN_GOALS[a]}" for a in ALL_ANTIGENS] + [ "", "### How it works", "Puncta antigens (CDA, PTAU, SYN): Gaussian denoise, " "smooth-background subtraction, " "Laplacian-of-Gaussian blob detection, rejection of " "low signal-to-noise candidates (background " "minimization), then a robust intensity-similarity " "filter that keeps puncta whose background-" "subtracted peak intensity lies within a tolerance " "band around the population median. SYN defaults " "use a smaller minimum punctum size and a more " "sensitive threshold to recover smaller, easily " "missed puncta.", "", "MBP: a coarse cell mask is segmented from a " "strongly smoothed copy of the image, and MBP-" "positive fluorescence is thresholded inside cells " "only, so background outside cells never " "contributes to the measurements. Fields with no " "clear cell region (for example negative controls) " "are reported as an error rather than a fabricated " "measurement, and a warning is shown when the cell " "mask may have collapsed onto the bright signal in " "confluent fields.", "", "### Calibration (training on your gold standard)", "Upload gold-standard ROI images with known manual " "counts; the app performs a supervised search over " "detector settings (punctum size, detection " "threshold, signal-to-noise cutoff, similarity " "tolerance) and keeps the combination that best " "reproduces your counts. Export the JSON to reuse " "the calibration later or on another deployment.", "", "### Input formats", "8/16-bit grayscale or RGB TIFF, PNG, JPG. " "Multi-channel stacks: choose the channel above. " "Z-stacks are reduced by maximum intensity " "projection.", ])) an_btn.click( cb_analyze, inputs=[an_file, antigen_dd, channel_dd, source_radio, pixel_size_num, calibrated_state] + sliders, outputs=[an_overlay, an_summary, an_table, an_csv, an_status]) b_btn.click( cb_batch, inputs=[b_files, antigen_dd, channel_dd, source_radio, pixel_size_num, calibrated_state] + sliders, outputs=[b_table, b_csv, b_status]) cal_btn.click( cb_calibrate, inputs=[cal_files, cal_counts, antigen_dd, channel_dd, calibrated_state], outputs=[calibrated_state, cal_report, cal_status]) exp_btn.click(cb_export_params, inputs=[calibrated_state], outputs=[exp_file]) imp_btn.click(cb_import_params, inputs=[imp_file, calibrated_state], outputs=[calibrated_state, imp_status]) antigen_dd.change(cb_antigen_change, inputs=[antigen_dd, calibrated_state], outputs=sliders + [goal_md]) if __name__ == "__main__": demo.launch()