[WIP] Octree ray tracing

[cphyc]

Implements a ray tracer for octree.
See the tests for how to use it.

The datasets are RAMSES datasets, with a single CPU and an AMR structure with levelmin=1, i.e. the octree is complete from the root node.

• ramses_shallow.tar.gz
• ramses_deep_noamr.tar.gz

TODO

• make it work for sparse octrees,
• handle multiple domains (need to switch the ray from one octree to another one),
• handle leaf cells that are not selected,
• make this work for rays sent in negative directions (-x, -y, -z),
• wire things together with the volume renderer.

[ngoldbaum]

Wow cool!

Note for reviewers that this depends on the octree ghost zones PR and that should be merged first.

Have any eye candy from a RAMSES output?

[cphyc]

In blue, the result of a YTRay and in red (with a small offset), the result from the ray tracing algorithm!

[cphyc]

Here's a somehow nicer example. This is a Kelvin-Helmoltz instability with RAMSES (only a single AMR domain for now). The shader is simply the maximum of the density field along the ray.

In this specific case, a simple mip projection could have been done with yt.(OffAxis)ProjectionPlot, but this demonstrates that the machinery works correctly.

---

Let's imagine you have a 3D data source (ds.all_data(), ds.sphere(...), ...) and let us call it ad.
The way the algorithm works is as follows:

1. Compute the mapping from oct and cell index (in-memory index) to index within data source (the order it appears in ad),
2. Cast a ray through the entire domain: this yields a list of oct and cell index, as well as t-values (coordinate along the line of sight),
3. Reconstruct the index of the data in the data source and extract data. We now have two vectors t and data along the line of sight (their size will be different for each ray),
4. Apply some shader to the vector(s)
5. Display the image.

In the example above, the shader is simply the function lambda data: data.max(). Note that step 1 need only be done once per domain and steps 2-4 are embarrassingly parallel. The computation time scales linearly with the number of pixels. I think the algorithm scales as N log(N), with N the number of octs in the whole domain. This could easily be reduced to n log(n) with n the number of octs in ad.

Regarding performance, the above picture (256² pixels) has been computed in 3s, and the dataset is made of 500,000 octs (~2,000,000 leaf cells).

[cphyc]

A bit more eye candy (in the extent of my artistic capacities): this is a volume rendering of a RAMSES cosmological simulation (64²). Red green and blue indicate the 10⁻²⁹, 10⁻²⁸ and 10⁻²⁷ g/cm³ contours (normalized so that the maximum intensity is 100% in each channel).

One main caveat is that the low-order of the integration makes the result look very boxy when going towards higher resolutions. Here's the same image, now with 512² pixels (and a slightly different integration).

The origin of the boxy look is probably the fact that the interpolation is only performed along the line of sight and does not make use of the neighbouring cell's values. This could in principle be solved using e.g. a trilinear interpolation with the vertices values obtained from ghost zones (see #2425).

[cphyc]

Closing this one in favour of #2610 .