A soft shadow¶
The share of the sun's disc a ball leaves a point, its soft shadow.
Seen from the point p, the sun is a disc of angular radius a about the direction l toward it, and a ball of radius r whose centre is d away a disc of angular radius b = asin(r / d) about the direction toward its centre, theta from l (taken by atan2 of the cross and dot products, so a small angle keeps its digits in binary32). The ball hides the two discs' overlap: nothing where theta >= a + b, all of the sun where the ball's disc holds it (theta <= b - a), the ball's own disc where the sun's holds it (theta <= a - b), and between them the lens of two circles theta apart; the share left is 1 - overlap / (pi a^2). The discs are taken as flat circles in angle (the sun's 0.27 degrees makes the error small; a ball close to the point is a large disc whose edge the sun's crosses nearly straight). A ball behind the point from the sun (its centre at 90 degrees or more from l), the ball the point is on among them, leaves all of it. The shares of several balls multiply, as lusine.render.occlusion's do. Source: the area of intersection of two circles (Weisstein, Circle-Circle Intersection, MathWorld); the sun's disc as an extended source over a sphere's angular disc (Pharr, Jakob and Humphreys, Physically Based Rendering, 4th ed., 12.2, area lights). Checked against a reference in numpy, which has its own tests.
- Source
- Weisstein, Circle-Circle Intersection (MathWorld).
- Proof
- reference, invariant, limit.
sun_share.trait
namespace lusine.render;
float sun_share(vector<float<m>, 3> centre, float<m> radius, vector<float<m>, 3> point, vector<float, 3> sun_direction,
float<rad> sun_radius)
{
let pi = 3.141592653589793;
let v = centre - point;
let toward = dot(v, sun_direction);
if (toward <= 0.0 m) return 1.0;
let d = sqrt(dot(v, v));
if (d <= radius) return 1.0;
let a = sun_radius;
let b = asin(radius / d);
let side = cross(v, sun_direction);
let theta = atan2(sqrt(dot(side, side)), toward);
if (theta >= a + b) return 1.0;
if (theta <= b - a) return 0.0;
if (theta <= a - b) return 1.0 - (b * b) / (a * a);
let lens = a * a * acos((theta * theta + a * a - b * b) / (2.0 * theta * a))
+ b * b * acos((theta * theta + b * b - a * a) / (2.0 * theta * b))
- 0.5 * sqrt((a + b - theta) * (theta + a - b) * (theta - a + b) * (theta + a + b));
return 1.0 - lens / (pi * a * a);
}
The share of the sun's disc a ball leaves a point, the overlap of the sun's angular disc and the ball's taken away. lusine.render.table multiplies every ball's share into a point's sunlight where the view gives the sun a radius, so a ball's shadow on the sand has a penumbra that widens with its distance.
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/sun_share.trait'))
// lusine.render.sun_share: the share of the sun's disc a ball leaves a point, its soft shadow. Seen from the point p,
// the sun is a disc of angular radius a about the direction l toward it, and a ball of radius r whose centre is d
// away a disc of angular radius b = asin(r / d) about the direction toward its centre, theta from l (taken by atan2 of
// the cross and dot products, so a small angle keeps its digits in binary32). The ball hides the two discs' overlap:
// nothing where theta >= a + b, all of the sun where the ball's disc holds it (theta <= b - a), the ball's own disc
// where the sun's holds it (theta <= a - b), and between them the lens of two circles theta apart; the share left is
// 1 - overlap / (pi a^2). The discs are taken as flat circles in angle (the sun's 0.27 degrees makes the error small;
// a ball close to the point is a large disc whose edge the sun's crosses nearly straight). A ball behind the point
// from the sun (its centre at 90 degrees or more from l), the ball the point is on among them, leaves all of it. The
// shares of several balls multiply, as lusine.render.occlusion's do.
// Source: the area of intersection of two circles (Weisstein, Circle-Circle Intersection, MathWorld); the sun's disc
// as an extended source over a sphere's angular disc (Pharr, Jakob and Humphreys, Physically Based Rendering, 4th ed.,
// 12.2, area lights).
// Checked against a reference in numpy, which has its own tests.
namespace lusine.render;
float sun_share(vector<float<m>, 3> centre, float<m> radius, vector<float<m>, 3> point, vector<float, 3> sun_direction,
float<rad> sun_radius)
{
let pi = 3.141592653589793;
let v = centre - point;
let toward = dot(v, sun_direction);
if (toward <= 0.0 m) return 1.0;
let d = sqrt(dot(v, v));
if (d <= radius) return 1.0;
let a = sun_radius;
let b = asin(radius / d);
let side = cross(v, sun_direction);
let theta = atan2(sqrt(dot(side, side)), toward);
if (theta >= a + b) return 1.0;
if (theta <= b - a) return 0.0;
if (theta <= a - b) return 1.0 - (b * b) / (a * a);
let lens = a * a * acos((theta * theta + a * a - b * b) / (2.0 * theta * a))
+ b * b * acos((theta * theta + b * b - a * a) / (2.0 * theta * b))
- 0.5 * sqrt((a + b - theta) * (theta + a - b) * (theta - a + b) * (theta + a + b));
return 1.0 - lens / (pi * a * a);
}
The program¶
print(trait_run.law_program('render/sun_share.trait', 'sun_share'))
from lusine.kernels.render.sun_share import sun_share 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> sun_direction;
host float<rad> sun_radius;
float result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.render.sun_share(items.centre, items.radius, items.point, items.sun_direction, items.sun_radius); }
Reference: across the penumbra¶
A ball 1.5 times the sun's angular size slid off the sun's disc, from dead ahead to clear of it: the Trait kernel within a small fraction of the share of the reference at every offset (binary32's sines and arc cosines near their limits move the edge by a few roundings of the angle).
def sun_share(centre, radius, point, sun_direction, sun_radius):
return trait_run.law('render/sun_share.trait', 'sun_share', centre=centre, radius=radius, point=point,
sun_direction=sun_direction, sun_radius=sun_radius)
a = float(F32(np.radians(0.2665)))
l = np.array([0.0, 0.6, 0.8])
side = np.cross(l, [1.0, 0.0, 0.0]); side /= np.linalg.norm(side)
offs = np.linspace(0.0, 4.0, 81)
centres = F32(10.0 * (np.cos(offs * a)[:, None] * l + np.sin(offs * a)[:, None] * side)).astype(float)
r = float(F32(10.0 * np.sin(1.5 * a)))
got = sun_share(centres, r, np.zeros(3), F32(l).astype(float), a)
want = reference(centres, r, np.zeros(3), F32(l).astype(float), a)
err = np.abs(got - want).max()
assert err < 2e-3, err
fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(offs, want, label='reference'); ax.plot(offs, got, '.', label='Trait, binary32')
ax.set_xlabel('the ball off the sun, in sun radii'); ax.set_ylabel('share of the sun left'); ax.legend()
print(f'across the penumbra within {err:.2e} of the reference')
across the penumbra within 5.59e-06 of the reference
Limit: a ball dead ahead¶
A ball whose disc holds the sun's leaves none of it; one half the sun's size dead ahead leaves a ring, 1 - 1/4.
one = lambda x: float(np.ravel(x)[0])
big = sun_share(10.0 * l, float(F32(10.0 * np.sin(2 * a))), np.zeros(3), l, a)
ring = sun_share(10.0 * l, float(F32(10.0 * np.sin(0.5 * a))), np.zeros(3), l, a)
assert one(big) == 0.0 and abs(one(ring) - 0.75) < 1e-3, (big, ring)
print('none, and a ring of', one(ring))
none, and a ring of 0.75
Invariant: behind and beside¶
A ball behind the point from the sun, or well clear of the sun's disc, leaves all of it, exactly.
one = lambda x: float(np.ravel(x)[0])
assert one(sun_share(-10.0 * l, 9.0, np.zeros(3), l, a)) == 1.0
clear = F32(10.0 * (np.cos(3 * a) * l + np.sin(3 * a) * side)).astype(float)
assert one(sun_share(clear, float(F32(10.0 * np.sin(a))), np.zeros(3), l, a)) == 1.0
print('behind and beside: all of the sun')
behind and beside: all of the sun