The sky's radiance¶
The radiance a ray sees where it meets nothing.
The sky's radiance runs from the horizon's at the horizon to the zenith's straight up, linearly in the ray's height (its direction's y), and is the horizon's below it; a ray within the sun's angular radius of the sun also sees the sun's disc, its irradiance spread evenly over the disc's solid angle, 4 pi sin^2(r / 2) (W/(m2 sr), linear RGB). A sky whose zenith and horizon are one radiance and a sun radius of 0 is the uniform sky, bit for bit. What lights a surface is lusine.render.sky_ambient, the same sky's cosine-weighted mean. Source: the disc's radiance, L = E / omega with omega = 2 pi (1 - cos r) = 4 pi sin^2(r / 2), the solid angle of a cone (Pharr, Jakob and Humphreys, Physically Based Rendering, 4th ed., 4.2.1); the gradient is this law's own. Checked against a reference in numpy, which has its own tests.
- Source
- the solid angle of a cone (Pharr, Jakob and Humphreys, Physically Based Rendering, 4th ed., 4.2.1); the gradient is the law's own.
- Proof
- reference, invariant, limit.
sky.trait
namespace lusine.render;
vector<float<W/m2/sr>, 3> sky(vector<float, 3> direction, vector<float<W/m2/sr>, 3> zenith,
vector<float<W/m2/sr>, 3> horizon, vector<float, 3> sun_direction,
vector<float<W/m2>, 3> sun, float<rad> sun_radius)
{
let pi = 3.141592653589793;
let height = direction.y > 0.0 ? direction.y : 0.0;
let gradient = horizon + (zenith - horizon) * height;
if (sun_radius <= 0.0 rad) return gradient;
if (dot(direction, sun_direction) < cos(sun_radius)) return gradient;
let half = sin(sun_radius / 2.0);
return gradient + sun / (4.0 * pi * half * half * 1.0 sr);
}
What a ray that meets nothing sees. The sky runs from the horizon's radiance at the horizon to the zenith's straight up, linearly in the ray's height, and is the horizon's below it; within the sun's angular radius of the sun a ray also sees the sun's disc, its irradiance spread over the disc's solid angle. What lights a surface is lusine.render.sky_ambient, the same sky's cosine-weighted mean. One of the kernels lusine.render.table draws a table's shapes with.
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
F32 = np.float32
U = float(np.finfo(F32).eps) / 2 # binary32's unit roundoff, 2^-24
rng = np.random.default_rng(2026)
The kernel¶
The law as the library states it.
print(trait_run.source('render/sky.trait'))
// lusine.render.sky: the radiance a ray sees where it meets nothing. The sky's radiance runs from the horizon's at the
// horizon to the zenith's straight up, linearly in the ray's height (its direction's y), and is the horizon's below
// it; a ray within the sun's angular radius of the sun also sees the sun's disc, its irradiance spread evenly over the
// disc's solid angle, 4 pi sin^2(r / 2) (W/(m2 sr), linear RGB). A sky whose zenith and horizon are one radiance and a
// sun radius of 0 is the uniform sky, bit for bit. What lights a surface
// is lusine.render.sky_ambient, the same sky's cosine-weighted mean.
// Source: the disc's radiance, L = E / omega with omega = 2 pi (1 - cos r) = 4 pi sin^2(r / 2), the solid angle of a
// cone (Pharr, Jakob and Humphreys, Physically Based Rendering, 4th ed., 4.2.1); the gradient is this law's own.
// Checked against a reference in numpy, which has its own tests.
namespace lusine.render;
vector<float<W/m2/sr>, 3> sky(vector<float, 3> direction, vector<float<W/m2/sr>, 3> zenith,
vector<float<W/m2/sr>, 3> horizon, vector<float, 3> sun_direction,
vector<float<W/m2>, 3> sun, float<rad> sun_radius)
{
let pi = 3.141592653589793;
let height = direction.y > 0.0 ? direction.y : 0.0;
let gradient = horizon + (zenith - horizon) * height;
if (sun_radius <= 0.0 rad) return gradient;
if (dot(direction, sun_direction) < cos(sun_radius)) return gradient;
let half = sin(sun_radius / 2.0);
return gradient + sun / (4.0 * pi * half * half * 1.0 sr);
}
The program¶
Each section runs the kernel over a lattice of items, one a case: the host sets every parameter, and main() binds the kernel over the items.
print(trait_run.law_program('render/sky.trait', 'sky'))
from lusine.kernels.render.sky import sky as reference
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host vector<float, 3> direction;
host vector<float<W/m2/sr>, 3> zenith;
host vector<float<W/m2/sr>, 3> horizon;
host vector<float, 3> sun_direction;
host vector<float<W/m2>, 3> sun;
host float<rad> sun_radius;
vector<float<W/m2/sr>, 3> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.render.sky(items.direction, items.zenith, items.horizon, items.sun_direction, items.sun, items.sun_radius); }
Reference: the gradient in every direction¶
Rays in every direction, a sky bluer at the zenith than at the horizon and no disc: each within a few of binary32's roundings of the reference.
def sky(direction, zenith, horizon, sun_direction, sun, sun_radius):
return trait_run.law('render/sky.trait', 'sky', direction=direction, zenith=zenith, horizon=horizon,
sun_direction=sun_direction, sun=sun, sun_radius=sun_radius)
v = rng.normal(size=(300, 3)); v /= np.linalg.norm(v, axis=1, keepdims=True)
v = F32(v).astype(float)
Z, H = F32([0.15, 0.25, 0.55]).astype(float), F32([0.5, 0.55, 0.6]).astype(float)
s = F32(np.array([-0.4, 0.85, -0.3]) / np.linalg.norm([-0.4, 0.85, -0.3])).astype(float)
E = F32([3.2, 3.0, 2.7]).astype(float)
got = sky(v, Z, H, s, E, 0.0)
want = reference(v, Z, H, s, E, 0.0)
err = np.abs(got - want) / np.abs(want)
assert err.max() < 8 * U, err.max()
print(f'every direction within {err.max() / U:.1f} roundings of the reference')
every direction within 1.3 roundings of the reference
Limit: a uniform sky¶
A zenith equal to the horizon and no disc is the uniform sky the table view drew before, bit for bit.
flat = F32([0.25, 0.3, 0.4]).astype(float)
assert np.array_equal(sky(v, flat, flat, s, E, 0.0), np.broadcast_to(flat, v.shape))
print('uniform, bit for bit')
uniform, bit for bit
Reference: the sun's disc¶
The sun 0.533 degrees across: rays inside its angular radius see the disc's radiance over the gradient, rays half a radius past its rim see none of it.
r = F32(np.radians(0.2665)).astype(float)
a = np.cross(s, [0.0, 0.0, 1.0]); a /= np.linalg.norm(a)
angles = np.array([0.0, 0.5, 0.9, 1.5, 3.0]) * r
rays = F32(np.cos(angles)[:, None] * s + np.sin(angles)[:, None] * a).astype(float)
got = sky(rays, Z, H, s, E, r)
want = reference(rays, Z, H, s, E, r)
disc = got - sky(rays, Z, H, s, E, 0.0)
assert np.all(disc[:3] > 1000) and np.all(disc[3:] == 0), disc
assert np.allclose(got, want, rtol=1e-5), np.abs(got / want - 1).max()
print('the disc: inside', disc[0], 'past the rim', disc[-1])
the disc: inside [47081.65323949 44139.05233029 39725.1462217 ] past the rim [0. 0. 0.]
Invariant: the horizon below¶
Every ray below the horizon sees the horizon's radiance, whatever its depth.
down = v[v[:, 1] <= 0]
assert np.array_equal(sky(down, Z, H, s, E, 0.0), np.broadcast_to(H, down.shape))
print(len(down), 'rays below the horizon, the horizon each')
135 rays below the horizon, the horizon each