A point in a cylinder¶
Whether a point, in a solid's own frame, lies in the cylinder of `radius` about the frame's y axis, from -half_height to half_height along it: p.x^2 + p.z^2 <= r^2 and |p.y| <= h, the surface inside, the radius compared squared.
A CSG primitive (lusine.render.solid); a ring on the ground is a cylinder less a narrower one, a disc a cylinder alone. 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_cylinder.trait
namespace lusine.render;
bool in_cylinder(vector<float<m>, 3> point, float<m> radius, float<m> half_height)
{
return point.x * point.x + point.z * point.z <= radius * radius && abs(point.y) <= half_height;
}
Whether a point in a solid's own frame lies within radius of its y axis and half_height of its origin along it, 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_cylinder.trait'))
// lusine.render.in_cylinder: whether a point, in a solid's own frame, lies in the cylinder of `radius` about the
// frame's y axis, from -half_height to half_height along it: p.x^2 + p.z^2 <= r^2 and |p.y| <= h, the surface inside,
// the radius compared squared. A CSG primitive (lusine.render.solid); a ring on the ground is a
// cylinder less a narrower one, a disc a cylinder alone.
// 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_cylinder(vector<float<m>, 3> point, float<m> radius, float<m> half_height)
{
return point.x * point.x + point.z * point.z <= radius * radius && abs(point.y) <= half_height;
}
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_cylinder.trait', 'in_cylinder'))
from lusine.kernels.render.in_cylinder import in_cylinder as reference
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host vector<float<m>, 3> point;
host float<m> radius;
host float<m> half_height;
bool result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.render.in_cylinder(items.point, items.radius, items.half_height); }
Reference: a throwing circle¶
Points about a boules throwing circle's hoop, 0.27 m in radius and a centimetre high: the kernel inside exactly where the reference is, apart from those within a few binary32 roundings of the surface.
p = rng.uniform(-1.0, 1.0, size=(3000, 3)) * np.array([0.3, 0.008, 0.3])
got = law('in_cylinder', p, radius=0.27, half_height=0.005)
want = reference(p, 0.27, 0.005)
near = (np.abs(np.hypot(p[:, 0], p[:, 2]) - 0.27) < 1e-6) | (np.abs(np.abs(p[:, 1]) - 0.005) < 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)')
1230 of 3000 inside, the same as the reference (0 on the surface left out)
Invariant: the share inside, and a turn about the axis¶
The share of a cube's points inside the cylinder is its share of the cube's volume, and a quarter turn about the axis leaves each point's answer as it was.
p = rng.uniform(-1.0, 1.0, size=(6000, 3))
inside = law('in_cylinder', p, radius=1.0, half_height=0.5)
turned = law('in_cylinder', p[:, [2, 1, 0]] * np.array([-1, 1, 1]), radius=1.0, half_height=0.5)
assert abs(inside.mean() - np.pi / 8) < 0.02 and np.array_equal(inside, turned)
print(f'{inside.mean():.4f} of the cube inside, pi / 8 = {np.pi / 8:.4f}; a quarter turn about y changes nothing')
0.4050 of the cube inside, pi / 8 = 0.3927; a quarter turn about y changes nothing
Limit: a disc and a ring¶
A disc on the ground is a cylinder alone; a ring is one less a narrower one, so the circle's centre is outside the ring and its rim inside.
pts = np.array([[0.0, 0.0, 0.0], [0.26, 0.0, 0.0], [0.0, 0.0, -0.26]])
outer = law('in_cylinder', pts, radius=0.27, half_height=0.005)
inner = law('in_cylinder', pts, radius=0.25, half_height=0.01)
assert outer.all() and list(outer & ~inner) == [False, True, True]
print('the disc holds all three; the ring, the cylinder less the narrower, holds the rim and not the centre')
the disc holds all three; the ring, the cylinder less the narrower, holds the rim and not the centre