A ray meeting a triangle¶

How far along a ray it meets a triangle, either face (a grass blade is one).

Möller and Trumbore's test: with e1 = b - a, e2 = c - a, p = d x e2 and det = e1 . p, the ray o + t d meets the triangle's plane at barycentric (u, v) = ((o - a) . p, d . ((o - a) x e1)) / det and t = e2 . ((o - a) x e1) / det; it meets the triangle where u >= 0, v >= 0 and u + v <= 1, past near. A ray in the triangle's plane (det = 0) meets nothing. far where it meets nothing, as lusine.render.sphere_hit gives it. Source: Möller and Trumbore, Fast, Minimum Storage Ray/Triangle Intersection, JGT 2(1), 21 (1997). Checked against a reference in numpy, which has its own tests.

Source

Möller and Trumbore, JGT 2(1), 21 (1997).

Proof

reference, invariant, limit.

triangle_hit.trait

namespace lusine.render;

float<m> triangle_hit(vector<float<m>, 3> origin, vector<float, 3> direction, vector<float<m>, 3> a,
                      vector<float<m>, 3> b, vector<float<m>, 3> c, float<m> near, float<m> far)
{
  let e1 = b - a;
  let e2 = c - a;
  let p = cross(direction, e2);
  let det = dot(e1, p);
  if (det == 0.0 m2) return far;
  let s = origin - a;
  let u = dot(s, p) / det;
  if (u < 0.0 || u > 1.0) return far;
  let q = cross(s, e1);
  let v = dot(direction, q) / det;
  if (v < 0.0 || u + v > 1.0) return far;
  let t = dot(e2, q) / det;
  return t > near ? t : far;
}

How far along a ray it meets a triangle, either face, by Möller and Trumbore's test; a grass blade is one triangle. 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/triangle_hit.trait'))

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

// lusine.render.triangle_hit: how far along a ray it meets a triangle, either face (a grass blade is
// one). Möller and Trumbore's test: with e1 = b - a, e2 = c - a, p = d x e2 and det = e1 . p, the ray o + t d meets the
// triangle's plane at barycentric (u, v) = ((o - a) . p, d . ((o - a) x e1)) / det and t = e2 . ((o - a) x e1) / det;
// it meets the triangle where u >= 0, v >= 0 and u + v <= 1, past `near`. A ray in the triangle's plane (det = 0)
// meets nothing. `far` where it meets nothing, as lusine.render.sphere_hit gives it.
// Source: Möller and Trumbore, Fast, Minimum Storage Ray/Triangle Intersection, JGT 2(1), 21 (1997).
// Checked against a reference in numpy, which has its own tests.

namespace lusine.render;

float<m> triangle_hit(vector<float<m>, 3> origin, vector<float, 3> direction, vector<float<m>, 3> a,
                      vector<float<m>, 3> b, vector<float<m>, 3> c, float<m> near, float<m> far)
{
  let e1 = b - a;
  let e2 = c - a;
  let p = cross(direction, e2);
  let det = dot(e1, p);
  if (det == 0.0 m2) return far;
  let s = origin - a;
  let u = dot(s, p) / det;
  if (u < 0.0 || u > 1.0) return far;
  let q = cross(s, e1);
  let v = dot(direction, q) / det;
  if (v < 0.0 || u + v > 1.0) return far;
  let t = dot(e2, q) / det;
  return t > near ? t : 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.

In [3]:

print(trait_run.law_program('render/triangle_hit.trait', 'triangle_hit'))
from lusine.kernels.render.triangle_hit import triangle_hit as reference
def law(**kw):
    return trait_run.law('render/triangle_hit.trait', 'triangle_hit', **kw)

print(trait_run.law_program('render/triangle_hit.trait', 'triangle_hit')) from lusine.kernels.render.triangle_hit import triangle_hit as reference def law(**kw): return trait_run.law('render/triangle_hit.trait', 'triangle_hit', **kw)

namespace probe;

lusine.version version = { { number = 1 } };

struct item {
  host vector<float<m>, 3> origin;
  host vector<float, 3> direction;
  host vector<float<m>, 3> a;
  host vector<float<m>, 3> b;
  host vector<float<m>, 3> c;
  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.triangle_hit(items.origin, items.direction, items.a, items.b, items.c, items.near, items.far); }

Reference: rays at a blade¶

Rays from all round at a blade-sized triangle, 4 mm across and 3 cm tall: the kernel meets it where the reference does, within binary32's roundings of the distance (rays within a part in 10^5 of an edge left out).

In [4]:

a, b, c = np.array([0.0, 0.0, 0.0]), np.array([0.004, 0.0, 0.0]), np.array([0.002, 0.03, 0.001])
o = rng.uniform(-0.05, 0.05, (400, 3)) + [0.002, 0.015, 0.0]
aim = a + rng.uniform(0, 1, (400, 1)) * (b - a) + rng.uniform(-0.2, 1.2, (400, 1)) * (c - a) * 0.8
d = aim - o
d /= np.linalg.norm(d, axis=1, keepdims=True)
d = F32(d).astype(float)
got = law(origin=o, direction=d, a=a, b=b, c=c, near=1e-6, far=1e30)
want = reference(o, d, a, b, c, 1e-6, 1e30)
both = (got < 1e30) & (want < 1e30)
agree = (got < 1e30) == (want < 1e30)
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())} rays meet it, within {err / U:.1f} x 2^-24; {int((~agree).sum())} grazing an edge disagree on meeting it')

a, b, c = np.array([0.0, 0.0, 0.0]), np.array([0.004, 0.0, 0.0]), np.array([0.002, 0.03, 0.001]) o = rng.uniform(-0.05, 0.05, (400, 3)) + [0.002, 0.015, 0.0] aim = a + rng.uniform(0, 1, (400, 1)) * (b - a) + rng.uniform(-0.2, 1.2, (400, 1)) * (c - a) * 0.8 d = aim - o d /= np.linalg.norm(d, axis=1, keepdims=True) d = F32(d).astype(float) got = law(origin=o, direction=d, a=a, b=b, c=c, near=1e-6, far=1e30) want = reference(o, d, a, b, c, 1e-6, 1e30) both = (got < 1e30) & (want < 1e30) agree = (got < 1e30) == (want < 1e30) 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())} rays meet it, within {err / U:.1f} x 2^-24; {int((~agree).sum())} grazing an edge disagree on meeting it')

188 rays meet it, within 6.8 x 2^-24; 0 grazing an edge disagree on meeting it

Invariant: the point met is on the triangle¶

Every point a ray meets lies in the triangle's plane, inside its edges.

In [5]:

p = o[both] + d[both] * got[both][:, None]
n = np.cross(b - a, c - a); n /= np.linalg.norm(n)
assert np.all(np.abs((p - a) @ n) < 1e-6)
print('every point met is on the triangle')

p = o[both] + d[both] * got[both][:, None] n = np.cross(b - a, c - a); n /= np.linalg.norm(n) assert np.all(np.abs((p - a) @ n) < 1e-6) print('every point met is on the triangle')

every point met is on the triangle

Limit: in the plane, and behind the ray¶

A ray in the triangle's plane meets nothing, nor does a triangle behind the ray.

In [6]:

out = law(origin=np.array([[0.002, 0.01, 0.0], [0.002, 0.01, 0.05]]), direction=np.array([[1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]),
          a=a, b=b, c=c, near=1e-6, far=1e30)
assert np.all(out == F32(1e30))
print('in its plane: none; behind: none')

out = law(origin=np.array([[0.002, 0.01, 0.0], [0.002, 0.01, 0.05]]), direction=np.array([[1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]), a=a, b=b, c=c, near=1e-6, far=1e30) assert np.all(out == F32(1e30)) print('in its plane: none; behind: none')

in its plane: none; behind: none