A ray meeting a sphere¶
How far along a ray it meets a sphere.
The ray o + t d (d a unit vector) meets the sphere |p - c| = r where t^2 + 2 b t + (|o - c|^2 - r^2) = 0, b = (o - c) . d: t = -b -+ sqrt(D), D = b^2 - (|o - c|^2 - r^2). The nearest root past near (the ray's own surface lies nearer), else the farther root (the ray starts inside), else far: no hit, as where D <= 0 (the ray passes by, or grazes). Every square root is IEEE's. near and far are the view's (the distance below which a hit is the ray's own surface, and the distance standing for none), never this law's. Source: Shirley and Marschner, Fundamentals of Computer Graphics, 3rd ed., section 4.4.1 (ray-sphere intersection). Checked against a reference in numpy, which has its own tests.
- Source
- Shirley and Marschner, Fundamentals of Computer Graphics, 3rd ed., section 4.4.1.
- Proof
- reference, invariant, limit.
sphere_hit.trait
namespace lusine.render;
float<m> sphere_hit(vector<float<m>, 3> centre, float<m> radius, vector<float<m>, 3> origin, vector<float, 3> direction,
float<m> near, float<m> far)
{
let oc = origin - centre;
let b = dot(oc, direction);
let disc = b * b - (dot(oc, oc) - radius * radius);
if (disc <= 0) return far;
let sq = sqrt(disc);
let t0 = -b - sq;
let t1 = -b + sq;
if (t0 > near) return t0;
if (t1 > near) return t1;
return far;
}
How far along a ray it meets a sphere, the nearest root of |o + t d - c| = r past near, else the farther one (a ray from inside), else far, the view's distance standing for none. 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/sphere_hit.trait'))
// lusine.render.sphere_hit: how far along a ray it meets a sphere. The ray o + t d (d a unit vector) meets the
// sphere |p - c| = r where t^2 + 2 b t + (|o - c|^2 - r^2) = 0, b = (o - c) . d: t = -b -+ sqrt(D), D = b^2 -
// (|o - c|^2 - r^2). The nearest root past `near` (the ray's own surface lies nearer), else the farther root (the ray
// starts inside), else `far`: no hit, as where D <= 0 (the ray passes by, or grazes). Every square root is IEEE's.
// `near` and `far` are the view's (the distance below which a hit is the ray's own surface, and the distance standing
// for none), never this law's.
// Source: Shirley and Marschner, Fundamentals of Computer Graphics, 3rd ed., section 4.4.1 (ray-sphere intersection).
// Checked against a reference in numpy, which has its own tests.
namespace lusine.render;
float<m> sphere_hit(vector<float<m>, 3> centre, float<m> radius, vector<float<m>, 3> origin, vector<float, 3> direction,
float<m> near, float<m> far)
{
let oc = origin - centre;
let b = dot(oc, direction);
let disc = b * b - (dot(oc, oc) - radius * radius);
if (disc <= 0) return far;
let sq = sqrt(disc);
let t0 = -b - sq;
let t1 = -b + sq;
if (t0 > near) return t0;
if (t1 > near) return t1;
return far;
}
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/sphere_hit.trait', 'sphere_hit'))
from lusine.kernels.render.sphere_hit import sphere_hit 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> origin;
host vector<float, 3> direction;
host float<m> near;
host float<m> far;
float<m> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.render.sphere_hit(items.centre, items.radius, items.origin, items.direction, items.near, items.far); }
Reference: rays aimed near a ball¶
Rays from random places aimed near a ball of 4 cm (a boule's): the kernel's distance within binary32's roundings of the float64 reference's as the formula conditions them, 16 x 2^-24 (|o - c|^2 / sqrt(D) + t): the discriminant D is a difference of two squares near |o - c|^2, and its square root magnifies their roundings near grazing. None where the reference finds none (rays within 5% of r^2 of grazing left out, where the two may disagree whether it is met).
def hit(centre, radius, origin, direction, near=1e-4, far=1e30):
return trait_run.law('render/sphere_hit.trait', 'sphere_hit', centre=centre, radius=radius, origin=origin,
direction=direction, near=near, far=far)
c, r = np.array([1.0, 0.04, 0.0]), 0.04
o = F32(rng.uniform(-1.0, 1.0, (500, 3)) + [0.0, 0.5, 0.0]).astype(float)
aim = c + rng.normal(size=(500, 3)) * 0.04 - o
d = F32(aim / np.linalg.norm(aim, axis=1, keepdims=True)).astype(float)
got = hit(c, r, o, d)
want = reference(c, r, o, d, 1e-4, 1e30)
oc = o - c; b = np.sum(oc * d, axis=1); disc = b * b - (np.sum(oc * oc, axis=1) - r * r)
clear = np.abs(disc) > 0.05 * r * r
met = clear & (want < 1e30)
assert np.array_equal(got[clear] < 1e30, want[clear] < 1e30)
bound = 16 * U * (np.sum(oc * oc, axis=1) / np.sqrt(np.abs(disc)) + np.abs(want))
worst = np.max(np.abs(got[met] - want[met]) / bound[met])
assert worst <= 1.0, worst
print(f'{int(met.sum())} hits within {worst:.2f} of their bound, {int((clear & ~met).sum())} misses agreed')
184 hits within 0.13 of their bound, 295 misses agreed
Invariant: a hit is on the sphere¶
Every point a ray meets lies on the sphere, |o + t d - c| = r, to within the distance's own bound above.
p = o + d * np.where(got < 1e30, got, 0.0)[:, None]
on = np.linalg.norm(p[met] - c, axis=1)
assert np.all(np.abs(on - r) <= bound[met] + 4 * U), np.max(np.abs(on - r) / bound[met])
print(f'every hit within {np.max(np.abs(on - r)):.2g} m of the surface')
every hit within 1.6e-05 m of the surface
Limit: from inside, behind, and by¶
From its centre a ray leaves a sphere at its radius; a sphere behind the ray, or one it passes by, is not met (far).
z = [0.0, 0.0, 1.0]
out = hit(np.array([[0.0, 0.0, 5.0]] * 3), 1.0, np.array([[0.0, 0.0, 5.0], [0.0, 0.0, 0.0], [0.0, 2.0, 0.0]]),
np.array([z, [0.0, 0.0, -1.0], z]))
assert out[0] == 1.0 and out[1] == 1e30 * F32(1) and out[2] == F32(1e30)
print('inside: the radius; behind or by: none')
inside: the radius; behind or by: none