""" Cochlear Neurofilament Tracer — HuggingFace Gradio app ====================================================== Traces auditory-nerve fibers (Neurofilament channel) in confocal z-stacks of the organ of Corti, uses the Myo7a hair-cell channel to separate IHC-innervating from OHC-innervating fibers, and reports per-region quantification (number of fibers, diameter, length, branch points, area covered, fibers-per-hair-cell and fiber direction) plus a black-background skeleton image, a trace-on-original overlay, and an Excel workbook. Accepts Zeiss **CZI** z-stacks and generic **TIFF** stacks. """ import os import tempfile import zipfile import traceback import numpy as np import pandas as pd import gradio as gr from skimage.io import imsave import processing as P OUT_DIR = tempfile.mkdtemp(prefix="neuron_tracer_") METRIC_COLS = [ "File", "Frequency region", "Region", "Number of fibers", "Hair cells (Myo7a)", "Fibers / hair cell", "Total length (um)", "Length / hair cell (um)", "Mean diameter (um)", "Median diameter (um)", "Branch points", "Radial fibers (IHC->OHC)", "Off-axis fibers", "Radial fibers (%)", "Area covered (um^2)", "FOV area (um^2)", "Area covered (% of FOV)", ] VIEW_OVERLAY = "Overlay (trace on original)" VIEW_SKELETON = "Traced skeleton (white on black)" VIEW_ORIGINAL = "Original (Neurofilament)" VIEW_CHOICES = [VIEW_OVERLAY, VIEW_SKELETON, VIEW_ORIGINAL] # --------------------------------------------------------------------------- # # Small helpers # --------------------------------------------------------------------------- # def _channel_choices(img: P.LoadedImage): choices = [] for i, ch in enumerate(img.channels): dye = P.channel_dye_label(ch) choices.append((f"Ch {i} — {ch.get('name', '?')} ({dye})", i)) return choices def _side_choices(axis: str): """OHC-side options phrased by IMAGE position, so they can't be misread. The stored value is the IHC side ('low'/'high' = smaller/larger coordinate) that the pipeline consumes; it is the OPPOSITE of the OHC position chosen here, so the labels read in OHC terms while the region math is unchanged. """ a = (axis or "Y").upper() if a == "X": return [("OHC is on the LEFT of the line", "high"), ("OHC is on the RIGHT of the line", "low")] if a == "Z": return [("OHC is the SHALLOW side (low Z)", "high"), ("OHC is the DEEP side (high Z)", "low")] return [("OHC is ABOVE the line (top of image)", "high")] def _boundary_preview(myo_u8, boundary_frac, axis="Y", angle=0.0, curvature=0.0, centroids=None, side=None): cents = centroids if centroids is not None else np.zeros((0, 2)) bfrac = None if boundary_frac is None else float(boundary_frac) return P.hair_cell_overlay(myo_u8, cents, bfrac, axis, side, float(angle), float(curvature)) def _save_png(arr, name): path = os.path.join(OUT_DIR, name) imsave(path, arr) return path def _manual(value): """Turn a manual hair-cell Number (0 = use detection) into an override.""" try: v = int(round(float(value))) except (TypeError, ValueError): return None return v if v > 0 else None def _build_excel(rows, path): """Write a tidy per-region sheet plus a frequency×region summary sheet.""" df = pd.DataFrame(rows, columns=METRIC_COLS) with pd.ExcelWriter(path, engine="openpyxl") as xl: df.to_excel(xl, sheet_name="Per region", index=False) # Summary: only IHC / OHC rows, pivoted by frequency region. sub = df[df["Region"].isin(["IHC region", "OHC region"])] if not sub.empty: for metric in ["Number of fibers", "Hair cells (Myo7a)", "Fibers / hair cell", "Total length (um)", "Length / hair cell (um)", "Mean diameter (um)", "Branch points", "Radial fibers (%)", "Area covered (um^2)"]: if metric not in sub: continue # Coerce to numeric first: blanks ("" for a non-applicable metric, # e.g. a region with no hair cells) become NaN so pivot_table's # mean works instead of erroring on an object-dtype column. vals = pd.to_numeric(sub[metric], errors="coerce") if vals.isna().all(): continue piv = sub.assign(_val=vals).pivot_table( index="Frequency region", columns="Region", values="_val", aggfunc="mean") sheet = metric.split(" (")[0][:28].replace("/", "per") piv.to_excel(xl, sheet_name=f"{sheet}") return df # --------------------------------------------------------------------------- # # Single-image interactive flow # --------------------------------------------------------------------------- # def load_and_preview(file_obj): empty = (None, gr.update(), gr.update(), gr.update(), gr.update(), None, None, None, None) if file_obj is None: return empty + ("Please upload a CZI or TIFF file.",) try: img = P.load_image(file_obj) except Exception as e: return empty + (f"Failed to load file:\n{e}\n{traceback.format_exc()}",) nf, myo = P.guess_channels(img) choices = _channel_choices(img) freq = P.detect_frequency(img.source_name) bfrac = P.suggest_boundary(img.data[myo]) nf_prev = P.channel_preview(img.data[nf]) dz_, dy_, dx_ = img.voxel status = (f"Loaded **{img.source_name}** — shape " f"{img.data.shape} (C,Z,Y,X)\n\n" f"Voxel size: dz={dz_:.3f}, dy={dy_:.4f}, dx={dx_:.4f} µm. " f"Detected frequency region: **{freq}**.\n\n" f"Auto-picked Neurofilament = Ch {nf}, Myo7a = Ch {myo}. " f"Pick an **analysis scope** to preview the Myo7a boundary, then " f"press **Run analysis**.") # Myo7a (MIP) + boundary preview is intentionally NOT shown on load — it # appears once the user chooses an analysis scope (see the scope radio). return (img, gr.update(choices=choices, value=nf), gr.update(choices=choices, value=myo), gr.update(value=freq), gr.update(value=round(bfrac, 3)), nf_prev, None, None, None, status) def update_side_labels(axis, current=None): """Re-label the OHC-side radio for the split axis, keeping the value valid. The Y axis offers a single option, so a value picked on X/Z ("low") must be clamped when switching back to Y — otherwise the Radio would carry a value that is no longer one of its choices and error out on the next run. """ choices = _side_choices(axis) valid = {v for _, v in choices} value = current if current in valid else choices[0][1] return gr.update(choices=choices, value=value) def switch_view(view, view_state): if not view_state: return None return view_state.get(view) def load_draw_background(img, myo_idx): """Load the Myo7a MIP into the drawing editor so the user can trace on it.""" if img is None or myo_idx is None: return None myo_u8 = P.channel_preview(img.data[int(myo_idx)]) rgb = np.stack([myo_u8] * 3, axis=-1) return {"background": rgb, "layers": [], "composite": rgb} def apply_drawn_boundary(editor_value, axis, mode, img, myo_idx, hc_cents, side): """Turn the freehand stroke into the active boundary and refresh the preview.""" prof = _extract_drawn_profile(editor_value, axis) if prof is None: return (gr.update(), gr.update(), "Draw the boundary across the image with the brush first, then " "press **Apply drawn boundary**. (The Z axis has no in-plane " "line to draw.)") myo_update = gr.update() if (mode or "").startswith("Split") and img is not None and myo_idx is not None: cents = hc_cents if hc_cents is not None else np.zeros((0, 2)) myo_update = _myo_preview_custom(img, myo_idx, prof["profile"], axis, cents, side) msg = (f"Custom **{prof['axis']}-axis** boundary applied — analysis and the " f"overlays will use your drawn curve instead of the sliders. Choose " f"**Split IHC vs OHC** to apply it; press **Clear drawn boundary** to " f"revert to the sliders.") return prof, myo_update, msg def clear_drawn_boundary(img, myo_idx, boundary, axis, hc_cents, side, angle, curvature, mode): """Discard the drawn curve; revert the preview to the slider boundary.""" prev = refresh_boundary_preview(img, myo_idx, boundary, axis, hc_cents, side, angle, curvature, mode, None) return None, prev, "Custom boundary cleared — using the sliders again." def _extract_drawn_profile(editor_value, axis): """Turn a freehand stroke drawn in the ImageEditor into a normalized boundary profile (values in [0,1]) for the given split axis. axis="Y" -> one boundary row per column (the usual, near-horizontal curve); axis="X" -> one boundary column per row. Returns ``{"axis": "Y"|"X", "profile": np.ndarray}`` or ``None`` if nothing was drawn (or the axis is Z, which has no in-plane line). """ a = (axis or "Y").upper() if not editor_value or a == "Z": return None layers = editor_value.get("layers") or [] paint = None for layer in layers: arr = np.asarray(layer) if arr.ndim == 3 and arr.shape[2] == 4: m = arr[..., 3] > 10 # painted where alpha > 0 elif arr.ndim == 3: m = arr.any(axis=-1) else: m = arr > 10 paint = m if paint is None else (paint | m) if paint is None or not paint.any(): # Fallback: difference between the composite and the background. comp, bg = editor_value.get("composite"), editor_value.get("background") if comp is not None and bg is not None: comp, bg = np.asarray(comp), np.asarray(bg) if comp.shape[:2] == bg.shape[:2] and comp.ndim == 3 and bg.ndim == 3: d = np.abs(comp[..., :3].astype(int) - bg[..., :3].astype(int)).sum(-1) paint = d > 30 if paint is None or not paint.any(): return None ny, nx = paint.shape if a == "Y": cols = np.where(paint.any(axis=0))[0] yvals = np.array([np.flatnonzero(paint[:, x]).mean() for x in cols]) yfrac = yvals / max(ny - 1, 1) prof = np.interp(np.linspace(0, 1, nx), cols / max(nx - 1, 1), yfrac) else: # a == "X" rows = np.where(paint.any(axis=1))[0] xvals = np.array([np.flatnonzero(paint[y, :]).mean() for y in rows]) xfrac = xvals / max(nx - 1, 1) prof = np.interp(np.linspace(0, 1, ny), rows / max(ny - 1, 1), xfrac) return {"axis": a, "profile": np.clip(prof, 0.0, 1.0)} def _draw_curve_on(rgb, profile, axis, color=(255, 238, 0), thickness=2): """Draw a normalized boundary profile as a curve on an RGB image (in place).""" ny, nx = rgb.shape[:2] half = max(thickness // 2, 0) if (axis or "Y").upper() == "Y": ys = np.clip(np.round(P._resample_profile(profile, nx) * ny), 0, ny - 1).astype(int) cols = np.arange(nx) for d in range(-half, half + 1): rgb[np.clip(ys + d, 0, ny - 1), cols] = color elif (axis or "Y").upper() == "X": xs = np.clip(np.round(P._resample_profile(profile, ny) * nx), 0, nx - 1).astype(int) rows = np.arange(ny) for d in range(-half, half + 1): rgb[rows, np.clip(xs + d, 0, nx - 1)] = color return rgb def _myo_preview_custom(img, myo_idx, profile, axis, cents, ihc_side): """Myo7a MIP with a freehand-drawn boundary curve (and hair cells coloured by the resulting region masks), matching hair_cell_overlay's look.""" myo_u8 = P.channel_preview(img.data[int(myo_idx)]) ny, nx = myo_u8.shape rgb = np.stack([myo_u8] * 3, axis=-1).copy() ihc_mask = None if ihc_side is not None and (axis or "Y").upper() in ("Y", "X"): ihc_mask, _ = P.make_region_masks((ny, nx), 0.5, ihc_side, axis, profile=profile) if cents is not None: for (y, x) in np.asarray(cents).astype(int): yc, xc = min(max(y, 0), ny - 1), min(max(x, 0), nx - 1) col = (0, 255, 0) if ihc_mask is not None: col = (0, 220, 255) if ihc_mask[yc, xc] else (255, 60, 200) rgb[max(0, y - 3):y + 4, max(0, x - 3):x + 4] = col _draw_curve_on(rgb, profile, axis) return rgb def refresh_boundary_preview(img, myo_idx, boundary, axis, hc_cents, side, angle, curvature, mode, custom=None): """Redraw the Myo7a preview — but only once an analysis scope is chosen. The preview is hidden until the user picks a scope (so it does not appear automatically on load). 'Split IHC vs OHC' draws the boundary and the IHC side; 'Inner hair cells only' draws no line (the whole field is IHC); 'Whole field only' draws a plain reference line. When a freehand curve has been applied for the current split axis, that curve is shown instead of the slider line. """ if img is None or myo_idx is None or not mode: return None myo_u8 = P.channel_preview(img.data[int(myo_idx)]) cents = hc_cents if hc_cents is not None else np.zeros((0, 2)) if (mode.startswith("Split") and custom and custom.get("axis") == (axis or "Y").upper()): return _myo_preview_custom(img, myo_idx, custom["profile"], axis, cents, side) if mode.startswith("Split"): return _boundary_preview(myo_u8, boundary, axis, angle, curvature, cents, side) if mode.startswith("Inner"): return _boundary_preview(myo_u8, None, axis, angle, curvature, cents, None) return _boundary_preview(myo_u8, boundary, axis, angle, curvature, cents, None) def detect_cells(img, myo_idx, axis, angle, curvature): """Run Cellpose-SAM on the Myo7a channel, overlay the detected hair cells, and propose an IHC/OHC boundary + side.""" if img is None or myo_idx is None: return (None, None, gr.update(), gr.update(), "Load an image and select the Myo7a channel first.") if axis.upper() == "Z": return (None, None, gr.update(), gr.update(), "Hair-cell auto-detection proposes an in-plane (Y/X) boundary; " "switch the split axis to Y or X to use it.") try: myo_idx = int(myo_idx) det = P.detect_hair_cells(img.data[myo_idx], img.voxel) reg = P.auto_regions_from_cells(det["centroids"], det["mip_shape"], um_per_px=img.voxel[2]) myo_u8 = P.channel_preview(img.data[myo_idx]) if reg["single_row"]: # Do not force a meaningless split — keep the current boundary and # just show the detected cells with a clear warning. ov = _boundary_preview(myo_u8, reg["boundary_frac"], axis, angle, curvature, det["centroids"], None) status = ( f"Detected **{det['count']}** hair cells " f"(engine: {det['engine']}).\n\n" f"**No IHC/OHC split suggested.** {reg['reason']}\n\n" f"The **Whole field** row is the reliable result here. If you " f"know the anatomy, set the boundary by hand; otherwise leave " f"it and read the whole-field metrics.") return (det["centroids"], ov, gr.update(), gr.update(), status) # Clamp the suggested IHC side to a value the radio still offers (the # "bottom" option may have been removed for the Y axis), so the overlay, # the status text, the radio, and what Run uses can never disagree. side_val = reg["ihc_side"] valid = {v for _, v in _side_choices(axis)} clamped = side_val not in valid if clamped: side_val = _side_choices(axis)[0][1] ov = _boundary_preview(myo_u8, reg["boundary_frac"], axis, angle, curvature, det["centroids"], side_val) # side_val is the IHC side; the radio reads in OHC terms, so describe # the OHC side here (the opposite) to stay consistent. side_word = "bottom/right" if side_val == "low" else "top/left" status = ( f"Detected **{det['count']}** hair cells " f"(engine: {det['engine']}).\n\n" f"Suggested boundary **{reg['boundary_frac']:.3f}**, " f"OHC on the **{side_word}** side, confidence: " f"**{reg['confidence']}** ({reg['reason']}).\n\n" f"Note: this is a *starting suggestion* — detection is often " f"incomplete on dense fields. **Check the overlay** (cyan = IHC, " f"magenta = OHC, yellow = boundary) and adjust the boundary / tilt / " f"curvature before running.") if clamped: status += ("\n\n*(Detection favoured the opposite side, but this " "axis only offers 'OHC is ABOVE the line', so that side " "is kept.)*") return (det["centroids"], ov, gr.update(value=round(reg["boundary_frac"], 3)), gr.update(value=side_val), status) except Exception as e: return (None, None, gr.update(), gr.update(), f"Hair-cell detection failed:\n{e}\n{traceback.format_exc()}") def run_single(img, nf_idx, myo_idx, freq, analyze_mode, axis, side, boundary, angle, curvature, sensitivity, capture, min_fiber, prune_um, min_diam, max_diam, man_ihc, man_ohc, view, hc_cents, custom=None): blank = (None, None, None, None, None, None) if img is None: return blank + ("Please load an image first.",) if nf_idx is None: return blank + ("Please select the Neurofilament channel.",) if not analyze_mode: return blank + ("Please choose an analysis scope (Split IHC vs OHC, " "Inner hair cells only, or Whole field only).",) # Step 3 runs only after Step 2 — hair-cell detection must have been run # (its detected cells feed the per-region hair-cell counts / normalization). # hc_cents is None until Step 2 runs; it becomes an array (possibly empty) # afterwards. This gate does not touch the neurofilament tracing/skeleton. # Exception: a Z (depth) split is set by hand — Step 2 is an in-plane (Y/X) # detector and refuses on Z, so requiring it there would deadlock. if hc_cents is None and (axis or "Y").upper() != "Z": return blank + ("Please run **Step 2 — Detect hair cells (Cellpose-SAM)** " "first — its results feed into Step 3.",) try: nf_idx, myo_idx = int(nf_idx), int(myo_idx) split = analyze_mode.startswith("Split") ihc_only = analyze_mode.startswith("Inner") # whole field is all IHC trace = P.trace_neurites(img.data[nf_idx], img.voxel, sensitivity=float(sensitivity), prune_um=float(prune_um), capture=float(capture)) shape3d = trace.skeleton.shape z_note = "" if split and axis.upper() == "Z" and shape3d[0] < 2: # A depth split is impossible on a single-plane image — fall back to # whole-field rather than producing an empty region. split = False z_note = ("\n\n**Note:** a Z (depth) split needs a multi-plane " "stack; this image has a single plane, so only the " "whole-field result is shown.") # Freehand-drawn boundary curve, applied only when it matches the split # axis (else fall back to the position/tilt/curvature sliders). use_profile = (custom.get("profile") if (custom and axis.upper() in ("Y", "X") and custom.get("axis") == axis.upper()) else None) cents = hc_cents if (hc_cents is not None and len(hc_cents)) else None man_i, man_o = _manual(man_ihc), _manual(man_ohc) man_whole = ((man_i or 0) + (man_o or 0)) if (man_i or man_o) else None common = dict(min_fiber_um=float(min_fiber), min_diameter_um=float(min_diam), max_diameter_um=float(max_diam), radial_axis_name=axis) m_ihc = m_ohc = None ihc_roi = ohc_roi = None if ihc_only: # Inner-hair-cell-only image: the whole field is IHC territory, so no # OHC region is created (nothing is mis-assigned as OHC). m_ihc = P.compute_metrics(trace, "IHC region", hair_cell_centroids=cents, manual_hair_cells=man_i, **common) metrics = [m_ihc] primary = m_ihc else: primary = P.compute_metrics(trace, "Whole field", hair_cell_centroids=cents, manual_hair_cells=man_whole, **common) metrics = [primary] if split: ihc_roi, ohc_roi = P.make_region_masks( shape3d, float(boundary), ihc_side=side, axis=axis, angle_deg=float(angle), curvature=float(curvature), profile=use_profile) m_ihc = P.compute_metrics(trace, "IHC region", ihc_roi, hair_cell_centroids=cents, manual_hair_cells=man_i, **common) m_ohc = P.compute_metrics(trace, "OHC region", ohc_roi, hair_cell_centroids=cents, manual_hair_cells=man_o, **common) metrics += [m_ihc, m_ohc] # Renders. nf_mip = P.channel_preview(img.data[nf_idx]) skel_img = P.skeleton_image(trace.skeleton) overlay_img = P.overlay_on_original(nf_mip, trace.skeleton) views = {VIEW_OVERLAY: overlay_img, VIEW_SKELETON: skel_img, VIEW_ORIGINAL: nf_mip} main_img = views.get(view, overlay_img) region_img = None if split and use_profile is not None: # Colour by the custom masks; draw the freehand curve (no slider line). region_img = P.region_overlay(trace.skeleton, ihc_roi, ohc_roi, None, axis) _draw_curve_on(region_img, use_profile, axis) elif split: region_img = P.region_overlay(trace.skeleton, ihc_roi, ohc_roi, float(boundary), axis, float(angle), float(curvature)) elif ihc_only: # Everything is IHC — show the whole trace as one (cyan) region. ny, nx = shape3d[1], shape3d[2] region_img = P.region_overlay( trace.skeleton, np.ones((ny, nx), bool), np.zeros((ny, nx), bool), None, axis) myo_u8 = P.channel_preview(img.data[myo_idx]) if split and use_profile is not None: myo_img = _myo_preview_custom( img, myo_idx, use_profile, axis, cents if cents is not None else np.zeros((0, 2)), side) elif split: myo_img = _boundary_preview(myo_u8, boundary, axis, angle, curvature, cents, side) else: # IHC-only draws no boundary; whole-field keeps the reference line. myo_img = _boundary_preview(myo_u8, None if ihc_only else boundary, axis, angle, curvature, cents, None) stem = os.path.splitext(img.source_name)[0] skel_path = _save_png(skel_img, f"{stem}_skeleton.png") overlay_path = _save_png(overlay_img, f"{stem}_overlay.png") rows = [{"File": img.source_name, "Frequency region": freq, **m.as_row()} for m in metrics] df = pd.DataFrame(rows, columns=METRIC_COLS) xl_path = os.path.join(OUT_DIR, f"{stem}_quantification.xlsx") _build_excel(rows, xl_path) # Status text. hc_note = "" if split and (cents is not None or man_i or man_o): hc_note = (f" Hair cells: IHC={m_ihc.n_hair_cells}, " f"OHC={m_ohc.n_hair_cells}.") elif ihc_only and (cents is not None or man_i): hc_note = f" Hair cells (IHC): {m_ihc.n_hair_cells}." dir_note = "" if primary.n_radial >= 0: dir_note = (f" Direction: {primary.n_radial} radial (IHC→OHC) / " f"{primary.n_offaxis} off-axis " f"({primary.pct_radial:.0f}% radial).") warn = "" if split: one_side_no_fibers = ((m_ihc.total_length_um == 0) != (m_ohc.total_length_um == 0)) one_side_no_cells = (cents is not None and not (man_i or man_o) and ((m_ihc.n_hair_cells == 0) != (m_ohc.n_hair_cells == 0))) if one_side_no_fibers or one_side_no_cells: warn = ("\n\n**Warning:** this boundary does not separate two " "hair-cell rows — one region has no fibers or no hair " "cells of its type. The field likely shows a single " "hair-cell row, so choose **Inner hair cells only** or " "trust the **Whole field** row, or place the boundary on " "the tunnel of Corti only if both rows are visible.") scope = ("Split IHC/OHC" if split else "Inner hair cells only" if ihc_only else "Whole field only") field_label = "Inner-hair-cell field" if ihc_only else "Whole field" status = (f"Done ({scope}). Traced {int(trace.skeleton.sum())} skeleton " f"voxels. {field_label}: {primary.n_fibers} fibers / " f"{primary.total_length_um:.0f} µm.") if split: status += (f" IHC={m_ihc.n_fibers} fibers / " f"{m_ihc.total_length_um:.0f} µm, OHC={m_ohc.n_fibers} " f"fibers / {m_ohc.total_length_um:.0f} µm.") status += hc_note + dir_note + warn + z_note files = [skel_path, overlay_path, xl_path] return (main_img, region_img, myo_img, df, files, views, status) except Exception as e: return blank + (f"Error:\n{e}\n{traceback.format_exc()}",) # --------------------------------------------------------------------------- # # Batch flow # --------------------------------------------------------------------------- # def switch_batch_view(view, view_items): """Re-render the batch gallery in the chosen view without reprocessing.""" if not view_items: return None v = view if view in VIEW_CHOICES else VIEW_OVERLAY return [(it[v], it["caption"]) for it in view_items] def run_batch(files, analyze_mode, axis, side, boundary, angle, curvature, sensitivity, capture, min_fiber, prune_um, min_diam, max_diam, detect_hc, gallery_view, progress=gr.Progress()): if not files: return None, None, None, "Please upload one or more files.", None split = analyze_mode.startswith("Split") ihc_only = analyze_mode.startswith("Inner") all_rows, view_items, img_paths = [], [], [] log = [] common = dict(min_fiber_um=float(min_fiber), min_diameter_um=float(min_diam), max_diameter_um=float(max_diam), radial_axis_name=axis) for f in progress.tqdm(files, desc="Processing"): path = f if isinstance(f, str) else f.name name = os.path.basename(path) try: img = P.load_image(path) nf, myo = P.guess_channels(img) freq = P.detect_frequency(name) trace = P.trace_neurites(img.data[nf], img.voxel, sensitivity=float(sensitivity), prune_um=float(prune_um), capture=float(capture)) shape3d = trace.skeleton.shape cents = None hc_tag = "" do_split = split bfrac, cur_side = float(boundary), side ang, curv = float(angle), float(curvature) if do_split and axis.upper() == "Z" and shape3d[0] < 2: do_split = False # depth split needs >= 2 planes hc_tag = " | single plane, whole-field only (Z split N/A)" if (split or ihc_only) and detect_hc and axis.upper() != "Z": det = P.detect_hair_cells(img.data[myo], img.voxel) cents = det["centroids"] if len(det["centroids"]) else None if do_split: reg = P.auto_regions_from_cells(det["centroids"], det["mip_shape"], um_per_px=img.voxel[2]) if reg["single_row"]: do_split = False hc_tag = (f" | {det['count']} hair cells " f"({det['engine']}) — single row, " "whole-field only") else: bfrac, cur_side = reg["boundary_frac"], reg["ihc_side"] ang, curv = 0.0, 0.0 # auto boundary is a straight line hc_tag = (f" | {det['count']} hair cells " f"({det['engine']}), auto-boundary " f"conf={reg['confidence']}") else: hc_tag = f" | {det['count']} hair cells ({det['engine']})" if ihc_only: # Whole field is IHC; no OHC region is created. metrics = [P.compute_metrics(trace, "IHC region", hair_cell_centroids=cents, **common)] else: metrics = [P.compute_metrics(trace, "Whole field", hair_cell_centroids=cents, **common)] if do_split: ihc_roi, ohc_roi = P.make_region_masks( shape3d, bfrac, ihc_side=cur_side, axis=axis, angle_deg=ang, curvature=curv) metrics += [ P.compute_metrics(trace, "IHC region", ihc_roi, hair_cell_centroids=cents, **common), P.compute_metrics(trace, "OHC region", ohc_roi, hair_cell_centroids=cents, **common)] for m in metrics: all_rows.append({"File": name, "Frequency region": freq, **m.as_row()}) # Save BOTH the black-on-white skeleton and the trace-on-original # overlay; the gallery shows whichever the user picked, and both go # into the ZIP so nothing is lost. nf_mip = P.channel_preview(img.data[nf]) skel_img = P.skeleton_image(trace.skeleton) overlay_img = P.overlay_on_original(nf_mip, trace.skeleton) stem = os.path.splitext(name)[0] sp = _save_png(skel_img, f"{stem}_skeleton.png") op = _save_png(overlay_img, f"{stem}_overlay.png") img_paths += [sp, op] view_items.append({VIEW_OVERLAY: overlay_img, VIEW_SKELETON: skel_img, VIEW_ORIGINAL: nf_mip, "caption": f"{name} ({freq})"}) log.append(f"{name}: {freq}{hc_tag}") except Exception as e: log.append(f"{name}: FAILED — {e}") if not all_rows: return None, None, [], "No files processed.\n" + "\n".join(log), [] xl_path = os.path.join(OUT_DIR, "batch_quantification.xlsx") df = _build_excel(all_rows, xl_path) zip_path = os.path.join(OUT_DIR, "batch_images.zip") with zipfile.ZipFile(zip_path, "w") as z: for p in img_paths: z.write(p, os.path.basename(p)) z.write(xl_path, os.path.basename(xl_path)) return (df, [xl_path, zip_path], switch_batch_view(gallery_view, view_items), "\n".join(log), view_items) # --------------------------------------------------------------------------- # # UI # --------------------------------------------------------------------------- # HEADER_HTML = """

Neuron Quantification using AI

Iman Sabir Ezzat, Randa K Ismail, Ayden Chavez, Marisa Zallocchi, PhD, Steven Fernandes, PhD
""" INTRO = """ Trace auditory-nerve fibers in confocal z-stacks and quantify them **per frequency region**, separating **IHC-innervating** from **OHC-innervating** fibers using the Myo7a hair-cell channel. **Channels expected:** *Neurofilament* (traces the neuron) and *Myo7a* (hair cells — reference to split IHC vs OHC). Many dyes are recognised automatically (405 / 488 / 514 / 555 / 568 / 594 / 633 / 647 / ATTO / Cy…); every channel in the file is selectable. IHCs form a single row, OHCs form three rows, so the Myo7a band is used to place the IHC/OHC boundary — which you can move, **tilt**, **curve**, or switch to an **X** or **Z** split axis. **Input:** Zeiss `.czi` z-stacks or generic `.tif/.tiff` stacks. """ # Professional light theme: white surfaces, near-black text, a restrained # slate/charcoal accent (no purple), clean system font. THEME = gr.themes.Base( primary_hue=gr.themes.colors.slate, secondary_hue=gr.themes.colors.slate, neutral_hue=gr.themes.colors.gray, font=["system-ui", "-apple-system", "Segoe UI", "Roboto", "Helvetica", "Arial", "sans-serif"], font_mono=["ui-monospace", "SFMono-Regular", "Menlo", "Consolas", "monospace"], ).set( body_background_fill="#ffffff", body_text_color="#1a1a1a", body_text_color_subdued="#4b5563", background_fill_primary="#ffffff", background_fill_secondary="#f7f7f8", block_background_fill="#ffffff", block_border_color="#e5e7eb", block_border_width="1px", block_label_text_color="#1a1a1a", block_title_text_color="#1a1a1a", border_color_primary="#e5e7eb", button_large_radius="8px", button_small_radius="8px", button_primary_background_fill="#1f2937", button_primary_background_fill_hover="#111827", button_primary_text_color="#ffffff", button_primary_border_color="#1f2937", button_secondary_background_fill="#ffffff", button_secondary_background_fill_hover="#f3f4f6", button_secondary_text_color="#1a1a1a", button_secondary_border_color="#d1d5db", input_background_fill="#ffffff", input_border_color="#c0c5cc", input_border_width="1px", # visible outline so dropdowns read as dropdowns input_border_color_focus="#1f2937", slider_color="#374151", checkbox_background_color="#ffffff", checkbox_background_color_selected="#1f2937", checkbox_border_width="2px", # visible ring so radios/checkboxes show checkbox_border_color="#9ca3af", checkbox_border_color_hover="#6b7280", checkbox_border_color_selected="#1f2937", checkbox_label_text_color="#1a1a1a", # Pin the *dark-mode* tokens to the same light values so the app keeps its # professional white look even when the browser / OS is in dark mode. body_background_fill_dark="#ffffff", body_text_color_dark="#1a1a1a", body_text_color_subdued_dark="#4b5563", background_fill_primary_dark="#ffffff", background_fill_secondary_dark="#f7f7f8", block_background_fill_dark="#ffffff", block_border_color_dark="#e5e7eb", block_label_background_fill_dark="#ffffff", block_label_text_color_dark="#1a1a1a", block_title_background_fill_dark="#ffffff", block_title_text_color_dark="#1a1a1a", panel_background_fill_dark="#ffffff", border_color_primary_dark="#e5e7eb", button_primary_background_fill_dark="#1f2937", button_primary_background_fill_hover_dark="#111827", button_primary_text_color_dark="#ffffff", button_primary_border_color_dark="#1f2937", button_secondary_background_fill_dark="#ffffff", button_secondary_background_fill_hover_dark="#f3f4f6", button_secondary_text_color_dark="#1a1a1a", button_secondary_border_color_dark="#d1d5db", input_background_fill_dark="#ffffff", input_border_color_dark="#c0c5cc", input_border_color_focus_dark="#1f2937", checkbox_background_color_dark="#ffffff", checkbox_background_color_selected_dark="#1f2937", checkbox_border_color_dark="#9ca3af", checkbox_border_color_hover_dark="#6b7280", checkbox_border_color_selected_dark="#1f2937", checkbox_label_text_color_dark="#1a1a1a", ) # Keep the layout tidy, and force the white professional look even if the # browser/OS is in dark mode (Gradio adds a `.dark` class we neutralise here). CSS = """ .gradio-container { max-width: 1200px !important; margin: 0 auto !important; } .prose h1 { font-weight: 650; letter-spacing: -0.01em; } :root, .dark { color-scheme: light; --body-background-fill: #ffffff; --background-fill-primary: #ffffff; --background-fill-secondary: #f7f7f8; --block-background-fill: #ffffff; --block-label-background-fill: #ffffff; --block-title-background-fill: #ffffff; --panel-background-fill: #ffffff; --body-text-color: #1a1a1a; --body-text-color-subdued: #4b5563; --block-label-text-color: #1a1a1a; --block-title-text-color: #1a1a1a; --block-info-text-color: #4b5563; --border-color-primary: #e5e7eb; --block-border-color: #e5e7eb; --input-background-fill: #ffffff; --input-border-color: #c0c5cc; --table-odd-background-fill: #ffffff; --table-even-background-fill: #f7f7f8; --code-background-fill: #f1f3f5; --neutral-950: #1a1a1a; } body, gradio-app, .gradio-container, .dark, .dark .gradio-container { background: #ffffff !important; color: #1a1a1a !important; } .prose code, .prose pre, .gradio-container code, .gradio-container pre { background: #f1f3f5 !important; color: #1a1a1a !important; } """ with gr.Blocks(title="Neuron Quantification using AI", theme=THEME, css=CSS) as demo: gr.HTML(HEADER_HTML) gr.Markdown(INTRO) with gr.Tab("Single image (interactive)"): img_state = gr.State() hc_state = gr.State() # detected hair-cell centroids (Nx2) view_state = gr.State() # {view name -> image} for the toggle custom_state = gr.State() # {"axis","profile"} freehand boundary with gr.Row(): with gr.Column(scale=1): file_in = gr.File(label="Upload CZI or TIFF", file_types=[".czi", ".tif", ".tiff"], type="filepath") load_btn = gr.Button("Step 1 — Load & preview", variant="primary") nf_dd = gr.Dropdown(label="Neurofilament channel", choices=[]) myo_dd = gr.Dropdown(label="Myo7a channel", choices=[]) freq_dd = gr.Dropdown(label="Frequency region", choices=P.FREQ_CHOICES, value="Other / unknown") mode_dd = gr.Radio( ["Split IHC vs OHC", "Inner hair cells only", "Whole field only"], value=None, label="Analysis scope — pick one to reveal the Myo7a (MIP) + " "boundary preview. Choose 'Inner hair cells only' when " "the image has no outer hair cells (no OHC region is " "created); the split controls below are then ignored") gr.Markdown("**IHC / OHC region split** (Myo7a-guided) — used " "only for 'Split IHC vs OHC'") axis_dd = gr.Radio(["Y", "X", "Z"], value="Y", label="Split axis (Y = radial, usual · " "Z = by depth)") side_dd = gr.Radio(_side_choices("Y"), value="high", label="Which side is OHC?") boundary_sl = gr.Slider(0.0, 1.0, value=0.5, step=0.005, label="Boundary position (fraction " "along split axis)") angle_sl = gr.Slider(-45.0, 45.0, value=0.0, step=1.0, label="Boundary tilt (degrees) — for a " "rotated organ") curve_sl = gr.Slider(-0.5, 0.5, value=0.0, step=0.02, label="Boundary curvature (bow, for a " "curving cochlea · 0 = straight)") with gr.Accordion("Draw a custom curved boundary (freehand)", open=False): gr.Markdown( "When the boundary isn't a straight or simply-bowed " "line: press **Load Myo7a here**, draw the boundary " "freehand across the whole width with the brush, then " "**Apply drawn boundary**. The analysis and overlays use " "your drawn curve (for the current split axis, Y or X). " "**Clear drawn boundary** reverts to the sliders " "above.") draw_editor = gr.ImageEditor( label="Draw boundary (freehand)", type="numpy", height=320, image_mode="RGB", brush=gr.Brush(default_size=5, colors=["#FFEE00"], color_mode="fixed")) with gr.Row(): draw_load_btn = gr.Button("Load Myo7a here") draw_apply_btn = gr.Button("Apply drawn boundary", variant="primary") draw_clear_btn = gr.Button("Clear drawn boundary") detect_btn = gr.Button( "Step 2 — Detect hair cells (Cellpose-SAM)", variant="primary") gr.Markdown( "Detection is a *visual assist*: it marks hair cells " "and proposes a starting boundary. Confirm/adjust against " "the Myo7a overlay — it does not replace your judgement. " "Detection can over-mark the tunnel of Corti, so you can " "override the counts below.") with gr.Accordion("Hair-cell count override (for normalization)", open=False): gr.Markdown( "Enter trusted hair-cell counts to normalize by " "(0 = use the detected count). Fibers-per-hair-cell and " "length-per-hair-cell are reported per region.") man_ihc = gr.Number(0, label="IHC hair cells (manual)", precision=0) man_ohc = gr.Number(0, label="OHC hair cells (manual)", precision=0) gr.Markdown("**Tracing**") sens_sl = gr.Slider(0.5, 1.5, value=1.0, step=0.05, label="Sensitivity (higher = capture more / " "thinner fibers)") grow_sl = gr.Slider(0.20, 1.0, value=0.40, step=0.05, label="Faint-fiber capture — recovers dim " "fibers (top-hat) to raise total length " "toward Imaris (lower = capture more · " "1.0 = legacy Otsu, off)") minfib_sl = gr.Slider(0.0, 20.0, value=5.0, step=0.5, label="Min fiber length to count (µm)") mindiam_sl = gr.Slider(0.0, 3.0, value=0.0, step=0.05, label="Min fiber diameter (µm, 0 = off)") maxdiam_sl = gr.Slider(0.0, 10.0, value=0.0, step=0.1, label="Max fiber diameter (µm, 0 = off)") prune_sl = gr.Slider(0.0, 15.0, value=3.0, step=0.5, label="Clean-up: remove fragments/spurs " "shorter than (µm) (0 = off)") run_btn = gr.Button("Step 3 — Run analysis", variant="primary") with gr.Column(scale=2): status = gr.Markdown() with gr.Row(): nf_prev = gr.Image(label="Neurofilament (MIP)", height=220) myo_prev = gr.Image(label="Myo7a (MIP) + boundary", height=220) view_dd = gr.Radio(VIEW_CHOICES, value=VIEW_OVERLAY, label="Main image view (toggle trace on/off " "over the original)") main_out = gr.Image(label="Traced neurons", height=320) region_out = gr.Image(label="Region overlay (cyan = IHC, " "magenta = OHC)", height=300) table = gr.Dataframe(label="Quantification", wrap=True) files_out = gr.Files(label="Downloads (skeleton + overlay PNG " "+ Excel)") load_btn.click(load_and_preview, [file_in], [img_state, nf_dd, myo_dd, freq_dd, boundary_sl, nf_prev, myo_prev, main_out, hc_state, status] ).then(lambda: (None, None), None, [view_state, custom_state]) # Adapt the IHC-side wording to the chosen split axis. axis_dd.change(update_side_labels, [axis_dd, side_dd], side_dd) # Live boundary preview: choosing a scope reveals it, and the split # controls keep it up to date (it stays hidden until a scope is picked). # A drawn curve, once applied, is shown in place of the slider line. for comp in (mode_dd, boundary_sl, myo_dd, axis_dd, side_dd, angle_sl, curve_sl): comp.change(refresh_boundary_preview, [img_state, myo_dd, boundary_sl, axis_dd, hc_state, side_dd, angle_sl, curve_sl, mode_dd, custom_state], myo_prev) # Freehand boundary drawing. draw_load_btn.click(load_draw_background, [img_state, myo_dd], draw_editor) draw_apply_btn.click(apply_drawn_boundary, [draw_editor, axis_dd, mode_dd, img_state, myo_dd, hc_state, side_dd], [custom_state, myo_prev, status]) draw_clear_btn.click(clear_drawn_boundary, [img_state, myo_dd, boundary_sl, axis_dd, hc_state, side_dd, angle_sl, curve_sl, mode_dd], [custom_state, myo_prev, status]) detect_btn.click(detect_cells, [img_state, myo_dd, axis_dd, angle_sl, curve_sl], [hc_state, myo_prev, boundary_sl, side_dd, status]) view_dd.change(switch_view, [view_dd, view_state], main_out) run_btn.click(run_single, [img_state, nf_dd, myo_dd, freq_dd, mode_dd, axis_dd, side_dd, boundary_sl, angle_sl, curve_sl, sens_sl, grow_sl, minfib_sl, prune_sl, mindiam_sl, maxdiam_sl, man_ihc, man_ohc, view_dd, hc_state, custom_state], [main_out, region_out, myo_prev, table, files_out, view_state, status]) with gr.Tab("Batch (multiple images)"): gr.Markdown( "Upload several z-stacks (e.g. all frequency regions of one " "cochlea). Each is auto-traced and results are combined into one " "Excel workbook organized by frequency region. Set the IHC/OHC " "boundary once below (used for every file), or let hair-cell " "detection place it per image.") with gr.Row(): with gr.Column(scale=1): batch_files = gr.File(label="Upload CZI/TIFF files", file_count="multiple", file_types=[".czi", ".tif", ".tiff"], type="filepath") b_mode = gr.Radio(["Split IHC vs OHC", "Inner hair cells only", "Whole field only"], value="Split IHC vs OHC", label="Analysis scope ('Inner hair cells " "only' = no OHC region)") b_axis = gr.Radio(["Y", "X", "Z"], value="Y", label="Split axis") b_side = gr.Radio(_side_choices("Y"), value="high", label="Which side is OHC?") b_boundary = gr.Slider(0.0, 1.0, value=0.5, step=0.005, label="Boundary position (used when " "detection is off)") b_angle = gr.Slider(-45.0, 45.0, value=0.0, step=1.0, label="Boundary tilt (degrees)") b_curve = gr.Slider(-0.5, 0.5, value=0.0, step=0.02, label="Boundary curvature (bow)") b_sens = gr.Slider(0.5, 1.5, value=1.0, step=0.05, label="Sensitivity") b_grow = gr.Slider(0.20, 1.0, value=0.40, step=0.05, label="Faint-fiber capture (lower = more " "length · 1.0 = legacy Otsu)") b_minfib = gr.Slider(0.0, 20.0, value=5.0, step=0.5, label="Min fiber length (µm)") b_mindiam = gr.Slider(0.0, 3.0, value=0.0, step=0.05, label="Min fiber diameter (µm, 0 = off)") b_maxdiam = gr.Slider(0.0, 10.0, value=0.0, step=0.1, label="Max fiber diameter (µm, 0 = off)") b_prune = gr.Slider(0.0, 15.0, value=3.0, step=0.5, label="Clean-up: remove fragments/spurs " "shorter than (µm)") b_detect = gr.Checkbox( value=False, label="Detect hair cells & auto-place boundary per image " "(Cellpose-SAM, adds ~15–20 s/image)") b_view = gr.Radio( VIEW_CHOICES, value=VIEW_OVERLAY, label="Gallery view (toggle trace on/off over the original)") batch_btn = gr.Button("Run batch", variant="primary") with gr.Column(scale=2): batch_log = gr.Textbox(label="Log", lines=6) batch_table = gr.Dataframe(label="Combined quantification", wrap=True) batch_files_out = gr.Files( label="Downloads (Excel + ZIP of skeleton & overlay PNGs)") batch_gallery = gr.Gallery(label="Traces", columns=3, height=400) batch_view_state = gr.State() b_axis.change(update_side_labels, [b_axis, b_side], b_side) batch_btn.click(run_batch, [batch_files, b_mode, b_axis, b_side, b_boundary, b_angle, b_curve, b_sens, b_grow, b_minfib, b_prune, b_mindiam, b_maxdiam, b_detect, b_view], [batch_table, batch_files_out, batch_gallery, batch_log, batch_view_state]) b_view.change(switch_batch_view, [b_view, batch_view_state], batch_gallery) gr.Markdown( "---\n*Method:* the Neurofilament channel is smoothed, thresholded " "(Otsu, scaled by the sensitivity control) and skeletonised in 3D. Each " "fiber is a connected skeleton component kept only if it is longer than " "the minimum length **and** its mean diameter falls in the chosen " "min/max diameter band. Short terminal spurs and small isolated " "fragments below the clean-up length are pruned (junction points " "preserved) so the fiber lines stay clean.\n\n" "*Diameter* is measured from the **3D Euclidean distance transform** of " "the segmented Neurofilament mask: at every skeleton voxel the distance " "to the nearest background voxel (in microns, using the real voxel " "spacing) is the local radius, so diameter = 2 × that distance. The " "reported mean/median diameter is taken over all kept skeleton voxels in " "the region.\n\n" "*Fiber direction* classifies each fiber by its principal (PCA) axis " "relative to the IHC→OHC (radial) axis: fibers aligned with it are " "**radial** (the expected innervation direction) and the rest are " "**off-axis** — the reviewers' 'majority correct / minority " "misdirected' split.\n\n" "*Normalization:* enter trusted IHC/OHC hair-cell counts (or use " "detection) and the app reports **fibers per hair cell** and **length " "per hair cell**.\n\n" "*Hair-cell detection* (Step 2) always runs **Cellpose-SAM** (the " "cpsam model shipped with cellpose 4.x; a fine-tuned model is used " "automatically if provided at models/hair_cell_cpsam) on the Myo7a " "max-projection. On dense fields this detection is often incomplete or " "over-marks the tunnel of Corti, so it is a **visual assist**: the " "numbers you trust come from the deterministic pipeline and the manual " "count override, and the boundary remains yours to set (position, tilt, " "curvature, or axis). Detection quality improves markedly on a GPU " "Space.") if __name__ == "__main__": demo.launch()