"""Preetham analytic sky: closed-form daytime sky radiance for a dynamic,
procedural environment cubemap (design D12).
Pure NumPy, no coefficient tables and no GPU: the Perez distribution
coefficients and the zenith luminance/chromaticity are the small linear /
cubic closed forms from Preetham et al. 1999, evaluated per texel. Both
renderer backends synthesize the same six float32 RGBA cube faces from this
one module (desktop feeds them to ``Engine.load_cubemap(faces=...)`` and thus
the existing skybox + IBL path; the web twin uploads them to its cubemap),
so a ``sky_mode="procedural"`` sky and its IBL ambient match across backends.
The faces are returned in Vulkan cube order ``[+X, -X, +Y, -Y, +Z, -Z]`` to
match :func:`simvx.graphics...cubemap_loader.gradient_cubemap_faces`.
"""
from __future__ import annotations
import numpy as np
__all__ = ["preetham_cubemap_faces", "sky_cache_key"]
# --- Perez distribution coefficients (linear in turbidity T) --------------
#
# Each row is (A, B, C, D, E) expressed as a·T + b, i.e. two constants per
# coefficient. From Preetham et al., "A Practical Analytic Model for Daylight".
# Y = luminance, x/y = CIE chromaticity.
_PEREZ_Y = np.array(
[
(0.1787, -1.4630), # A
(-0.3554, 0.4275), # B
(-0.0227, 5.3251), # C
(0.1206, -2.5771), # D
(-0.0670, 0.3703), # E
],
dtype=np.float64,
)
_PEREZ_X = np.array(
[
(-0.0193, -0.2592),
(-0.0665, 0.0008),
(-0.0004, 0.2125),
(-0.0641, -0.8989),
(-0.0033, 0.0452),
],
dtype=np.float64,
)
_PEREZ_Y_CHROMA = np.array(
[
(-0.0167, -0.2608),
(-0.0950, 0.0092),
(-0.0079, 0.2102),
(-0.0441, -1.6537),
(-0.0109, 0.0529),
],
dtype=np.float64,
)
# Zenith chromaticity polynomial matrices (multiply [T^2, T, 1] by a cubic in
# the sun zenith angle theta_s). From the same paper.
_ZENITH_X = np.array(
[
[0.00166, -0.00375, 0.00209, 0.0],
[-0.02903, 0.06377, -0.03202, 0.00394],
[0.11693, -0.21196, 0.06052, 0.25886],
],
dtype=np.float64,
)
_ZENITH_Y = np.array(
[
[0.00275, -0.00610, 0.00317, 0.0],
[-0.04214, 0.08970, -0.04153, 0.00516],
[0.15346, -0.26756, 0.06670, 0.26688],
],
dtype=np.float64,
)
# CIE XYZ -> linear sRGB (D65).
_XYZ_TO_RGB = np.array(
[
[3.2406, -1.5372, -0.4986],
[-0.9689, 1.8758, 0.0415],
[0.0557, -0.2040, 1.0570],
],
dtype=np.float64,
)
# Vulkan cube-face basis: for each face, (right, up, forward) so a texel at
# uv in [-1, 1]^2 maps to direction normalize(forward + u*right + v*up). These
# reproduce the exact directions of ``gradient_cubemap_faces`` +X..-Z.
_FACE_BASIS = [
((0.0, 0.0, -1.0), (0.0, -1.0, 0.0), (1.0, 0.0, 0.0)), # +X
((0.0, 0.0, 1.0), (0.0, -1.0, 0.0), (-1.0, 0.0, 0.0)), # -X
((1.0, 0.0, 0.0), (0.0, 0.0, 1.0), (0.0, 1.0, 0.0)), # +Y
((1.0, 0.0, 0.0), (0.0, 0.0, -1.0), (0.0, -1.0, 0.0)), # -Y
((1.0, 0.0, 0.0), (0.0, -1.0, 0.0), (0.0, 0.0, 1.0)), # +Z
((-1.0, 0.0, 0.0), (0.0, -1.0, 0.0), (0.0, 0.0, -1.0)), # -Z
]
# Below-horizon Perez blows up (cos(theta) crosses zero), so clamp the zenith
# cosine to a small floor and shade genuinely-downward rays as ground bounce.
_COS_THETA_FLOOR = 0.045
# Night colour the sky fades toward once the sun drops below the horizon.
_NIGHT_RGB = np.array([0.008, 0.015, 0.035], dtype=np.float64)
def _perez(cos_theta: np.ndarray, gamma: np.ndarray, cos_gamma: np.ndarray, coeff: np.ndarray) -> np.ndarray:
a, b, c, d, e = coeff
return (1.0 + a * np.exp(b / cos_theta)) * (1.0 + c * np.exp(d * gamma) + e * cos_gamma * cos_gamma)
def _coeff(row_ab: np.ndarray, turbidity: float) -> np.ndarray:
"""Resolve the five (A..E) Perez coefficients at *turbidity*."""
return row_ab[:, 0] * turbidity + row_ab[:, 1]
def _zenith(matrix: np.ndarray, turbidity: float, theta_s: float) -> float:
tvec = np.array([turbidity * turbidity, turbidity, 1.0], dtype=np.float64)
svec = np.array([theta_s**3, theta_s**2, theta_s, 1.0], dtype=np.float64)
return float(tvec @ matrix @ svec)
[docs]
def sky_cache_key(
sun_direction: tuple[float, float, float],
turbidity: float,
ground_albedo: tuple[float, float, float],
size: int,
) -> tuple:
"""Quantised identity for a synthesized sky.
Quantises the sun direction to ~1.5 degrees and the scalar params to a
coarse grid so a slowly-sweeping sun (the day-night example) re-bakes the
cubemap + IBL only in discrete steps rather than every frame. A static sky
hashes to one key and is synthesized exactly once (zero re-cost).
"""
d = np.asarray(sun_direction, dtype=np.float64)
n = float(np.linalg.norm(d))
if n > 1e-8:
d = d / n
q = tuple(int(round(float(x) * 40.0)) for x in d)
return (q, round(float(turbidity), 2), tuple(round(float(c), 3) for c in ground_albedo[:3]), int(size))
[docs]
def preetham_cubemap_faces(
sun_direction: tuple[float, float, float],
turbidity: float = 2.5,
ground_albedo: tuple[float, float, float] = (0.15, 0.15, 0.16),
size: int = 64,
exposure: float = 0.18,
) -> list[np.ndarray]:
"""Synthesize six Preetham sky cube faces.
Args:
sun_direction: unit(-ish) vector pointing TOWARD the sun in world space
(i.e. the negation of a ``DirectionalLight3D`` travel direction).
turbidity: atmospheric haze, ~2 (clear) .. ~10 (hazy). Clamped to
``[1.7, 10]`` where the analytic model stays well-behaved.
ground_albedo: RGB reflectance shading genuinely-downward rays as a
flat ground bounce (Preetham models only the sky dome).
size: per-face resolution (64-128).
exposure: master linear scale applied to the analytic radiance before
it is written into the HDR cube faces.
Returns:
Six ``(size, size, 4)`` float32 RGBA arrays in Vulkan face order,
non-negative linear radiance (HDR: values may exceed 1 near the sun).
"""
sun = np.asarray(sun_direction, dtype=np.float64)
n = float(np.linalg.norm(sun))
sun = sun / n if n > 1e-8 else np.array([0.0, 1.0, 0.0])
turbidity = float(np.clip(turbidity, 1.7, 10.0))
ground = np.asarray(ground_albedo, dtype=np.float64)[:3]
# Sun zenith angle + a smooth day factor that fades the whole dome to night
# as the sun sinks below the horizon (sun.y crossing 0).
sun_y = float(np.clip(sun[1], -1.0, 1.0))
theta_s = float(np.arccos(np.clip(sun_y, _COS_THETA_FLOOR, 1.0)))
day = float(np.clip((sun_y + 0.12) / 0.24, 0.0, 1.0))
day = day * day * (3.0 - 2.0 * day) # smoothstep
cy = _coeff(_PEREZ_Y, turbidity)
cx = _coeff(_PEREZ_X, turbidity)
cyc = _coeff(_PEREZ_Y_CHROMA, turbidity)
xz = _zenith(_ZENITH_X, turbidity, theta_s)
yz = _zenith(_ZENITH_Y, turbidity, theta_s)
cos_ts = np.array(float(np.cos(theta_s)))
gamma_s = np.array(theta_s)
denom_x = _perez(np.array(1.0), gamma_s, cos_ts, cx)
denom_y = _perez(np.array(1.0), gamma_s, cos_ts, cyc)
# Master luminance target: zenith reads ~1 in linear HDR at exposure 0.18;
# turbidity/elevation modulate through the Perez shape + chromaticity.
grid = (np.arange(size, dtype=np.float64) + 0.5) / size * 2.0 - 1.0
u, v = np.meshgrid(grid, grid)
faces: list[np.ndarray] = []
for right, up, fwd in _FACE_BASIS:
r = np.asarray(right)
p = np.asarray(up)
f = np.asarray(fwd)
d = f[None, None, :] + u[..., None] * r[None, None, :] + v[..., None] * p[None, None, :]
d = d / np.linalg.norm(d, axis=-1, keepdims=True)
dy = d[..., 1]
cos_theta = np.clip(dy, _COS_THETA_FLOOR, 1.0)
cos_gamma = np.clip(np.sum(d * sun[None, None, :], axis=-1), -1.0, 1.0)
gamma = np.arccos(cos_gamma)
rel_x = _perez(cos_theta, gamma, cos_gamma, cx) / denom_x
rel_y = _perez(cos_theta, gamma, cos_gamma, cyc) / denom_y
rel_lum = _perez(cos_theta, gamma, cos_gamma, cy) / _perez(np.array(1.0), gamma_s, cos_ts, cy)
chroma_x = np.clip(xz * rel_x, 1e-3, 0.75)
chroma_y = np.clip(yz * rel_y, 1e-3, 0.75)
lum = np.clip(rel_lum, 0.0, None)
# xyY -> XYZ -> linear sRGB.
big_x = (chroma_x / chroma_y) * lum
big_z = ((1.0 - chroma_x - chroma_y) / chroma_y) * lum
xyz = np.stack([big_x, lum, big_z], axis=-1)
rgb = xyz @ _XYZ_TO_RGB.T
rgb = np.clip(rgb, 0.0, None) * exposure
# Sun disc + layered halo. The sharp disc is a real HDR highlight; the
# broader falloff terms guarantee a visible warm glow even when the disc
# centre falls between texels of a 64px face (and survive the IBL
# prefilter), so the sun reads regardless of its exact direction.
disc = np.clip(cos_gamma, 0.0, 1.0)
halo = np.power(disc, 1500.0) * 16.0 + np.power(disc, 200.0) * 2.4 + np.power(disc, 32.0) * 0.6
sun_rgb = np.array([1.0, 0.86, 0.68])
rgb = rgb + halo[..., None] * sun_rgb[None, None, :] * day
# Ground bounce for downward rays (Preetham only models the dome).
below = np.clip(-dy, 0.0, 1.0)[..., None]
horizon_lum = float(np.clip(yz * 0.35, 0.05, 0.5)) * exposure * 6.0
ground_rgb = ground[None, None, :] * horizon_lum * max(day, 0.15)
rgb = rgb * (1.0 - below) + ground_rgb * below
# Night fade: mix the whole dome toward a dark night colour.
rgb = _NIGHT_RGB[None, None, :] * (1.0 - day) + rgb * day
out = np.empty((size, size, 4), dtype=np.float32)
out[..., :3] = np.clip(rgb, 0.0, 65504.0).astype(np.float32)
out[..., 3] = 1.0
faces.append(np.ascontiguousarray(out))
return faces