ENH: use Lagrangian tesselation for smoothing N-body plotting

[cphyc]

One solution to plot smooth-looking N-body data is to approximate them as SPH particles.

This PR explores an alternative (see https://github.com/yipihey/LagrangianPhaseSpaceSheetEasy, credits to @yipihey) using "Lagrangian tesselation".

In the initial conditions, DM particles reside on a grid, from which we can compute a tesselation made of tetrahedron only. As the simulation progresses, these tetrahedrons will be distorted. The method implemented here samples positions within each individual tetrahedron to reconstruct a smoother version of the density field.

Comparison

This is obtained for a 64^3 tiny dataset.
The cherry on the cake: this PR is actually much faster than SPH interpolation because there's no need to build the KD-tree necessary for the SPH interpolation.

Method	Result
This PR ("exact")
This PR (1 sample / tetrahedron)
This PR (10 samples / tetrahedron)
This PR (100 samples / tetrahedron)
NGP interp
SPH interp

TODO:

• add documentation
• add tests
• support cases where particles do not form a uniform grid (e.g. when plotting high-resolution DM particles).

[neutrinoceros]

Too soon for a review I guess but just wanted to say this looks amazing.

Cherry on the cake: this PR is actually much faster than SPH interpolation because there's no need to build the KD-tree necessary for the SPH interpolation.

Can you clarify what you compared ? (since the quality can be adjusted I assume increasing it has a performance cost)

[cphyc]

Can you clarify what you compared? (since the quality can be adjusted I assume increasing it has a performance cost)

I've updated the PR description, so hopefully, that should be better. The most crucial bit is that the SPH-interpolation scales as N log(N) while this implementation scales as N (times the number of samples).

The reason there needs to be some sampling is that we need to compute the volume of the intersection of each tetrahedron with each pixel. In principle, this could be computed exactly but the lazy way is just to (randomly) sample the tetrahedron to estimate the volume.

[neutrinoceros]

The reason there needs to be some sampling is that we need to compute the volume of the intersection of each tetrahedron with each pixel. In principle, this could be computed exactly but the lazy way is just to (randomly) sample the tetrahedron to estimate the volume.

Interesting. I'm guessing there should be a point (critical number of samples) at which the lazy evaluation becomes as expensive as the exact evaluation. Is there an estimate for it ? It's not critical, but if it's known it should probably be documented.

[cphyc]

TBH I'd eventually want to use https://github.com/devonmpowell/r3d/ to compute the intersection between voxels and tetrahedron, but that may be too optimistic.

[yipihey]

Another option that would work well in your cythonic code would be to use the recursive algorithm I first did in IDL and we describe in this paper: Noiseless Gravitational Lensing Simulations https://arxiv.org/abs/1309.1161 There you split tetrahedra recursively until they are smaller than a voxel and then just deposit them as nearest grid point deposit.  Exciting to see this capability in yt! Tom

[cphyc]

After comments from @yipihey, I went ahead and implemented an "exact" integrator.
The good thing is that it is now "exact" (see comments below) and doesn't require any fudge parameter.
The less good news is that it is fairly slow (about ~30s to project 64^3 particles on a 800x800 grid).

The resulting picture looks like this 🍿

It works by noticing that the projection of a uniformly-sampled tetrahedron onto the plane is a simple geometrical figure (see below). There are only two cases to treat: either the projection of the tetrahedron is a triangle or a quadrilateral.

Triangular case	Quadrilateral case

My implementation works as follows:

1. Project the four vertices defining each tetrahedron onto a plane,
2. Find the convex hull,
3. If the hull is made of three points, the maximum density is reached at the location of the fourth point (point 2 on the left figure above).
   If the hull is a quadrilateral, split it into two triangles with the density maximum at the intersection of the two diagonals.
4. For each triangle, iterate over pixels within its bounding box.
5. If the pixel centre is within the triangle, linearly interpolate the projected density at this location.
6. If no pixel was found in the triangle, add the value of the tetrahedron to the pixel closest to the centre of mass and return

At this stage, we have an estimate of the density of all the pixels with centre in the hull. Unfortunately, if we sum these values, we don't get exactly the right integral due to aliasing issues.

8. Repeat 5-6 correcting the density so that it sums to the right value.

The method is actually fairly accurate: since the projected density is a linear function, the mean value in a pixel is exactly the value at the pixel centre, except around the edges.

[yipihey]

Very nice! There is a fast way we developed for projections that is particularly efficient on GPUs that we implemented here (see section 4.4) A Novel Approach to Visualizing Dark Matter Simulations https://arxiv.org/pdf/1208.3206.pdf You first find the back and the front faces of tetrahedra by dotting the normal of each face with the unit vector defining the direction of your view. Then you add the density times the depth of the face (for each image pixel ) of the front facing faces into one image and for the back facing one’s in another image. When you are done with all tetrahedra you subtract the two images and you get back the integrated density along the line of sight. Voila exact projection of the linear elements with few if statements. Hope this approach may help to speed also your current implementation.  Tom

[chummels]

This looks awesome! Excited to see this functionality. Thanks for the work on this, @cphyc !