diff --git a/docs/changelog.rst b/docs/changelog.rst index 01a82cf3..887848d1 100644 --- a/docs/changelog.rst +++ b/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 ------ diff --git a/src/hofmann/model/projection.py b/src/hofmann/model/projection.py index be25f632..1401646d 100644 --- a/src/hofmann/model/projection.py +++ b/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)) diff --git a/src/hofmann/model/view_state.py b/src/hofmann/model/view_state.py index f0f54758..b0894443 100644 --- a/src/hofmann/model/view_state.py +++ b/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. diff --git a/tests/test_model/test_view_state.py b/tests/test_model/test_view_state.py index 9233321b..4786ea66 100644 --- a/tests/test_model/test_view_state.py +++ b/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 diff --git a/tests/test_rendering/test_bond_geometry.py b/tests/test_rendering/test_bond_geometry.py index ce6e094e..37db9274 100644 --- a/tests/test_rendering/test_bond_geometry.py +++ b/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()