# Experiment 07 (plan): image-in, image-out ray tracer on the GPU

Status: **plan, nothing built yet.**

## 1. Idea

Describe the scene only as images, and trace rays directly on the pixels instead of on building polygons:

| Input image (uint8, same grid) | Meaning | Proposed encoding |
|---|---|---|
| `H` height | obstacle top above ground | 0.5 m per level → 0–127.5 m (0 = open ground) |
| `R` reflectivity | share of power that reflects specularly | reflection loss in dB, 0.1 dB per level → 0–25.5 dB (255 = absorber) |
| `S` scattering | share of the hit power that leaves as diffuse scatter instead of a mirror bounce | 0–255 → 0–1 |
| Tx | transmitter position and height | a parameter for now (later a 4th image) |

Output: one float image, received power per pixel in dB (and optionally a per-bounce breakdown: LOS / 1 bounce / 2+).

Terms. What the brief calls "dispersed" is **diffuse scattering**: a rough surface sends part of the power in
many directions instead of one mirror direction (Sionna's and ITU's *scattering coefficient* S). "Splitting a
beam into 2–3" is **ray splitting** (branching), one way to simulate it. "Dispersion" in radio means a
frequency-dependent effect, so I avoid that word.

Why bother, given we already have `rcm_ml/raytrace.py` and Sionna:

- **Speed and scale.** On polygons the cost is rays × walls (65 s for the 2 km mosaic, no spatial index). On a raster, a ray step costs O(1) whatever the number of buildings, and all rays step together on the GPU. Target: 10⁵–10⁶ rays per map in about a second.
- **Same format as the ML models.** The U-Nets already take images. A tracer that reads the same images can make unlimited training labels for the parallel "learn ray movement" experiment, and the material images become something a model can predict or edit.
- **Optionally differentiable.** Written in PyTorch, the material images (`R`, `S`) can be fitted by gradient descent to WinProp/Sionna maps (phase 4).

## 2. Physics, and the choices that matter

**Ray marching.** Every ray steps through the grid with a DDA (Amanatides–Woo): exact cell crossings, no fixed
step. 2D mode: a hit is entering a cell with `H > 0`. 2.5D mode: the ray has a height too, and a hit is
`z_ray < H(cell)`; rays pass over low buildings, which is what a drone link needs.

**Spreading loss: careful in 2D.** Ray density from a point source falls as 1/d in 2D (rays spread in a circle,
not a sphere). Our polygon tracer fixes this by carrying an extra 1/d^(n−1) weight along each ray (n = 2.6 was
fitted to WinProp). Same here: 2D mode carries the weight, 2.5D/3D mode gets 1/d² from ray density alone.

**Energy per pixel: track-length estimator** (as in `raytrace.py`): each ray adds `power × length of its path
inside the cell` to the cell. Unbiased with any ray count, and smooth with far fewer rays than counting
ray hits per pixel. On the GPU it is one `scatter_add` per step.

**Wall normals from a raster.** The hard part. A wall drawn in pixels is a staircase, and the per-pixel normal
from it is noisy, which spreads mirror reflections into a fan. Plan: precompute a normal map from the gradient
of a lightly blurred occupancy image (Sobel on Gaussian σ ≈ 1 px), and test the error with a single tilted wall
(§4, test 2). If it is not good enough: sub-pixel walls from a signed-distance field.

**Reflection and scattering at a hit**, with a constant number of rays (GPU-friendly):

- with probability `1 − S`: mirror bounce, power × 10^(−R_dB/10);
- with probability `S`: new direction drawn from a Lambertian lobe around the normal, same power factor.

This "one child per hit, chosen at random" is unbiased on average and never grows the ray count. The literal
"split into k children" (power / k each) is also implemented as a switch for small bounce counts, because it is
less noisy per ray but grows as k^bounces.

**Diffraction around corners.** Without it, everything behind a building is dark, and the RadioMapSeer
experiments showed that is most of the street-level map. Options, in order:
1. none (v0: shows the gap);
2. corner pixels found from the raster (convex corners of the occupancy mask) emit a fan of secondary rays with the calibrated loss from `raytrace.py` (6 dB + 20 dB per 90°);
3. in 2.5D, also roof-edge diffraction (knife-edge over the building top) for drones.

**Not modelled at first:** phase/coherent sums (powers add incoherently, like our tracer), transmission through
buildings, ground reflection, antenna patterns.

## 3. Implementation

```
rcm_ml/beam/
  scene.py     polygons / RadioMapSeer / mosaics → H, R, S PNGs; normal map
  march.py     PyTorch wavefront tracer (MPS): all rays step together, scatter_add deposit
  ref_numba.py same algorithm on the CPU with numba, for checking the GPU version
  eval.py      tests and maps vs our polygon tracer, IRT2, Sionna; timing
```

Wavefront loop: state tensors (x, y, z, dx, dy, dz, power, bounces, alive) of length N_rays; each iteration does
one DDA step for all rays, deposits, handles hits, kills rays below a power floor (−160 dB) or past max
bounces. Fixed tensor sizes, so no recompaction until a large share is dead (then compact once).

Rough cost: 256 px map, 10⁵ rays, ~400 steps per bounce, 3 bounces ≈ 10⁸ cell visits ≈ 1,200 loop iterations
of ~15 small kernels. Probably 0.3–1 s on the M2 Max in PyTorch. If the per-kernel overhead dominates, move the
inner loop into one Metal kernel (Taichi or `torch.mps` custom kernel) and keep the same interface.

## 4. Tests and comparisons

| # | Test | Pass criterion |
|---|---|---|
| 1 | Empty map | matches the analytic distance law to < 0.5 dB mean, < 1.5 dB p99 |
| 2 | One wall at 0°, 30°, 45° (image-source solution) | reflected power within 1 dB; tells us if raster normals are good enough |
| 3 | `S` = 0 vs 0.3 vs 1 on one street canyon | scatter spreads energy into side streets; total energy conserved |
| 4 | RadioMapSeer, 20 test maps × 2 Tx (same as 03 §7), 2D, same mechanisms and losses as our polygon tracer | within ~1 dB RMSE of our polygon tracer; vs IRT2 around its 4.3 dB |
| 5 | Ray count 10³ → 10⁶ | RMSE vs a 10⁷-ray run; where the curve flattens = rays needed |
| 6 | 512 m / 1 km / 2 km mosaics | time per map vs the polygon tracer's 1.4 / 7.3 / 65 s; Sionna far-field agreement (06) |
| 7 | 2.5D, drone at 30 / 60 / 120 m | vs `sionna_drone.py` maps (seed-averaged), cells with signal and RMSE |

## 5. Phases

1. **Spike (½ day):** 2D marcher on MPS, no hits handled, empty map + test 1, rays/second. Decides PyTorch vs a custom kernel.
2. **v0 2D:** reflection + scattering, tests 2–5. Material images: `R` = 10 dB everywhere (the calibrated value), `S` swept.
3. **Diffraction from raster corners**, test 4 again, test 6. This is where it becomes a usable fast label generator.
4. **Optional:** differentiable version (soft hits); fit `R`/`S` images, or a few global values, to IRT2.
5. **2.5D heights**, test 7.

## 6. Open questions

- 2D first (RadioMapSeer has flat 25 m buildings and WinProp truth) or heights from day one (drone links)?
- Do we want the per-pixel `R` and `S` to come from real data (OSM building material tags, land use), or are they free parameters to fit?
- Pixel size: 1 m (RadioMapSeer) for street level; 2–4 m is likely enough at drone altitudes.
