Add soft IoU validation to test_fire_2d_angle.py

Adds class TestSoftIoU with three test tiers that measure spatial
accuracy of the FIRE fiber extraction pipeline by comparing 1-pixel-thick
skeleton images.

Soft IoU metric
---------------
Both skeletons (ground truth and extracted) are first smoothed with a
Gaussian at sigma=5 (≈10-px FWHM) and rescaled to [0, 1].  The IoU is
then computed as:

    IoU = (m1 * m2).sum() / (m1² + m2² - m1 * m2).sum()

Operating at ~10-px FWHM measures agreement at fiber scale rather than
pixel scale, tolerating the sub-pixel registration differences that are
normal for a distance-transform medial-axis algorithm.

Thresholds: 0.70 for synthetic images (known GT), 0.60 for MATLAB
comparison (real images, higher variation expected).

Test tiers
----------
A. test_soft_iou_synthetic (seed=42, always runs)
   - Generates 8 straight-line Gaussian-profile fibers on a low background
   - Records the exact 1-px skeleton as ground truth at generation time
   - Asserts soft IoU > 0.70, total extracted length > 20% of GT length,
     and angle std > 0.3 rad (confirms angle computation is non-degenerate)
   - Uses thresh_im=0.2 (fractional) to cleanly separate fiber signal
     (peak ~180) from background (mean ~10); thresh_im2=5 (absolute) would
     include virtually all background pixels and produce hundreds of spurious
     short zigzag fibers

B. test_soft_iou_multiple_seeds (seeds 0, 7, 99)
   - Same structure as A; guards against a lucky pass on a single seed

C. test_soft_iou_matlab_vs_python (@pytest.mark.matlab, skips if absent)
   - Rasterizes MATLAB Xa/Fa as ground truth, compares to Python Xf/Ff
   - Also checks total length within ±40% and mean |angle| within 10°

Supporting additions
--------------------
- _smooth_mask, _soft_iou, _rasterize_fibers, _make_synthetic_fiber_image
  helper functions (pattern adapted from tme_quant/tests/test_ctfire.py)
- _default_fire_params() reads JSON params for real-image tests
- _synthetic_fire_params() overrides thresh_im=0.2 for synthetic tests
- __main__ block replaced with standalone demo (saves overlay PNG via
  plot_fiber_overlay from src/ctfire_py/test_fire_2d.py)
- pyproject.toml: register 'matlab' pytest marker

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Author: yuminguw

pyproject.toml

@@ -52,3 +52,6 @@ napari_curvealign = ["napari.yaml"]
 testpaths = [
   "tests",
 ]
+markers = [
+  "matlab: tests requiring MATLAB reference .mat files (skip if absent)",
+]

tests/test_fire_2d_angle.py

@@ -186,6 +186,150 @@ def load_matlab_reference(mat_file_path):
     pytest.skip("Neither h5py nor scipy.io available for loading MATLAB files")
 
 
+# ============================================================================
+# Soft IoU helpers
+# ============================================================================
+
+# Gaussian sigma for smoothing 1-px skeletons before computing soft IoU.
+# sigma=5 → ~10-px FWHM, which measures agreement at fiber scale rather than
+# pixel scale.  Adapted from tme_quant/tests/test_ctfire.py.
+_SMOOTH_SIGMA = 5.0
+
+SOFT_IOU_THRESHOLD_SYNTHETIC = 0.70  # extracted skeleton vs known GT skeleton
+SOFT_IOU_THRESHOLD_MATLAB    = 0.60  # Python centerlines vs MATLAB centerlines
+
+
+def _smooth_mask(mask, smooth_sigma=_SMOOTH_SIGMA):
+    """Gaussian-smooth a binary mask and rescale to [0, 1]."""
+    from skimage import exposure, filters
+    mask = mask.astype(np.float32)
+    mask = exposure.rescale_intensity(mask, out_range=(0.0, 1.0))
+    density = filters.gaussian(mask, sigma=smooth_sigma, preserve_range=False)
+    return exposure.rescale_intensity(density, out_range=(0.0, 1.0)).astype(np.float32)
+
+
+def _soft_iou(mask_1, mask_2, beta=1e-3):
+    """
+    Soft IoU on two float masks already smoothed to [0, 1].
+
+    Formula: (m1·m2).sum() / (m1²+m2²-m1·m2).sum()
+    beta prevents 0/0 when both masks are empty.
+    """
+    intersection = mask_1 * mask_2
+    union = mask_1**2 + mask_2**2 - mask_1 * mask_2
+    return float((intersection.sum() + beta) / (union.sum() + beta))
+
+
+def _rasterize_fibers(X, F, image_shape):
+    """
+    Rasterize fire_2d_angle fiber output to a 1-px skeleton image.
+
+    Uses Bresenham line drawing between consecutive fiber vertices then
+    morphological skeletonize to guarantee 1-px thickness.  Vertex indices
+    follow the 0-based convention (X[v] = coordinates for vertex v) as
+    established in CLAUDE_CTFIRE.md.
+    """
+    from skimage.draw import line as draw_line
+    from skimage.morphology import skeletonize
+    canvas = np.zeros(image_shape, dtype=np.uint8)
+    H, W = image_shape
+    X_arr = np.asarray(X)
+    for fiber in F:
+        v_list = fiber['v'] if isinstance(fiber, dict) else list(fiber)
+        for seg in range(len(v_list) - 1):
+            v0, v1 = v_list[seg], v_list[seg + 1]
+            if v0 >= len(X_arr) or v1 >= len(X_arr):
+                continue
+            r0 = int(np.clip(round(float(X_arr[v0, 0])), 0, H - 1))
+            c0 = int(np.clip(round(float(X_arr[v0, 1])), 0, W - 1))
+            r1 = int(np.clip(round(float(X_arr[v1, 0])), 0, H - 1))
+            c1 = int(np.clip(round(float(X_arr[v1, 1])), 0, W - 1))
+            rr, cc = draw_line(r0, c0, r1, c1)
+            canvas[rr, cc] = 1
+    return skeletonize(canvas > 0)
+
+
+def _make_synthetic_fiber_image(shape=(256, 256), n_fibers=8, fiber_sigma=2.5, rng_seed=42):
+    """
+    Generate a synthetic fiber image with a known ground-truth skeleton.
+
+    Straight-line fibers are drawn with Gaussian cross-section profiles
+    (signal amplitude 180, background mean 10 with Gaussian noise σ=4)
+    so the FIRE distance-transform pipeline can find them at thresh_im2=5.
+    Returns the image and the skeletonized 1-px ground-truth centerline mask.
+
+    Why synthetic images:
+      - Exact centerline positions are known at generation time — soft IoU
+        measures genuine spatial recovery, not just "did any fibers come out".
+      - Deterministic RNG seeds make CI results reproducible.
+      - No MATLAB .mat reference files required.
+      - Straight-line Gaussian fibers are the canonical FIRE input; regressions
+        in extend_xlink, trimxfv, or filtering stages show up as IoU drops.
+    """
+    from scipy.ndimage import distance_transform_edt
+    from skimage.draw import line as draw_line
+    from skimage.morphology import skeletonize
+    rng = np.random.default_rng(rng_seed)
+    H, W = shape
+    margin = 20
+    skeleton = np.zeros(shape, dtype=bool)
+    generated = 0
+    while generated < n_fibers:
+        r0 = int(rng.integers(margin, H - margin))
+        c0 = int(rng.integers(margin, W - margin))
+        angle = rng.uniform(0, np.pi)
+        length = int(rng.integers(60, min(H, W) - 2 * margin))
+        r1 = int(np.clip(r0 + length * np.sin(angle), margin, H - margin))
+        c1 = int(np.clip(c0 + length * np.cos(angle), margin, W - margin))
+        if np.hypot(r1 - r0, c1 - c0) < 50:
+            continue
+        rr, cc = draw_line(r0, c0, r1, c1)
+        skeleton[rr, cc] = True
+        generated += 1
+    dist = distance_transform_edt(~skeleton).astype(np.float32)
+    signal = 180.0 * np.exp(-dist**2 / (2.0 * fiber_sigma**2))
+    bg = rng.normal(10.0, 4.0, size=shape).astype(np.float32)
+    image = np.clip(bg + signal, 0.0, 255.0).astype(np.float32)
+    return image, skeletonize(skeleton)
+
+
+def _default_fire_params():
+    """
+    Return the standard FIRE algorithm parameters from test_cases_fire_2d.json.
+
+    Uses the 'real1_ctfire_params' test case (the only one whose image is
+    available; 2B_D9_ROI1.tif is not present in the repo).
+    Intended for use with real biological images.
+    """
+    config_path = (
+        Path(__file__).parent
+        / "test_results"
+        / "fire_2d_test_files"
+        / "test_cases_fire_2d.json"
+    )
+    with open(config_path, "r") as f:
+        cfg = json.load(f)
+    case = next(c for c in cfg["test_cases"] if c["name"] == "real1_ctfire_params")
+    return dict(case["params"])
+
+
+def _synthetic_fire_params():
+    """
+    Return FIRE algorithm parameters tuned for synthetic test images.
+
+    Synthetic images have Gaussian-profile fibers (peak ~180, background ~10).
+    The JSON default of thresh_im2=5 (absolute) includes nearly all background
+    pixels at that level, producing hundreds of spurious short zigzag fibers.
+    Using a fractional threshold (thresh_im=0.2, i.e. 20% of image maximum)
+    raises the effective cutoff to ~36, cleanly separating fiber signal from
+    background and suppressing spurious detections.
+    """
+    p = _default_fire_params()
+    p["thresh_im"]  = 0.2   # fractional: keep pixels > 20% of max (~36 for peak 180)
+    p["thresh_im2"] = []    # disable absolute threshold when thresh_im is set
+    return p
+
+
 # ============================================================================
 # Basic Functionality Tests
 # ============================================================================
@@ -445,6 +589,162 @@ def test_fire_2d_matches_matlab_angles(test_name, test_case):
             print(f"\nAngle distribution correlation: {correlation:.3f}")
 
 
+# ============================================================================
+# Soft IoU Validation
+# ============================================================================
+
+
+class TestSoftIoU:
+    """
+    Validate that fire_2d_angle recovers fiber centerlines at fiber scale.
+
+    Soft IoU at sigma=5 (≈10-px FWHM) measures spatial overlap of 1-px
+    skeleton images after Gaussian smoothing — a score > 0.30 means the
+    extracted skeleton substantially overlaps the ground truth at fiber
+    granularity, not pixel precision.
+    """
+
+    def test_soft_iou_synthetic(self):
+        """
+        Synthetic image with known GT skeleton — always runs, no .mat needed.
+
+        Validates:
+        1. Soft IoU of extracted centerlines vs ground-truth skeleton > 0.30
+        2. Total extracted fiber length is at least 20% of GT skeleton length
+        3. Extracted fiber angles span a reasonable range (std > 0.3 rad),
+           confirming the angle computation is not degenerate
+        """
+        if not CPP_AVAILABLE:
+            pytest.skip("C++ backend not available")
+
+        image, gt_skeleton = _make_synthetic_fiber_image(rng_seed=42)
+        p = _synthetic_fire_params()
+        data = fire_2d_angle(p=p, im=image, plotflag=0)
+
+        assert len(data["Ff"]) > 0, "no filtered fibers extracted"
+
+        # 1. Soft IoU
+        pred = _rasterize_fibers(data["Xf"], data["Ff"], image.shape)
+        iou = _soft_iou(_smooth_mask(gt_skeleton.astype(np.float32)),
+                        _smooth_mask(pred.astype(np.float32)))
+        assert iou > SOFT_IOU_THRESHOLD_SYNTHETIC, (
+            f"soft IoU {iou:.3f} < {SOFT_IOU_THRESHOLD_SYNTHETIC} — "
+            "extracted centerlines do not overlap ground truth at fiber scale"
+        )
+
+        # 2. Total extracted length vs GT skeleton pixel count (arc-length proxy)
+        tot_L = float(data["M"]["totL"])
+        assert tot_L > 0, "total extracted fiber length is zero"
+        gt_length = float(gt_skeleton.sum())
+        assert tot_L > 0.2 * gt_length, (
+            f"extracted total length {tot_L:.1f} px < 20% of GT skeleton "
+            f"length {gt_length:.1f} px — pipeline may be discarding too many fibers"
+        )
+
+        # 3. Angle range sanity check
+        angles = np.asarray(data["M"].get("angle_xy", []))
+        assert len(angles) > 0, "no fiber angles in M['angle_xy']"
+        assert np.all(np.abs(angles) <= np.pi + 1e-6), "angle value outside [-π, π]"
+        assert angles.std() > 0.3, (
+            f"angle std {angles.std():.3f} rad unexpectedly small — "
+            "all fibers have nearly the same orientation on a random image"
+        )
+
+    @pytest.mark.parametrize("seed", [0, 7, 99])
+    def test_soft_iou_multiple_seeds(self, seed):
+        """
+        Soft IoU > 0.30 across multiple random synthetic configurations.
+
+        Guards against a lucky pass on a single seed by checking that the
+        pipeline recovers fibers consistently across different random images.
+        """
+        if not CPP_AVAILABLE:
+            pytest.skip("C++ backend not available")
+
+        image, gt_skeleton = _make_synthetic_fiber_image(rng_seed=seed)
+        p = _synthetic_fire_params()
+        data = fire_2d_angle(p=p, im=image, plotflag=0)
+
+        assert len(data["Ff"]) > 0, f"no fibers extracted (seed={seed})"
+        pred = _rasterize_fibers(data["Xf"], data["Ff"], image.shape)
+        iou = _soft_iou(_smooth_mask(gt_skeleton.astype(np.float32)),
+                        _smooth_mask(pred.astype(np.float32)))
+        assert iou > SOFT_IOU_THRESHOLD_SYNTHETIC, (
+            f"soft IoU {iou:.3f} < {SOFT_IOU_THRESHOLD_SYNTHETIC} (seed={seed})"
+        )
+
+    @pytest.mark.matlab
+    @pytest.mark.parametrize(
+        "test_name,test_case",
+        load_test_cases(matlab_only=True),
+        ids=[name for name, _ in load_test_cases(matlab_only=True)],
+    )
+    def test_soft_iou_matlab_vs_python(self, test_name, test_case):
+        """
+        Compare Python centerlines against MATLAB Xa/Fa via soft IoU.
+
+        Skips gracefully if the .mat reference file is absent.
+        Also checks total length (within ±40%) and mean |angle| (within 10°).
+        """
+        if not CPP_AVAILABLE:
+            pytest.skip("C++ backend not available")
+
+        # Skip if test image is not present (e.g. 2B_D9_ROI1.tif is not in the repo)
+        img_path = Path(__file__).parent / "test_images" / test_case["image"]
+        if not img_path.exists():
+            pytest.skip(f"Test image not found: {test_case['image']}")
+
+        # Load image and run Python extraction
+        img = load_test_image(test_case["image"])
+        if img.ndim == 2:
+            img = img[np.newaxis, :, :]
+
+        data_py = fire_2d_angle(p=test_case["params"], im=img, plotflag=0)
+
+        # Load MATLAB reference (skips if file absent)
+        mat_path = (
+            Path(__file__).parent
+            / "test_results"
+            / "fire_2d_test_files"
+            / test_case["matlab_reference_mat"]
+        )
+        data_mat = load_matlab_reference(mat_path)
+
+        image_2d = img[0] if img.ndim == 3 else img
+        H, W = image_2d.shape
+
+        # Soft IoU — rasterize MATLAB Xa/Fa vs Python Xf/Ff
+        if data_mat["Xa"] is not None and data_mat.get("Fa") is not None:
+            mat_skel = _rasterize_fibers(data_mat["Xa"], data_mat["Fa"], (H, W))
+            py_skel  = _rasterize_fibers(data_py["Xf"], data_py["Ff"], (H, W))
+            iou = _soft_iou(_smooth_mask(mat_skel.astype(np.float32)),
+                            _smooth_mask(py_skel.astype(np.float32)))
+            assert iou > SOFT_IOU_THRESHOLD_MATLAB, (
+                f"soft IoU {iou:.3f} < {SOFT_IOU_THRESHOLD_MATLAB} "
+                f"({test_name}): Python and MATLAB centerlines diverge spatially"
+            )
+            print(f"\n{test_name} soft IoU = {iou:.3f}")
+
+        # Total length within 40% of MATLAB reference
+        mat_totL = data_mat["M"].get("totL", 0)
+        py_totL  = float(data_py["M"]["totL"])
+        if mat_totL > 0 and py_totL > 0:
+            assert 0.6 * mat_totL <= py_totL <= 1.4 * mat_totL, (
+                f"total length {py_totL:.1f} not within 40% of MATLAB {mat_totL:.1f}"
+            )
+
+        # Mean |angle| within 10° of MATLAB reference
+        mat_angles = np.asarray(data_mat["M"].get("angle_xy", []))
+        py_angles  = np.asarray(data_py["M"].get("angle_xy", []))
+        if len(mat_angles) > 0 and len(py_angles) > 0:
+            delta_deg = np.degrees(
+                abs(np.mean(np.abs(py_angles)) - np.mean(np.abs(mat_angles)))
+            )
+            assert delta_deg < 10.0, (
+                f"mean |angle| differs by {delta_deg:.1f}° > 10° ({test_name})"
+            )
+
+
 # ============================================================================
 # Utility Tests
 # ============================================================================
@@ -514,5 +814,32 @@ def test_implementation_status_documented():
 
 
 if __name__ == "__main__":
-    # Run tests with verbose output
-    pytest.main([__file__, "-v", "-s", "--tb=short"])
+    # Standalone demo: generate synthetic image, run FIRE, print metrics, save overlay.
+    # Run with:  python tests/test_fire_2d_angle.py
+    from ctfire_py.test_fire_2d import plot_fiber_overlay
+
+    image, gt_skeleton = _make_synthetic_fiber_image(rng_seed=42)
+    p = _synthetic_fire_params()
+    data = fire_2d_angle(p=p, im=image, plotflag=0)
+
+    pred = _rasterize_fibers(data["Xf"], data["Ff"], image.shape)
+    iou = _soft_iou(_smooth_mask(gt_skeleton.astype(np.float32)),
+                    _smooth_mask(pred.astype(np.float32)))
+
+    angles = np.asarray(data["M"].get("angle_xy", []))
+    mean_angle_deg = float(np.degrees(np.mean(np.abs(angles)))) if len(angles) > 0 else float("nan")
+
+    print(
+        f"Soft IoU    = {iou:.4f}  (threshold {SOFT_IOU_THRESHOLD_SYNTHETIC})\n"
+        f"Fibers      = {len(data['Ff'])}\n"
+        f"Total length= {data['M']['totL']:.1f} px\n"
+        f"Mean |angle|= {mean_angle_deg:.1f}°"
+    )
+
+    plot_fiber_overlay(
+        image,
+        data["Xf"],
+        data["Ff"],
+        title=f"FIRE overlay  (soft IoU={iou:.3f})",
+        save_path="fire_2d_soft_iou_overlay.png",
+    )