Source code for spindoctor.nav_model.nav_model_body_simulated

"""Simulated-body NavModel.

Renders a body from operator-supplied geometric parameters (centre, axes,
rotation, lighting) rather than from SPICE.  Used by the simulated-image
GUI to compose synthetic test scenes; the rendered body becomes a
``BODY_DISC`` ``NavFeature`` that the standard pipeline can navigate
against.  A well-resolved low-phase body also emits a ``LIMB_ARC``, and a
body at appreciable phase emits a ``TERMINATOR_ARC`` -- both matching the
SPICE-backed ``NavModelBody``'s geometry and gating semantics (polyline
conventions, emission gates, reliability inputs), so a sim scene exercises
the same techniques a real frame would.  The per-vertex sigma model is the
deliberate exception: the real model derives per-vertex sigmas from the
PSF, limb softness, and albedo terms, while the noise-free sim render
carries fixed per-vertex values (see ``_TERMINATOR_SIGMA_NORMAL_PX``).
"""

from __future__ import annotations

import math
from typing import TYPE_CHECKING, Any

import numpy as np
from oops import Observation

from spindoctor.annotation import Annotations
from spindoctor.config import Config
from spindoctor.feature.feature import NavFeature, NavReliabilityBreakdown
from spindoctor.feature.feature_type import NavFeatureType
from spindoctor.feature.flags import BodyDiscFlags, LimbArcFlags, TerminatorArcFlags
from spindoctor.feature.geometry import BodyDiscGeometry, LimbPolyline, TerminatorPolyline
from spindoctor.nav_model.body_shape import load_body_shape
from spindoctor.nav_model.nav_model import NavModel
from spindoctor.nav_model.nav_model_body import (
    LIMB_ARC_MIN_VERTICES,
    TERMINATOR_MIN_PHASE_FACTOR,
    TERMINATOR_MIN_VERTICES,
    limb_reliability,
    shape_features_suppressed,
    terminator_reliability,
)
from spindoctor.nav_model.nav_model_body_base import BODY_BLOB_MIN_DIAMETER_PX, NavModelBodyBase
from spindoctor.nav_model.sim_body import create_simulated_body
from spindoctor.sim.ellipsoid_geometry import DARK_SIDE_ILLUM_STRENGTH
from spindoctor.sim.mesh_geometry import mesh_spec_from_params, render_mesh_body_image
from spindoctor.support.filters import NavFilterKind, NavFilterSpec
from spindoctor.support.image import shift_array
from spindoctor.support.time import now_dt
from spindoctor.support.types import NDArrayBoolType, NDArrayFloatType

if TYPE_CHECKING:  # pragma: no cover - typing-only import
    from spindoctor.nav_orchestrator.nav_context import NavContext

__all__ = ['NavModelBodySimulated']

# Minimum limb-polyline vertices to emit a LIMB_ARC.  Matches BodyLimbNav's
# ``min_arc_vertices`` feasibility floor: a shorter arc cannot constrain the
# fit.  Imported from the catalog body model so the two emission gates cannot
# desync.
_MIN_LIMB_ARC_VERTICES: int = LIMB_ARC_MIN_VERTICES
# Minimum silhouette diameter to emit a LIMB_ARC.  The catalog body model gates
# the limb on its ellipsoid-fit uncertainty (``<= 3 px``) plus the shared
# vertex-count floor; the sim has no km scale, so it gates on resolution
# directly.  Below this the limb fit is imprecise and -- because it is an
# LM-refined distance-transform fit -- it injects cross-process jitter into the
# fused offset of any body scene, so only a well-resolved body emits a limb.
_MIN_LIMB_DIAMETER_PX: float = 100.0
# Per-vertex limb uncertainty for a simulated body.  The rendered silhouette edge
# is sharp and noise-free, so the predicted limb sits within ~1 px of the image
# edge; the tangent sigma reflects the one-pixel polyline sampling resolution.
_LIMB_SIGMA_NORMAL_PX: float = 1.0
_LIMB_SIGMA_TANGENT_PX: float = 0.5
# Above this phase the lit limb is a thin crescent and most of the silhouette
# boundary is terminator (a soft brightness gradient, not a sharp image edge), so
# a limb fit is unreliable -- the body is then navigated by the blob centroid.
# Gating emission here keeps LIMB_ARC off the high-phase scenes and matches the
# catalog body model's limb/blob handoff.
_LIMB_MAX_PHASE_DEG: float = 60.0
# Per-vertex terminator uncertainty for a simulated body.  The terminator is a
# soft photometric boundary (a shading ramp crossing zero) rather than a sharp
# silhouette edge, so its normal sigma is one pixel wider than the limb's; the
# tangent sigma reflects the one-pixel polyline sampling resolution.  The
# catalog body model derives these from the ellipsoid residual and albedo; the
# noise-free sim render has neither, so fixed values stand in.
_TERMINATOR_SIGMA_NORMAL_PX: float = 2.0
_TERMINATOR_SIGMA_TANGENT_PX: float = 0.5


def _limb_polyline_from_mask(
    limb_mask: NDArrayBoolType,
    body_mask: NDArrayBoolType,
) -> tuple[NDArrayFloatType, NDArrayFloatType]:
    """Extract limb vertices and outward normals from a 1-pixel limb mask.

    Each ``True`` pixel of ``limb_mask`` (a body pixel adjacent to sky) becomes
    one vertex.  The outward normal points away from the body interior: it is the
    normalized sum of the unit directions toward each non-body 4-neighbour, so a
    vertex on the sunward limb gets a normal pointing into the sky.

    Parameters:
        limb_mask: Extfov-shape boolean limb mask (the silhouette boundary).
        body_mask: Extfov-shape boolean body silhouette mask.

    Returns:
        ``(vertices_vu, normals_vu)`` each shaped ``(N, 2)``; empty when the mask
        has no pixels.
    """
    if not limb_mask.any():
        empty: NDArrayFloatType = np.empty((0, 2), dtype=np.float64)
        return empty, empty
    vs, us = np.where(limb_mask)
    vertices_vu = np.stack([vs.astype(np.float64), us.astype(np.float64)], axis=1)
    rows, cols = body_mask.shape
    normals_vu = np.zeros_like(vertices_vu)
    for i, (v, u) in enumerate(zip(vs, us, strict=True)):
        dv = 0.0
        du = 0.0
        if v > 0 and not body_mask[v - 1, u]:
            dv -= 1.0
        if v < rows - 1 and not body_mask[v + 1, u]:
            dv += 1.0
        if u > 0 and not body_mask[v, u - 1]:
            du -= 1.0
        if u < cols - 1 and not body_mask[v, u + 1]:
            du += 1.0
        norm = float(np.hypot(dv, du)) or 1.0
        normals_vu[i, 0] = dv / norm
        normals_vu[i, 1] = du / norm
    return vertices_vu, normals_vu


def _terminator_ridge_mask(
    model_img: NDArrayFloatType,
) -> tuple[NDArrayBoolType, NDArrayBoolType]:
    """Interior lit/unlit boundary ridge of a rendered body image.

    The shading floors the visible-but-unlit hemisphere at
    ``DARK_SIDE_ILLUM_STRENGTH``, so brightness ``> 0`` is the whole visible
    disc while brightness above the floor is the lit region; the unlit disc is
    their difference.  A ridge pixel is a lit pixel with an unlit *interior*
    disc pixel as a 4-neighbour -- the interior restriction drops the
    anti-aliased limb ring (unlit only because its edge brightness has ramped
    below the floor), which keeps the ridge off the silhouette everywhere the
    lit and unlit regions are separated by more than a pixel (the cusps of a
    very thin crescent are the exception; see ``_terminator_polyline``).

    Parameters:
        model_img: A rendered body image (any canvas).

    Returns:
        ``(ridge_mask, lit_mask)`` boolean arrays of the same shape.
    """
    disc = model_img > 0.0
    lit = model_img > DARK_SIDE_ILLUM_STRENGTH
    dark_disc = disc & ~lit
    touches_sky = (
        shift_array(~disc, (-1, 0))
        | shift_array(~disc, (1, 0))
        | shift_array(~disc, (0, -1))
        | shift_array(~disc, (0, 1))
    )
    dark_disc_interior = dark_disc & ~touches_sky
    neighbor_dark = (
        shift_array(dark_disc_interior, (-1, 0))
        | shift_array(dark_disc_interior, (1, 0))
        | shift_array(dark_disc_interior, (0, -1))
        | shift_array(dark_disc_interior, (0, 1))
    )
    ridge = lit & neighbor_dark
    return ridge, lit


def _silhouette_diameter_px(body_mask: NDArrayBoolType) -> float:
    """Return the longer pixel extent of a rendered body silhouette.

    Measures the bounding-box span of the lit silhouette in each axis and
    returns the larger of the two.  Returns ``0.0`` for an empty mask.
    """
    if not body_mask.any():
        return 0.0
    rows = np.any(body_mask, axis=1)
    cols = np.any(body_mask, axis=0)
    v_indices = np.where(rows)[0]
    u_indices = np.where(cols)[0]
    v_extent = float(v_indices[-1] - v_indices[0] + 1)
    u_extent = float(u_indices[-1] - u_indices[0] + 1)
    return max(v_extent, u_extent)


def _tight_bbox_extfov(
    body_mask: NDArrayBoolType,
    *,
    ext_margin_vu: tuple[int, int],
    extfov_shape: tuple[int, int],
    diameter_px: float,
) -> tuple[int, int, int, int]:
    """Return a tight extfov-coord bbox around a data-coord body silhouette.

    The bbox is the silhouette's bounding box inflated by a slop margin
    (so the body stays inside it under a modest pointing error) and
    clamped to the extfov shape.  Returns the full extfov frame when the
    mask is empty (a degenerate render the caller will not emit features
    for).

    Parameters:
        body_mask: Data-shape boolean silhouette mask.
        ext_margin_vu: Extfov margins ``(v, u)`` to shift data coords into
            extfov coords.
        extfov_shape: Extfov array shape ``(h, w)`` to clamp against.
        diameter_px: Predicted silhouette diameter, used to size the slop.

    Returns:
        ``(v_min, u_min, v_max, u_max)`` half-open bbox in extfov coords.
    """
    ext_margin_v, ext_margin_u = ext_margin_vu
    h, w = extfov_shape
    if not body_mask.any():
        return (0, 0, int(h), int(w))
    rows = np.where(np.any(body_mask, axis=1))[0]
    cols = np.where(np.any(body_mask, axis=0))[0]
    slop = max(round(0.1 * diameter_px), 4)
    v_min = int(rows[0]) + ext_margin_v - slop
    v_max = int(rows[-1]) + ext_margin_v + slop + 1
    u_min = int(cols[0]) + ext_margin_u - slop
    u_max = int(cols[-1]) + ext_margin_u + slop + 1
    return (
        max(0, v_min),
        max(0, u_min),
        min(int(h), v_max),
        min(int(w), u_max),
    )