pose3d

3D pose from the 2D: the [pose3d], [smoothing] and [correction] stages.

triangulation and pictorial implement the [pose3d] methods, smoothing is the ensemble Kalman smoother, and correction holds the anatomical corrections.

Robust multi-view triangulation over a rig: the inverse of Rig.project.

A NaN 2D observation means the view did not see the point, so visibility travels as NaNs, not a mask. pts2d is view-leading, (V, *pts, 2), and the points come back as (*pts, 3). triangulate_ransac runs its hypothesis search inside one jitted jax.lax.scan; the GT-honoring solve at the end is host-side numpy on one point at a time (ADR-0068).

consensus_support

consensus_support(
    inliers: Bool[ndarray, V * pts],
    conf: Float[ndarray, V * pts] | None = None,
) -> Float[np.ndarray, "*pts"]

How much a consensus rests on: the summed vote of the views it kept.

Each kept view counts its precision_from_conf, the vote of triangulate_ransac's weighted_vote (ADR-0129), so six sure views give 6 and two half-sure ones 0.5; without conf every kept view counts 1. Below LOW_SUPPORT the point is one to distrust (ADR-0132).

Source code in src/deeperfly/pose3d/triangulation.py

def consensus_support(
    inliers: Bool[np.ndarray, "V *pts"],
    conf: Float[np.ndarray, "V *pts"] | None = None,
) -> Float[np.ndarray, "*pts"]:
    """How much a consensus rests on: the summed vote of the views it kept.

    Each kept view counts its [precision_from_conf][], the vote of
    [triangulate_ransac][]'s `weighted_vote` (ADR-0129), so six sure views give 6 and
    two half-sure ones 0.5; without `conf` every kept view counts 1. Below
    `LOW_SUPPORT` the point is one to distrust (ADR-0132).
    """
    inliers = np.asarray(inliers, dtype=bool)
    votes = np.ones(inliers.shape) if conf is None else precision_from_conf(conf)
    return np.where(inliers, votes, 0.0).sum(axis=0)

depth_on_ray

depth_on_ray(
    rig: Rig,
    gt_view: int,
    gt_xy,
    stab_obs: Float[ndarray, "V 2"],
    *,
    huber_px: float = _STAB_SIGMA_PX,
    max_iter: int = 12,
) -> Float[np.ndarray, 3] | None

The point on the single GT view's ray whose depth the stabilizers agree on.

Only the depth is refit, by Huber IRLS at huber_px over max_iter steps, so the point reprojects onto the GT pixel (ADR-0068). stab_obs is (V, 2), NaN where none; None when no stabilizer fixes a depth.

Source code in src/deeperfly/pose3d/triangulation.py

def depth_on_ray(
    rig: Rig,
    gt_view: int,
    gt_xy,
    stab_obs: Float[np.ndarray, "V 2"],
    *,
    huber_px: float = _STAB_SIGMA_PX,
    max_iter: int = 12,
) -> Float[np.ndarray, "3"] | None:
    """The point on the single GT view's ray whose depth the stabilizers agree on.

    Only the depth is refit, by Huber IRLS at `huber_px` over `max_iter` steps, so the
    point reprojects onto the GT pixel (ADR-0068). `stab_obs` is `(V, 2)`, NaN where
    none; `None` when no stabilizer fixes a depth.
    """
    stab_obs = np.asarray(stab_obs, dtype=float)
    usable = np.nonzero(np.isfinite(stab_obs).all(axis=-1))[0]
    if usable.size == 0:
        return None
    (origin, direction), *_ = _rays(rig, [gt_view], [gt_xy])
    stab_rays = _rays(rig, usable, stab_obs[usable])
    c = float(huber_px)
    lam = _depth_from_rays(origin, direction, stab_rays)
    if lam is None:
        return None
    for _ in range(max_iter):
        x = origin + lam * direction
        err = _point_px_err(rig, x, usable, stab_obs)
        w = np.where(err <= c, 1.0, c / np.maximum(err, 1e-12))
        nxt = _depth_from_rays(origin, direction, stab_rays, weights=w)
        if nxt is None:
            break
        if abs(nxt - lam) <= 1e-9 * max(1.0, abs(lam)):
            lam = nxt
            break
        lam = nxt
    x = origin + lam * direction
    return x if np.all(np.isfinite(x)) else None

placed_in_depth

placed_in_depth(
    rig: Rig,
    pts3d: Float[ndarray, "*pts 3"],
    used: Bool[ndarray, V * pts],
) -> Bool[np.ndarray, "*pts"]

Whether some two of each point's used views see it MIN_PARALLAX_DEG apart.

And as far from opposite: two views looking along one line fix a point across it but not along it, so a point resting on such views alone has an arbitrary depth (ADR-0119, ADR-0130). False where the point is NaN.

Source code in src/deeperfly/pose3d/triangulation.py

def placed_in_depth(
    rig: Rig, pts3d: Float[np.ndarray, "*pts 3"], used: Bool[np.ndarray, "V *pts"]
) -> Bool[np.ndarray, "*pts"]:
    """Whether some two of each point's `used` views see it `MIN_PARALLAX_DEG` apart.

    And as far from opposite: two views looking along one line fix a point across
    it but not along it, so a point resting on such views alone has an arbitrary depth
    (ADR-0119, ADR-0130). False where the point is NaN.
    """
    pts3d = np.asarray(pts3d, dtype=float)
    used = np.asarray(used, dtype=bool)
    centers = view_centers(rig)
    rays = pts3d[None] - centers.reshape(-1, *([1] * (pts3d.ndim - 1)), 3)
    with np.errstate(invalid="ignore", divide="ignore"):
        rays = rays / np.linalg.norm(rays, axis=-1, keepdims=True)
    min_sin = np.sin(np.radians(MIN_PARALLAX_DEG))
    a, b = np.triu_indices(len(centers), k=1)  # every pair of views, at once
    sin = np.linalg.norm(np.cross(rays[a], rays[b]), axis=-1)
    return (used[a] & used[b] & (sin >= min_sin)).any(axis=0)

precision_from_conf

precision_from_conf(
    conf: Float[ndarray, V * pts],
    *,
    exponent: float = CONF_PRECISION_EXPONENT,
) -> Float[np.ndarray, V * pts]

Detector confidence as a triangulation weight, conf ** exponent (ADR-0067).

A non-finite confidence becomes 0, which drops the view.

Source code in src/deeperfly/pose3d/triangulation.py

def precision_from_conf(
    conf: Float[np.ndarray, "V *pts"], *, exponent: float = CONF_PRECISION_EXPONENT
) -> Float[np.ndarray, "V *pts"]:
    """Detector confidence as a triangulation weight, `conf ** exponent` (ADR-0067).

    A non-finite confidence becomes 0, which drops the view.
    """
    return np.nan_to_num(np.asarray(conf, dtype=float), nan=0.0) ** exponent

reprojection_error

reprojection_error(
    rig: Rig,
    pts3d: Float[ndarray, "*pts 3"],
    pts2d: Float[ndarray, "V *pts 2"],
) -> Float[np.ndarray, V * pts]

Per-(view, point) reprojection error in pixels, NaN where either side is NaN.

Source code in src/deeperfly/pose3d/triangulation.py

def reprojection_error(
    rig: Rig,
    pts3d: Float[np.ndarray, "*pts 3"],
    pts2d: Float[np.ndarray, "V *pts 2"],
) -> Float[np.ndarray, "V *pts"]:
    """Per-(view, point) reprojection error in pixels, NaN where either side is NaN."""
    err = _reproj_err(
        np.asarray(pts3d),
        rig.rvecs,
        rig.tvecs,
        rig.intrs,
        rig.dists,
        rig.mirrored,
        np.asarray(pts2d),
    )
    return np.asarray(err)

triangulate

triangulate(
    rig: Rig,
    pts2d: Float[ndarray, "V *pts 2"],
    weights: Float[ndarray, V * pts] | None = None,
    *,
    method: str = "fast",
) -> Float[np.ndarray, "*pts 3"]

Triangulate 3D points from 2D observations by NaN-aware DLT.

Fewer than two views give NaN; triangulate_ransac rejects gross outliers. weights scale each view's rows by sqrt(weight). method="fast" solves by inverse iteration and re-solves with eigh the rare points too ill-conditioned for it, matching "eigh" everywhere the triangulation is well posed (ADR-0043).

Source code in src/deeperfly/pose3d/triangulation.py

def triangulate(
    rig: Rig,
    pts2d: Float[np.ndarray, "V *pts 2"],
    weights: Float[np.ndarray, "V *pts"] | None = None,
    *,
    method: str = "fast",
) -> Float[np.ndarray, "*pts 3"]:
    """Triangulate 3D points from 2D observations by NaN-aware DLT.

    Fewer than two views give NaN; [triangulate_ransac][] rejects gross outliers.
    `weights` scale each view's rows by `sqrt(weight)`. `method="fast"` solves by
    inverse iteration and re-solves with `eigh` the rare points too ill-conditioned for
    it, matching `"eigh"` everywhere the triangulation is well posed (ADR-0043).
    """
    if method not in ("fast", "eigh"):
        raise ValueError(f"method must be 'fast' or 'eigh', got {method!r}")
    pts2d = np.asarray(pts2d)
    w = None if weights is None else np.asarray(weights)
    proj = pmats(rig.rvecs, rig.tvecs, rig.intrs, rig.mirrored)
    if method == "eigh":
        return np.asarray(triangulate_dlt(pts2d, proj, w, method="eigh"))

    pts3d, illcond = triangulate_dlt_flagged(pts2d, proj, w)
    pts3d, illcond = np.asarray(pts3d), np.asarray(illcond)
    if illcond.any():
        pts3d = _eigh_refine(pts3d, illcond, pts2d, proj, w)
    return pts3d

triangulate_ransac

triangulate_ransac(
    rig: Rig,
    pts2d: Float[ndarray, "V *pts 2"],
    *,
    threshold: float = 15.0,
    min_views: int = 2,
    weights: Float[ndarray, V * pts] | None = None,
    method: str = "fast",
    weighted_vote: bool = False,
    view_precision: Float[ndarray, " V"] | None = None,
) -> tuple[
    Float[np.ndarray, "*pts 3"], Bool[np.ndarray, V * pts]
]

Robustly triangulate 3D points, rejecting gross 2D outliers (RANSAC).

A rig has a handful of views, so every C(V, 2) two-view hypothesis is tried (the deterministic limit of RANSAC). Each candidate counts the views reprojecting within threshold px; the largest consensus wins, ties to the smaller inlier error, and the point is refit from its inliers. Fewer than min_views agreeing views give NaN. weights enter the candidate fits and the refit (ADR-0067); with weighted_vote they are the vote too, the consensus being the hypothesis whose agreeing views weigh the most rather than the one with the most views (ADR-0129). view_precision, (V,), scales each view's weight in the fits but not in the vote (view_precision). A pair of views within MIN_PARALLAX_DEG of parallel or opposite at its point forms no hypothesis, since it leaves the depth arbitrary (ADR-0130). A consensus can still end on such views alone, when the hypothesis' other view falls outside the gate, and its refit then has an arbitrary depth: a caller that needs the depth checks placed_in_depth on the inliers, as [pose3d] does (ADR-0135). method is the DLT solver's, as for triangulate.

Returns:

Type	Description
Float[ndarray, '*pts 3']	A tuple (pts3d, inliers): the points, NaN below min_views, and the
Bool[ndarray, V * pts]	(V, *pts) mask of the views kept per point.

Source code in src/deeperfly/pose3d/triangulation.py

def triangulate_ransac(
    rig: Rig,
    pts2d: Float[np.ndarray, "V *pts 2"],
    *,
    threshold: float = 15.0,
    min_views: int = 2,
    weights: Float[np.ndarray, "V *pts"] | None = None,
    method: str = "fast",
    weighted_vote: bool = False,
    view_precision: Float[np.ndarray, " V"] | None = None,
) -> tuple[Float[np.ndarray, "*pts 3"], Bool[np.ndarray, "V *pts"]]:
    """Robustly triangulate 3D points, rejecting gross 2D outliers (RANSAC).

    A rig has a handful of views, so every `C(V, 2)` two-view hypothesis is tried (the
    deterministic limit of RANSAC). Each candidate counts the views reprojecting
    within `threshold` px; the largest consensus wins, ties to the smaller inlier
    error, and the point is refit from its inliers. Fewer than `min_views` agreeing
    views give NaN. `weights` enter the candidate fits and the refit (ADR-0067); with
    `weighted_vote` they are the vote too, the consensus being the hypothesis whose
    agreeing views weigh the most rather than the one with the most views (ADR-0129).
    `view_precision`, `(V,)`, scales each view's weight in the fits but not in the vote
    ([view_precision][deeperfly.pose3d.triangulation.view_precision]). A pair of views
    within `MIN_PARALLAX_DEG` of parallel or opposite at its point forms no hypothesis,
    since it leaves the depth arbitrary (ADR-0130). A consensus can still end on such
    views alone, when the hypothesis' other view falls outside the gate, and its refit
    then has an arbitrary depth: a caller that needs the depth checks
    [placed_in_depth][] on the inliers, as `[pose3d]` does (ADR-0135). `method` is the
    DLT solver's, as for [triangulate][].

    Returns:
        A tuple `(pts3d, inliers)`: the points, NaN below `min_views`, and the
        `(V, *pts)` mask of the views kept per point.
    """
    if min_views < 2:
        raise ValueError(f"min_views must be >= 2, got {min_views}")
    if method not in ("fast", "eigh"):
        raise ValueError(f"method must be 'fast' or 'eigh', got {method!r}")
    pts2d = np.asarray(pts2d, dtype=float)
    w = None if weights is None else np.asarray(weights, dtype=float)
    votes = None
    if weighted_vote:
        votes = np.ones(pts2d.shape[:-1]) if w is None else w
    if view_precision is not None:
        per_view = np.asarray(view_precision, dtype=float)
        per_view = per_view.reshape(-1, *([1] * (pts2d.ndim - 2)))
        w = (np.ones(pts2d.shape[:-1]) if w is None else w) * per_view
    pts3d, inliers = _ransac_kernel(
        pts2d,
        view_centers(rig),
        rig.rvecs,
        rig.tvecs,
        rig.intrs,
        rig.dists,
        rig.mirrored,
        w,
        votes,
        float(threshold),
        int(min_views),
        "invit" if method == "fast" else "eigh",
    )
    return np.asarray(pts3d), np.asarray(inliers)

triangulate_stabilized

triangulate_stabilized(
    rig: Rig,
    x0: Float[ndarray, 3],
    gt_obs: Float[ndarray, "V 2"],
    stab_obs: Float[ndarray, "V 2"],
    *,
    gt_sigma_px: float = _GT_SIGMA_PX,
    huber_px: float = _STAB_SIGMA_PX,
    max_iter: int = 12,
) -> Float[np.ndarray, 3]

Two or more GT views, with the stabilizer supplying only what GT cannot see.

Gauss-Newton, for a fixed max_iter steps, on

sum_{GT}  |proj(x) - g|^2 / sigma_gt^2
  + sum_{stab} rho_c(|proj(x) - s|) / huber_px^2

with rho_c Huber at huber_px (ADR-0068). x0 is the GT-only solution, the start and the anchor, returned unchanged when no stabilizer is usable; gt_obs and stab_obs are (V, 2) raw pixels, NaN where none.

Source code in src/deeperfly/pose3d/triangulation.py

def triangulate_stabilized(
    rig: Rig,
    x0: Float[np.ndarray, "3"],
    gt_obs: Float[np.ndarray, "V 2"],
    stab_obs: Float[np.ndarray, "V 2"],
    *,
    gt_sigma_px: float = _GT_SIGMA_PX,
    huber_px: float = _STAB_SIGMA_PX,
    max_iter: int = 12,
) -> Float[np.ndarray, "3"]:
    """Two or more GT views, with the stabilizer supplying only what GT cannot see.

    Gauss-Newton, for a fixed `max_iter` steps, on

        sum_{GT}  |proj(x) - g|^2 / sigma_gt^2
          + sum_{stab} rho_c(|proj(x) - s|) / huber_px^2

    with `rho_c` Huber at `huber_px` (ADR-0068). `x0` is the GT-only solution, the
    start and the anchor, returned unchanged when no stabilizer is usable; `gt_obs` and
    `stab_obs` are `(V, 2)` raw pixels, NaN where none.
    """
    gt_obs = np.asarray(gt_obs, dtype=float)
    stab_obs = np.asarray(stab_obs, dtype=float)
    x = np.asarray(x0, dtype=float)
    gt_views = np.nonzero(np.isfinite(gt_obs).all(axis=-1))[0]
    stab_views = np.nonzero(np.isfinite(stab_obs).all(axis=-1))[0]
    if gt_views.size == 0 or stab_views.size == 0 or not np.all(np.isfinite(x)):
        return x

    c = float(huber_px)
    w_gt = 1.0 / max(float(gt_sigma_px), 1e-6) ** 2
    w_stab = 1.0 / max(c, 1e-6) ** 2
    eye = np.eye(3)
    for _ in range(max_iter):
        proj, jac = _proj_and_jac(rig, x)
        r_gt = proj[gt_views] - gt_obs[gt_views]
        j_gt = jac[gt_views]
        normal = w_gt * np.einsum("vij,vik->jk", j_gt, j_gt)
        grad = w_gt * np.einsum("vij,vi->j", j_gt, r_gt)
        r_stab = proj[stab_views] - stab_obs[stab_views]
        j_stab = jac[stab_views]
        err = np.linalg.norm(r_stab, axis=-1)
        # Huber IRLS: a stabilizer that agrees keeps its full vote, a grossly wrong one
        # keeps c/|r| of it -- down-weighted on its merits, never hard-rejected.
        wt = w_stab * np.where(err <= c, 1.0, c / np.maximum(err, 1e-12))
        normal += np.einsum("v,vij,vik->jk", wt, j_stab, j_stab)
        grad += np.einsum("v,vij,vi->j", wt, j_stab, r_stab)
        # One GT plus one stabilizer can leave the normal matrix singular; a
        # trace-scaled ridge regularizes that without biasing a healthy solve.
        ridge = 1e-12 * max(float(np.trace(normal)), 1e-30)
        step = -np.linalg.solve(normal + ridge * eye, grad)
        if not np.all(np.isfinite(step)):
            return np.asarray(x0, dtype=float)
        x = x + step
        if np.max(np.abs(step)) <= 1e-12:
            break
    return x if np.all(np.isfinite(x)) else np.asarray(x0, dtype=float)

view_centers

view_centers(rig: Rig) -> Float[np.ndarray, 'V 3']

Each view's center in world coordinates, in the rig's order.

Source code in src/deeperfly/pose3d/triangulation.py

def view_centers(rig: Rig) -> Float[np.ndarray, "V 3"]:
    """Each view's center in world coordinates, in the rig's order."""
    return np.stack([rig.views[name].position for name in rig.names])

view_precision

view_precision(
    rig: Rig,
    pts2d: Float[ndarray, "V *pts 2"],
    weights: Float[ndarray, V * pts] | None = None,
    *,
    threshold: float = 15.0,
    max_points: int = 20000,
) -> Float[np.ndarray, " V"]

Each view's precision relative to the steadiest view, (V,) in (0, 1].

(sigma_min / sigma_v) ** 2, with sigma_v the median residual of view v's inliers to the RANSAC consensus of pts2d (at most max_points points, spread evenly, weights as for triangulate_ransac). The noise level is measured, so a rig whose views differ in quality has its steadier views pull harder without naming them (ADR-0129). A view with no inlier takes the lowest precision measured; with nothing measured every view is 1.

Source code in src/deeperfly/pose3d/triangulation.py

def view_precision(
    rig: Rig,
    pts2d: Float[np.ndarray, "V *pts 2"],
    weights: Float[np.ndarray, "V *pts"] | None = None,
    *,
    threshold: float = 15.0,
    max_points: int = 20_000,
) -> Float[np.ndarray, " V"]:
    """Each view's precision relative to the steadiest view, `(V,)` in `(0, 1]`.

    `(sigma_min / sigma_v) ** 2`, with `sigma_v` the median residual of view `v`'s
    inliers to the RANSAC consensus of `pts2d` (at most `max_points` points, spread
    evenly, `weights` as for [triangulate_ransac][]). The noise level is measured, so a
    rig whose views differ in quality has its steadier views pull harder without
    naming them (ADR-0129). A view with no inlier takes the lowest precision measured;
    with nothing measured every view is 1.
    """
    pts2d = np.asarray(pts2d, dtype=float)
    n_views = pts2d.shape[0]
    flat = pts2d.reshape(n_views, -1, 2)
    w = None if weights is None else np.asarray(weights, dtype=float)
    w = None if w is None else w.reshape(n_views, -1)
    if flat.shape[1] > max_points:
        take = np.unique(np.linspace(0, flat.shape[1] - 1, max_points).round())
        take = take.astype(int)
        flat = flat[:, take]
        w = None if w is None else w[:, take]
    x, inliers = triangulate_ransac(rig, flat, threshold=threshold, weights=w)
    err = reprojection_error(rig, x, flat)
    sigma = np.array(
        [
            np.median(err[v][inliers[v]]) if inliers[v].any() else np.nan
            for v in range(n_views)
        ]
    )
    if not np.isfinite(sigma).any():
        return np.ones(n_views)
    sigma = np.maximum(np.where(np.isfinite(sigma), sigma, np.nanmax(sigma)), 1e-6)
    return (sigma.min() / sigma) ** 2

Pictorial structures: one 3D per keypoint, chosen among the candidate peaks.

When a view's strongest peak is the wrong one (a crossed leg, a left/right swap), the right answer is often its second or third, and a selection by multi-view consistency recovers what a RANSAC veto cannot (Gunel et al., DeepFly3D, 2019; ADR-0038). solve triangulates every candidate pair of every view pair (700 hypotheses at eight views and five candidates, 11.4 ms of a frame's 25.8), scores each by its views' Gaussian-weighted support, prunes and refits them, and picks one per keypoint by exact max-product over the skeleton's forest with a Huber bone-length prior. elect is the cheaper leave-one-out rule. A veto keeps a weakly calibrated view out of the geometry while its 2D is still re-selected.

PictorialResult dataclass

What the selection committed to, per frame.

Attributes:

Name	Type	Description
points	Float[ndarray, 'T K 3']	The chosen 3D, (T, K, 3); NaN where no hypothesis reached the cell.
points2d	Float[ndarray, 'V T K 2']	Per view the candidate that supported it, (V, T, K, 2): NaN where the view did not, the decoded 2D where the keypoint has no 3D.
conf2d	Float[ndarray, 'V T K']	That candidate's score, (V, T, K); 0 where unsupported.
inliers	Bool[ndarray, 'V T K']	Which views supported the chosen hypothesis, (V, T, K).

Source code in src/deeperfly/pose3d/pictorial.py

@dataclass(frozen=True)
class PictorialResult:
    """What the selection committed to, per frame.

    Attributes:
        points: The chosen 3D, `(T, K, 3)`; NaN where no hypothesis reached the cell.
        points2d: Per view the candidate that supported it, `(V, T, K, 2)`: NaN where
            the view did not, the decoded 2D where the keypoint has no 3D.
        conf2d: That candidate's score, `(V, T, K)`; 0 where unsupported.
        inliers: Which views supported the chosen hypothesis, `(V, T, K)`.
    """

    points: Float[np.ndarray, "T K 3"]
    points2d: Float[np.ndarray, "V T K 2"]
    conf2d: Float[np.ndarray, "V T K"]
    inliers: Bool[np.ndarray, "V T K"]

bone_targets

bone_targets(
    rig: Rig,
    pts2d: Float[ndarray, "V T K 2"],
    edges: Int[ndarray, "E 2"],
    *,
    frames: int = 200,
    threshold_px: float = 15.0,
    min_views: int = 2,
) -> tuple[Float[np.ndarray, " E"], float]

The prior: each edge's median 3D length over pts2d's frames, and the scale.

Up to frames frames of the decoded 2D, triangulated robustly (threshold_px, min_views); NaN for an edge never triangulated. The scale is the median of the medians, the Huber residual's and the suppression radius's unit, or a tenth of the pose's median extent with no edge measured. Returns (targets, scale).

Source code in src/deeperfly/pose3d/pictorial.py

def bone_targets(
    rig: Rig,
    pts2d: Float[np.ndarray, "V T K 2"],
    edges: Int[np.ndarray, "E 2"],
    *,
    frames: int = 200,
    threshold_px: float = 15.0,
    min_views: int = 2,
) -> tuple[Float[np.ndarray, " E"], float]:
    """The prior: each edge's median 3D length over `pts2d`'s frames, and the scale.

    Up to `frames` frames of the decoded 2D, triangulated robustly (`threshold_px`,
    `min_views`); NaN for an edge never triangulated. The scale is the median of the
    medians, the Huber residual's and the suppression radius's unit, or a tenth of the
    pose's median extent with no edge measured. Returns `(targets, scale)`.
    """
    pts2d = np.asarray(pts2d, dtype=float)
    edges = np.asarray(edges, dtype=int).reshape(-1, 2)
    t = pts2d.shape[1]
    sel = np.unique(np.linspace(0, t - 1, min(frames, t)).round().astype(int))
    x, _ = triangulate_ransac(
        rig, pts2d[:, sel], threshold=threshold_px, min_views=min_views
    )
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", RuntimeWarning)  # all-NaN medians are answers
        if len(edges):
            lengths = np.linalg.norm(x[:, edges[:, 0]] - x[:, edges[:, 1]], axis=-1)
            targets = np.nanmedian(lengths, axis=0)
        else:
            targets = np.zeros(0)
        scale = float(np.nanmedian(targets)) if np.isfinite(targets).any() else np.nan
        if not np.isfinite(scale):
            extent = np.nanmax(x, axis=1) - np.nanmin(x, axis=1)
            scale = 0.1 * float(np.nanmedian(np.linalg.norm(extent, axis=-1)))
    return targets, scale if np.isfinite(scale) and scale > 0 else 1.0

edge_forest

edge_forest(
    edges: Int[ndarray, "E 2"], n_points: int
) -> tuple[Int[np.ndarray, " K"], list[int]]

Root the skeleton's edges: each keypoint's parent (-1 at a root) and an order.

Each component is rooted at its lowest index and walked breadth first, parents before children, which keeps the max-product pass exact.

Raises:

Type	Description
ValueError	If the edges close a cycle; the exact pass needs a forest.

Source code in src/deeperfly/pose3d/pictorial.py

def edge_forest(
    edges: Int[np.ndarray, "E 2"], n_points: int
) -> tuple[Int[np.ndarray, " K"], list[int]]:
    """Root the skeleton's edges: each keypoint's parent (-1 at a root) and an order.

    Each component is rooted at its lowest index and walked breadth first, parents
    before children, which keeps the max-product pass exact.

    Raises:
        ValueError: If the edges close a cycle; the exact pass needs a forest.
    """
    edges = np.asarray(edges, dtype=int).reshape(-1, 2)
    adjacency: list[list[int]] = [[] for _ in range(n_points)]
    for a, b in edges:
        adjacency[int(a)].append(int(b))
        adjacency[int(b)].append(int(a))
    parents = np.full(n_points, -1, dtype=int)
    order: list[int] = []
    seen = np.zeros(n_points, dtype=bool)
    components = 0
    for root in range(n_points):
        if seen[root]:
            continue
        components += 1
        seen[root] = True
        queue = [root]
        while queue:
            node = queue.pop(0)
            order.append(node)
            for nxt in adjacency[node]:
                if not seen[nxt]:
                    seen[nxt] = True
                    parents[nxt] = node
                    queue.append(nxt)
    if len(edges) != n_points - components:
        raise ValueError(
            "the skeleton's edges close a cycle (or repeat an edge); the pictorial "
            "selection needs a forest of chains"
        )
    return parents, order

elect

elect(
    rig: Rig,
    cand_xy: Float[ndarray, "V T K C 2"],
    cand_score: Float[ndarray, "V T K C"],
    *,
    fallback_xy: Float[ndarray, "V T K 2"],
    fallback_conf: Float[ndarray, "V T K"],
    threshold_px: float = 15.0,
    min_views: int = 2,
    reach: float = ELECTION_REACH,
    veto=None,
) -> PictorialResult

Re-elect each view's candidate by the other views' decode: leave-one-out.

For view v the other voting views' decode is triangulated (RANSAC at threshold_px, min_views inliers) and reprojected into v, and each keypoint takes its nearest candidate within reach * threshold_px, a bound on the geometry's error, not on the move (the swaps worth recovering sit 100 to 300 px from the decode). An incumbent that is nearest stays, a cell the others cannot place keeps its decode, and the result is triangulated once more. No hypothesis pool, so flat in the candidate count, but the higher-variance estimator: on the old repo's held-out animals it held 14 to 69% of the available gain where solve held 53 to 69%, so it is chosen for cost, never the default. A vetoed view is still re-elected.

Raises:

Type	Description
ValueError	On min_views below 2, or a veto leaving fewer than two views.

Source code in src/deeperfly/pose3d/pictorial.py

def elect(
    rig: Rig,
    cand_xy: Float[np.ndarray, "V T K C 2"],
    cand_score: Float[np.ndarray, "V T K C"],
    *,
    fallback_xy: Float[np.ndarray, "V T K 2"],
    fallback_conf: Float[np.ndarray, "V T K"],
    threshold_px: float = 15.0,
    min_views: int = 2,
    reach: float = ELECTION_REACH,
    veto=None,
) -> PictorialResult:
    """Re-elect each view's candidate by the other views' decode: leave-one-out.

    For view `v` the other voting views' decode is triangulated (RANSAC at
    `threshold_px`, `min_views` inliers) and reprojected into `v`, and each keypoint
    takes its nearest candidate within `reach * threshold_px`, a bound on the
    geometry's error, not on the move (the swaps worth recovering sit 100 to 300 px
    from the decode). An incumbent that is nearest stays, a cell the others cannot
    place keeps its decode, and the result is triangulated once more. No hypothesis
    pool, so flat in the candidate count, but the higher-variance estimator: on the old
    repo's held-out animals it held 14 to 69% of the available gain where [solve][]
    held 53 to 69%, so it is chosen for cost, never the default. A vetoed view is still
    re-elected.

    Raises:
        ValueError: On `min_views` below 2, or a veto leaving fewer than two views.
    """
    if min_views < 2:
        raise ValueError(f"min_views must be >= 2, got {min_views}")
    cand_xy = np.asarray(cand_xy, dtype=float)
    cand_score = np.asarray(cand_score, dtype=float)
    points2d = np.array(fallback_xy, dtype=float)
    conf2d = np.array(fallback_conf, dtype=float)
    n_v = cand_xy.shape[0]
    voting = ~_veto_mask(veto, n_v)
    if voting.sum() < 2:
        raise ValueError("veto leaves fewer than two views to triangulate with")
    seen = np.where(voting[:, None, None, None], points2d, np.nan)
    limit = reach * threshold_px
    for v in range(n_v):
        held = seen.copy()
        held[v] = np.nan  # a view must not vote for itself
        x3, inl = triangulate_ransac(
            rig, held, threshold=threshold_px, min_views=min_views
        )
        proj = np.asarray(rig.views[rig.names[v]].project(np.asarray(x3)))
        d = np.linalg.norm(cand_xy[v] - proj[..., None, :], axis=-1)  # (T, K, C)
        d = np.where(np.isfinite(d), d, np.inf)
        best = d.argmin(-1)
        ok = (
            (np.asarray(inl).sum(0) >= min_views)
            & np.isfinite(proj).all(-1)
            & (d.min(-1) < limit)
            & (best != 0)
        )
        ti, ki = np.nonzero(ok)
        points2d[v, ti, ki] = cand_xy[v, ti, ki, best[ok]]
        conf2d[v, ti, ki] = cand_score[v, ti, ki, best[ok]]
    pts3d, inliers = triangulate_ransac(
        rig,
        np.where(voting[:, None, None, None], points2d, np.nan),
        threshold=threshold_px,
        min_views=min_views,
    )
    inliers = np.asarray(inliers) & voting[:, None, None]
    return PictorialResult(np.asarray(pts3d), points2d, conf2d, inliers)

solve

solve(
    rig: Rig,
    cand_xy: Float[ndarray, "V T K C 2"],
    cand_score: Float[ndarray, "V T K C"],
    edges: Int[ndarray, "E 2"],
    *,
    fallback_xy: Float[ndarray, "V T K 2"],
    fallback_conf: Float[ndarray, "V T K"],
    threshold_px: float = 15.0,
    min_views: int = 2,
    weighted: bool = False,
    bone_prior: float = 0.0,
    targets: Float[ndarray, " E"] | None = None,
    scale: float | None = None,
    max_hypotheses: int = DEFAULT_MAX_HYPOTHESES,
    veto=None,
) -> PictorialResult

Choose one 3D per keypoint and frame among the candidates' triangulations.

cand_xy (V, T, K, C, 2) and cand_score (V, T, K, C) are the candidates, slot 0 the decode; fallback_xy and fallback_conf are the decoded 2D, what the prior is fitted from and what an unreached keypoint keeps. A view supports a hypothesis within threshold_px; a hypothesis needs min_views support; weighted scales a vote by the candidate's score; bone_prior weighs the Huber bone cost against one unit of evidence per view, over targets and scale (fitted when omitted); at most max_hypotheses are kept per keypoint. A veto view (a mask or indices) seeds, supports and refits nothing but is committed to its nearest candidate.

Raises:

Type	Description
ValueError	On min_views below 2, edges that close a cycle, or a veto leaving fewer than two views.

Source code in src/deeperfly/pose3d/pictorial.py

def solve(
    rig: Rig,
    cand_xy: Float[np.ndarray, "V T K C 2"],
    cand_score: Float[np.ndarray, "V T K C"],
    edges: Int[np.ndarray, "E 2"],
    *,
    fallback_xy: Float[np.ndarray, "V T K 2"],
    fallback_conf: Float[np.ndarray, "V T K"],
    threshold_px: float = 15.0,
    min_views: int = 2,
    weighted: bool = False,
    bone_prior: float = 0.0,
    targets: Float[np.ndarray, " E"] | None = None,
    scale: float | None = None,
    max_hypotheses: int = DEFAULT_MAX_HYPOTHESES,
    veto=None,
) -> PictorialResult:
    """Choose one 3D per keypoint and frame among the candidates' triangulations.

    `cand_xy` `(V, T, K, C, 2)` and `cand_score` `(V, T, K, C)` are the candidates, slot
    0 the decode; `fallback_xy` and `fallback_conf` are the decoded 2D, what the prior
    is fitted from and what an unreached keypoint keeps. A view supports a hypothesis
    within `threshold_px`; a hypothesis needs `min_views` support; `weighted` scales a
    vote by the candidate's score; `bone_prior` weighs the Huber bone cost against one
    unit of evidence per view, over `targets` and `scale` (fitted when omitted); at
    most `max_hypotheses` are kept per keypoint. A `veto` view (a mask or indices)
    seeds, supports and refits nothing but is committed to its nearest candidate.

    Raises:
        ValueError: On `min_views` below 2, edges that close a cycle, or a veto
            leaving fewer than two views.
    """
    if min_views < 2:
        raise ValueError(f"min_views must be >= 2, got {min_views}")
    cand_xy = np.asarray(cand_xy, dtype=float)
    cand_score = np.asarray(cand_score, dtype=float)
    fallback_xy = np.asarray(fallback_xy, dtype=float)
    fallback_conf = np.asarray(fallback_conf, dtype=float)
    edges = np.asarray(edges, dtype=int).reshape(-1, 2)
    n_v, n_t, n_k, n_c, _ = cand_xy.shape
    voting = ~_veto_mask(veto, n_v)
    if voting.sum() < 2:
        raise ValueError("veto leaves fewer than two views to triangulate with")
    seen = np.where(voting[:, None, None, None], fallback_xy, np.nan)
    if targets is None or scale is None:
        fit_targets, fit_scale = bone_targets(
            rig, seen, edges, threshold_px=threshold_px, min_views=min_views
        )
        targets = fit_targets if targets is None else np.asarray(targets, dtype=float)
        scale = fit_scale if scale is None else float(scale)
    parents, order = edge_forest(edges, n_k)
    edge_target = np.full(n_k, np.nan)
    for (a, b), length in zip(edges, targets, strict=True):
        child = int(b) if parents[b] == a else int(a)
        edge_target[child] = length

    proj = np.asarray(pmats(rig.rvecs, rig.tvecs, rig.intrs, rig.mirrored))
    pairs = list(combinations(np.flatnonzero(voting).tolist(), 2))
    n_m = len(pairs) * n_c * n_c
    block = max(1, int(BLOCK_BUDGET // max(1, n_v * n_k * n_m * n_c)))

    points = np.full((n_t, n_k, 3), np.nan)
    points2d = np.full((n_v, n_t, n_k, 2), np.nan)
    conf2d = np.zeros((n_v, n_t, n_k))
    inliers = np.zeros((n_v, n_t, n_k), dtype=bool)
    for f0 in range(0, n_t, block):
        f1 = min(n_t, f0 + block)
        cxy, csc = cand_xy[:, f0:f1], cand_score[:, f0:f1]
        hyp = _hypotheses(proj, pairs, cxy)
        nearest, close, evidence, n_in = _support(
            rig, hyp, cxy, csc, threshold_px, weighted, voting
        )
        inlier = close & voting[:, None, None, None]
        kept = _prune(
            hyp,
            evidence,
            n_in,
            min_views=min_views,
            radius=NMS_FRACTION * scale,
            max_hyp=max_hypotheses,
        )
        n_f = f1 - f0
        fi = np.arange(n_f)[:, None, None]
        ki = np.arange(n_k)[None, :, None]
        valid = kept >= 0
        kept0 = np.where(valid, kept, 0)
        obs = np.full((n_v, n_f, n_k, kept.shape[-1], 2), np.nan)
        for v in range(n_v):
            supports = inlier[v][fi, ki, kept0] & valid
            near = nearest[v][fi, ki, kept0]
            obs[v] = np.where(supports[..., None], cxy[v][fi, ki, near], np.nan)
        refit = triangulate(rig, obs)
        nearest, close, evidence, n_in = _support(
            rig, refit, cxy, csc, threshold_px, weighted, voting
        )
        valid &= (n_in >= min_views) & np.isfinite(refit).all(-1)
        state = _select(
            refit,
            -evidence,
            valid,
            parents,
            order,
            edge_target,
            lam=bone_prior,
            scale=scale,
        )
        chosen = state >= 0
        st = np.where(chosen, state, 0)[..., None]
        points[f0:f1] = np.where(
            chosen[..., None], refit[fi, ki, st][..., 0, :], np.nan
        )
        for v in range(n_v):
            # A view commits to its candidate nearest the chosen 3D whether or not
            # it voted; only a voting view counts as an inlier.
            committed = close[v][fi, ki, st][..., 0] & chosen
            near = nearest[v][fi, ki, st][..., 0]
            xy = cxy[v][fi[..., 0], ki[..., 0], near]
            points2d[v, f0:f1] = np.where(committed[..., None], xy, np.nan)
            conf2d[v, f0:f1] = np.where(
                committed, csc[v][fi[..., 0], ki[..., 0], near], 0.0
            )
            inliers[v, f0:f1] = committed & voting[v]
    unsolved = ~np.isfinite(points).all(-1)  # (T, K)
    points2d[:, unsolved] = fallback_xy[:, unsolved]
    conf2d[:, unsolved] = fallback_conf[:, unsolved]
    return PictorialResult(points, points2d, conf2d, inliers)

The ensemble Kalman smoother: multi-view 2D in, one smoothed 3D track per keypoint.

The nonlinear, geometric EKS of Lightning Pose 3D (Aharon, Whiteway et al., 2026), reimplemented on project_full_one (ADR-0062). The latent state is a 3D point observed through the rig's projection, and a cross-view variance inflation down-weights the views that disagree.

EksResult dataclass

What the smoother produced, in deeperfly's array conventions.

Attributes:

Name	Type	Description
pts3d	Float[ndarray, 'T P 3']	The smoothed 3D; NaN for a keypoint two views never saw at once.
pts2d	Float[ndarray, 'V T P 2']	pts3d reprojected; NaN in the unobserved cells unless filled.
posterior_var	Float[ndarray, 'T P 3']	The per-axis posterior variance, in world units squared.
smooth_param	Float[ndarray, P]	The process-noise scale per keypoint; small smooths hard.
n_inflated	int	Observations the variance inflation down-weighted.
n_testable	int	Observations the inflation could test at all.

Source code in src/deeperfly/pose3d/smoothing/smoother.py

@dataclass(frozen=True)
class EksResult:
    """What the smoother produced, in deeperfly's array conventions.

    Attributes:
        pts3d: The smoothed 3D; NaN for a keypoint two views never saw at once.
        pts2d: `pts3d` reprojected; NaN in the unobserved cells unless filled.
        posterior_var: The per-axis posterior variance, in world units squared.
        smooth_param: The process-noise scale per keypoint; small smooths hard.
        n_inflated: Observations the variance inflation down-weighted.
        n_testable: Observations the inflation could test at all.
    """

    pts3d: Float[np.ndarray, "T P 3"]
    pts2d: Float[np.ndarray, "V T P 2"]
    posterior_var: Float[np.ndarray, "T P 3"]
    smooth_param: Float[np.ndarray, "P"]
    n_inflated: int = 0
    n_testable: int = 0

ensemble_statistics

ensemble_statistics(
    pts2d: Float[ndarray, "M V T P 2"],
    conf: Float[ndarray, "M V T P"] | None = None,
    *,
    avg_mode: str = "median",
    var_mode: str = "confidence_weighted_var",
) -> tuple[
    Float[np.ndarray, "V T P 2"],
    Float[np.ndarray, "V T P 2"],
]

Collapse an ensemble of detectors into a center and a per-axis variance in px^2.

With one member the variance is 1 / confidence, a prior rather than a measurement (ADR-0062); conf=None treats every observation as equally confident. avg_mode is "median" (robust to one member blowing up) or "mean"; var_mode is "var", or "confidence_weighted_var" to divide the variance by the mean confidence.

Source code in src/deeperfly/pose3d/smoothing/ensemble.py

def ensemble_statistics(
    pts2d: Float[np.ndarray, "M V T P 2"],
    conf: Float[np.ndarray, "M V T P"] | None = None,
    *,
    avg_mode: str = "median",
    var_mode: str = "confidence_weighted_var",
) -> tuple[Float[np.ndarray, "V T P 2"], Float[np.ndarray, "V T P 2"]]:
    """Collapse an ensemble of detectors into a center and a per-axis variance in px^2.

    With one member the variance is `1 / confidence`, a prior rather than a
    measurement (ADR-0062); `conf=None` treats every observation as equally
    confident. `avg_mode` is `"median"` (robust to one member blowing up) or
    `"mean"`; `var_mode` is `"var"`, or `"confidence_weighted_var"` to divide the
    variance by the mean confidence.
    """
    pts2d = np.asarray(pts2d, dtype=float)
    if pts2d.ndim != 5:
        raise ValueError(f"pts2d must be (M, V, T, P, 2); got {pts2d.shape}")
    if avg_mode not in ("median", "mean"):
        raise ValueError(f"avg_mode must be 'median' or 'mean'; got {avg_mode!r}")
    if var_mode not in ("var", "confidence_weighted_var"):
        raise ValueError(
            f"var_mode must be 'var' or 'confidence_weighted_var'; got {var_mode!r}"
        )
    n_models = pts2d.shape[0]

    with warnings.catch_warnings():  # all-NaN cells are expected (unobserved)
        warnings.simplefilter("ignore", RuntimeWarning)
        center = (np.nanmedian if avg_mode == "median" else np.nanmean)(pts2d, axis=0)
        if conf is None:
            mean_conf = np.ones(pts2d.shape[1:4], dtype=float)
        else:
            mean_conf = np.asarray(conf, dtype=float).mean(axis=0)
        if n_models == 1:
            single = 1.0 / np.maximum(mean_conf, _MIN_CONF)
            var = np.repeat(single[..., None], 2, axis=-1)
        else:
            var = np.nanvar(pts2d, axis=0)
            if var_mode == "confidence_weighted_var":
                var = var / np.maximum(mean_conf, _MIN_CONF)[..., None]

    var = np.where(np.isfinite(var), var, _UNDEFINED_VAR)
    return center, var

inflate_variances

inflate_variances(
    center: Float[ndarray, "V T P 2"],
    var: Float[ndarray, "V T P 2"],
    mask: Bool[ndarray, "V T P"],
    pts3d: Float[ndarray, "T P 3"],
    rig: Rig,
    *,
    threshold: float = 5.0,
    factor: float = 10.0,
    n_rounds: int = 8,
) -> tuple[Float[np.ndarray, "V T P 2"], int]

Inflate the variance of each view that disagrees with the cross-view consensus.

pts3d is the linearization point; a cell without a finite one, or with fewer than two observed views, is not tested. A view over the Mahalanobis threshold has its variance multiplied by factor, for up to n_rounds rounds (defaults: ADR-0062). Returns the variances and how many (view, frame, point) observations were inflated.

Source code in src/deeperfly/pose3d/smoothing/ensemble.py

def inflate_variances(
    center: Float[np.ndarray, "V T P 2"],
    var: Float[np.ndarray, "V T P 2"],
    mask: Bool[np.ndarray, "V T P"],
    pts3d: Float[np.ndarray, "T P 3"],
    rig: Rig,
    *,
    threshold: float = 5.0,
    factor: float = 10.0,
    n_rounds: int = 8,
) -> tuple[Float[np.ndarray, "V T P 2"], int]:
    """Inflate the variance of each view that disagrees with the cross-view consensus.

    `pts3d` is the linearization point; a cell without a finite one, or with fewer
    than two observed views, is not tested. A view over the Mahalanobis `threshold`
    has its variance multiplied by `factor`, for up to `n_rounds` rounds (defaults:
    ADR-0062). Returns the variances and how many `(view, frame, point)`
    observations were inflated.
    """
    n_views, n_frames, n_points = mask.shape
    finite3d = np.isfinite(np.asarray(pts3d, dtype=float)).all(axis=-1)  # (T, P)
    # An untestable cell gets a placeholder point, so no NaN enters the kernel.
    with warnings.catch_warnings():  # a never-triangulated point is all-NaN
        warnings.simplefilter("ignore", RuntimeWarning)
        fallback = np.nanmedian(np.asarray(pts3d, dtype=float), axis=0)  # (P, 3)
    fallback = np.where(np.isfinite(fallback), fallback, 0.0)
    safe3d = np.where(finite3d[..., None], pts3d, fallback[None])

    x = stack_views(center)  # (T, P, 2V)
    v = stack_views(var)
    m = np.repeat(np.transpose(mask, (1, 2, 0)), 2, axis=-1)  # (T, P, 2V)
    # Two observed views is the minimum at which "the others" is a real opinion.
    testable = finite3d & (mask.sum(axis=0) >= 2)  # (T, P)

    rig_arrays = (
        jnp.asarray(rig.rvecs),
        jnp.asarray(rig.tvecs),
        jnp.asarray(rig.intrs),
        jnp.asarray(rig.dists),
        jnp.asarray(rig.mirrored),
    )
    kernel = jax.jit(
        jax.vmap(
            _cell_inflation,
            in_axes=(0, 0, 0, 0, 0, None, None, None, None),
        ),
        static_argnums=(6, 7, 8),
    )
    flat = (n_frames * n_points,)
    out = kernel(
        jnp.asarray(x.reshape(*flat, 2 * n_views)),
        jnp.asarray(v.reshape(*flat, 2 * n_views)),
        jnp.asarray(m.reshape(*flat, 2 * n_views)),
        jnp.asarray(safe3d.reshape(*flat, 3)),
        jnp.asarray(testable.reshape(flat)),
        rig_arrays,
        float(threshold),
        float(factor),
        int(n_rounds),
    )
    inflated = np.asarray(out).reshape(n_frames, n_points, 2 * n_views)
    # Count views, not dimensions: both of a view's axes inflate together.
    n_inflated = int(np.count_nonzero(inflated[..., ::2] > v[..., ::2] * (1.0 + 1e-9)))
    return unstack_views(inflated, n_views), n_inflated

smooth

smooth(
    rig: Rig,
    pts2d: Float[ndarray, "V T P 2"]
    | Float[ndarray, "M V T P 2"],
    conf: Float[ndarray, "V T P"]
    | Float[ndarray, "M V T P"]
    | None = None,
    *,
    init3d: Float[ndarray, "T P 3"] | None = None,
    smooth_param: float | Float[ndarray, P] | None = None,
    avg_mode: str = "median",
    var: str = "confidence_weighted_var",
    inflate: bool = True,
    inflate_threshold: float = 5.0,
    inflate_factor: float = 10.0,
    fit_frames: int = 2000,
    fit_iterations: int = 24,
    fill_unobserved: bool = False,
) -> EksResult

Smooth a run's multi-view 2D into one latent 3D trajectory per keypoint.

pts2d is one detector's 2D or an ensemble's, with a leading member axis. The rig's projection is the observation model, so pass the rig the run uses (bundle-adjusted when that stage is on). init3d defaults to a DLT triangulation of the ensemble center. smooth_param=None fits one scale per keypoint on the first fit_frames frames (0 = all); smaller smooths harder. avg_mode and var are the modes of ensemble_statistics. fill_unobserved also reports the reprojection in the cells nothing observed (ADR-0062).

Source code in src/deeperfly/pose3d/smoothing/smoother.py

def smooth(
    rig: Rig,
    pts2d: Float[np.ndarray, "V T P 2"] | Float[np.ndarray, "M V T P 2"],
    conf: Float[np.ndarray, "V T P"] | Float[np.ndarray, "M V T P"] | None = None,
    *,
    init3d: Float[np.ndarray, "T P 3"] | None = None,
    smooth_param: float | Float[np.ndarray, "P"] | None = None,
    avg_mode: str = "median",
    var: str = "confidence_weighted_var",
    inflate: bool = True,
    inflate_threshold: float = 5.0,
    inflate_factor: float = 10.0,
    fit_frames: int = 2000,
    fit_iterations: int = 24,
    fill_unobserved: bool = False,
) -> EksResult:
    """Smooth a run's multi-view 2D into one latent 3D trajectory per keypoint.

    `pts2d` is one detector's 2D or an ensemble's, with a leading member axis. The
    rig's projection is the observation model, so pass the rig the run uses
    (bundle-adjusted when that stage is on). `init3d` defaults to a DLT triangulation
    of the ensemble center. `smooth_param=None` fits one scale per keypoint on the
    first `fit_frames` frames (0 = all); smaller smooths harder. `avg_mode` and `var`
    are the modes of [ensemble_statistics][]. `fill_unobserved` also reports the
    reprojection in the cells nothing observed (ADR-0062).
    """
    pts2d = np.asarray(pts2d, dtype=float)
    if pts2d.ndim == 4:
        pts2d = pts2d[None]
    if pts2d.ndim != 5 or pts2d.shape[-1] != 2:
        raise ValueError(
            f"pts2d must be (V, T, P, 2) or (M, V, T, P, 2); got {pts2d.shape}"
        )
    _, n_views, n_frames, n_points, _ = pts2d.shape
    if n_views != len(rig.views):
        raise ValueError(
            f"pts2d has {n_views} views but the rig has {len(rig.views)} views"
        )
    if n_frames < 2:
        raise ValueError("the smoother needs at least two frames")
    if conf is not None:
        conf = np.asarray(conf, dtype=float)
        if conf.ndim == 3:
            conf = conf[None]
        if conf.shape != pts2d.shape[:4]:
            raise ValueError(f"conf must be {pts2d.shape[:4]}; got {conf.shape}")

    center, var2d = ensemble_statistics(pts2d, conf, avg_mode=avg_mode, var_mode=var)
    mask: Bool[np.ndarray, "V T P"] = np.isfinite(center).all(axis=-1)
    if init3d is None:
        init3d = triangulate(rig, np.where(mask[..., None], center, np.nan))
    init3d = np.asarray(init3d, dtype=float)
    if init3d.shape != (n_frames, n_points, 3):
        raise ValueError(
            f"init3d must be ({n_frames}, {n_points}, 3); got {init3d.shape}"
        )

    # A keypoint two views never see at once is unobservable and ends as NaN; one
    # seen so in some frames is kept and interpolated (ADR-0062).
    determined = np.isfinite(init3d).all(axis=-1).any(axis=0)  # (P,)
    if not determined.all():
        # RANSAC can drop a keypoint that plain DLT resolves; rescue it before the
        # inflation, which skips cells without a finite 3D.
        fallback = np.asarray(
            triangulate(rig, np.where(mask[..., None], center, np.nan))
        )
        rescued = np.isfinite(fallback).all(axis=-1).any(axis=0) & ~determined
        if rescued.any():
            init3d = init3d.copy()
            init3d[:, rescued] = fallback[:, rescued]
            determined = determined | rescued
        if not determined.all():
            log.info(
                "eks: %d keypoint(s) are never seen by two views at once and stay NaN",
                int((~determined).sum()),
            )

    # The gate inflate_variances applies, so the reported rate is of cells it tested.
    finite3d = np.isfinite(init3d).all(axis=-1)  # (T, P)
    n_testable = int(
        np.count_nonzero(mask & (finite3d & (mask.sum(axis=0) >= 2))[None])
    )
    n_inflated = 0
    if inflate:
        var2d, n_inflated = inflate_variances(
            center,
            var2d,
            mask,
            init3d,
            rig,
            threshold=inflate_threshold,
            factor=inflate_factor,
        )
        log.info(
            "eks: variance inflation down-weighted %d of %d testable observations "
            "(%.2f%%)",
            n_inflated,
            n_testable,
            100.0 * n_inflated / max(n_testable, 1),
        )

    # Zero the unobserved cells so no NaN reaches the kernels; the mask excludes them.
    obs = _per_keypoint(np.where(mask[..., None], center, 0.0))  # (P, T, 2V)
    obs_var = _per_keypoint(var2d)
    obs_mask = _per_keypoint(np.repeat(mask[..., None], 2, axis=-1)).astype(bool)
    m0, prior_cov, q_shape = _initial_state(init3d)
    seen = obs_mask.any(axis=(1, 2)) & determined  # (P,) -- keypoints worth solving

    rig_arrays = (
        jnp.asarray(rig.rvecs),
        jnp.asarray(rig.tvecs),
        jnp.asarray(rig.intrs),
        jnp.asarray(rig.dists),
        jnp.asarray(rig.mirrored),
    )

    if smooth_param is None:
        n_fit = n_frames if not fit_frames else min(int(fit_frames), n_frames)
        log.info(
            "eks: fitting the smoothing parameter for %d keypoint(s) on %d frame(s)",
            int(seen.sum()),
            n_fit,
        )
        fit = jax.jit(
            jax.vmap(
                lambda y, m, v, mu, s0, q: fit_smooth_param(
                    y, m, v, mu, s0, q, rig_arrays, n_iter=fit_iterations
                )
            )
        )
        log_s = np.asarray(
            fit(
                jnp.asarray(obs[:, :n_fit]),
                jnp.asarray(obs_mask[:, :n_fit]),
                jnp.asarray(obs_var[:, :n_fit]),
                jnp.asarray(m0),
                jnp.asarray(prior_cov),
                jnp.asarray(q_shape),
            )
        )
        params = np.where(seen, np.exp(log_s), np.nan)
    else:
        params = np.broadcast_to(
            np.asarray(smooth_param, dtype=float), (n_points,)
        ).astype(float)
        if not np.all(np.isfinite(params) & (params > 0)):
            raise ValueError("smooth_param must be positive and finite")
        # NaN for the unsolved keypoints, as the fitted branch reports.
        params = np.where(seen, params, np.nan)

    log.info("eks: smoothing %d frames x %d keypoints", n_frames, n_points)
    run = jax.jit(
        jax.vmap(
            lambda y, m, v, mu, s0, q, s: smooth_keypoint(
                y, m, v, mu, s0, q, s, rig_arrays
            )
        )
    )
    means, covs = run(
        jnp.asarray(obs),
        jnp.asarray(obs_mask),
        jnp.asarray(obs_var),
        jnp.asarray(m0),
        jnp.asarray(prior_cov),
        jnp.asarray(q_shape),
        jnp.asarray(np.where(seen, params, 1.0)),
    )
    # np.asarray of a JAX array is read-only; copy before masking.
    pts3d = np.array(np.swapaxes(np.asarray(means), 0, 1))  # (T, P, 3)
    posterior_var = np.array(
        np.swapaxes(np.diagonal(np.asarray(covs), axis1=-2, axis2=-1), 0, 1)
    )  # (T, P, 3)
    unseen = ~seen
    pts3d[:, unseen] = np.nan
    posterior_var[:, unseen] = np.nan

    reprojected = np.asarray(rig.project(pts3d))  # (V, T, P, 2)
    if not fill_unobserved:
        reprojected = np.where(mask[..., None], reprojected, np.nan)

    return EksResult(
        pts3d=pts3d,
        pts2d=reprojected,
        posterior_var=posterior_var,
        smooth_param=params,
        n_inflated=n_inflated,
        n_testable=n_testable,
    )

Anatomical corrections over a 3D track: hold a point still, mirror a pair.

A tethered fly's thorax-coxa joints and neck do not move in its frame, so their jitter is noise and one position per track the better estimate; and the fly is bilaterally symmetric, a left point and its partner reflections through one sagittal plane. Pure functions over arrays and point indices (ADR-0039); the [correction] stage alone knows which points a config named.

freeze

freeze(
    pts3d: Float[ndarray, "T P 3"],
    points: Sequence[int],
    *,
    method: str = "median",
    trim: float = 0.1,
) -> Float[np.ndarray, "T P 3"]

Replace each of points with one position, held over every frame.

Occluded frames are filled in, and a point no frame observed stays NaN; method is one of STATIC_METHODS, trim the tail fraction of "trimmed_mean". Returns a copy.

Raises:

Type	Description
ValueError	If method is unknown, or trim is outside [0, 0.5).

Source code in src/deeperfly/pose3d/correction.py

def freeze(
    pts3d: Float[np.ndarray, "T P 3"],
    points: Sequence[int],
    *,
    method: str = "median",
    trim: float = 0.1,
) -> Float[np.ndarray, "T P 3"]:
    """Replace each of `points` with one position, held over every frame.

    Occluded frames are filled in, and a point no frame observed stays NaN; `method`
    is one of `STATIC_METHODS`, `trim` the tail fraction of `"trimmed_mean"`. Returns a
    copy.

    Raises:
        ValueError: If `method` is unknown, or `trim` is outside `[0, 0.5)`.
    """
    out = np.array(pts3d, dtype=float)
    cols = list(points)
    if not cols:
        return out
    block = out[:, cols, :]  # (T, R, 3)
    if method == "geometric_median":
        # Per point, over time: the coordinate axis is a vector, not a cell.
        center = np.stack([_geometric_median(block[:, i]) for i in range(len(cols))])
    else:
        center = _temporal_center(block, method=method, trim=trim)
    out[:, cols, :] = center
    return out

sagittal_plane

sagittal_plane(
    pts3d: Float[ndarray, "*t P 3"],
    pairs: Sequence[tuple[int, int]],
) -> tuple[Float[np.ndarray, " 3"] | None, float | None]

Fit the animal's mirror plane from left/right pairs: (normal, offset).

A pair's midpoint lies on the plane and its difference is normal to it; the normal is the leading singular vector of the differences (their signs are arbitrary, so a mean could cancel) and the offset puts the plane through the midpoints' centroid. (None, None) below two usable pairs, the under-determined case.

Source code in src/deeperfly/pose3d/correction.py

def sagittal_plane(
    pts3d: Float[np.ndarray, "*t P 3"], pairs: Sequence[tuple[int, int]]
) -> tuple[Float[np.ndarray, " 3"] | None, float | None]:
    """Fit the animal's mirror plane from left/right pairs: `(normal, offset)`.

    A pair's midpoint lies on the plane and its difference is normal to it; the normal
    is the leading singular vector of the differences (their signs are arbitrary, so a
    mean could cancel) and the offset puts the plane through the midpoints' centroid.
    `(None, None)` below two usable pairs, the under-determined case.
    """
    if not pairs:
        return None, None
    left = pts3d[..., [a for a, _ in pairs], :]
    right = pts3d[..., [b for _, b in pairs], :]
    ok = np.isfinite(left).all(-1) & np.isfinite(right).all(-1)
    if ok.sum() < 2:
        return None, None
    diff = (left - right)[ok]  # (N, 3), normal directions with arbitrary sign
    mid = ((left + right) / 2.0)[ok]  # (N, 3), points on the plane
    _, _, vt = np.linalg.svd(diff, full_matrices=False)
    normal = vt[0] / np.linalg.norm(vt[0])
    return normal, float(normal @ mid.mean(axis=0))

symmetrize

symmetrize(
    pts3d: Float[ndarray, "T P 3"],
    pairs: Sequence[tuple[int, int]],
    midline: Sequence[int] = (),
    *,
    per_frame: bool = False,
    strength: float = 1.0,
) -> tuple[
    Float[np.ndarray, "T P 3"],
    Float[np.ndarray, " 3"] | None,
    float | None,
]

Make pairs mirror-symmetric and put midline on the plane.

Each side moves half-way toward the other's mirror (times strength), which privileges neither; a midline point moves onto the plane along its normal. The plane is fitted from the pairs alone, since the midline is about to be moved onto it. With per_frame=False one plane serves the track, a fixed affine map, so a point freeze made static stays static. Returns (pts3d, normal, offset), normal None (and the track unchanged) when no plane fits.

Source code in src/deeperfly/pose3d/correction.py

def symmetrize(
    pts3d: Float[np.ndarray, "T P 3"],
    pairs: Sequence[tuple[int, int]],
    midline: Sequence[int] = (),
    *,
    per_frame: bool = False,
    strength: float = 1.0,
) -> tuple[Float[np.ndarray, "T P 3"], Float[np.ndarray, " 3"] | None, float | None]:
    """Make `pairs` mirror-symmetric and put `midline` on the plane.

    Each side moves half-way toward the other's mirror (times `strength`), which
    privileges neither; a midline point moves onto the plane along its normal. The
    plane is fitted from the pairs alone, since the midline is about to be moved onto
    it. With `per_frame=False` one plane serves the track, a fixed affine map, so a
    point [freeze][deeperfly.pose3d.correction.freeze] made static stays static. Returns
    `(pts3d, normal, offset)`, `normal` `None` (and the track unchanged) when no plane
    fits.
    """
    out = np.array(pts3d, dtype=float)
    if strength <= 0.0 or (not pairs and not midline):
        return out, None, None
    if per_frame:
        normal = offset = None
        for t in range(out.shape[0]):
            normal, offset = sagittal_plane(out[t : t + 1], pairs)
            _reflect_onto(out[t], pairs, midline, normal, offset, strength)
        return out, normal, offset
    normal, offset = sagittal_plane(out, pairs)
    _reflect_onto(out, pairs, midline, normal, offset, strength)
    return out, normal, offset