A ball's contact shadow¶
The share of the sky's light a ball leaves a point, its contact shadow.
A sphere of radius r whose centre is d from the point p, at theta from p's normal n, covers about (r/d)^2 cos theta of the cosine-weighted hemisphere above p (the solid angle of a small sphere, pi r^2 / d^2, weighted by cos theta over the hemisphere's pi); scaled by the view's strength (0 to 1) and held to 0..1, the share left is 1 - clamp(strength (r/d)^2 max(cos theta, 0), 0, 1), cos theta = n . (c - p) / d. A sphere behind the surface (cos theta <= 0) leaves all of it, as does the sphere the point is on (its centre is behind its own surface) and one whose centre is the point. The shares of several balls multiply. Source: the view factor from a differential area to a small sphere, (r/d)^2 cos theta (Modest, Radiative Heat Transfer, 3rd ed., chapter 4, view factors), as Quilez's analytic sphere occlusion applies it (iquilezles.org). Checked against a reference in numpy, which has its own tests.
- Source
- the view factor from a differential area to a small sphere, (r/d)^2 cos theta (Modest, Radiative Heat Transfer, 3rd ed., chapter 4).
- Proof
- reference, invariant, limit.
occlusion.trait
namespace lusine.render;
float occlusion(vector<float<m>, 3> centre, float<m> radius, vector<float<m>, 3> point, vector<float, 3> normal,
float strength)
{
let v = centre - point;
let d2 = dot(v, v);
if (d2 <= 0) return 1.0;
let facing = dot(normal, v);
let c = (facing > 0 ? facing : 0) / sqrt(d2);
return 1.0 - clamp(strength * (radius * radius / d2) * c, 0.0, 1.0);
}
The share of the sky's light a ball leaves a surface point, 1 - clamp(strength (r/d)^2 max(cos theta, 0), 0, 1): the view factor of a small sphere, so a ball darkens the ground under it and its neighbours softly. 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/occlusion.trait'))
// lusine.render.occlusion: the share of the sky's light a ball leaves a point, its contact shadow. A sphere of radius
// r whose centre is d from the point p, at theta from p's normal n, covers about (r/d)^2 cos theta of the
// cosine-weighted hemisphere above p (the solid angle of a small sphere, pi r^2 / d^2, weighted by cos theta over the
// hemisphere's pi); scaled by the view's `strength` (0 to 1) and held to 0..1, the share left is
// 1 - clamp(strength (r/d)^2 max(cos theta, 0), 0, 1), cos theta = n . (c - p) / d.
// A sphere behind the surface (cos theta <= 0) leaves all of it, as does the sphere the point is on (its centre is
// behind its own surface) and one whose centre is the point. The shares of several balls multiply.
// Source: the view factor from a differential area to a small sphere, (r/d)^2 cos theta (Modest, Radiative Heat
// Transfer, 3rd ed., chapter 4, view factors), as Quilez's analytic sphere occlusion applies it (iquilezles.org).
// Checked against a reference in numpy, which has its own tests.
namespace lusine.render;
float occlusion(vector<float<m>, 3> centre, float<m> radius, vector<float<m>, 3> point, vector<float, 3> normal,
float strength)
{
let v = centre - point;
let d2 = dot(v, v);
if (d2 <= 0) return 1.0;
let facing = dot(normal, v);
let c = (facing > 0 ? facing : 0) / sqrt(d2);
return 1.0 - clamp(strength * (radius * radius / d2) * c, 0.0, 1.0);
}
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/occlusion.trait', 'occlusion'))
from lusine.kernels.render.occlusion import occlusion as reference
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host vector<float<m>, 3> centre;
host float<m> radius;
host vector<float<m>, 3> point;
host vector<float, 3> normal;
host float strength;
float result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.render.occlusion(items.centre, items.radius, items.point, items.normal, items.strength); }
Reference: points around a ball¶
Points on the ground around a 4 cm ball resting on it: the kernel's share within binary32's roundings of the float64 reference's.
def occ(centre, radius, point, normal, strength=1.0):
return trait_run.law('render/occlusion.trait', 'occlusion', centre=centre, radius=radius, point=point,
normal=normal, strength=strength)
c, r = np.array([0.0, 0.04, 0.0]), 0.04
x = np.linspace(-0.4, 0.4, 161)
pts = np.stack([x, np.zeros_like(x), np.zeros_like(x)], axis=1)
up = np.tile([0.0, 1.0, 0.0], (len(x), 1))
got = occ(c, r, pts, up)
want = reference(c, r, pts, up, 1.0)
err = np.max(np.abs(got - want))
assert err <= 16 * U, err / U
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(x, got)
ax.set_xlabel('distance from the contact (m)'); ax.set_ylabel('sky left')
plt.show()
print(f'within {err / U:.1f} x 2^-24 of the reference')
within 0.7 x 2^-24 of the reference
Invariant: (r/d)^2, the view factor¶
A ball straight overhead at a height d covers (r/d)^2 of the cosine-weighted sky, the closed form of a spherical cap's view factor; twice as far, a quarter as much.
hs = np.array([0.2, 0.4, 0.8, 1.6])
taken = 1.0 - occ(np.stack([np.zeros(4), hs, np.zeros(4)], axis=1), 0.04, np.zeros((4, 3)), np.tile([0.0, 1.0, 0.0], (4, 1)))
assert np.allclose(taken, (float(F32(0.04)) / hs) ** 2, rtol=0, atol=4 * U)
print(taken)
[0.04000002 0.00999999 0.0025 0.00062501]
Limit: the contact, and a ball beside or below¶
The point a ball rests on is wholly covered; a ball level with the surface or under it takes nothing, as does strength 0.
out = occ(np.array([[0, 0.04, 0], [0.5, 0.0, 0], [0, -0.5, 0], [0, 0.3, 0]]), 0.04, np.zeros((4, 3)),
np.tile([0.0, 1.0, 0.0], (4, 1)), np.array([1.0, 1.0, 1.0, 0.0]))
assert out[0] == 0.0 and np.all(out[1:] == 1.0)
print('covered at the contact; nothing beside, below, or at strength 0')
covered at the contact; nothing beside, below, or at strength 0