A point in a sphere¶
Whether a point, in a solid's own frame, lies in a sphere of `radius` about the frame's origin: |p|^2 <= r^2, the boundary inside, compared squared so no square root rounds it.
A CSG primitive (lusine.render.solid), the frame's placement the instance's and the step's. Source: Requicha, Representations for Rigid Solids (1980), point membership classification. Checked against a reference in numpy, which has its own tests.
- Source
- Requicha, Representations for Rigid Solids (1980), point membership classification.
- Proof
- reference, invariant, limit.
in_sphere.trait
namespace lusine.render;
bool in_sphere(vector<float<m>, 3> point, float<m> radius)
{
return dot(point, point) <= radius * radius;
}
Whether a point in a solid's own frame lies within radius of its origin, the surface inside. A primitive of the library's constructive solid geometry (lusine.render.solid), which a shape's steps place and combine and lusine.render.table's paints use to mark the ground.
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
rng = np.random.default_rng(2026)
def law(name, point, **kw):
# the kernel over many points, a few hundred a run (a step's host queue holds that many)
return np.concatenate([trait_run.law(f'render/{name}.trait', name, point=point[k:k + 500], **kw)
for k in range(0, len(point), 500)])
The kernel¶
The law as the library states it.
print(trait_run.source('render/in_sphere.trait'))
// lusine.render.in_sphere: whether a point, in a solid's own frame, lies in a sphere of `radius` about the frame's
// origin: |p|^2 <= r^2, the boundary inside, compared squared so no square root rounds it. A CSG primitive
// (lusine.render.solid), the frame's placement the instance's and the step's.
// Source: Requicha, Representations for Rigid Solids (1980), point membership classification.
// Checked against a reference in numpy, which has its own tests.
namespace lusine.render;
bool in_sphere(vector<float<m>, 3> point, float<m> radius)
{
return dot(point, point) <= radius * radius;
}
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/in_sphere.trait', 'in_sphere'))
from lusine.kernels.render.in_sphere import in_sphere as reference
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host vector<float<m>, 3> point;
host float<m> radius;
bool result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.render.in_sphere(items.point, items.radius); }
Reference: points about a ball¶
Points scattered about a boules ball of 37.5 mm: the kernel inside exactly where the reference is, apart from those within a few binary32 roundings of the surface.
p = rng.uniform(-0.05, 0.05, size=(3000, 3))
got = law('in_sphere', p, radius=0.0375)
want = reference(p, 0.0375)
near = np.abs(np.linalg.norm(p, axis=1) - 0.0375) < 1e-8
assert np.array_equal(got[~near], want[~near])
print(f'{int(got.sum())} of {len(p)} inside, the same as the reference ({int(near.sum())} on the surface left out)')
674 of 3000 inside, the same as the reference (0 on the surface left out)
Invariant: the share inside¶
The share of a cube's points inside its inscribed sphere is pi / 6.
p = rng.uniform(-1.0, 1.0, size=(6000, 3))
share = law('in_sphere', p, radius=1.0).mean()
assert abs(share - np.pi / 6) < 0.02, share
print(f'{share:.4f} of the cube inside; pi / 6 is {np.pi / 6:.4f}')
0.5245 of the cube inside; pi / 6 is 0.5236
Limit: the surface and a sphere of no radius¶
A point on the surface is inside, one just past it outside; a sphere of no radius holds its centre alone.
out = law('in_sphere', np.array([[0.0, 0.5, 0.0], [0.0, 0.50001, 0.0]]), radius=0.5)
assert list(out) == [True, False]
zero = law('in_sphere', np.array([[0.0, 0.0, 0.0], [1e-6, 0.0, 0.0]]), radius=0.0)
assert list(zero) == [True, False]
print('on the surface inside, past it outside; no radius, the centre alone')
on the surface inside, past it outside; no radius, the centre alone