Upload 317 files

README.md

@@ -16,15 +16,21 @@ A variant of the neuron tracer that applies the **MedCLIPSeg** approach
 (Koleilat et al., *Probabilistic Vision–Language Adaptation for Data-Efficient
 and Generalizable Medical Image Segmentation*, CVPR 2026): a **text prompt**
 describes the target ("nerve fibers …") and a frozen **CLIP** backbone produces a
-**fiber-probability map**, which is then traced into a skeleton.
+**fiber-probability map** that acts as a *semantic region prior* — telling the
+tracer **where** the nerve fibers are.
 
 ## What the app does
 - **Neurofilament (MIP):** the maximum-intensity projection of the input channel.
-- **Tracked skeleton (white on black):** the MedCLIP fiber-probability map is
-  thresholded, skeletonised and pruned into a fiber skeleton.
-- **Total trace length (µm):** the spacing-aware total length of that MedCLIP
-  skeleton, reported above the images. The text prompt and probability threshold
-  are editable.
+- **Tracked skeleton (white on black):** MedCLIP-gated tubular tracing. The
+  MedCLIP probability map (WHERE) seeds a **Sato tubeness** ridge filter run on
+  the raw neurofilament signal (the thin fiber GEOMETRY); seeds grow along the
+  raw ridge across the whole field, so the long radial fibers below the band are
+  recovered, and a width gate removes any residual medial-axis mesh. This traces
+  clean thin fibers instead of skeletonising the coarse probability blob (which
+  produced a "cracked-mud" mesh).
+- **Total trace length (µm):** the spacing-aware total length of that skeleton,
+  reported above the images. The text prompt and the MedCLIP **semantic gate**
+  (region-prior threshold) are editable.
 
 ## Few-shot training (included)
 A runnable, few-shot version of the MedCLIPSeg approach — **frozen CLIP encoders +

app.py

@@ -26,10 +26,11 @@ def analyze(file_obj, fg_prompt, threshold):
         if fg_prompt and fg_prompt.strip():
             MC.FG_PROMPTS = [p.strip() for p in fg_prompt.split("|") if p.strip()]
 
-        # --- MedCLIP text-prompted fiber-probability map -> tracked skeleton ---
+        # --- MedCLIP text-prompted region prior -> raw-signal fiber geometry ---
         prob, _ = MC.segment_best(nf_mip)
-        trace = P.trace_from_probability(prob, img.voxel,
-                                         threshold=float(threshold), prune_um=3.0)
+        trace = P.trace_fibers_medclip_gated(nf_mip, prob, img.voxel,
+                                             gate_threshold=float(threshold),
+                                             prune_um=3.0)
         m = P.compute_metrics(trace, "MedCLIP", min_fiber_um=5.0)
         skel_img = P.skeleton_image(trace.skeleton, dilate=1)   # white on black
 
@@ -37,8 +38,10 @@ def analyze(file_obj, fg_prompt, threshold):
                  else "zero-shot CLIP prior")
         status = (f"**MedCLIP total trace length: "
                   f"{m.total_length_um:,.1f} µm.**\n\n"
-                  f"Tracked from the MedCLIP fiber-probability map "
-                  f"({model}, text-prompted) at threshold {float(threshold):.2f}. "
+                  f"Fiber geometry traced from the raw neurofilament signal "
+                  f"(Sato tubeness), gated by the MedCLIP region prior "
+                  f"({model}, text-prompted) at semantic gate "
+                  f"{float(threshold):.2f}. "
                   f"Prompt(s): *{', '.join(MC.FG_PROMPTS)}*.")
         return nf_mip, skel_img, status
     except Exception as e:  # noqa: BLE001
@@ -101,8 +104,9 @@ with gr.Blocks(title="Neuron Quantification using AI — MedCLIPSeg",
                 label="Text prompt(s) for the target (separate with | )",
                 value="a fluorescence microscopy image of nerve fibers | "
                       "neurofilament nerve fibers and axons")
-            thresh = gr.Slider(0.1, 0.9, value=0.5, step=0.05,
-                               label="Fiber-probability threshold")
+            thresh = gr.Slider(0.2, 0.8, value=0.4, step=0.05,
+                               label="MedCLIP semantic gate (region prior) — "
+                                     "lower = include more of the band")
             btn = gr.Button("Analyze", variant="primary")
         with gr.Column(scale=2):
             status = gr.Markdown()

processing.py

@@ -21,8 +21,9 @@ from typing import Optional
 
 import numpy as np
 from scipy import ndimage as ndi
-from skimage.filters import gaussian, threshold_otsu
-from skimage.morphology import remove_small_objects, skeletonize
+from skimage.filters import gaussian, threshold_otsu, sato
+from skimage.morphology import (remove_small_objects, remove_small_holes,
+                                skeletonize)
 from skan import Skeleton, summarize
 
 
@@ -598,6 +599,130 @@ def trace_from_probability(prob: np.ndarray, voxel: tuple,
     return TraceResult(mask[None], skel, dist[None], voxel)
 
 
+def _safe_otsu(values: np.ndarray, fallback: float) -> float:
+    """threshold_otsu that never raises on constant / tiny inputs."""
+    v = np.asarray(values).ravel()
+    v = v[np.isfinite(v)]
+    if v.size < 8 or np.unique(v).size < 2:
+        return float(fallback)
+    try:
+        return float(threshold_otsu(v))
+    except Exception:
+        return float(fallback)
+
+
+def trace_fibers_medclip_gated(nf_mip: np.ndarray, prob: np.ndarray, voxel: tuple,
+                               gate_threshold: float = 0.4,
+                               fiber_radii_um: tuple = (0.2, 0.35, 0.6),
+                               gate_dilate_um: float = 4.0,
+                               max_fiber_half_width_um: float = 1.5,
+                               bright_floor: float = 0.5,
+                               min_object_px: int = 32,
+                               prune_um: float = 3.0) -> TraceResult:
+    """MedCLIP-gated tubular fiber tracer.
+
+    MedCLIP ``prob`` is used ONLY as a semantic region prior (WHERE fibers are).
+    The thin fiber GEOMETRY comes from a multiscale Sato tubeness response on the
+    RAW neurofilament MIP, so ``skeletonize`` acts on thin ridges (clean lines),
+    never on the coarse probability slab (which would yield a medial-axis
+    honeycomb).  Seeds are placed where strong tubeness meets the MedCLIP region;
+    those seeds are then grown along the raw ridge over the WHOLE image, so the
+    long radial (vertical) spiral-ganglion fibers descending below the band -
+    where ``prob`` is ~0 - are recovered.  A final EDT half-width gate deletes any
+    residual wide medial-axis ridge.  Returns a Z=1 ``TraceResult`` so
+    ``compute_metrics`` / ``skeleton_image`` / the length path are unchanged.
+
+    Parameters
+    ----------
+    nf_mip : 2D uint8/float RAW neurofilament MIP (``channel_preview`` output).
+    prob   : 2D MedCLIP fiber-probability in [0, 1] at MIP resolution.
+    voxel  : (dz, dy, dx) in microns.
+    gate_threshold : MedCLIP prob level defining the semantic region (UI slider).
+    """
+    dz, dy, dx = voxel
+    px = float(np.mean([dy, dx]))                          # microns / pixel
+    raw = np.asarray(nf_mip, np.float32)
+    prob = np.asarray(prob, np.float32)
+    if prob.shape != raw.shape:                            # robust to any mismatch
+        from skimage.transform import resize
+        prob = resize(prob, raw.shape, order=1, preserve_range=True).astype(np.float32)
+
+    def _empty():
+        z = np.zeros((1,) + raw.shape, bool)
+        return TraceResult(z, z.copy(), np.zeros(z.shape, np.float32), voxel)
+
+    # 1. Raw intensity, robustly normalised (reuse existing helper).
+    I = _norm(raw)
+
+    # 2. Multiscale tubeness on the RAW signal -> thin, orientation-invariant
+    #    ridge response (present for vertical radials as much as band fibers).
+    sigmas = sorted({max(1.0, round(float(r) / max(dx, 1e-6), 2))
+                     for r in fiber_radii_um})
+    T = sato(I, sigmas=sigmas, black_ridges=False).astype(np.float32)
+    T = np.nan_to_num(T, nan=0.0, posinf=0.0, neginf=0.0)
+    thi = np.percentile(T, 99.5)
+    T = np.clip(T / (thi + 1e-6), 0.0, 1.0)
+
+    # 3. MedCLIP semantic region (WHERE) -> generous, dilated gate.
+    region = prob >= float(gate_threshold)
+    if not region.any():
+        region = prob >= np.percentile(prob, 90.0)         # never empty
+    if region.any():
+        region = remove_small_objects(region, int(min_object_px))
+        region = ndi.binary_closing(region, structure=np.ones((3, 3), bool))
+    if not region.any():
+        return _empty()
+    iters = max(1, int(round(float(gate_dilate_um) / px)))
+    gate = ndi.binary_dilation(region, iterations=iters)
+
+    # 4. Weak brightness floor (rejects pure-noise ridges; keeps dim radials).
+    ipos = I[I > 0]
+    bright = I >= float(bright_floor) * _safe_otsu(ipos, 0.15)
+
+    # 5. Seed = strong tubeness INSIDE the semantic gate.
+    thr_hi = _safe_otsu(T[gate & (T > 0)], 0.3)
+    thr_hi = float(np.clip(thr_hi, 0.02, 0.99))
+    thr_lo = 0.5 * thr_hi
+    core = (T >= thr_hi) & gate & bright
+
+    # 6. Growth domain = thin raw ridges over the WHOLE image (radial recovery).
+    permissive = (T >= thr_lo) & bright
+
+    if core.any():
+        # 8-connected reconstruction: follow each in-band seed down its fiber,
+        # out of the band, recovering the long vertical radials.
+        fiber = ndi.binary_propagation(core, mask=permissive,
+                                       structure=np.ones((3, 3), bool))
+    else:
+        # Empty-seed fallback: thin tubeness inside the gate.  NEVER skeletonise
+        # the prob blob (that is the original honeycomb bug).
+        fiber = (T >= thr_hi) & gate & bright
+
+    fiber = remove_small_objects(fiber, int(min_object_px))
+    if not fiber.any():
+        return _empty()
+    # Bridge sub-pixel gaps without fusing adjacent parallel fibers.
+    fiber = remove_small_holes(fiber, area_threshold=4)
+
+    # 7. Skeletonise the THIN mask -> clean single-pixel centrelines.
+    dist = ndi.distance_transform_edt(fiber, sampling=(dy, dx)).astype(np.float32)
+    skel = skeletonize(fiber)
+
+    # 8. Artifact (filled-slab) rejection: drop skeleton pixels whose local
+    #    half-width exceeds a real fiber -> removes any residual medial-axis mesh
+    #    in the dense band, keyed on WIDTH (not orientation), so thin verticals
+    #    are never touched.
+    max_hw = max(px, float(max_fiber_half_width_um))
+    skel = skel & (dist <= max_hw)
+    if not skel.any():
+        return _empty()
+
+    skel = skel[None]
+    if prune_um and prune_um > 0:
+        skel = prune_skeleton(skel, voxel, spur_um=float(prune_um))
+    return TraceResult(fiber[None], skel, dist[None], voxel)
+
+
 # --------------------------------------------------------------------------- #
 #  Region definition (IHC vs OHC) from Myo7a
 # --------------------------------------------------------------------------- #