stereotaxic.py

"""Utilities for calculating tilted/rotated stereotaxic target coordinates.

This module ports and cleans up the original MATLAB/Python ``target3d`` routine
by Tom Davidson (Stanford University, 2010-2016).

Coordinate convention
---------------------
Atlas coordinates are ``[ML, AP, DV]`` with ``+ML`` to the right, ``+AP`` toward
the nose, and ``+DV`` dorsal (up). Ventral targets therefore have negative DV.

Arm configurations
------------------
- ``LR``: the manipulator's ML axis tilts together with the DV axis (typical
  left-right tilting arm). The arm's lean direction must include a left-right
  component, so ``rotation_degrees in {0, 180}`` is rejected as singular.
- ``FB``: the manipulator's ML axis stays horizontal under tilt (typical
  front-back tilting arm). The arm's lean direction must include a
  nose-tail component, so ``rotation_degrees in {90, 270}`` is rejected as
  singular.

In the ``FB`` configuration the basis is *not* orthonormal: the manipulator's
AP stage moves along the atlas AP direction regardless of tilt or rotation
(this matches typical stereotax mechanics, where the AP stage sits below the
tilt/rotate mechanism). The returned ``[ML, AP, DV]`` coordinates are still
the correct dial settings to move the manipulator to the target.

Routines
--------
target3d
    Top-level entry point. Computes manipulator dial settings (and
    optionally the skull-entry point) for a tilted/rotated approach.
stereotaxic_basis
    Builds the 3x3 basis matrix whose columns are the manipulator's
    ML, AP, and DV axes in atlas coordinates.
plane_intersection
    Intersects an arbitrary approach path with a horizontal DV plane.
spherical_to_cartesian
    Low-level spherical-to-Cartesian conversion using the stereotaxic
    angle convention.
format_coordinates
    Human-readable formatter for ``[ML, AP, DV]`` triples.

Example
-------
>>> import numpy as np
>>> from stereotaxic import target3d
>>> result = target3d([0.25, 0.7, -6.6], "LR", 20, 90)
>>> np.round(result.target_coordinates, 2)
array([ 2.49,  0.7 , -6.12])
"""

from __future__ import annotations

import math
import warnings
from typing import Literal, NamedTuple, TypeAlias, cast

import numpy as np
from numpy.typing import ArrayLike, NDArray

FloatArray: TypeAlias = NDArray[np.float64]
ArmConfig: TypeAlias = Literal["LR", "FB", "lr", "fb"]
NormalizedArmConfig: TypeAlias = Literal["lr", "fb"]

# Atlas AP unit vector, used as the manipulator's AP axis in every basis.
# Defined once at module load to avoid reallocation on every call. The
# read-only flag prevents accidental mutation by callers.
_ATLAS_AP_AXIS: FloatArray = np.array([0.0, 1.0, 0.0], dtype=np.float64)
_ATLAS_AP_AXIS.setflags(write=False)

__all__ = [
    "ArmConfig",
    "Target3DResult",
    "format_coordinates",
    "plane_intersection",
    "spherical_to_cartesian",
    "stereotaxic_basis",
    "target3d",
]


class Target3DResult(NamedTuple):
    """Output from :func:`target3d`.

    Attributes
    ----------
    target_coordinates : ndarray of shape (3,)
        Target coordinates ``[ML, AP, DV]`` in the tilted/rotated stereotaxic
        arm coordinate system.
    plane_intersection : ndarray of shape (3,) or None
        Coordinates ``[ML, AP, DV]`` of the approach path's intersection with
        the requested horizontal DV plane, in the original/unrotated atlas
        coordinate system. ``None`` when no ``dv_intersect`` value was supplied.
    """

    target_coordinates: FloatArray
    plane_intersection: FloatArray | None


def _as_3vector(value: ArrayLike, name: str = "vector") -> FloatArray:
    """Return ``value`` as a finite ``float64`` vector of shape ``(3,)``.

    Parameters
    ----------
    value : array_like of shape (3,)
        Coordinates ``[ML, AP, DV]``.
    name : str, default="vector"
        Argument name used in error messages.

    Returns
    -------
    vector : ndarray of shape (3,)
        Validated 3-vector with ``float64`` dtype.

    Raises
    ------
    ValueError
        If ``value`` is not one-dimensional with length 3, or contains NaN or
        infinite values.
    """
    vector = np.asarray(value, dtype=np.float64)

    if vector.shape != (3,):
        raise ValueError(
            f"{name} must be a one-dimensional array-like with shape (3,) "
            f"representing [ML, AP, DV]."
        )
    if not np.all(np.isfinite(vector)):
        raise ValueError(f"{name} must contain only finite numeric values.")

    return vector


def _normalize_arm_config(arm_config: str) -> NormalizedArmConfig:
    """Normalize and validate the stereotaxic arm configuration.

    Parameters
    ----------
    arm_config : {'LR', 'FB', 'lr', 'fb'}
        ``'LR'`` for left-right tilt or ``'FB'`` for front-back tilt.

    Returns
    -------
    normalized_config : {'lr', 'fb'}
        Lowercase arm configuration.

    Raises
    ------
    ValueError
        If ``arm_config`` is not ``'LR'`` or ``'FB'``.
    """
    normalized = arm_config.lower()
    if normalized not in {"lr", "fb"}:
        raise ValueError("arm_config must be either 'LR' or 'FB'.")
    return cast(NormalizedArmConfig, normalized)


def _validate_tilt(tilt_degrees: float) -> float:
    """Validate a tilt angle.

    Parameters
    ----------
    tilt_degrees : float
        Tilt angle in degrees. Valid range is ``0 <= tilt_degrees <= 90``.

    Returns
    -------
    tilt_degrees : float
        Validated tilt angle as a ``float``.

    Raises
    ------
    ValueError
        If the tilt is outside the valid range or is not finite.
    """
    tilt = float(tilt_degrees)
    if not math.isfinite(tilt):
        raise ValueError("tilt_degrees must be finite.")
    if not 0.0 <= tilt <= 90.0:
        raise ValueError("tilt_degrees must be between 0 and 90 degrees.")
    return tilt


def _normalize_rotation(rotation_degrees: float) -> float:
    """Normalize rotation to the interval ``[0, 360)`` degrees.

    Parameters
    ----------
    rotation_degrees : float
        Clockwise rotation angle in degrees, where ``0`` points toward nose,
        ``90`` points right/lateral, and ``180`` points posterior/tail.

    Returns
    -------
    normalized_rotation : float
        Rotation angle in the interval ``[0, 360)``.

    Raises
    ------
    ValueError
        If ``rotation_degrees`` is not finite.
    """
    rotation = float(rotation_degrees)
    if not math.isfinite(rotation):
        raise ValueError("rotation_degrees must be finite.")
    return rotation % 360.0


def _effective_tilt_and_rotation(
    tilt_degrees: float,
    rotation_degrees: float,
) -> tuple[float, float]:
    """Return the signed tilt and adjusted rotation used internally.

    Parameters
    ----------
    tilt_degrees : float
        Validated tilt angle in degrees.
    rotation_degrees : float
        Normalized clockwise rotation angle in degrees, in ``[0, 360)``.

    Returns
    -------
    effective_tilt_degrees : float
        Signed tilt angle in degrees. Posterior approaches use negative tilt.
    effective_rotation_degrees : float
        Adjusted rotation angle in degrees.
    """
    if 90.0 < rotation_degrees < 270.0:
        return -tilt_degrees, (rotation_degrees - 180.0) % 360.0
    return tilt_degrees, rotation_degrees


def spherical_to_cartesian(
    theta_radians: float,
    phi_radians: float,
    radius: float = 1.0,
) -> FloatArray:
    """Convert spherical coordinates to a Cartesian ``[ML, AP, DV]`` vector.

    Parameters
    ----------
    theta_radians : float
        Azimuth in the ML-AP plane, in radians. ``0`` points toward ``+ML``
        (right) and ``pi/2`` points toward ``+AP`` (nose).
    phi_radians : float
        Elevation above the horizontal ML-AP plane, in radians. ``0`` is
        horizontal and ``pi/2`` is vertical (``+DV``, up).
    radius : float, default=1.0
        Vector length.

    Returns
    -------
    vector : ndarray of shape (3,)
        Cartesian coordinates ``[ML, AP, DV]``.

    Raises
    ------
    ValueError
        If any input is not finite.

    See Also
    --------
    stereotaxic_basis : Builds the rotated basis using this function.

    Notes
    -----
    The azimuth convention here is not the math/physics default.
    ``theta=0`` is aligned with the ``+ML`` axis (not ``+AP``) so that
    ``theta = pi/2 - rot`` maps cleanly to the stereotaxic rotation
    parameter, where ``rot=0`` points toward the nose.

    The components are::

        ML = r * cos(phi) * cos(theta)
        AP = r * cos(phi) * sin(theta)
        DV = r * sin(phi)

    Examples
    --------
    >>> spherical_to_cartesian(0.0, 0.0)
    array([1., 0., 0.])
    >>> np.allclose(spherical_to_cartesian(0.0, np.pi / 2), [0, 0, 1])
    True
    >>> np.allclose(spherical_to_cartesian(np.pi / 2, 0.0), [0, 1, 0])
    True
    """
    theta = float(theta_radians)
    phi = float(phi_radians)
    r = float(radius)

    if not (math.isfinite(theta) and math.isfinite(phi) and math.isfinite(r)):
        raise ValueError("theta_radians, phi_radians, and radius must be finite.")

    cos_phi = math.cos(phi)
    sin_phi = math.sin(phi)
    cos_theta = math.cos(theta)
    sin_theta = math.sin(theta)
    return np.array(
        [r * cos_phi * cos_theta, r * cos_phi * sin_theta, r * sin_phi],
        dtype=np.float64,
    )


def stereotaxic_basis(
    arm_config: ArmConfig | str,
    tilt_degrees: float,
    rotation_degrees: float,
) -> FloatArray:
    """Construct the tilted/rotated stereotaxic basis matrix.

    Parameters
    ----------
    arm_config : {'LR', 'FB', 'lr', 'fb'}
        ``'LR'`` for a left-right tilt mechanism or ``'FB'`` for a front-back
        tilt mechanism. Case-insensitive.
    tilt_degrees : float
        Tilt angle from vertical, in degrees. Must satisfy
        ``0 <= tilt_degrees <= 90``.
    rotation_degrees : float
        Clockwise rotation in degrees, normalized modulo 360. ``0`` points
        toward the nose, ``90`` toward the animal's right, ``180`` toward
        the tail.

    Returns
    -------
    basis : ndarray of shape (3, 3)
        Basis matrix whose columns are the manipulator's ``ML``, ``AP``,
        and ``DV`` axes, respectively, expressed in the atlas coordinate
        system. A target ``t`` in atlas coordinates can be recovered as
        ``t = basis @ c``, where ``c`` is the target in manipulator
        coordinates.

    Raises
    ------
    ValueError
        If inputs are invalid, or if the requested ``(arm_config,
        rotation_degrees)`` combination is singular. ``LR`` cannot lean
        toward nose/tail (``rot in {0, 180}``); ``FB`` cannot lean toward
        left/right (``rot in {90, 270}``).

    See Also
    --------
    target3d : Solves ``basis @ c == target`` for the manipulator dial
        settings ``c``.
    spherical_to_cartesian : Computes each axis from spherical angles.

    Notes
    -----
    The ``FB`` basis is intentionally *not* orthonormal: the manipulator's
    AP axis is fixed at the atlas ``[0, 1, 0]`` regardless of tilt or
    rotation, because the AP stage on typical front-back arms sits below
    the tilt/rotate mechanism. The basis is still invertible for non-
    singular ``(arm, rotation)`` pairs, and ``solve(basis, target)``
    returns the correct manipulator dial settings.

    The ``LR`` basis is orthonormal at ``rotation_degrees in {90, 270}``
    (the lateral approaches) and non-orthonormal elsewhere, for the same
    reason.

    Examples
    --------
    >>> B = stereotaxic_basis("LR", 20, 90)
    >>> np.round(B, 4)
    array([[ 0.9397,  0.    ,  0.342 ],
           [ 0.    ,  1.    ,  0.    ],
           [-0.342 ,  0.    ,  0.9397]])

    >>> stereotaxic_basis("LR", 30, 0)
    Traceback (most recent call last):
        ...
    ValueError: LR arm cannot lean toward nose/tail ...
    """
    normalized_config = _normalize_arm_config(str(arm_config))
    tilt = _validate_tilt(tilt_degrees)
    rotation = _normalize_rotation(rotation_degrees)

    # Reject configurations where the arm mechanism cannot lean in the
    # requested direction. These produce a rank-2 basis at every tilt;
    # the original MATLAB/Python code silently returned ~1e16-magnitude
    # garbage from explicit matrix inversion.
    if normalized_config == "lr" and (
        np.isclose(rotation, 0.0) or np.isclose(rotation, 180.0)
    ):
        raise ValueError(
            f"LR arm cannot lean toward nose/tail "
            f"(rotation_degrees={rotation_degrees}, normalized to {rotation}). "
            f"Use the 'FB' arm config instead."
        )
    if normalized_config == "fb" and (
        np.isclose(rotation, 90.0) or np.isclose(rotation, 270.0)
    ):
        raise ValueError(
            f"FB arm cannot lean toward left/right "
            f"(rotation_degrees={rotation_degrees}, normalized to {rotation}). "
            f"Use the 'LR' arm config instead."
        )

    effective_tilt, effective_rotation = _effective_tilt_and_rotation(tilt, rotation)

    theta = np.deg2rad(-effective_rotation + 90.0)
    phi = np.deg2rad(-effective_tilt + 90.0)

    dv_axis = spherical_to_cartesian(theta, phi)
    if normalized_config == "lr":
        ml_axis = spherical_to_cartesian(theta, phi - np.pi / 2.0)
    else:
        ml_axis = spherical_to_cartesian(theta - np.pi / 2.0, 0.0)

    return np.column_stack((ml_axis, _ATLAS_AP_AXIS, dv_axis))


def plane_intersection(
    target: ArrayLike,
    direction: ArrayLike,
    dv_intersect: float,
) -> FloatArray:
    """Find where an approach path intersects a horizontal DV plane.

    Parameters
    ----------
    target : array_like of shape (3,)
        Atlas coordinates ``[ML, AP, DV]`` of a point on the approach path
        (typically the actual target deep in the brain).
    direction : array_like of shape (3,)
        Direction vector ``[ML, AP, DV]`` of the approach path in atlas
        space. Does not need to be unit-normalized, but its DV component
        must be non-trivial relative to its overall magnitude.
    dv_intersect : float
        DV value of the horizontal plane to intersect (for example, skull
        surface or a plane slightly below bregma).

    Returns
    -------
    intersection : ndarray of shape (3,)
        Atlas coordinates ``[ML, AP, DV]`` of the intersection point. By
        construction, ``intersection[2] == dv_intersect``.

    Raises
    ------
    ValueError
        If ``dv_intersect`` is not finite, the direction vector has a
        negligible DV component (path is essentially horizontal), or
        ``target`` / ``direction`` fail shape/finiteness validation.

    See Also
    --------
    target3d : Calls this with the approach path's DV axis to compute
        skull-entry points.

    Examples
    --------
    >>> # Vertical approach through (1, 1, -5), intersect DV=0 plane.
    >>> plane_intersection([1.0, 1.0, -5.0], [0.0, 0.0, 1.0], 0.0)
    array([1., 1., 0.])

    >>> # 45-degree approach intersecting DV=0 plane.
    >>> plane_intersection([0.0, 0.0, -10.0], [1.0, 0.0, 1.0], 0.0)
    array([10.,  0.,  0.])
    """
    target_vector = _as_3vector(target, "target")
    direction_vector = _as_3vector(direction, "direction")
    dv_plane = float(dv_intersect)

    if not math.isfinite(dv_plane):
        raise ValueError("dv_intersect must be finite.")

    # Use a relative tolerance against the direction-vector norm so the
    # check is meaningful regardless of how the caller scales `direction`.
    direction_norm = float(np.linalg.norm(direction_vector))
    if abs(direction_vector[2]) < 1e-9 * max(direction_norm, 1.0):
        raise ValueError(
            "Cannot compute plane intersection because the approach path has "
            "a negligible DV component. This usually occurs at 90-degree "
            "horizontal tilt."
        )

    scale = (dv_plane - target_vector[2]) / direction_vector[2]
    return target_vector + scale * direction_vector


def format_coordinates(label: str, coordinates: ArrayLike) -> str:
    """Format a 3-vector of stereotaxic coordinates for display.

    Parameters
    ----------
    label : str
        Label printed before the coordinates (e.g. ``'Target'``).
    coordinates : array_like of shape (3,)
        Coordinates ``[ML, AP, DV]``, in millimeters.

    Returns
    -------
    line : str
        A single-line, signed, two-decimal-place representation matching
        the format printed by stereotaxic digital readouts.

    Raises
    ------
    ValueError
        If ``coordinates`` does not have shape ``(3,)`` or contains
        non-finite values.

    Examples
    --------
    >>> format_coordinates("Target", [1.0, -2.5, 3.14])
    'Target: ML:+1.00 / AP:-2.50 / DV:+3.14'

    >>> format_coordinates("Entry", np.array([-0.05, 1.234, -0.3]))
    'Entry: ML:-0.05 / AP:+1.23 / DV:-0.30'
    """
    vector = _as_3vector(coordinates, "coordinates")
    return f"{label}: ML:{vector[0]:+0.2f} / AP:{vector[1]:+0.2f} / DV:{vector[2]:+0.2f}"


def target3d(
    target: ArrayLike,
    arm_config: ArmConfig | str,
    tilt_degrees: float,
    rotation_degrees: float,
    dv_intersect: float | None = None,
    *,
    verbose: bool = False,
) -> Target3DResult:
    """Calculate dial settings for a tilted and rotated stereotaxic approach.

    Parameters
    ----------
    target : array_like of shape (3,)
        Target atlas coordinates ``[ML, AP, DV]``, in millimeters.
        ``+ML`` is right, ``+AP`` is anterior (toward the nose), and
        ``+DV`` is dorsal (up).
    arm_config : {'LR', 'FB', 'lr', 'fb'}
        Stereotaxic arm mechanism. ``'LR'`` is a left-right tilting arm,
        ``'FB'`` is a front-back tilting arm. Case-insensitive.
    tilt_degrees : float
        Tilt angle from vertical, in degrees. Must satisfy
        ``0 <= tilt_degrees <= 90``.
    rotation_degrees : float
        Clockwise rotation in degrees, normalized modulo 360. ``0`` points
        toward the nose, ``90`` toward the animal's right, ``180`` toward
        the tail.
    dv_intersect : float or None, default=None
        DV value of a horizontal plane to intersect with the approach
        path, in millimeters (typically the skull surface or a plane
        slightly below bregma). If ``None``, no intersection is computed.
    verbose : bool, default=False
        If ``True``, print formatted coordinates to ``stdout``.

    Returns
    -------
    result : Target3DResult
        Named tuple with fields:

        ``target_coordinates`` : ndarray of shape (3,)
            Manipulator dial settings ``[ML, AP, DV]`` to reach the target
            in the tilted/rotated arm coordinate system.
        ``plane_intersection`` : ndarray of shape (3,) or None
            Atlas coordinates of the approach path's intersection with the
            ``dv_intersect`` plane, or ``None`` if ``dv_intersect`` was not
            supplied.

    Raises
    ------
    ValueError
        If inputs are invalid, or if the requested ``(arm_config,
        rotation_degrees)`` pair is singular. ``LR`` cannot lean toward
        nose/tail (``rot in {0, 180}``); ``FB`` cannot lean toward
        left/right (``rot in {90, 270}``).
    RuntimeWarning
        Issued (not raised) if ``target[2] > 0``, since ventral
        coordinates are conventionally negative.

    See Also
    --------
    stereotaxic_basis : Builds the basis matrix solved here.
    plane_intersection : Computes the skull-entry point when
        ``dv_intersect`` is supplied.
    format_coordinates : Formats the output for printed display.

    Notes
    -----
    The returned ``target_coordinates`` are the dial settings to move the
    manipulator's tip from origin to ``target`` along the manipulator's
    own axes. For ``FB`` configurations, the manipulator's AP axis stays
    aligned with the atlas AP axis regardless of tilt or rotation (see
    :func:`stereotaxic_basis`).

    Examples
    --------
    Lateral approach with the LR arm tilted 20 degrees to the right:

    >>> result = target3d([0.25, 0.7, -6.6], "LR", 20, 90)
    >>> np.round(result.target_coordinates, 2)
    array([ 2.49,  0.7 , -6.12])

    FB arm tilted 30 degrees forward and rotated 20 degrees counter-
    clockwise, with a skull-surface plane at DV = -0.3:

    >>> result = target3d([-2, -2, -3.1], "FB", 30, -20, dv_intersect=-0.3)
    >>> np.round(result.target_coordinates, 2)
    array([-2.78,  0.63, -3.58])
    >>> np.round(result.plane_intersection, 2)
    array([-2.55, -0.48, -0.3 ])
    """
    target_vector = _as_3vector(target, "target")
    if target_vector[2] > 0:
        # Warning rather than error because the original MATLAB/Python behavior
        # allowed positive-DV targets, but they are unusual in this convention.
        warnings.warn(
            "target DV is positive. Ventral coordinates are usually represented "
            "with negative values.",
            RuntimeWarning,
            stacklevel=2,
        )

    # `stereotaxic_basis` rejects the known-singular (arm, rotation) pairs
    # explicitly. The try/except below is a defensive safety net for any
    # other singular geometry that slips through, and exposes the underlying
    # LinAlgError for debugging via exception chaining.
    basis = stereotaxic_basis(arm_config, tilt_degrees, rotation_degrees)
    try:
        target_coordinates = np.linalg.solve(basis, target_vector)
    except np.linalg.LinAlgError as exc:
        raise ValueError(
            "Computed coordinate basis is singular. This should not occur "
            "for any valid input; please report it as a bug."
        ) from exc

    intersection = None
    if dv_intersect is not None:
        dv_axis = basis[:, 2]
        intersection = plane_intersection(target_vector, dv_axis, dv_intersect)

    result = Target3DResult(target_coordinates, intersection)

    if verbose:
        print("\n--- Calculated Stereotaxic Coordinates ---")
        print(format_coordinates("Target (new coords)", result.target_coordinates))
        if result.plane_intersection is not None:
            print(format_coordinates("Plane intersection", result.plane_intersection))

    return result


if __name__ == "__main__":
    example = target3d([1.3, 1.2, -5.3], "LR", 10, 90, verbose=True)
    print(example)

test_stereotaxic.py

"""Tests for the ``stereotaxic`` module.

Run with ``pytest tests/test_stereotaxic.py``.
"""

from __future__ import annotations

import numpy as np
import pytest

from stereotaxic import (
    Target3DResult,
    format_coordinates,
    plane_intersection,
    spherical_to_cartesian,
    stereotaxic_basis,
    target3d,
)


# ---------------------------------------------------------------------------
# Canonical examples from the original notebook
# ---------------------------------------------------------------------------
class TestDocstringExamples:
    """The canonical ``(target, arm, tilt, rot)`` examples documented in
    the original MATLAB/notebook implementation. Tolerances match the
    2-decimal display the digital readout uses."""

    def test_lr_arm_tilted_right(self):
        # "Stereotax arm tilted 20 degrees to the right."
        result = target3d([0.25, 0.7, -6.6], "LR", 20, 90)
        np.testing.assert_allclose(
            result.target_coordinates, [2.49, 0.70, -6.12], atol=5e-3
        )
        assert result.plane_intersection is None

    def test_fb_arm_forward_and_left(self):
        # "Arm leaning forward and towards the animal's left."
        result = target3d([-2, -2, -3.1], "FB", 30, -20)
        np.testing.assert_allclose(
            result.target_coordinates, [-2.78, 0.63, -3.58], atol=5e-3
        )

    def test_fb_arm_forward_and_left_with_skull_plane(self):
        # Same approach as above plus a skull-surface plane at DV=-0.3.
        result = target3d([-2, -2, -3.1], "FB", 30, -20, dv_intersect=-0.3)
        np.testing.assert_allclose(
            result.target_coordinates, [-2.78, 0.63, -3.58], atol=5e-3
        )
        np.testing.assert_allclose(
            result.plane_intersection, [-2.55, -0.48, -0.30], atol=5e-3
        )

    def test_fb_arm_tilted_backward_and_right(self):
        # rot=150 triggers the "backward approach" sign-flip code path.
        result = target3d([1.2, -1.3, -4.4], "FB", 25, 150)
        np.testing.assert_allclose(
            result.target_coordinates, [2.57, -4.36, -4.85], atol=5e-3
        )


# ---------------------------------------------------------------------------
# Singular configurations
# ---------------------------------------------------------------------------
class TestSingularConfigurations:
    """LR cannot lean nose/tail; FB cannot lean left/right.

    The original implementation silently returned ~1e16-magnitude garbage for
    these via explicit matrix inversion. The refactor must reject them.
    """

    @pytest.mark.parametrize("rotation_degrees", [0.0, 180.0, 360.0, -180.0])
    def test_lr_rejected_when_lean_is_nose_or_tail(self, rotation_degrees):
        with pytest.raises(ValueError, match="LR arm cannot lean"):
            target3d([1.0, 2.0, -3.0], "LR", 30, rotation_degrees)

    @pytest.mark.parametrize("rotation_degrees", [90.0, 270.0, -90.0, 450.0])
    def test_fb_rejected_when_lean_is_left_or_right(self, rotation_degrees):
        with pytest.raises(ValueError, match="FB arm cannot lean"):
            target3d([1.0, 2.0, -3.0], "FB", 30, rotation_degrees)

    def test_singular_check_fires_at_zero_tilt(self):
        # Singularity is a property of (arm, rotation), not tilt.
        with pytest.raises(ValueError, match="LR arm cannot lean"):
            target3d([1.0, 2.0, -3.0], "LR", 0, 0)
        with pytest.raises(ValueError, match="FB arm cannot lean"):
            target3d([1.0, 2.0, -3.0], "FB", 0, 90)

    def test_singular_error_message_suggests_alternate_arm(self):
        with pytest.raises(ValueError, match="Use the 'FB' arm config"):
            stereotaxic_basis("LR", 30, 0)
        with pytest.raises(ValueError, match="Use the 'LR' arm config"):
            stereotaxic_basis("FB", 30, 90)


# ---------------------------------------------------------------------------
# tilt=0: manipulator-frame rotation
# ---------------------------------------------------------------------------
class TestZeroTiltRotation:
    """At ``tilt=0`` the approach is vertical, but the manipulator's
    ML/AP axes are still rotated by ``rotation_degrees`` about DV.
    Non-singular rotations must therefore produce a correspondingly
    rotated coordinate."""

    def test_lr_natural_lateral_is_identity(self):
        # rot=90 is the LR arm's natural lateral orientation.
        result = target3d([1.0, 2.0, -3.0], "LR", 0, 90)
        np.testing.assert_allclose(result.target_coordinates, [1.0, 2.0, -3.0])

    def test_lr_45_degree_rotation_rotates_ml(self):
        # 45-degree base rotation: 1.414 along the rotated ML axis,
        # 1.0 along the (unrotated) AP stage.
        result = target3d([1.0, 2.0, -3.0], "LR", 0, 45)
        np.testing.assert_allclose(
            result.target_coordinates, [np.sqrt(2), 1.0, -3.0], atol=1e-12
        )

    def test_lr_lateral_left_negates_ml(self):
        result = target3d([1.0, 2.0, -3.0], "LR", 0, 270)
        np.testing.assert_allclose(result.target_coordinates, [-1.0, 2.0, -3.0])

    def test_fb_natural_forward_is_identity(self):
        result = target3d([1.0, 2.0, -3.0], "FB", 0, 0)
        np.testing.assert_allclose(result.target_coordinates, [1.0, 2.0, -3.0])


# ---------------------------------------------------------------------------
# Input validation
# ---------------------------------------------------------------------------
class TestInputValidation:
    """Inputs that the original silently accepted now raise informatively."""

    @pytest.mark.parametrize(
        "bad_target",
        [[1.0, 2.0], [1.0, 2.0, 3.0, 4.0], [[1.0, 2.0, 3.0]]],
    )
    def test_target_wrong_shape_rejected(self, bad_target):
        with pytest.raises(ValueError, match="shape"):
            target3d(bad_target, "LR", 20, 90)

    def test_target_nan_rejected(self):
        with pytest.raises(ValueError, match="finite"):
            target3d([np.nan, 2.0, -3.0], "LR", 20, 90)

    def test_target_inf_rejected(self):
        with pytest.raises(ValueError, match="finite"):
            target3d([np.inf, 2.0, -3.0], "LR", 20, 90)

    @pytest.mark.parametrize("tilt_degrees", [-1.0, 91.0, -0.001, 90.001])
    def test_tilt_out_of_range_rejected(self, tilt_degrees):
        with pytest.raises(ValueError, match="between 0 and 90"):
            target3d([1.0, 2.0, -3.0], "LR", tilt_degrees, 90)

    @pytest.mark.parametrize("bad_tilt", [np.nan, np.inf, -np.inf])
    def test_tilt_non_finite_rejected(self, bad_tilt):
        with pytest.raises(ValueError, match="tilt_degrees must be finite"):
            target3d([1.0, 2.0, -3.0], "LR", bad_tilt, 90)

    @pytest.mark.parametrize("bad_rotation", [np.nan, np.inf, -np.inf])
    def test_rotation_non_finite_rejected(self, bad_rotation):
        with pytest.raises(ValueError, match="rotation_degrees must be finite"):
            target3d([1.0, 2.0, -3.0], "LR", 20, bad_rotation)

    @pytest.mark.parametrize("bad_arm", ["XY", "lL", "", "FBL", "LRR"])
    def test_invalid_arm_config_rejected(self, bad_arm):
        with pytest.raises(ValueError, match="arm_config"):
            target3d([1.0, 2.0, -3.0], bad_arm, 20, 90)

    def test_positive_dv_warns_but_still_computes(self):
        with pytest.warns(RuntimeWarning, match="DV is positive"):
            result = target3d([1.0, 2.0, 1.0], "LR", 20, 90)
        # The warning is advisory; the function still returns a valid result.
        assert result.target_coordinates.shape == (3,)
        assert np.all(np.isfinite(result.target_coordinates))

    def test_verbose_is_keyword_only(self):
        # `verbose` must be keyword-only so it can't be confused with
        # `dv_intersect` if a caller passes too many positional args.
        with pytest.raises(TypeError):
            target3d([1.0, 2.0, -3.0], "LR", 20, 90, None, True)


# ---------------------------------------------------------------------------
# arm_config case insensitivity
# ---------------------------------------------------------------------------
class TestArmConfigCasing:
    """The ``arm_config`` argument must accept any letter casing."""

    @pytest.mark.parametrize("arm", ["LR", "lr", "Lr", "lR"])
    def test_lr_variants_produce_identical_output(self, arm):
        reference = target3d([0.25, 0.7, -6.6], "LR", 20, 90)
        result = target3d([0.25, 0.7, -6.6], arm, 20, 90)
        np.testing.assert_array_equal(
            reference.target_coordinates, result.target_coordinates
        )

    @pytest.mark.parametrize("arm", ["FB", "fb", "Fb", "fB"])
    def test_fb_variants_produce_identical_output(self, arm):
        reference = target3d([-2, -2, -3.1], "FB", 30, -20)
        result = target3d([-2, -2, -3.1], arm, 30, -20)
        np.testing.assert_array_equal(
            reference.target_coordinates, result.target_coordinates
        )


# ---------------------------------------------------------------------------
# plane_intersection
# ---------------------------------------------------------------------------
class TestPlaneIntersection:
    def test_intersection_lies_on_requested_plane(self):
        result = target3d([-2, -2, -3.1], "FB", 30, -20, dv_intersect=-0.3)
        assert result.plane_intersection is not None
        assert result.plane_intersection[2] == pytest.approx(-0.3)

    def test_horizontal_approach_path_rejected(self):
        # tilt=90 means the DV axis is horizontal — nothing to intersect.
        with pytest.raises(ValueError, match="DV component"):
            target3d([1.0, 0.0, -3.0], "LR", 90, 90, dv_intersect=-0.3)

    def test_vertical_path_preserves_xy(self):
        result = plane_intersection([1.0, 1.0, -5.0], [0.0, 0.0, 1.0], 0.0)
        np.testing.assert_allclose(result, [1.0, 1.0, 0.0])

    def test_45_degree_path(self):
        # Path from (0, 0, -10) along (1, 0, 1): reaches DV=0 at ML=10.
        result = plane_intersection([0.0, 0.0, -10.0], [1.0, 0.0, 1.0], 0.0)
        np.testing.assert_allclose(result, [10.0, 0.0, 0.0])

    @pytest.mark.parametrize("bad_value", [np.nan, np.inf, -np.inf])
    def test_dv_intersect_non_finite_rejected(self, bad_value):
        with pytest.raises(ValueError, match="dv_intersect must be finite"):
            plane_intersection([1.0, 1.0, -5.0], [0.0, 0.0, 1.0], bad_value)

    def test_zero_magnitude_direction_rejected(self):
        with pytest.raises(ValueError, match="DV component"):
            plane_intersection([1.0, 1.0, -5.0], [0.0, 0.0, 0.0], 0.0)

    def test_relatively_negligible_dv_component_rejected(self):
        # [1e3, 1e3, 1e-7] is essentially horizontal even though abs(DV) > 0.
        # The new relative-tolerance check must catch this.
        with pytest.raises(ValueError, match="DV component"):
            plane_intersection([1.0, 1.0, -5.0], [1e3, 1e3, 1e-7], 0.0)

    def test_invariant_under_direction_scaling(self):
        # The intersection point should not depend on how the caller
        # scales the direction vector.
        result_small = plane_intersection([0.0, 0.0, -10.0], [1.0, 0.0, 1.0], 0.0)
        result_large = plane_intersection([0.0, 0.0, -10.0], [42.0, 0.0, 42.0], 0.0)
        np.testing.assert_allclose(result_small, result_large)


# ---------------------------------------------------------------------------
# stereotaxic_basis
# ---------------------------------------------------------------------------
class TestStereotaxicBasis:
    def test_lr_lateral_right_is_orthonormal(self):
        basis = stereotaxic_basis("LR", 20, 90)
        np.testing.assert_allclose(basis.T @ basis, np.eye(3), atol=1e-12)

    def test_lr_lateral_left_is_orthonormal(self):
        basis = stereotaxic_basis("LR", 30, 270)
        np.testing.assert_allclose(basis.T @ basis, np.eye(3), atol=1e-12)

    def test_fb_basis_not_orthonormal_under_tilt(self):
        # FB convention: AP axis stays as atlas [0, 1, 0] regardless of tilt
        # or rotation, so DV . AP = sin(tilt)*cos(rot) is generally nonzero.
        basis = stereotaxic_basis("FB", 20, 0)
        gram = basis.T @ basis
        assert abs(gram[1, 2]) == pytest.approx(np.sin(np.deg2rad(20)), abs=1e-12)

    def test_fb_off_axis_basis_not_orthonormal(self):
        basis = stereotaxic_basis("FB", 20, 45)
        assert not np.allclose(basis.T @ basis, np.eye(3), atol=1e-3)

    def test_column_order_is_ml_ap_dv(self):
        # Multiplying the basis by a unit column-selector should return
        # the corresponding column.
        basis = stereotaxic_basis("LR", 20, 90)
        np.testing.assert_allclose(basis @ np.array([1.0, 0.0, 0.0]), basis[:, 0])
        np.testing.assert_allclose(basis @ np.array([0.0, 1.0, 0.0]), basis[:, 1])
        np.testing.assert_allclose(basis @ np.array([0.0, 0.0, 1.0]), basis[:, 2])

    @pytest.mark.parametrize("arm", ["LR", "FB"])
    @pytest.mark.parametrize("tilt_degrees", [10, 30, 60])
    @pytest.mark.parametrize("rotation_degrees", [30, 45, 120, 200, 240])
    def test_ap_axis_is_atlas_ap_for_all_valid_configs(
        self, arm, tilt_degrees, rotation_degrees
    ):
        # By design the manipulator's AP stage moves along atlas AP,
        # regardless of arm config, tilt, or rotation.
        basis = stereotaxic_basis(arm, tilt_degrees, rotation_degrees)
        np.testing.assert_allclose(basis[:, 1], [0.0, 1.0, 0.0])

    def test_dv_axis_for_lr_lateral_matches_expected(self):
        # LR tilted 30° right (rot=90): DV axis tips into +ML with the
        # original cos(30°) component remaining along atlas DV.
        basis = stereotaxic_basis("LR", 30, 90)
        expected_dv = np.array([np.sin(np.deg2rad(30)), 0.0, np.cos(np.deg2rad(30))])
        np.testing.assert_allclose(basis[:, 2], expected_dv, atol=1e-12)

    def test_forward_lean_dv_points_anterior(self):
        # rot=30 is a lean toward the nose. The manipulator's +DV axis
        # should have positive AP component.
        basis = stereotaxic_basis("FB", 30, 30)
        assert basis[1, 2] > 0

    def test_backward_lean_dv_points_posterior(self):
        # rot=150 is a lean toward the tail.
        basis = stereotaxic_basis("FB", 30, 150)
        assert basis[1, 2] < 0


# ---------------------------------------------------------------------------
# Return type and backward compatibility
# ---------------------------------------------------------------------------
class TestReturnType:
    def test_namedtuple_unpacks_as_two_tuple(self):
        # Notebook style: `R, B = target3d(...)` must still work.
        target_coords, intersection = target3d([0.25, 0.7, -6.6], "LR", 20, 90)
        assert target_coords.shape == (3,)
        assert intersection is None

    def test_result_is_target3d_result_instance(self):
        result = target3d([0.25, 0.7, -6.6], "LR", 20, 90, dv_intersect=-0.3)
        assert isinstance(result, Target3DResult)
        assert isinstance(result.target_coordinates, np.ndarray)
        assert isinstance(result.plane_intersection, np.ndarray)

    def test_plane_intersection_none_without_dv_intersect(self):
        result = target3d([0.25, 0.7, -6.6], "LR", 20, 90)
        assert result.plane_intersection is None

    def test_target_coordinates_dtype_is_float64(self):
        result = target3d([0.25, 0.7, -6.6], "LR", 20, 90)
        assert result.target_coordinates.dtype == np.float64


# ---------------------------------------------------------------------------
# Verbose printing
# ---------------------------------------------------------------------------
class TestVerboseOutput:
    def test_silent_by_default(self, capsys):
        target3d([0.25, 0.7, -6.6], "LR", 20, 90)
        captured = capsys.readouterr()
        assert captured.out == ""

    def test_verbose_prints_target_coordinates(self, capsys):
        target3d([0.25, 0.7, -6.6], "LR", 20, 90, verbose=True)
        captured = capsys.readouterr()
        assert "Target" in captured.out
        assert "+2.49" in captured.out  # ML component
        assert "-6.12" in captured.out  # DV component

    def test_verbose_prints_intersection_when_supplied(self, capsys):
        target3d([-2, -2, -3.1], "FB", 30, -20, dv_intersect=-0.3, verbose=True)
        captured = capsys.readouterr()
        assert "Plane intersection" in captured.out
        assert "-0.30" in captured.out  # DV value of the intersection

    def test_verbose_omits_intersection_when_not_supplied(self, capsys):
        target3d([0.25, 0.7, -6.6], "LR", 20, 90, verbose=True)
        captured = capsys.readouterr()
        assert "Plane intersection" not in captured.out


# ---------------------------------------------------------------------------
# spherical_to_cartesian
# ---------------------------------------------------------------------------
class TestSphericalToCartesian:
    def test_zero_angles_return_ml_unit_vector(self):
        v = spherical_to_cartesian(0.0, 0.0)
        np.testing.assert_allclose(v, [1.0, 0.0, 0.0], atol=1e-12)

    def test_phi_pi_over_2_returns_dv_unit_vector(self):
        v = spherical_to_cartesian(0.0, np.pi / 2)
        np.testing.assert_allclose(v, [0.0, 0.0, 1.0], atol=1e-12)

    def test_theta_pi_over_2_returns_ap_unit_vector(self):
        v = spherical_to_cartesian(np.pi / 2, 0.0)
        np.testing.assert_allclose(v, [0.0, 1.0, 0.0], atol=1e-12)

    def test_radius_scales_output(self):
        v = spherical_to_cartesian(0.0, 0.0, radius=2.5)
        np.testing.assert_allclose(v, [2.5, 0.0, 0.0], atol=1e-12)

    def test_returns_float64_array(self):
        v = spherical_to_cartesian(0.3, 0.4)
        assert v.dtype == np.float64
        assert v.shape == (3,)

    @pytest.mark.parametrize(
        "theta,phi,radius",
        [
            (np.nan, 0.0, 1.0),
            (np.inf, 0.0, 1.0),
            (0.0, np.nan, 1.0),
            (0.0, np.inf, 1.0),
            (0.0, 0.0, np.nan),
            (0.0, 0.0, np.inf),
        ],
    )
    def test_non_finite_inputs_rejected(self, theta, phi, radius):
        with pytest.raises(ValueError, match="must be finite"):
            spherical_to_cartesian(theta, phi, radius)


# ---------------------------------------------------------------------------
# format_coordinates
# ---------------------------------------------------------------------------
class TestFormatCoordinates:
    def test_signed_two_decimal_format(self):
        s = format_coordinates("Target", [1.0, -2.5, 3.14])
        assert s == "Target: ML:+1.00 / AP:-2.50 / DV:+3.14"

    def test_label_appears_first(self):
        s = format_coordinates("Entry", [0.0, 0.0, 0.0])
        assert s.startswith("Entry:")

    def test_zero_renders_with_explicit_sign(self):
        s = format_coordinates("X", [0.0, 0.0, 0.0])
        assert "+0.00" in s

    def test_rounds_to_two_decimals(self):
        s = format_coordinates("X", [1.234, -5.678, 9.0])
        assert "+1.23" in s
        assert "-5.68" in s
        assert "+9.00" in s

    def test_accepts_numpy_array(self):
        s = format_coordinates("X", np.array([1.0, 2.0, 3.0]))
        assert "+1.00" in s and "+2.00" in s and "+3.00" in s

    def test_accepts_tuple(self):
        s = format_coordinates("X", (1.0, 2.0, 3.0))
        assert "+1.00" in s

    def test_rejects_wrong_shape(self):
        with pytest.raises(ValueError, match="coordinates"):
            format_coordinates("X", [1.0, 2.0])

    def test_rejects_non_finite(self):
        with pytest.raises(ValueError, match="finite"):
            format_coordinates("X", [1.0, np.nan, 3.0])


# ---------------------------------------------------------------------------
# Defensive error handling
# ---------------------------------------------------------------------------
class TestDefensiveErrorHandling:
    """The ``target3d`` body wraps :func:`numpy.linalg.solve` in a
    defensive ``try/except`` so that any singularity not caught by
    :func:`stereotaxic_basis` still surfaces as a ``ValueError`` rather
    than a raw ``LinAlgError``. This branch should not fire in practice,
    but the conversion path must work."""

    def test_singular_basis_converted_to_value_error(self, monkeypatch):
        import stereotaxic

        # Force ``stereotaxic_basis`` to return a singular (all-zero) matrix
        # so that ``np.linalg.solve`` raises ``LinAlgError``.
        monkeypatch.setattr(
            stereotaxic,
            "stereotaxic_basis",
            lambda *args, **kwargs: np.zeros((3, 3), dtype=np.float64),
        )
        with pytest.raises(ValueError, match="should not occur"):
            stereotaxic.target3d([1.0, 2.0, -3.0], "LR", 20, 90)

    def test_underlying_linalg_error_chained(self, monkeypatch):
        # The LAPACK message should be preserved via ``raise ... from exc``
        # so that the original error survives in ``__cause__``.
        import stereotaxic

        monkeypatch.setattr(
            stereotaxic,
            "stereotaxic_basis",
            lambda *args, **kwargs: np.zeros((3, 3), dtype=np.float64),
        )
        with pytest.raises(ValueError) as exc_info:
            stereotaxic.target3d([1.0, 2.0, -3.0], "LR", 20, 90)
        assert isinstance(exc_info.value.__cause__, np.linalg.LinAlgError)


# ---------------------------------------------------------------------------
# Round-trip via the basis (sanity)
# ---------------------------------------------------------------------------
class TestRoundTripWithBasis:
    """``basis @ result.target_coordinates`` should reproduce the input
    target to machine precision for any non-singular geometry."""

    @pytest.mark.parametrize(
        "arm,tilt_degrees,rotation_degrees",
        [
            ("LR", 20, 90),
            ("LR", 45, 45),
            ("LR", 80, 89),  # near-singular but well-defined
            ("FB", 30, -20),
            ("FB", 25, 150),
            ("FB", 60, 200),
        ],
    )
    def test_roundtrip(self, arm, tilt_degrees, rotation_degrees):
        target = np.array([1.5, -2.3, -4.7])
        result = target3d(target, arm, tilt_degrees, rotation_degrees)
        basis = stereotaxic_basis(arm, tilt_degrees, rotation_degrees)
        recovered = basis @ result.target_coordinates
        np.testing.assert_allclose(recovered, target, atol=1e-12)