Foreshorten to the true perspective eye at view_distance / strength

docs/changelog.rst

@@ -24,6 +24,17 @@ Unreleased
   ``P`` returns to an orthographic view at the bottom of the ladder,
   as before.
 
+- Under perspective, a sphere's silhouette now converges to its
+  orthographic size as the perspective strength falls, measured from
+  the eye at ``view_distance / strength``.  Previously the apparent
+  size did not approach the orthographic size correctly.  A sphere
+  large enough to contain the eye now warns rather than ballooning
+  silently.
+
+- Under perspective, bond end caps now foreshorten to the same
+  corrected eye as the sphere silhouettes, so their drawn shape changes
+  at perspective strengths other than one.
+
 0.20.0
 ------
 

src/hofmann/model/projection.py

@@ -3,6 +3,7 @@
 from __future__ import annotations
 
 import math
+import warnings
 from abc import ABC, abstractmethod
 from dataclasses import dataclass
 
@@ -14,6 +15,17 @@
 _PARALLEL_EYE_DISTANCE = 1e6
 
 
+def _sqrt_difference_of_squares(a: np.ndarray, b: np.ndarray) -> np.ndarray:
+    """``sqrt(a**2 - b**2)`` for non-negative *a*, *b*, without overflow.
+
+    Evaluated as ``sqrt(a - b) * sqrt(a + b)``: forming ``a**2 - b**2``
+    first overflows to ``inf`` for *a* beyond about 1.3e154, and the
+    caller's division by that ``inf`` then drives the silhouette to zero.
+    ``a < b`` clamps to zero (the eye is inside the sphere; the caller warns).
+    """
+    return np.sqrt(np.maximum(a - b, 0.0)) * np.sqrt(a + b)
+
+
 class Projection(ABC):
     """A camera projection mode.
 
@@ -119,11 +131,20 @@ def silhouette_radius(
         self, depth: np.ndarray, radii: np.ndarray,
     ) -> np.ndarray:
         d = self.view_distance - depth * self.strength
-        # Silhouette radius r*D/sqrt(d^2 - r^2), approximate: the eye is
-        # at D/strength, for which the exact form carries (r*strength)^2.
-        # Bond end caps use the same reference distance D (bond_geometry),
-        # so the two share an eye.
-        denom = np.sqrt(np.maximum(d**2 - radii**2, 1e-12))
+        rs = radii * self.strength
+        if np.any(np.abs(d) <= rs):
+            warnings.warn(
+                "a sphere contains the perspective eye and is drawn "
+                "with an unbounded silhouette; increase view_distance "
+                "or reduce strength.",
+                UserWarning,
+                stacklevel=2,
+            )
+        # Silhouette radius r*D/sqrt(d^2 - (r*s)^2), the eye at D/s.
+        # abs(d) keeps a real denominator for atoms behind the eye
+        # (d < 0); the 1e-6 floor gives a huge finite radius when the
+        # eye is inside a sphere.
+        denom = np.maximum(_sqrt_difference_of_squares(np.abs(d), rs), 1e-6)
         return radii * self.view_distance / denom
 
     def max_magnification(self, worst_depth: float) -> float:
@@ -132,9 +153,7 @@ def max_magnification(self, worst_depth: float) -> float:
 
     @property
     def eye_distance(self) -> float:
-        # Exact only at full strength: the eye actually sits at
-        # view_distance / strength.
-        return self.view_distance
+        return self.view_distance / self.strength
 
     def reaches_eye_plane(self, depth: np.ndarray) -> bool:
         return bool(np.any(self.view_distance - depth * self.strength <= 0))

src/hofmann/model/view_state.py

@@ -107,20 +107,16 @@ def project_camera(self, camera: np.ndarray) -> np.ndarray:
             *xy* of shape ``(n, 2)`` — screen positions with zoom applied.
         """
         camera = np.asarray(camera, dtype=float)
-        proj = self.projection
-        if proj.reaches_eye_plane(camera[:, 2]):
-            assert isinstance(proj, Perspective)  # only Perspective reaches it
+        if self.projection.reaches_eye_plane(camera[:, 2]):
             warnings.warn(
-                "one or more points lie at or behind the "
-                f"perspective eye plane (view_distance="
-                f"{proj.view_distance:g}, strength={proj.strength:g}); "
-                "they are drawn mirrored through the origin and "
-                "sorted as if nearest the viewer.  Increase "
+                "one or more points lie at or behind the perspective "
+                "eye plane; they are drawn mirrored through the origin "
+                "and sorted as if nearest the viewer.  Increase "
                 "view_distance or reduce strength.",
                 UserWarning,
                 stacklevel=2,
             )
-        return proj.to_screen(camera) * self.zoom
+        return self.projection.to_screen(camera) * self.zoom
 
     def set_orthographic(self) -> ViewState:
         """Draw without perspective foreshortening.

tests/test_model/test_view_state.py

@@ -456,3 +456,81 @@ def test_perspective_to_screen_not_clamped_at_or_behind_eye(self):
         assert not np.all(np.isfinite(on_plane))
         behind = persp.to_screen(np.array([[1.0, 0.0, 11.0]]))  # depth > D
         assert behind[0, 0] < 0  # mirrored through the origin
+
+    def test_silhouette_converges_to_orthographic_as_strength_falls(self):
+        depth = np.array([2.0])
+        radii = np.array([1.5])
+        ortho = Orthographic().silhouette_radius(depth, radii)  # == radii
+        r_strong = Perspective(1.0, 10.0).silhouette_radius(depth, radii)
+        r_weak = Perspective(0.01, 10.0).silhouette_radius(depth, radii)
+        assert abs(r_weak[0] - ortho[0]) < abs(r_strong[0] - ortho[0])
+        np.testing.assert_allclose(
+            Perspective(1e-6, 10.0).silhouette_radius(depth, radii),
+            ortho, rtol=1e-3,
+        )
+
+    def test_silhouette_pins_effective_eye_below_full_strength(self):
+        # eye at D/s = 20; d = D - depth*s = 8, rs = r*s = 0.75, so the
+        # radius is 15 / sqrt(64 - 0.5625) = 1.8833.
+        np.testing.assert_allclose(
+            Perspective(0.5, 10.0).silhouette_radius(
+                np.array([4.0]), np.array([1.5])
+            ),
+            [1.8833], rtol=1e-4,
+        )
+
+    def test_silhouette_finite_at_huge_view_distance(self):
+        r = Perspective(1.0, 1e200).silhouette_radius(
+            np.array([0.0]), np.array([1.5])
+        )
+        assert np.all(np.isfinite(r)) and np.all(r > 0)
+
+    def test_silhouette_finite_for_atoms_behind_the_eye(self):
+        r = Perspective(1.0, 10.0).silhouette_radius(
+            np.array([20.0]), np.array([1.5])
+        )
+        assert np.all(np.isfinite(r))
+
+    def test_sphere_containing_the_eye_warns(self):
+        with pytest.warns(UserWarning, match="contains the perspective eye"):
+            r = Perspective(1.0, 10.0).silhouette_radius(
+                np.array([9.5]), np.array([1.0])
+            )
+        # the 1e-6 denominator floor keeps the radius huge but finite
+        assert np.all(np.isfinite(r)) and np.all(r >= 1e6)
+
+    def test_sphere_warning_dedups_across_calls(self):
+        """One warning line for many eye-containing spheres, not one each."""
+        with warnings.catch_warnings(record=True) as caught:
+            warnings.simplefilter("default")
+            for view_distance in (10.0, 11.0, 12.0, 13.0, 14.0):
+                Perspective(1.0, view_distance).silhouette_radius(
+                    np.array([view_distance - 0.5]), np.array([1.0])
+                )
+        contained = [
+            w for w in caught if "contains the perspective eye" in str(w.message)
+        ]
+        assert len(contained) == 1
+
+    def test_ordinary_atom_does_not_warn(self):
+        with warnings.catch_warnings():
+            warnings.simplefilter("error")
+            Perspective(1.0, 10.0).silhouette_radius(
+                np.array([0.0]), np.array([1.0])
+            )
+
+    def test_eye_distance_is_the_effective_pinhole(self):
+        assert Perspective(1.0, 10.0).eye_distance == 10.0
+        assert Perspective(0.5, 5.0).eye_distance == 10.0
+        assert Perspective(0.5, 10.0).eye_distance == 20.0
+
+    def test_eye_plane_warning_dedups_as_camera_varies(self):
+        """A drag that changes view_distance each frame emits one warning
+        line, not one per frame."""
+        behind = np.array([[1.0, 0.0, 100.0]])  # behind the eye at every step
+        with warnings.catch_warnings(record=True) as caught:
+            warnings.simplefilter("default")
+            for view_distance in (10.0, 11.0, 12.0, 13.0, 14.0):
+                ViewState(projection=Perspective(1.0, view_distance)).project(behind)
+        eye = [w for w in caught if "eye plane" in str(w.message)]
+        assert len(eye) == 1

tests/test_rendering/test_bond_geometry.py

@@ -2,7 +2,7 @@
 
 import numpy as np
 
-from hofmann.model import BondSpec, ViewState
+from hofmann.model import BondSpec, Perspective, ViewState
 from hofmann.rendering.bond_geometry import (
     _bond_polygon,
     _bond_polygons_batch,
@@ -237,6 +237,35 @@ def test_matches_scalar_rotated_view(self):
             s_verts, s_start, s_end = scalar
             np.testing.assert_allclose(full_verts[i], s_verts, atol=1e-12)
 
+    def test_matches_scalar_perspective_view(self):
+        """Batch matches scalar under perspective: the effective-eye
+        foreshortening must be applied identically on both paths."""
+        rot = _rotation_y(0.7) @ _rotation_x(0.3)
+        view = ViewState(rotation=rot, projection=Perspective(0.6, 12.0))
+        d = self._ch4_scene_data(view=view)
+        atom_scale = d["atom_scale"]
+
+        full_verts, start_2d, end_2d, *_, valid = _bond_polygons_batch(
+            d["rotated"], d["xy"],
+            d["radii_3d"] * atom_scale, d["screen_radii"],
+            d["bond_ia"], d["bond_ib"], d["bond_radii"],
+            view,
+        )
+        for i, bond in enumerate(d["bonds"]):
+            ia, ib = bond.index_a, bond.index_b
+            scalar = _bond_polygon(
+                d["rotated"][ia], d["rotated"][ib],
+                d["radii_3d"][ia] * atom_scale,
+                d["radii_3d"][ib] * atom_scale,
+                bond.spec.radius,
+                d["screen_radii"][ia], d["screen_radii"][ib],
+                view,
+            )
+            assert scalar is not None
+            assert valid[i]
+            s_verts, s_start, s_end = scalar
+            np.testing.assert_allclose(full_verts[i], s_verts, atol=1e-9)
+
     def test_empty_bonds(self):
         """No bonds produces empty arrays."""
         view = ViewState()