simam3d_core.py

"""Pure NumPy geometry utilities for Simam3D.

This module intentionally has no Torch, Gradio, or CUDA dependency. It is the
testable contract between model adapters and the export/UI layer.
"""

from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path

import numpy as np


@dataclass(frozen=True)
class FusionResult:
    points: np.ndarray
    colors: np.ndarray
    weights: np.ndarray
    source_views: np.ndarray
    source_masks: np.ndarray


def normalize_map(values: np.ndarray) -> np.ndarray:
    values = np.asarray(values, dtype=np.float32)
    finite = np.isfinite(values)
    if not finite.any():
        return np.zeros_like(values, dtype=np.float32)
    lo = float(values[finite].min())
    hi = float(values[finite].max())
    return np.clip(np.nan_to_num((values - lo) / max(hi - lo, 1e-6)), 0.0, 1.0)


def analyze_light_field(colors: np.ndarray) -> dict[str, float | str]:
    """Estimate image-space illumination cues from an RGB image.

    This is an auditable cue detector, not physical inverse rendering. It uses
    luminance, contrast, and the centroid of bright pixels to provide a stable
    lighting hypothesis for the UI and future learned relighting work.
    """
    rgb = np.asarray(colors, dtype=np.float32)
    if rgb.ndim != 3 or rgb.shape[2] < 3 or min(rgb.shape[:2]) < 2:
        raise ValueError("colors must have shape HxWx3 or HxWx4")
    rgb = np.clip(rgb[..., :3], 0.0, 255.0) / 255.0
    luminance = 0.2126 * rgb[..., 0] + 0.7152 * rgb[..., 1] + 0.0722 * rgb[..., 2]
    mean = float(luminance.mean())
    contrast = float(luminance.std())
    bright_threshold = float(np.percentile(luminance, 90))
    bright = luminance >= max(bright_threshold, 0.65)
    yy, xx = np.mgrid[0:luminance.shape[0], 0:luminance.shape[1]]
    mass = luminance * bright
    total = float(mass.sum())
    if total > 1e-6:
        cx = float((xx * mass).sum() / total / max(luminance.shape[1] - 1, 1) * 2.0 - 1.0)
        cy = float(1.0 - (yy * mass).sum() / total / max(luminance.shape[0] - 1, 1) * 2.0)
    else:
        cx = cy = 0.0
    highlight_fraction = float(bright.mean())
    shadow_fraction = float((luminance <= np.percentile(luminance, 10)).mean())
    azimuth = float(np.degrees(np.arctan2(cx, max(cy, 1e-6))))
    elevation = float(np.clip(35.0 + 45.0 * cy, 5.0, 85.0))
    if contrast < 0.08:
        label = "soft / low-contrast light"
    elif highlight_fraction < 0.03 and shadow_fraction > 0.18:
        label = "directional light with strong shadow cues"
    else:
        label = "directional light with diffuse fill"
    confidence = float(np.clip(contrast * 3.0 + min(highlight_fraction * 2.0, 0.25), 0.0, 1.0))
    return {
        "lighting_label": label,
        "mean_luminance": mean,
        "contrast": contrast,
        "highlight_fraction": highlight_fraction,
        "shadow_fraction": shadow_fraction,
        "bright_centroid_x": cx,
        "bright_centroid_y": cy,
        "estimated_azimuth_deg": azimuth,
        "estimated_elevation_deg": elevation,
        "detector_confidence": confidence,
    }


def estimate_depth_light(depth: np.ndarray, colors: np.ndarray) -> dict[str, float | str]:
    """Estimate a camera-space light direction from depth-derived normals.

    This is a low-cost photometric heuristic. It is useful for ranking light
    hypotheses and checking generated views, but it is not physical inverse
    rendering and cannot disambiguate albedo from illumination by itself.
    """
    depth = np.asarray(depth, dtype=np.float32)
    rgb = np.asarray(colors, dtype=np.float32)
    if depth.ndim != 2 or rgb.ndim != 3 or rgb.shape[:2] != depth.shape or rgb.shape[2] < 3:
        raise ValueError("depth must be HxW and colors must be matching HxWx3")
    h, w = depth.shape
    z = 1.0 - normalize_map(depth)
    yy, xx = np.mgrid[0:h, 0:w].astype(np.float32)
    x = (xx / max(w - 1, 1) - 0.5) * 2.0
    y = (0.5 - yy / max(h - 1, 1)) * 2.0
    points = np.stack([x, y, z], axis=-1)
    dx = np.gradient(points, axis=1)
    dy = np.gradient(points, axis=0)
    normals = np.cross(dx, dy)
    lengths = np.linalg.norm(normals, axis=-1, keepdims=True)
    normals = normals / np.maximum(lengths, 1e-6)
    normals = np.where(normals[..., 2:3] < 0.0, -normals, normals)
    luminance = (0.2126 * rgb[..., 0] + 0.7152 * rgb[..., 1] + 0.0722 * rgb[..., 2]) / 255.0
    signal = np.clip(luminance - float(luminance.mean()), 0.0, None)
    total = float(signal.sum())
    if total <= 1e-6:
        direction = np.array([0.0, 0.0, 1.0], dtype=np.float32)
        confidence = 0.0
    else:
        direction = (normals * signal[..., None]).sum(axis=(0, 1)) / total
        direction = direction / max(float(np.linalg.norm(direction)), 1e-6)
        confidence = float(np.clip(np.linalg.norm((normals * signal[..., None]).sum(axis=(0, 1))) / total, 0.0, 1.0))
    azimuth = float(np.degrees(np.arctan2(float(direction[0]), float(direction[2]))))
    elevation = float(np.degrees(np.arcsin(np.clip(float(direction[1]), -1.0, 1.0))))
    return {
        "light_direction_x": float(direction[0]),
        "light_direction_y": float(direction[1]),
        "light_direction_z": float(direction[2]),
        "estimated_azimuth_deg": azimuth,
        "estimated_elevation_deg": elevation,
        "photometric_confidence": confidence,
        "status": "depth-normal photometric heuristic; albedo and exposure are confounds",
    }


def deterministic_demo_depth(colors: np.ndarray) -> np.ndarray:
    """Create an explicitly non-AI depth-like map for export-pipeline smoke tests."""
    rgb = np.asarray(colors, dtype=np.float32)
    if rgb.ndim != 3 or rgb.shape[2] < 3:
        raise ValueError("colors must have shape HxWx3 or HxWx4")
    luminance = 0.2126 * rgb[..., 0] + 0.7152 * rgb[..., 1] + 0.0722 * rgb[..., 2]
    # Invert brightness only to create stable synthetic geometry; this has no
    # semantic or metric meaning and must remain labelled as demo output.
    return (255.0 - luminance).astype(np.float32)


def classify_scene_hypothesis(text: str | None) -> dict[str, str | list[str]]:
    """Turn the user's reveal guess into explicit, inspectable semantic cues."""
    raw = (text or "").strip()
    normalized = " ".join(raw.split())
    lower = normalized.lower()
    cue_groups = {
        "outdoor": ("sky", "street", "garden", "forest", "mountain", "building", "road"),
        "indoor": ("room", "wall", "kitchen", "office", "hallway", "inside"),
        "person": ("person", "face", "body", "portrait", "human"),
        "object": ("object", "chair", "car", "shoe", "bottle", "product", "box"),
        "nature": ("tree", "plant", "garden", "water", "animal", "cloud", "rock"),
    }
    cues = [name for name, words in cue_groups.items() if any(word in lower for word in words)]
    return {
        "user_guess": normalized or "No prior guess supplied",
        "semantic_cues": cues or ["unspecified"],
        "status": "user hypothesis only; not verified hidden content",
    }


def reveal_uncertainty(depth: np.ndarray, confidence: np.ndarray | None = None) -> np.ndarray:
    """Build a visible-space uncertainty map for deciding where to inspect next.

    Low confidence and strong depth discontinuities are useful acquisition cues,
    but this map does not assert what exists outside the camera's field of view.
    """
    depth = np.asarray(depth, dtype=np.float32)
    if depth.ndim != 2:
        raise ValueError("depth must be a 2-D map")
    normalized = normalize_map(depth)
    gy, gx = np.gradient(normalized)
    edges = np.sqrt(gx * gx + gy * gy)
    edges = normalize_map(edges)
    if confidence is None:
        conf = np.ones_like(normalized)
    else:
        conf = np.asarray(confidence, dtype=np.float32)
        if conf.shape != depth.shape:
            conf = np.ones_like(normalized)
        conf = np.nan_to_num(conf, nan=0.0, posinf=0.0, neginf=0.0)
        conf = np.clip(conf, 0.0, None)
        peak = float(conf.max()) if conf.size else 0.0
        conf = conf / peak if peak > 1e-6 else np.zeros_like(normalized)
    return np.clip(0.7 * (1.0 - conf) + 0.3 * edges, 0.0, 1.0).astype(np.float32)


def next_view_plan(view_count: int, poses: np.ndarray | None = None) -> list[dict[str, float | str]]:
    """Suggest yaw targets that fill the largest observed camera gaps.

    When camera poses are available, use their camera centers instead of
    assuming that submitted views are evenly spaced. This remains a 2-D
    azimuth heuristic and is not a full next-best-view optimizer.
    """
    if view_count < 1:
        raise ValueError("view_count must be positive")
    existing = None
    if poses is not None:
        raw_poses = np.asarray(poses, dtype=np.float32)
        if raw_poses.ndim == 3 and raw_poses.shape[0] == view_count:
            centers = []
            for raw_pose in raw_poses:
                pose = ensure_pose(raw_pose)
                center = np.linalg.inv(pose)[:3, 3]
                if float(np.linalg.norm(center[[0, 2]])) > 1e-6:
                    centers.append(float(np.degrees(np.arctan2(center[0], center[2])) % 360.0))
            if centers:
                existing = np.asarray(centers, dtype=np.float32)
    if existing is None:
        existing = np.linspace(0.0, 360.0, view_count, endpoint=False)
    targets = [180.0, 90.0, 270.0, 45.0, 315.0]
    suggestions = []
    for yaw in targets:
        distance = float(np.min(np.abs(((existing - yaw + 180.0) % 360.0) - 180.0)))
        suggestions.append({"yaw_deg": yaw, "gap_from_existing_deg": distance, "reason": "largest observed azimuth gap"})
    return sorted(suggestions, key=lambda item: (-float(item["gap_from_existing_deg"]), float(item["yaw_deg"])))


def orbit_project_points(points: np.ndarray, yaw_deg: float) -> tuple[np.ndarray, np.ndarray]:
    """Project reconstructed points to a normalized orbit view for previews."""
    points = np.asarray(points, dtype=np.float32)
    if points.ndim != 2 or points.shape[1] != 3:
        raise ValueError("points must have shape Nx3")
    angle = np.deg2rad(float(yaw_deg))
    cosine, sine = np.cos(angle), np.sin(angle)
    rotation = np.array([[cosine, 0.0, -sine], [0.0, 1.0, 0.0], [sine, 0.0, cosine]], dtype=np.float32)
    view = points @ rotation.T
    extent = max(float(np.max(np.abs(view[:, :2]))) if len(view) else 0.0, 1e-6)
    xy = np.clip(view[:, :2] / extent, -1.0, 1.0)
    return xy.astype(np.float32), view[:, 2].astype(np.float32)


def default_intrinsics(width: int, height: int) -> np.ndarray:
    focal = float(max(width, height))
    return np.array(
        [[focal, 0.0, width / 2.0], [0.0, focal, height / 2.0], [0.0, 0.0, 1.0]],
        dtype=np.float32,
    )


def prepare_intrinsics(matrix: np.ndarray | None, width: int, height: int) -> np.ndarray:
    """Validate intrinsics and convert normalized matrices to pixel units."""
    if matrix is None:
        return default_intrinsics(width, height)
    k = np.asarray(matrix, dtype=np.float32).copy()
    if k.shape != (3, 3) or not np.isfinite(k).all() or abs(float(k[2, 2])) < 1e-6:
        return default_intrinsics(width, height)
    # The V2 placeholder is identity, which means unknown rather than a real
    # camera. Normalized camera matrices commonly use fx/fy and cx/cy in [0, 1].
    if np.allclose(k, np.eye(3, dtype=np.float32)):
        return default_intrinsics(width, height)
    if abs(float(k[0, 0])) <= 2.0 and abs(float(k[1, 1])) <= 2.0:
        k[0, :] *= width
        k[1, :] *= height
    if abs(float(k[0, 0])) < 1e-6 or abs(float(k[1, 1])) < 1e-6:
        return default_intrinsics(width, height)
    return k


def ensure_pose(matrix: np.ndarray | None) -> np.ndarray:
    if matrix is None:
        return np.eye(4, dtype=np.float32)
    pose = np.asarray(matrix, dtype=np.float32)
    if pose.shape == (3, 4):
        pose = np.vstack([pose, [0.0, 0.0, 0.0, 1.0]])
    if pose.shape != (4, 4) or not np.isfinite(pose).all():
        return np.eye(4, dtype=np.float32)
    return pose


def synthetic_orbit_poses(view_count: int, radius: float = 0.35) -> np.ndarray:
    """Return a deterministic orbit prior for views without camera metadata.

    These are world-to-camera poses for visualization/fusion only. They are not
    a substitute for calibrated camera poses and are labeled as synthetic by
    the application manifest.
    """
    if view_count < 1:
        return np.empty((0, 4, 4), dtype=np.float32)
    if view_count == 1:
        return np.eye(4, dtype=np.float32)[None, ...]
    poses = []
    up = np.array([0.0, 1.0, 0.0], dtype=np.float32)
    for index in range(view_count):
        angle = 2.0 * np.pi * index / max(view_count, 1)
        position = np.array([radius * np.sin(angle), 0.0, radius * np.cos(angle)], dtype=np.float32)
        forward = -position / max(float(np.linalg.norm(position)), 1e-6)
        right = np.cross(up, forward)
        right /= max(float(np.linalg.norm(right)), 1e-6)
        true_up = np.cross(forward, right)
        camera_to_world = np.eye(4, dtype=np.float32)
        camera_to_world[:3, :3] = np.stack([right, true_up, forward], axis=1)
        camera_to_world[:3, 3] = position
        poses.append(np.linalg.inv(camera_to_world).astype(np.float32))
    return np.stack(poses, axis=0)


def normalize_prediction_contract(prediction, view_count: int):
    """Validate and normalize the official DA3 Prediction output contract."""
    depth = np.asarray(getattr(prediction, "depth", None), dtype=np.float32)
    if depth.ndim == 2 and view_count == 1:
        depth = depth[None, ...]
    if depth.ndim != 3 or depth.shape[0] != view_count or not np.isfinite(depth).all():
        raise ValueError(f"DA3 depth must have shape ({view_count}, H, W), got {depth.shape}")
    raw_confidence = getattr(prediction, "conf", None)
    if raw_confidence is None:
        confidence = np.ones_like(depth, dtype=np.float32)
        confidence_source = "baseline_uniform"
    else:
        confidence = np.asarray(raw_confidence, dtype=np.float32)
        if confidence.ndim == 2 and view_count == 1:
            confidence = confidence[None, ...]
        if confidence.shape != depth.shape or not np.isfinite(confidence).all():
            raise ValueError(f"DA3 confidence must match depth shape {depth.shape}, got {confidence.shape}")
        confidence_source = "depth_anything_3"
    raw_extrinsics = getattr(prediction, "extrinsics", None)
    if raw_extrinsics is not None:
        extrinsics = np.asarray(raw_extrinsics, dtype=np.float32)
        valid_pose_shapes = {(view_count, 3, 4), (view_count, 4, 4)}
        if extrinsics.shape not in valid_pose_shapes or not np.isfinite(extrinsics).all():
            raise ValueError(f"DA3 extrinsics have invalid shape {extrinsics.shape}")
        try:
            for pose in extrinsics:
                np.linalg.inv(ensure_pose(pose))
        except np.linalg.LinAlgError:
            extrinsics = synthetic_orbit_poses(view_count)
            pose_source = "synthetic_orbit_prior"
        else:
            pose_source = "depth_anything_3"
    else:
        extrinsics = synthetic_orbit_poses(view_count)
        pose_source = "synthetic_orbit_prior"
    raw_intrinsics = getattr(prediction, "intrinsics", None)
    if raw_intrinsics is not None:
        intrinsics = np.asarray(raw_intrinsics, dtype=np.float32)
        if intrinsics.shape != (view_count, 3, 3) or not np.isfinite(intrinsics).all():
            raise ValueError(f"DA3 intrinsics have invalid shape {intrinsics.shape}")
    else:
        intrinsics = np.repeat(np.eye(3, dtype=np.float32)[None, ...], view_count, axis=0)
    return depth, confidence, extrinsics, intrinsics, confidence_source, pose_source


def project_world_points(
    points: np.ndarray,
    intrinsics: np.ndarray | None,
    extrinsics: np.ndarray | None,
    width: int,
    height: int,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Project world-space points into a view and return pixels, depth, validity."""
    points = np.asarray(points, dtype=np.float32)
    if points.ndim != 2 or points.shape[1] != 3:
        raise ValueError("points must have shape Nx3")
    pose = ensure_pose(extrinsics)
    camera = (pose @ np.c_[points, np.ones(len(points), dtype=np.float32)].T).T[:, :3]
    k = prepare_intrinsics(intrinsics, width, height)
    positive = camera[:, 2] > 1e-6
    pixels = np.full((len(points), 2), np.nan, dtype=np.float32)
    pixels[positive, 0] = k[0, 0] * camera[positive, 0] / camera[positive, 2] + k[0, 2]
    pixels[positive, 1] = k[1, 1] * camera[positive, 1] / camera[positive, 2] + k[1, 2]
    inside = positive & (pixels[:, 0] >= 0.0) & (pixels[:, 0] < width) & (pixels[:, 1] >= 0.0) & (pixels[:, 1] < height)
    return pixels, camera[:, 2].astype(np.float32), inside


def reprojection_consistency(
    points: np.ndarray,
    depth: np.ndarray,
    intrinsics: np.ndarray | None,
    extrinsics: np.ndarray | None,
) -> dict[str, float]:
    """Compare fused-point depth against a view's relative depth prediction."""
    depth = np.asarray(depth, dtype=np.float32)
    if depth.ndim != 2:
        raise ValueError("depth must be a 2-D map")
    height, width = depth.shape
    pixels, point_depth, inside = project_world_points(
        points, intrinsics, extrinsics, width, height
    )
    if not inside.any():
        return {"coverage": 0.0, "mean_absolute_error": 0.0, "median_absolute_error": 0.0}
    raster = np.full((height, width), np.nan, dtype=np.float32)
    pixel_x = np.clip(np.rint(pixels[inside, 0]).astype(np.int64), 0, width - 1)
    pixel_y = np.clip(np.rint(pixels[inside, 1]).astype(np.int64), 0, height - 1)
    for x, y, value in zip(pixel_x, pixel_y, point_depth[inside]):
        if not np.isfinite(raster[y, x]) or value < raster[y, x]:
            raster[y, x] = value
    valid = np.isfinite(raster)
    if not valid.any():
        return {"coverage": 0.0, "mean_absolute_error": 0.0, "median_absolute_error": 0.0}
    target = 1.0 - normalize_map(depth)
    error = np.abs(raster[valid] - target[valid])
    return {
        "coverage": float(valid.mean()),
        "mean_absolute_error": float(error.mean()),
        "median_absolute_error": float(np.median(error)),
    }


def depth_to_camera_points(
    depth: np.ndarray,
    intrinsics: np.ndarray | None = None,
    confidence: np.ndarray | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Project a depth map into camera-space points and return point weights."""
    depth = np.asarray(depth, dtype=np.float32)
    if depth.ndim != 2:
        raise ValueError(f"depth must be HxW, got {depth.shape}")
    height, width = depth.shape
    k = prepare_intrinsics(intrinsics, width, height)
    yy, xx = np.mgrid[0:height, 0:width]
    # Depth Anything outputs relative depth; the scaffold uses inverse-depth
    # convention consistently with the GLB/PLY exporters (higher prediction
    # values place a point closer to the camera). Metric scale is not claimed.
    z = 1.0 - normalize_map(depth)
    z = np.clip(z, 0.02, 1.0)
    points = np.stack(
        [(xx - k[0, 2]) / max(float(k[0, 0]), 1e-6) * z,
         (yy - k[1, 2]) / max(float(k[1, 1]), 1e-6) * z,
         z],
        axis=-1,
    ).reshape(-1, 3)
    if confidence is None:
        weights = np.ones((height, width), dtype=np.float32)
    else:
        weights = np.nan_to_num(np.asarray(confidence, dtype=np.float32), nan=0.0)
        if weights.shape != (height, width):
            weights = np.ones((height, width), dtype=np.float32)
        weights = np.clip(weights, 0.0, None)
        peak = float(weights.max()) if weights.size else 0.0
        # Preserve absolute confidence ratios across views: min-max scaling
        # would turn a valid 0.5 confidence into zero when it is the minimum.
        weights = weights / peak if peak > 1e-6 else np.zeros_like(weights)
    return points, weights.reshape(-1)


def fuse_point_views(
    views: list[tuple[np.ndarray, np.ndarray, np.ndarray | None, np.ndarray | None, np.ndarray]],
    voxel_size: float = 0.02,
    max_points: int = 150_000,
) -> FusionResult:
    """Fuse `(depth, confidence, intrinsics, extrinsics, colors)` views.

    Extrinsics are interpreted as world-to-camera, matching DA3's API. Fusion
    is deterministic: points are grouped into voxels and confidence-weighted
    centroids/colors are retained, then capped by voxel confidence.
    """
    points_all, colors_all, weights_all, sources_all, masks_all = [], [], [], [], []
    for source_view, (depth, confidence, intrinsics, extrinsics, colors) in enumerate(views):
        depth = np.asarray(depth, dtype=np.float32)
        colors = np.asarray(colors, dtype=np.uint8)
        if colors.shape[:2] != depth.shape or colors.shape[-1] < 3:
            raise ValueError("colors must have shape HxWx3 or HxWx4 matching depth")
        camera_points, weights = depth_to_camera_points(depth, intrinsics, confidence)
        pose = ensure_pose(extrinsics)
        world_points = (np.linalg.inv(pose) @ np.c_[camera_points, np.ones(len(camera_points))].T).T[:, :3]
        colors_all.append(colors[..., :3].reshape(-1, 3))
        points_all.append(world_points)
        weights_all.append(weights)
        sources_all.append(np.full(len(weights), source_view, dtype=np.int32))
        masks_all.append(np.full(len(weights), 1 << source_view, dtype=np.int32))
    if not points_all:
        empty = np.empty((0, 3), dtype=np.float32)
        return FusionResult(
            empty, np.empty((0, 3), dtype=np.uint8), np.empty((0,), dtype=np.float32),
            np.empty((0,), dtype=np.int32), np.empty((0,), dtype=np.int32)
        )
    if len(views) > 30:
        raise ValueError("At most 30 views are supported by the signed int32 provenance bitmask")
    points = np.concatenate(points_all).astype(np.float32)
    colors = np.concatenate(colors_all).astype(np.uint8)
    weights = np.concatenate(weights_all).astype(np.float32)
    source_views = np.concatenate(sources_all).astype(np.int32)
    source_masks = np.concatenate(masks_all).astype(np.int32)
    valid = np.isfinite(points).all(axis=1) & np.isfinite(weights) & (weights > 0)
    points, colors, weights = points[valid], colors[valid], weights[valid]
    source_views, source_masks = source_views[valid], source_masks[valid]
    if voxel_size > 0 and len(points):
        voxels = np.floor(points / voxel_size).astype(np.int64)
        unique_voxels, inverse = np.unique(voxels, axis=0, return_inverse=True)
        grouped_points = np.zeros((len(unique_voxels), 3), dtype=np.float64)
        grouped_colors = np.zeros((len(unique_voxels), 3), dtype=np.float64)
        grouped_weight = np.zeros(len(unique_voxels), dtype=np.float64)
        np.add.at(grouped_points, inverse, points * weights[:, None])
        np.add.at(grouped_colors, inverse, colors * weights[:, None])
        np.add.at(grouped_weight, inverse, weights)
        counts = np.bincount(inverse, minlength=len(unique_voxels)).astype(np.float64)
        # Pick the strongest individual observation as an auditable provenance
        # pointer for each voxel; full multi-source lineage remains future work.
        raw_order = np.argsort(-weights, kind="stable")
        _, first_source = np.unique(inverse[raw_order], return_index=True)
        source_by_voxel = np.zeros(len(unique_voxels), dtype=np.int32)
        source_by_voxel[inverse[raw_order[first_source]]] = source_views[raw_order[first_source]]
        mask_by_voxel = np.zeros(len(unique_voxels), dtype=np.int32)
        np.bitwise_or.at(mask_by_voxel, inverse, source_masks)
        points = (grouped_points / np.maximum(grouped_weight[:, None], 1e-6)).astype(np.float32)
        colors = np.clip(
            grouped_colors / np.maximum(grouped_weight[:, None], 1e-6), 0.0, 255.0
        ).astype(np.uint8)
        weights = (grouped_weight / np.maximum(counts, 1.0)).astype(np.float32)
        source_views = source_by_voxel
        source_masks = mask_by_voxel
    if len(points) > max_points:
        keep = np.argsort(-weights, kind="stable")[:max_points]
        points, colors, weights = points[keep], colors[keep], weights[keep]
        source_views, source_masks = source_views[keep], source_masks[keep]
    return FusionResult(points, colors, weights, source_views, source_masks)


def fusion_metrics(result: FusionResult) -> dict[str, float | int]:
    if not len(result.points):
        return {
            "point_count": 0,
            "confidence_mean": 0.0,
            "confidence_p95": 0.0,
            "dominant_source_view_count": 0,
            "supported_source_view_count": 0,
            "multi_view_voxel_fraction": 0.0,
        }
    supported_mask = int(np.bitwise_or.reduce(result.source_masks))
    multi_view = np.fromiter(
        (int(mask).bit_count() > 1 for mask in result.source_masks),
        dtype=bool,
        count=len(result.source_masks),
    )
    return {
        "point_count": int(len(result.points)),
        "confidence_mean": float(result.weights.mean()),
        "confidence_p95": float(np.percentile(result.weights, 95)),
        "dominant_source_view_count": int(len(np.unique(result.source_views))),
        "supported_source_view_count": supported_mask.bit_count(),
        "multi_view_voxel_fraction": float(multi_view.mean()),
    }


def classify_fusion_support(
    metrics: dict[str, float | int],
    view_count: int,
    unique_view_count: int | None = None,
    pose_source: str | None = None,
    near_duplicate_pair_count: int = 0,
) -> str:
    """Classify support without overstating duplicate or synthetic-pose evidence."""
    if int(view_count) < 2:
        return "single_view_baseline"
    if unique_view_count is not None and int(unique_view_count) < int(view_count):
        return "duplicate_input_support"
    if int(near_duplicate_pair_count) > 0:
        return "near_duplicate_input_support"
    if pose_source is not None and pose_source != "depth_anything_3":
        return "synthetic_prior_support"
    fraction = float(metrics.get("multi_view_voxel_fraction", 0.0))
    if fraction >= 0.5:
        return "strong_multi_view_support"
    if fraction >= 0.1:
        return "partial_multi_view_support"
    return "low_multi_view_support"


def filter_fusion_by_support(result: FusionResult, minimum_views: int) -> FusionResult:
    """Retain only voxels supported by at least ``minimum_views`` source views."""
    minimum_views = int(minimum_views)
    if minimum_views <= 1 or not len(result.points):
        return result
    keep = np.fromiter(
        (int(mask).bit_count() >= minimum_views for mask in result.source_masks),
        dtype=bool,
        count=len(result.source_masks),
    )
    return FusionResult(
        result.points[keep], result.colors[keep], result.weights[keep],
        result.source_views[keep], result.source_masks[keep],
    )


def input_diversity_metrics(
    input_hashes: list[str] | tuple[str, ...], near_duplicate_pair_count: int = 0
) -> dict[str, int | bool | str]:
    """Describe whether submitted views contain duplicate normalized images.

    A repeated image can still exercise batching and export code, but it is not
    evidence of additional visual coverage. Keeping this diagnostic beside the
    fusion metrics prevents a multi-view run from being over-interpreted.
    """
    hashes = [str(value) for value in input_hashes]
    unique_count = len(set(hashes))
    duplicate_count = max(0, len(hashes) - unique_count)
    near_duplicate_pair_count = max(0, int(near_duplicate_pair_count))
    if duplicate_count:
        status = "duplicate_inputs"
    elif near_duplicate_pair_count:
        status = "near_duplicate_inputs"
    else:
        status = "distinct_inputs"
    return {
        "input_view_count": len(hashes),
        "unique_input_view_count": unique_count,
        "duplicate_input_view_count": duplicate_count,
        "near_duplicate_pair_count": near_duplicate_pair_count,
        "has_duplicate_input_views": duplicate_count > 0,
        "input_diversity_status": status,
    }


def gaussian_initialization(weights: np.ndarray | None, point_count: int) -> tuple[np.ndarray, np.ndarray]:
    """Map confidence weights to initial opacity and isotropic log-scales.

    This is deliberately a lightweight initialization heuristic. It does not
    claim to replace a learned Gaussian optimizer.
    """
    if weights is None or len(weights) != point_count:
        normalized = np.ones(point_count, dtype=np.float32)
    else:
        normalized = np.clip(np.asarray(weights, dtype=np.float32), 0.0, 1.0)
    opacity = 0.05 + 0.90 * normalized
    log_scale = -4.605 + 0.8 * (1.0 - normalized)
    return opacity, log_scale


def write_confidence_ply(
    path: str | Path,
    points: np.ndarray,
    colors: np.ndarray,
    confidence: np.ndarray,
    source_views: np.ndarray | None = None,
    source_masks: np.ndarray | None = None,
) -> None:
    """Write an ASCII PLY retaining confidence and optional source provenance."""
    points = np.asarray(points, dtype=np.float32)
    colors = np.asarray(colors, dtype=np.uint8)[..., :3]
    confidence = np.asarray(confidence, dtype=np.float32)
    if points.ndim != 2 or points.shape[1] != 3 or colors.shape != (len(points), 3):
        raise ValueError("points and colors must contain matching Nx3 arrays")
    if confidence.shape != (len(points),):
        raise ValueError("confidence must contain one value per point")
    if source_views is None:
        source_views = np.full(len(points), -1, dtype=np.int32)
    source_views = np.asarray(source_views, dtype=np.int32)
    if source_views.shape != (len(points),):
        raise ValueError("source_views must contain one value per point")
    if source_masks is None:
        source_masks = np.where(source_views >= 0, 1 << source_views, 0).astype(np.int32)
    source_masks = np.asarray(source_masks, dtype=np.int32)
    if source_masks.shape != (len(points),):
        raise ValueError("source_masks must contain one value per point")
    header = [
        "ply", "format ascii 1.0", f"element vertex {len(points)}",
        "property float x", "property float y", "property float z",
        "property uchar red", "property uchar green", "property uchar blue",
        "property float confidence", "property int source_view", "property int source_mask", "end_header",
    ]
    with Path(path).open("w", encoding="utf-8") as handle:
        handle.write("\n".join(header) + "\n")
        for point, color, score, source_view, source_mask in zip(points, colors, confidence, source_views, source_masks):
            handle.write(
                "{:.7g} {:.7g} {:.7g} {} {} {} {:.7g} {} {}\n".format(
                    *point, *color, float(score), int(source_view), int(source_mask)
                )
            )


def write_gaussian_ply(
    path: str | Path,
    points: np.ndarray,
    colors: np.ndarray,
    weights: np.ndarray | None = None,
    source_views: np.ndarray | None = None,
    source_masks: np.ndarray | None = None,
    scales: np.ndarray | None = None,
    rotations: np.ndarray | None = None,
    opacities: np.ndarray | None = None,
) -> None:
    """Write a 3DGS PLY with evidence fields.

    When direct Gaussian parameters are supplied, preserve their scales,
    rotations, and opacities. Otherwise initialize those fields from the
    confidence weights for the portable depth-to-Gaussian baseline.
    """
    points = np.asarray(points, dtype=np.float32)
    colors = np.asarray(colors, dtype=np.uint8)[..., :3]
    if points.ndim != 2 or points.shape[1] != 3 or colors.shape != (len(points), 3):
        raise ValueError("points and colors must contain matching Nx3 arrays")
    if scales is None and rotations is None and opacities is None:
        opacity, log_scale = gaussian_initialization(weights, len(points))
        log_scale = np.repeat(log_scale[:, None], 3, axis=1)
        rotations = np.zeros((len(points), 4), dtype=np.float32)
        rotations[:, 0] = 1.0
    else:
        if scales is None or rotations is None or opacities is None:
            raise ValueError("scales, rotations, and opacities must be supplied together")
        scales = np.asarray(scales, dtype=np.float32)
        rotations = np.asarray(rotations, dtype=np.float32)
        opacity = np.asarray(opacities, dtype=np.float32)
        if scales.shape != (len(points), 3) or rotations.shape != (len(points), 4) or opacity.shape != (len(points),):
            raise ValueError("direct Gaussian parameters must have shapes Nx3, Nx4, and N")
        if not np.isfinite(scales).all() or not np.isfinite(rotations).all() or not np.isfinite(opacity).all():
            raise ValueError("direct Gaussian parameters must be finite")
        # DA3 exposes standard deviations; the PLY convention stores log scales.
        log_scale = np.log(np.clip(np.abs(scales), 1e-6, None))
        opacity = np.clip(opacity, 0.0, 1.0)
        rotation_norm = np.linalg.norm(rotations, axis=1, keepdims=True)
        rotations = rotations / np.clip(rotation_norm, 1e-8, None)
    if weights is None or len(weights) != len(points):
        exported_confidence = np.ones(len(points), dtype=np.float32)
    else:
        exported_confidence = np.clip(np.asarray(weights, dtype=np.float32), 0.0, 1.0)
    if source_views is None or len(source_views) != len(points):
        source_views = np.full(len(points), -1, dtype=np.int32)
    else:
        source_views = np.asarray(source_views, dtype=np.int32)
    if source_masks is None:
        source_masks = np.where(source_views >= 0, 1 << source_views, 0).astype(np.int32)
    source_masks = np.asarray(source_masks, dtype=np.int32)
    if source_masks.shape != (len(points),):
        raise ValueError("source_masks must contain one value per point")
    dc = (colors.astype(np.float32) / 255.0 - 0.5) / 0.2820947918
    header = [
        "ply", "format ascii 1.0", f"element vertex {len(points)}",
        "property float x", "property float y", "property float z",
        "property float nx", "property float ny", "property float nz",
        "property float f_dc_0", "property float f_dc_1", "property float f_dc_2",
        "property float opacity", "property float scale_0", "property float scale_1",
        "property float scale_2", "property float rot_0", "property float rot_1",
        "property float rot_2", "property float rot_3",
        "property float confidence", "property int source_view", "property int source_mask", "end_header",
    ]
    with Path(path).open("w", encoding="utf-8") as handle:
        handle.write("\n".join(header) + "\n")
        for point, color, alpha, scale, rotation, score, source_view, source_mask in zip(
            points, dc, opacity, log_scale, rotations, exported_confidence, source_views, source_masks
        ):
            numeric = [*point, 0.0, 0.0, 0.0, *color, float(alpha), *scale, *rotation, float(score)]
            fields = [f"{float(value):.7g}" for value in numeric]
            fields.extend((str(int(source_view)), str(int(source_mask))))
            handle.write(" ".join(fields) + "\n")