A ray meeting a region of the ground¶

How far along a ray it meets a region of the ground.

The ground is the plane y = 0 (y is up), and a region of it the rectangle x0 <= x <= x1, z0 <= z <= z1, or the whole plane past it where it reaches the horizon. A ray o + t d going down (d.y < -1e-12; level or rising, it never meets the plane) meets the plane at t = -o.y / d.y; the region where t is past near and the point o + t d is in the rectangle (or the region reaches the horizon), else far: no hit. near and far are the view's, as lusine.render.sphere_hit takes them. Source: Shirley and Marschner, Fundamentals of Computer Graphics, 3rd ed., section 4.4.2 (ray-plane 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.2.

Proof

reference, invariant, limit.

ground_hit.trait

namespace lusine.render;

float<m> ground_hit(vector<float<m>, 3> origin, vector<float, 3> direction, float<m> x0, float<m> z0, float<m> x1,
                    float<m> z1, bool horizon, float<m> near, float<m> far)
{
  if (direction.y >= -0.000000000001) return far;
  let t = -origin.y / direction.y;
  if (t <= near) return far;
  let x = origin.x + direction.x * t;
  let z = origin.z + direction.z * t;
  let inside = x >= x0 && x <= x1 && z >= z0 && z <= z1;
  return horizon || inside ? t : far;
}

How far along a ray it meets the plane y = 0 within a region's rectangle (or anywhere, for a region reaching the horizon), past near; else far. One of the kernels lusine.render.table draws a table's shapes with.

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

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.

In [2]:

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

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

// lusine.render.ground_hit: how far along a ray it meets a region of the ground. The ground is the plane y = 0 (y is
// up), and a region of it the rectangle x0 <= x <= x1, z0 <= z <= z1, or the whole plane past it where it reaches the
// horizon. A ray o + t d going down (d.y < -1e-12; level or rising, it never meets the plane) meets the plane at
// t = -o.y / d.y; the region where t is past `near` and the point o + t d is in the rectangle (or the region reaches
// the horizon), else `far`: no hit. `near` and `far` are the view's, as lusine.render.sphere_hit takes them.
// Source: Shirley and Marschner, Fundamentals of Computer Graphics, 3rd ed., section 4.4.2 (ray-plane intersection).
// Checked against a reference in numpy, which has its own tests.

namespace lusine.render;

float<m> ground_hit(vector<float<m>, 3> origin, vector<float, 3> direction, float<m> x0, float<m> z0, float<m> x1,
                    float<m> z1, bool horizon, float<m> near, float<m> far)
{
  if return far;
  let t = -origin.y / direction.y;
  if (t <= near) return far;
  let x = origin.x + direction.x * t;
  let z = origin.z + direction.z * t;
  let inside = x >= x0 && x <= x1 && z >= z0 && z <= z1;
  return horizon || inside ? 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/ground_hit.trait', 'ground_hit'))
from lusine.kernels.render.ground_hit import ground_hit as reference

print(trait_run.law_program('render/ground_hit.trait', 'ground_hit')) from lusine.kernels.render.ground_hit import ground_hit as reference

namespace probe;

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

struct item {
  host vector<float<m>, 3> origin;
  host vector<float, 3> direction;
  host float<m> x0;
  host float<m> z0;
  host float<m> x1;
  host float<m> z1;
  host bool horizon;
  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.ground_hit(items.origin, items.direction, items.x0, items.z0, items.x1, items.z1, items.horizon, items.near, items.far); }

Reference: rays over a pitch¶

Rays from a camera's height in every direction over a 1 m by 0.9 m pitch: the kernel's distance within binary32's roundings of the float64 reference's, and the same rays met (those within a part in ten thousand of the rectangle's edge left out).

In [4]:

def ground(origin, direction, x0=0.6, z0=-0.5, x1=1.6, z1=0.4, horizon=False, near=1e-4, far=1e30):
    return trait_run.law('render/ground_hit.trait', 'ground_hit', origin=origin, direction=direction, x0=x0,
                         z0=z0, x1=x1, z1=z1, horizon=horizon, near=near, far=far)

o = np.tile([0.2, 0.35, -0.55], (600, 1))
v = rng.normal(size=(600, 3)); v /= np.linalg.norm(v, axis=1, keepdims=True)
d = F32(v).astype(float)
got, want = ground(o, d), reference(o, d, 0.6, -0.5, 1.6, 0.4, False, 1e-4, 1e30)
t = np.where(d[:, 1] < 0, -o[:, 1] / np.where(d[:, 1] < 0, d[:, 1], -1), 0)
x, zz = o[:, 0] + d[:, 0] * t, o[:, 2] + d[:, 2] * t
edge = np.minimum.reduce([np.abs(x - 0.6), np.abs(x - 1.6), np.abs(zz + 0.5), np.abs(zz - 0.4)]) < 1e-4 * (1 + t)
keep = ~edge
assert np.array_equal(got[keep] < 1e30, want[keep] < 1e30)
met = keep & (want < 1e30)
err = np.max(np.abs(got[met] - want[met]) / want[met])
assert err <= 8 * U, err / U
print(f'{int(met.sum())} rays meet the pitch, within {err / U:.1f} x 2^-24 of the reference; the rest agree none')

def ground(origin, direction, x0=0.6, z0=-0.5, x1=1.6, z1=0.4, horizon=False, near=1e-4, far=1e30): return trait_run.law('render/ground_hit.trait', 'ground_hit', origin=origin, direction=direction, x0=x0, z0=z0, x1=x1, z1=z1, horizon=horizon, near=near, far=far) o = np.tile([0.2, 0.35, -0.55], (600, 1)) v = rng.normal(size=(600, 3)); v /= np.linalg.norm(v, axis=1, keepdims=True) d = F32(v).astype(float) got, want = ground(o, d), reference(o, d, 0.6, -0.5, 1.6, 0.4, False, 1e-4, 1e30) t = np.where(d[:, 1] < 0, -o[:, 1] / np.where(d[:, 1] < 0, d[:, 1], -1), 0) x, zz = o[:, 0] + d[:, 0] * t, o[:, 2] + d[:, 2] * t edge = np.minimum.reduce([np.abs(x - 0.6), np.abs(x - 1.6), np.abs(zz + 0.5), np.abs(zz - 0.4)]) < 1e-4 * (1 + t) keep = ~edge assert np.array_equal(got[keep] < 1e30, want[keep] < 1e30) met = keep & (want < 1e30) err = np.max(np.abs(got[met] - want[met]) / want[met]) assert err <= 8 * U, err / U print(f'{int(met.sum())} rays meet the pitch, within {err / U:.1f} x 2^-24 of the reference; the rest agree none')

22 rays meet the pitch, within 1.1 x 2^-24 of the reference; the rest agree none

Invariant: the point met is on the ground, in the rectangle¶

Every point a ray meets has y = 0 to a rounding and lies in the rectangle; a region reaching the horizon meets every ray going down, at the same distance.

In [5]:

p = o + d * np.where(got < 1e30, got, 0)[:, None]
m = got < 1e30
assert np.all(np.abs(p[m, 1]) <= 1e-6) and np.all((p[m, 0] >= 0.6 - 1e-6) & (p[m, 0] <= 1.6 + 1e-6))
far_ = ground(o, d, horizon=True)
assert np.array_equal(far_ < 1e30, d[:, 1] < -1e-12) and np.array_equal(far_[m], got[m])
print('on the ground and in the rectangle; to the horizon, every downward ray')

p = o + d * np.where(got < 1e30, got, 0)[:, None] m = got < 1e30 assert np.all(np.abs(p[m, 1]) <= 1e-6) and np.all((p[m, 0] >= 0.6 - 1e-6) & (p[m, 0] <= 1.6 + 1e-6)) far_ = ground(o, d, horizon=True) assert np.array_equal(far_ < 1e30, d[:, 1] < -1e-12) and np.array_equal(far_[m], got[m]) print('on the ground and in the rectangle; to the horizon, every downward ray')

on the ground and in the rectangle; to the horizon, every downward ray

Limit: level, rising and grazing rays¶

A level or rising ray never meets the ground, nor one starting on it.

In [6]:

out = ground(np.array([[1.0, 1.0, 0.0], [1.0, 1.0, 0.0], [1.0, 0.0, 0.0]]),
             np.array([[1.0, 0.0, 0.0], [0.0, 0.6, 0.8], [0.0, -1.0, 0.0]]), horizon=True)
assert np.all(out == F32(1e30))
print('level, rising and from the ground: none')

out = ground(np.array([[1.0, 1.0, 0.0], [1.0, 1.0, 0.0], [1.0, 0.0, 0.0]]), np.array([[1.0, 0.0, 0.0], [0.0, 0.6, 0.8], [0.0, -1.0, 0.0]]), horizon=True) assert np.all(out == F32(1e30)) print('level, rising and from the ground: none')

level, rising and from the ground: none