Metadata-Version: 2.4
Name: isoray
Version: 0.2.0
Summary: Robust iso-surface rendering via interval-arithmetic ray casting for arbitrary implicit geometries
Project-URL: Repository, https://github.com/narnia-ai-mason/iso-surface-rendering
Author: Mason Seo
License-Expression: MIT
License-File: LICENSE
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: Visualization
Requires-Python: >=3.11
Requires-Dist: numpy>=1.24.0
Requires-Dist: pillow>=11.0.0
Provides-Extra: mesh
Requires-Dist: point-cloud-utils>=0.34.0; extra == 'mesh'
Provides-Extra: torch
Requires-Dist: torch>=2.0.0; extra == 'torch'
Description-Content-Type: text/markdown

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

```bash
# Using uv (recommended)
uv pip install isoray

# Or using pip
pip install isoray
```

### Dependencies

- `numpy >= 1.24.0`
- `pillow >= 11.0.0`

**Optional extras:**

```bash
pip install isoray[mesh]    # Mesh SDF support (point-cloud-utils)
pip install isoray[torch]   # PyTorch differentiable rendering
```

## Quick Start

```python
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:

```python
# 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:

```python
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:

```bash
uv run python examples/render_simple.py
```

### Low-Level Pipeline

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

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

```python
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), lipschitz=2.5)

# Depth is differentiable w.r.t. radius
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)

```python
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)  # estimate Lipschitz constant
```

### Analytic function (NIA/AA recipes)

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

```python
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

```bash
# Run all tests (191 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          |

## 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](LICENSE) for details.
