isoray 0.3.0

Robust iso-surface rendering via interval-arithmetic ray casting for arbitrary implicit geometries

isoray

Robust iso-surface rendering via interval-arithmetic ray casting for arbitrary implicit geometries.

isoray renders the zero-level-set of any implicit function f: R^3 -> R using interval-arithmetic bisection. The primary target is black-box functions where only point evaluations are available — no symbolic formula access required. When the analytic formula is known, Natural Interval Arithmetic (NIA) and Affine Arithmetic (AA) provide tighter bounds as an optimization path.

Key Features

• Black-box compatible — works with any function f(point) -> scalar. No symbolic access needed.
• Guaranteed robustness — interval arithmetic ensures no surface features are missed, even for thin shells, cusps, high-frequency surfaces, and near-tangent intersections.
• Four inclusion strategies — GradientInclusion (default, uses gradient for tighter bounds), SamplingInclusion (black-box baseline), NaturalIntervalInclusion (analytic), AffineInclusion (analytic, tightest bounds).
• Newton refinement — optional post-bisection refinement for sub-pixel hit-point accuracy.
• Mesh SDF support — render triangle meshes via signed distance fields (MeshSDF).
• CSG operations — Union, Intersection, Difference compositing of implicit functions.
• Batch-vectorized — all evaluations operate on (N, 3) arrays, no per-point loops.
• Differentiable depth rendering — PyTorch integration with exact gradients via implicit function theorem (IFT). No need to differentiate through the bisection.

Installation

Requires Python >= 3.11.

# Using uv (recommended)
uv add isoray

# Or using pip
pip install isoray

Dependencies

• numpy >= 1.24.0
• pillow >= 11.0.0

Optional extras:

pip install isoray[mesh]    # Mesh SDF support (point-cloud-utils)
pip install isoray[torch]   # or: uv add isoray[torch]

Quick Start

import numpy as np
from isoray import render

# Define any implicit function: f(pos) -> (values, gradients_or_None)
def my_sphere(pos):
    return np.sum(pos**2, axis=1) - 1.0, None

# Render with one call — just provide the function and bounding box
result = render(my_sphere, bounds=([-1.5]*3, [1.5]*3), save_to="sphere.png")

Camera direction and distance are auto-calculated. Override as needed:

# Custom camera direction, explicit distance, high resolution
result = render(my_sphere, bounds=([-1.5]*3, [1.5]*3),
                direction=(0, 0, 1), distance=5.0,
                width=1024, height=1024, save_to="sphere_hd.png")

# Depth map instead of shaded image
result = render(my_sphere, bounds=([-1.5]*3, [1.5]*3),
                mode="depth", save_to="sphere_depth.png")

# Orthographic (parallel) projection
result = render(my_sphere, bounds=([-1.5]*3, [1.5]*3),
                projection="parallel", save_to="sphere_ortho.png")

Using a built-in geometry:

from isoray import Sphere, render

sphere = Sphere(radius=1.0)
result = render(sphere, bounds=sphere.bounds(), save_to="sphere.png")

Run the built-in example:

uv run python examples/render_simple.py

Low-Level Pipeline

For full control, use the six-step pipeline directly:

from isoray import (
    SamplingInclusion, generate_rays, ray_aabb_intersect,
    bisect_trace, newton_refine, shade_lambertian, save_image,
)

origins, directions = generate_rays(256, 256, fov_deg=60.0, eye=[0, 0, 3])
t_min, t_max, hit = ray_aabb_intersect(origins, directions, [-1.5]*3, [1.5]*3)
inclusion = SamplingInclusion(my_sphere, lipschitz=2.5)
result = bisect_trace(origins, directions, t_min, t_max, hit,
                      inclusion, my_sphere, max_depth=40)
result.t = newton_refine(result.t, origins, directions, result.hit, my_sphere,
                         t_lo=result.t_lo, t_hi=result.t_hi)
result.position = origins + result.t[:, None] * directions
# ... compute normals, shade, save

Rendering Pipeline

generate_rays -> ray_aabb_intersect -> bisect_trace -> newton_refine -> shade_lambertian -> save_image

1. Ray generation — Pinhole camera model emits rays through each pixel.
2. AABB intersection — Slab test clips rays to the bounding box, discarding misses early.
3. Bisection trace — Level-synchronized interval bisection walks along each ray. The inclusion function bounds f over axis-aligned boxes; intervals that cannot contain zero are pruned.
4. Newton refinement — Optional iterative refinement of hit points using directional derivatives.
5. Shading — Lambertian shading from surface normals (analytic gradient or centered finite difference).
6. Output — PNG images for both shaded color and depth maps.

Architecture

Core Protocols

Protocol	Module	Signature
ImplicitFunction	geometry/interface.py	__call__(pos: (N,3)) -> (val: (N,), grad: (N,3)|None) + bounds()
InclusionFunction	arithmetic/interface.py	evaluate(lo: (N,3), hi: (N,3)) -> (f_lo: (N,), f_hi: (N,))

Inclusion Strategies

Strategy	Module	Use Case
GradientInclusion	arithmetic/interval.py	Default in render(). Tightens SamplingInclusion via gradient sign test + Mean Value Form. Falls back to SamplingInclusion when gradient is unavailable
SamplingInclusion	arithmetic/interval.py	Black-box baseline: samples 8 corners + center, adds Lipschitz margin
NaturalIntervalInclusion	arithmetic/interval.py	Analytic: evaluates formula with interval arithmetic ops
AffineInclusion	arithmetic/affine.py	Analytic: tightest bounds via noise-symbol correlation tracking

Module Layout

isoray/
  geometry/       # ImplicitFunction implementations
    interface.py  #   Protocol definition
    primitives.py #   Sphere, Torus, ThinShell, HighFrequencySurface, WhitneyUmbrella, CuspSurface
    csg.py        #   Union, Intersection, Difference
    mesh_sdf.py   #   MeshSDF (triangle mesh via point_cloud_utils)
  arithmetic/     # Interval & affine arithmetic
    interface.py  #   InclusionFunction protocol
    interval.py   #   ia_* ops, NIA recipes, SamplingInclusion, GradientInclusion, NaturalIntervalInclusion
    affine.py     #   AffineForm dataclass, aa_* ops, AA recipes, AffineInclusion
  tracer/         # Ray casting core
    bisection.py  #   Level-synchronized bisection (the core algorithm)
    refine.py     #   Newton refinement
    ray.py        #   Ray-AABB slab intersection
  renderer/       # Image output
    render.py     #   High-level render() API
    camera.py     #   Pinhole + orthographic ray generation
    shader.py     #   Lambertian shading
    image.py      #   PNG output (color + depth map)
  torch/          # PyTorch differentiable rendering (optional)
    render.py     #   render_depth() — IFT-based differentiable depth
  utils.py        # centered_difference_normals (6N-point finite diff)

Differentiable Rendering (PyTorch)

isoray.torch provides differentiable depth rendering for shape optimization, inverse problems, and learning-based geometry. The forward pass uses the robust NumPy bisection pipeline; the backward pass computes exact depth gradients via the implicit function theorem (IFT) — no need to differentiate through the bisection.

import torch
from isoray.torch import render_depth

# Parameters you want to optimize
radius = torch.tensor(1.0, requires_grad=True)

# Define your implicit function in PyTorch (can be arbitrarily complex)
def my_sdf(pos):
    return (pos ** 2).sum(dim=1) - radius ** 2

result = render_depth(my_sdf, bounds=([-1.5]*3, [1.5]*3))

# Depth is differentiable w.r.t. radius (Lipschitz constant auto-estimated)
loss = result.depth[result.hit].mean()
loss.backward()
print(radius.grad)  # d(mean_depth) / d(radius)

The function f can use any PyTorch operations — neural networks, complex math, etc. As long as f is differentiable w.r.t. its parameters, gradients will flow through the depth map.

Built-in Geometries

Geometry	Description
Sphere	|p - c|^2 - r^2
Torus	(sqrt(x^2+y^2) - R)^2 + z^2 - r^2
ThinShell	|x^2+y^2+z^2 - r^2| - epsilon
HighFrequencySurface	sin(wx)sin(wy)sin(wz) - c
WhitneyUmbrella	x^2 - y^2*z (self-intersecting)
CuspSurface	x^3 - y^2 (degenerate gradient)
MeshSDF	Signed distance field from triangle mesh (STL/OBJ)
Union / Intersection / Difference	CSG compositing

Adding Custom Geometry

Black-box function (no formula access)

import numpy as np
from isoray import SamplingInclusion, bisect_trace

class MyFunction:
    def __call__(self, pos: np.ndarray) -> tuple[np.ndarray, None]:
        # pos: (N, 3) -> val: (N,), grad: None (unknown)
        x, y, z = pos[:, 0], pos[:, 1], pos[:, 2]
        val = x**2 + y**2 + z**2 - 1.0  # your function here
        return val, None  # gradient not available

    def bounds(self) -> tuple[np.ndarray, np.ndarray]:
        return np.array([-2.0, -2.0, -2.0]), np.array([2.0, 2.0, 2.0])

f = MyFunction()
inclusion = SamplingInclusion(f, lipschitz=3.0)  # or use estimate_lipschitz(f, f.bounds())

Analytic function (NIA/AA recipes)

For tighter bounds, provide interval arithmetic or affine arithmetic recipes:

from isoray import NaturalIntervalInclusion, AffineInclusion
from isoray.arithmetic.interval import ia_mul, ia_add, ia_sub

# NIA recipe: pure function using ia_* interval ops
def my_nia(x_lo, x_hi, y_lo, y_hi, z_lo, z_hi, params):
    xx = ia_mul(x_lo, x_hi, x_lo, x_hi)
    yy = ia_mul(y_lo, y_hi, y_lo, y_hi)
    zz = ia_mul(z_lo, z_hi, z_lo, z_hi)
    sum_lo, sum_hi = ia_add(*ia_add(xx, yy), zz)
    return ia_sub(sum_lo, sum_hi, np.float64(1.0), np.float64(1.0))

inclusion = NaturalIntervalInclusion(my_nia, params=None)

Development

# Run all tests (201 tests)
uv run pytest

# Run specific test module
uv run pytest tests/test_tracer/test_bisection.py

# Lint and format
uv run ruff check isoray/ tests/
uv run ruff format isoray/ tests/

# Type check
uv run ty check

# Render all robustness test cases
uv run python examples/render_robustness.py

Test Organization

Directory	Coverage
tests/test_arithmetic/	Interval and affine arithmetic operations, inclusion properties
tests/test_geometry/	Primitive function correctness
tests/test_tracer/	Bisection convergence, Newton refinement accuracy
tests/test_robustness/	Challenging cases: thin shell, cusp, grazing ray, near-tangent, high-frequency, self-intersecting, mesh SDF
tests/test_renderer/	End-to-end rendering pipeline
tests/test_torch/	Differentiable depth rendering: forward correctness, IFT gradient vs finite difference, edge cases
tests/test_utils/	Lipschitz estimation (analytic gradient, centered difference, auto mode)

How It Works

The core algorithm is level-synchronized interval bisection. For each ray:

1. The ray's parameter range [t_min, t_max] (from the AABB intersection) is represented as an interval.
2. At each bisection step, the inclusion function computes guaranteed bounds [f_lo, f_hi] on f over the axis-aligned box swept by the ray segment.
3. If f_lo > 0 or f_hi < 0 (the interval doesn't contain zero), the segment is pruned — no surface can exist there.
4. Otherwise, the segment is split and both halves are recursively examined.
5. When the interval width falls below a threshold, a surface hit is recorded.
6. Optional Newton refinement then converges the hit point to machine precision.

This approach guarantees no surface features are missed, unlike fixed-step ray marching which can skip thin features or high-frequency detail.

License

MIT License. See LICENSE for details.