The nearest grass blade along a ray¶

The nearest grass blade along a ray through the lawn's layer.

The ray enters the layer of blades (from y = top down to the ground) and crosses its cells in order, Amanatides and Woo's grid traversal: from the cell it enters by, it steps each time into the cell across the nearer of the two cell boundaries ahead of it, and in each cell it crosses meets that cell's blade (lusine.render.blade, render. triangle_hit). The walk is an ordered fold over a fixed table of steps, 48 cells, its state the cell, the distances to the next boundaries, and the nearest blade met so far and its cell; past stop (where the ray meets the ground, or the walk's own limit) a step does nothing. A blade leaning out of its own cell is met only from the cells the ray crosses into, so a ray clipping a neighbour's tip misses it; the reference does the same. Source: Amanatides and Woo, A Fast Voxel Traversal Algorithm for Ray Tracing, Eurographics 1987. Checked against a reference in numpy, which has its own tests.

Source

Amanatides and Woo, Eurographics 1987.

Proof

reference, invariant, limit.

blade_walk.trait

namespace lusine.render;

// the walk's steps: one a cell crossed. A stored table, not a const one: a fold over a const table is unrolled, and
// unrolled the 48 steps in each of lusine.render.table's rays made its render's shader 10 MB of GLSL (4.7 MB looped)
struct blade_step { integer index; };
table<blade_step> blade_steps(capacity = 48, key = index) = {
  { index = 0 },
  { index = 1 },
  { index = 2 },
  { index = 3 },
  { index = 4 },
  { index = 5 },
  { index = 6 },
  { index = 7 },
  { index = 8 },
  { index = 9 },
  { index = 10 },
  { index = 11 },
  { index = 12 },
  { index = 13 },
  { index = 14 },
  { index = 15 },
  { index = 16 },
  { index = 17 },
  { index = 18 },
  { index = 19 },
  { index = 20 },
  { index = 21 },
  { index = 22 },
  { index = 23 },
  { index = 24 },
  { index = 25 },
  { index = 26 },
  { index = 27 },
  { index = 28 },
  { index = 29 },
  { index = 30 },
  { index = 31 },
  { index = 32 },
  { index = 33 },
  { index = 34 },
  { index = 35 },
  { index = 36 },
  { index = 37 },
  { index = 38 },
  { index = 39 },
  { index = 40 },
  { index = 41 },
  { index = 42 },
  { index = 43 },
  { index = 44 },
  { index = 45 },
  { index = 46 },
  { index = 47 },
};

// the walk so far: the cell it is in and the distance it entered it at, the distances to the next x and z boundaries,
// and the nearest blade met and its cell (best = far: none)
struct blade_walk_state { integer ix; integer iz; float<m> at; float<m> next_x; float<m> next_z; float<m> best;
                          integer bx; integer bz; };

// one step of the walk (`step`, its place in blade_steps): the cell's blade met, then into the next cell
blade_walk_state blade_walk_step(blade_walk_state w, vector<float<m>, 3> origin, vector<float, 3> direction,
                                 float<m> spacing, float<m> height, float spread, float<m> width, float<rad> lean,
                                 float<rad> lean_spread, float<rad> grain, float grain_share, integer sx, integer sz,
                                 float<m> dx, float<m> dz, float<m> stop, float<m> near, integer step)
{
  if (w.at > stop) return w;
  let shape = blade(ix = w.ix, iz = w.iz, spacing = spacing, height = height, spread = spread, width = width,
                    lean = lean, lean_spread = lean_spread, grain = grain, grain_share = grain_share);
  let t = triangle_hit(origin = origin, direction = direction, a = shape.a, b = shape.b, c = shape.tip, near = near,
                       far = w.best);
  let closer = t < w.best;
  let across_x = w.next_x < w.next_z;
  return blade_walk_state { ix = across_x ? w.ix + sx : w.ix, iz = across_x ? w.iz : w.iz + sz,
                            at = across_x ? w.next_x : w.next_z,
                            next_x = across_x ? w.next_x + dx : w.next_x, next_z = across_x ? w.next_z : w.next_z + dz,
                            best = closer ? t : w.best, bx = closer ? w.ix : w.bx, bz = closer ? w.iz : w.bz };
}

// the walk: where the ray enters the layer (at y = top, or `near` where it starts inside), its cell, and the 48 steps;
// the nearest blade's distance (far: none) and its cell
blade_walk_state blade_walk(vector<float<m>, 3> origin, vector<float, 3> direction, float<m> spacing, float<m> height,
                            float spread, float<m> width, float<rad> lean, float<rad> lean_spread, float<rad> grain,
                            float grain_share, float<m> stop, float<m> near, float<m> far)
{
  let top = height * (1.0 + spread);
  let inside = origin.y <= top;
  let down = direction.y < 0.0;
  let enter = inside ? near : (down ? (origin.y - top) / -direction.y : far);
  let start = enter < near ? near : enter;
  let p = origin + direction * start;
  let ix = floor(p.x / spacing);
  let iz = floor(p.z / spacing);
  let sx = direction.x >= 0.0 ? 1 : -1;
  let sz = direction.z >= 0.0 ? 1 : -1;
  let next_x = direction.x > 0.0 ? ((ix + 1) * spacing - origin.x) / direction.x
      : (direction.x < 0.0 ? (ix * spacing - origin.x) / direction.x : far);
  let next_z = direction.z > 0.0 ? ((iz + 1) * spacing - origin.z) / direction.z
      : (direction.z < 0.0 ? (iz * spacing - origin.z) / direction.z : far);
  let dx = direction.x != 0.0 ? spacing / abs(direction.x) : far;
  let dz = direction.z != 0.0 ? spacing / abs(direction.z) : far;
  let first = blade_walk_state { ix = ix, iz = iz, at = start, next_x = next_x, next_z = next_z, best = far,
                                 bx = ix, bz = iz };
  return fold(w = first) blade_walk_step(w = w, origin = origin, direction = direction, spacing = spacing,
                                         height = height, spread = spread, width = width, lean = lean,
                                         lean_spread = lean_spread, grain = grain, grain_share = grain_share, sx = sx,
                                         sz = sz, dx = dx, dz = dz, stop = start > stop ? start - 1.0 m : stop,
                                         near = near, step = blade_steps.index);
}

The nearest blade a ray meets in the lawn's layer, by stepping cell by cell (Amanatides and Woo's traversal) as an ordered fold over 48 steps, meeting each cell's blade. Part of the table view's grass.

In [1]:

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
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)

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 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.

In [2]:

print(trait_run.source('render/blade_walk.trait'))

print(trait_run.source('render/blade_walk.trait'))

// lusine.render.blade_walk: the nearest grass blade along a ray through the lawn's layer. The ray
// enters the layer of blades (from y = top down to the ground) and crosses its cells in order, Amanatides and Woo's
// grid traversal: from the cell it enters by, it steps each time into the cell across the nearer of the two cell
// boundaries ahead of it, and in each cell it crosses meets that cell's blade (lusine.render.blade, render.
// triangle_hit). The walk is an ordered fold over a fixed table of steps, 48 cells, its state the cell, the distances to
// the next boundaries, and the nearest blade met so far and its cell; past `stop` (where the ray meets the ground, or
// the walk's own limit) a step does nothing. A blade leaning out of its own cell is met only from the cells the ray
// crosses into, so a ray clipping a neighbour's tip misses it; the reference does the same.
// Source: Amanatides and Woo, A Fast Voxel Traversal Algorithm for Ray Tracing, Eurographics 1987.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.render;

// the walk's steps: one a cell crossed. A stored table, not a const one: a fold over a const table is unrolled, and
// unrolled the 48 steps in each of lusine.render.table's rays made its render's shader 10 MB of GLSL (4.7 MB looped)
struct blade_step { integer index; };
table<blade_step> blade_steps(capacity = 48, key = index) = {
  { index = 0 },
  { index = 1 },
  { index = 2 },
  { index = 3 },
  { index = 4 },
  { index = 5 },
  { index = 6 },
  { index = 7 },
  { index = 8 },
  { index = 9 },
  { index = 10 },
  { index = 11 },
  { index = 12 },
  { index = 13 },
  { index = 14 },
  { index = 15 },
  { index = 16 },
  { index = 17 },
  { index = 18 },
  { index = 19 },
  { index = 20 },
  { index = 21 },
  { index = 22 },
  { index = 23 },
  { index = 24 },
  { index = 25 },
  { index = 26 },
  { index = 27 },
  { index = 28 },
  { index = 29 },
  { index = 30 },
  { index = 31 },
  { index = 32 },
  { index = 33 },
  { index = 34 },
  { index = 35 },
  { index = 36 },
  { index = 37 },
  { index = 38 },
  { index = 39 },
  { index = 40 },
  { index = 41 },
  { index = 42 },
  { index = 43 },
  { index = 44 },
  { index = 45 },
  { index = 46 },
  { index = 47 },
};

// the walk so far: the cell it is in and the distance it entered it at, the distances to the next x and z boundaries,
// and the nearest blade met and its cell (best = far: none)
struct blade_walk_state { integer ix; integer iz; float<m> at; float<m> next_x; float<m> next_z; float<m> best;
                          integer bx; integer bz; };

// one step of the walk (`step`, its place in blade_steps): the cell's blade met, then into the next cell
blade_walk_state blade_walk_step(blade_walk_state w, vector<float<m>, 3> origin, vector<float, 3> direction,
                                 float<m> spacing, float<m> height, float spread, float<m> width, float<rad> lean,
                                 float<rad> lean_spread, float<rad> grain, float grain_share, integer sx, integer sz,
                                 float<m> dx, float<m> dz, float<m> stop, float<m> near, integer step)
{
  if (w.at > stop) return w;
  let shape = blade(ix = w.ix, iz = w.iz, spacing = spacing, height = height, spread = spread, width = width,
                    lean = lean, lean_spread = lean_spread, grain = grain, grain_share = grain_share);
  let t = triangle_hit(origin = origin, direction = direction, a = shape.a, b = shape.b, c = shape.tip, near = near,
                       far = w.best);
  let closer = t < w.best;
  let across_x = w.next_x < w.next_z;
  return blade_walk_state { ix = across_x ? w.ix + sx : w.ix, iz = across_x ? w.iz : w.iz + sz,
                            at = across_x ? w.next_x : w.next_z,
                            next_x = across_x ? w.next_x + dx : w.next_x, next_z = across_x ? w.next_z : w.next_z + dz,
                            best = closer ? t : w.best, bx = closer ? w.ix : w.bx, bz = closer ? w.iz : w.bz };
}

// the walk: where the ray enters the layer (at y = top, or `near` where it starts inside), its cell, and the 48 steps;
// the nearest blade's distance (far: none) and its cell
blade_walk_state blade_walk(vector<float<m>, 3> origin, vector<float, 3> direction, float<m> spacing, float<m> height,
                            float spread, float<m> width, float<rad> lean, float<rad> lean_spread, float<rad> grain,
                            float grain_share, float<m> stop, float<m> near, float<m> far)
{
  let top = height * (1.0 + spread);
  let inside = origin.y <= top;
  let down = direction.y < 0.0;
  let enter = inside ? near : (down ? (origin.y - top) / -direction.y : far);
  let start = enter < near ? near : enter;
  let p = origin + direction * start;
  let ix = floor(p.x / spacing);
  let iz = floor(p.z / spacing);
  let sx = direction.x >= 0.0 ? 1 : -1;
  let sz = direction.z >= 0.0 ? 1 : -1;
  let next_x = direction.x > 0.0 ? ((ix + 1) * spacing - origin.x) / direction.x
      : (direction.x < 0.0 ? (ix * spacing - origin.x) / direction.x : far);
  let next_z = direction.z > 0.0 ? ((iz + 1) * spacing - origin.z) / direction.z
      : (direction.z < 0.0 ? (iz * spacing - origin.z) / direction.z : far);
  let dx = direction.x != 0.0 ? spacing / abs(direction.x) : far;
  let dz = direction.z != 0.0 ? spacing / abs(direction.z) : far;
  let first = blade_walk_state { ix = ix, iz = iz, at = start, next_x = next_x, next_z = next_z, best = far,
                                 bx = ix, bz = iz };
  return fold(w = first) blade_walk_step(w = w, origin = origin, direction = direction, spacing = spacing,
                                         height = height, spread = spread, width = width, lean = lean,
                                         lean_spread = lean_spread, grain = grain, grain_share = grain_share, sx = sx,
                                         sz = sz, dx = dx, dz = dz, stop = start > stop ? start - 1.0 m : stop,
                                         near = near, step = blade_steps.index);
}

The program¶

The kernel returns the walk's state, a struct; the program binds a small kernel reading its distance over a lattice of rays.

In [3]:

from lusine.kernels.render.blade_walk import blade_walk as reference, GRASS
TEXT = '''namespace probe;
lusine.version version = { { number = 1 } };
struct item { host vector<float<m>, 3> o; host vector<float, 3> d; host float<m> stop; float<m> best; };
lattice<item, 1> items(spacing = 1 m, boundary = wall);
float<m> best_of(vector<float<m>, 3> o, vector<float, 3> d, float<m> stop) {
  return lusine.render.blade_walk(origin = o, direction = d, spacing = 0.014 m, height = 0.03 m, spread = 0.3,
      width = 0.004 m, lean = 0.26 rad, lean_spread = 0.17 rad, grain = 0.5 rad, grain_share = 0.5,
      stop = stop, near = 0.0001 m, far = 1.0e30 m).best;
}
void main() { items.best = best_of(items.o, items.d, items.stop); }
'''
print(TEXT)
def walk(o, d, stop):
    out = []
    for k in range(0, len(o), 200):
        run = trait_run.Program([], TEXT, {'probe.items': [len(o[k:k + 200])]})
        run.store('probe.items', 'o', o[k:k + 200]); run.store('probe.items', 'd', d[k:k + 200])
        run.store('probe.items', 'stop', stop[k:k + 200])
        out.append(np.asarray(run.steps('probe.main', 1, [('probe.items', 'best')])[0][0]))
    return np.concatenate(out)

from lusine.kernels.render.blade_walk import blade_walk as reference, GRASS TEXT = '''namespace probe; lusine.version version = { { number = 1 } }; struct item { host vector<float<m>, 3> o; host vector<float, 3> d; host float<m> stop; float<m> best; }; lattice<item, 1> items(spacing = 1 m, boundary = wall); float<m> best_of(vector<float<m>, 3> o, vector<float, 3> d, float<m> stop) { return lusine.render.blade_walk(origin = o, direction = d, spacing = 0.014 m, height = 0.03 m, spread = 0.3, width = 0.004 m, lean = 0.26 rad, lean_spread = 0.17 rad, grain = 0.5 rad, grain_share = 0.5, stop = stop, near = 0.0001 m, far = 1.0e30 m).best; } void main() { items.best = best_of(items.o, items.d, items.stop); } ''' print(TEXT) def walk(o, d, stop): out = [] for k in range(0, len(o), 200): run = trait_run.Program([], TEXT, {'probe.items': [len(o[k:k + 200])]}) run.store('probe.items', 'o', o[k:k + 200]); run.store('probe.items', 'd', d[k:k + 200]) run.store('probe.items', 'stop', stop[k:k + 200]) out.append(np.asarray(run.steps('probe.main', 1, [('probe.items', 'best')])[0][0])) return np.concatenate(out)

namespace probe;
lusine.version version = { { number = 1 } };
struct item { host vector<float<m>, 3> o; host vector<float, 3> d; host float<m> stop; float<m> best; };
lattice<item, 1> items(spacing = 1 m, boundary = wall);
float<m> best_of(vector<float<m>, 3> o, vector<float, 3> d, float<m> stop) {
  return lusine.render.blade_walk(origin = o, direction = d, spacing = 0.014 m, height = 0.03 m, spread = 0.3,
      width = 0.004 m, lean = 0.26 rad, lean_spread = 0.17 rad, grain = 0.5 rad, grain_share = 0.5,
      stop = stop, near = 0.0001 m, far = 1.0e30 m).best;
}
void main() { items.best = best_of(items.o, items.d, items.stop); }

Reference: rays down into a lawn¶

Rays from a few centimetres to 20 cm above the lawn, down at every angle to the ground: the kernel meets a blade where the reference does, at the same distance within binary32's roundings (a ray grazing a blade's edge may disagree).

In [4]:

n = 600
o = np.tile([0.3, 0.0, 0.1], (n, 1)) + rng.uniform(-0.05, 0.05, (n, 3)) * [1, 0, 1] + np.stack([np.zeros(n), rng.uniform(0.02, 0.2, n), np.zeros(n)], 1)
d = np.stack([rng.uniform(-1, 1, n), -rng.uniform(0.2, 1, n), rng.uniform(-1, 1, n)], 1)
d /= np.linalg.norm(d, axis=1, keepdims=True)
d = F32(d).astype(float)
stop = o[:, 1] / -d[:, 1]
got = walk(o, d, stop)
want = reference(o, d, stop=stop, near=1e-4, far=1e30, **GRASS)[0]
agree = (got < 1e29) == (want < 1e29)
both = (got < 1e29) & (want < 1e29)
assert agree.mean() >= 0.99, agree.mean()
err = np.max(np.abs(got[both] - want[both]) / want[both])
assert err <= 64 * U, err / U
print(f'{int(both.sum())} of {n} rays meet a blade, within {err / U:.1f} x 2^-24; {int((~agree).sum())} grazing an edge disagree')

n = 600 o = np.tile([0.3, 0.0, 0.1], (n, 1)) + rng.uniform(-0.05, 0.05, (n, 3)) * [1, 0, 1] + np.stack([np.zeros(n), rng.uniform(0.02, 0.2, n), np.zeros(n)], 1) d = np.stack([rng.uniform(-1, 1, n), -rng.uniform(0.2, 1, n), rng.uniform(-1, 1, n)], 1) d /= np.linalg.norm(d, axis=1, keepdims=True) d = F32(d).astype(float) stop = o[:, 1] / -d[:, 1] got = walk(o, d, stop) want = reference(o, d, stop=stop, near=1e-4, far=1e30, **GRASS)[0] agree = (got < 1e29) == (want < 1e29) both = (got < 1e29) & (want < 1e29) assert agree.mean() >= 0.99, agree.mean() err = np.max(np.abs(got[both] - want[both]) / want[both]) assert err <= 64 * U, err / U print(f'{int(both.sum())} of {n} rays meet a blade, within {err / U:.1f} x 2^-24; {int((~agree).sum())} grazing an edge disagree')

138 of 600 rays meet a blade, within 27.3 x 2^-24; 0 grazing an edge disagree

Invariant: before the ground, within the layer¶

Every blade met is nearer than the ground and within the layer's height.

In [5]:

p = o[both] + d[both] * got[both][:, None]
assert np.all(got[both] <= stop[both] * (1 + 1e-6)) and np.all((p[:, 1] >= -1e-6) & (p[:, 1] <= 0.039 + 1e-6))
print('every blade met is before the ground, inside the layer')

p = o[both] + d[both] * got[both][:, None] assert np.all(got[both] <= stop[both] * (1 + 1e-6)) and np.all((p[:, 1] >= -1e-6) & (p[:, 1] <= 0.039 + 1e-6)) print('every blade met is before the ground, inside the layer')

every blade met is before the ground, inside the layer

Limit: stopped, and level above¶

A walk stopped before the layer meets nothing, nor does a level ray above it.

In [6]:

none = walk(o[:4], d[:4], np.zeros(4))
level = walk(np.array([[0.3, 0.1, 0.1]]), np.array([[1.0, 0.0, 0.0]]), np.array([1e30]))
assert np.all(none >= F32(1e30)) and level[0] >= F32(1e30)
print('stopped: none; level above the blades: none')

none = walk(o[:4], d[:4], np.zeros(4)) level = walk(np.array([[0.3, 0.1, 0.1]]), np.array([[1.0, 0.0, 0.0]]), np.array([1e30])) assert np.all(none >= F32(1e30)) and level[0] >= F32(1e30) print('stopped: none; level above the blades: none')

stopped: none; level above the blades: none