Volume Rendering datasets in spherical coordinates

[chrishavlin]

So this is a verion of #64 that seems to work pretty well with a few caveats.

This version does 2 things differently:

1. It calculates the cartesian bounding boxes of every block initially (using code from ENH: cartesian cutting planes through non-cartesian geometries  yt#4847), and passes that info down to the shaders as additional vertex attributes.
2. And then those cartesian bounding boxes are used in the ray tracing shader for calculating ray intersections. At each sample point along the ray, the spherical coordinates are calculated and used to sample the texture if it falls in the actual spherical element.

The main caveat is that you need a fairly large number of sample points for the ray tracing since intersection with the cartesian bounding box does not guarantee intersection with the enclosed spherical element.

EDIT: cartsian bbox calculation now works with big elements (via a recursive division when its needed). Also, at least one thing that should be fixed before this goes in is that the cartesian bounding box calculation was designed for fairly small spherical elements and breaks down when blocks span whole quadrants of the spherical coordinate system. e.g., the new example script uses load_uniform_grid with nproc=128, if nproc is smaller (less than 8 or so), the blocks get too big and the bounding box calculation breaks down and results start to get wonky (or nonexistant). So, should fix that...

Finally -- I decided to go with a new PR rather than modifying #64 because I'm abandoning the complicated intersection checks with the underlying spherical element bounds. While I think that approach could work (and could be added on here...), I think it's worth getting the simpler version here merged first and wanted to preserve #64 in its current state to make that attempt easier in the future.

[chrishavlin]

Assuming tests still pass, I think this is ready!

I just updated the docs. and also added an additional pre-processor directive, NONCARTESIAN_GEOM that is used for directive checks in the rendering pipeline that would be common to any non-cartesian coordinate system. So the existing SPHERICAL_GEOM one is now only used where the coordiante conversion is defined in the fragment shader. This makes it easier to add support for other coordinate systems.

[matthewturk]

This is amazing! I've left a few comments, including one about possible artifacts. But I think it's generally good to go!

[matthewturk]

I did not realize you could return tuples with acceptable performance! This is great.

[matthewturk]

Why 0.2 rad / 11 deg?

[chrishavlin]

I... don't remember. Best guess looking at it now is that when playing with what to use, 10 degrees worked well without adding too much extra computation, but I wanted it to be a nice number in radians?

[matthewturk]

Great catch -- this is the kind of thing I usually miss doing.

[matthewturk]

I'm not 100% sure but I think that this may be the source of the (mild!) artifacts we see when we viz ("index", "ones") with projective integration. I attached just-off-axis and almost-on-axis to demonstrate.

To be clear, these artifacts are not deal breakers! But I'm wondering if this is a source of the artifacts, and if we could use more precise values to compute them.

[chrishavlin]

it could be! and i've banged my head on many things looking for the source of those artifacts but it didn't occur to me that it might a precision issue.

[chrishavlin]

so hopefully you're right! will look a bit more closely.

[chrishavlin]

I'm not 100% sure but I think that this may be the source of the (mild!) artifacts we see when we viz

You were 99% right! it was dx_cart (a few lines after this), see #159 (comment)

[chrishavlin]

@matthewturk I've been looking at that artifact issue on and off. and I'm not positive but I think it might be a limitation of the current implementation that uses the intersection of the ray with the cartesian bounding box of each block as the ray start/end: if you imagine discretizing a ray, then the number of points along the ray that fall within a sphere will depend on the curvature of the sphere and the stepsize along the ray and there are situations that can lead to artifacts. I played around with the idea in 2D over in this notebook, here's a figure that plots the integrated distance along a 2D ray that falls within a circle for rays that are parallel to the y axis along with the derivative of that distance with respect to the x position of the ray:

I think those little bumps may be the artifact we're seeing: there are spots where the integrated distance increases depending on how the circle's radius lines up with the discretized ray. Not sure there's a fix for this beyond working out the full intersection calculations.

[matthewturk]

Let's imagine we did do the full intersection calculations. Wouldn't we be able to do this only for the regions that actually intersect a ray -- i.e., only in the fragment shader, and so potentially with a reduced number of calculations? The bounding boxes are how we compute which fragments to compute, but potentially couldn't we inside the fragment shader also use that to discard fragments, or compute the actual intersections?

Also, if there are computations (potentially, for instance, inverse trig functions, square roots, etc) that we can bundle we could do them inside a compute shader step, or otherwise pre-compute and store them.

[chrishavlin]

Let's imagine we did do the full intersection calculations. Wouldn't we be able to do this only for the regions that actually intersect a ray -- i.e., only in the fragment shader, and so potentially with a reduced number of calculations? The bounding boxes are how we compute which fragments to compute, but potentially couldn't we inside the fragment shader also use that to discard fragments, or compute the actual intersections?

yup! and that's what i've done in my attempts -- use the cartesian box check to initially discard fragments before doing the full intersection calculations. let me dig up my last attempt which I think was almost working and remind myself where things were going wrong...

[chrishavlin]

well. that was a big fix!

The latest push calculates an approximate "dx", which in the cartesian branch of the shader is just the grid spacing in x, y,z for a block. For the spherical branch, what I just pushed estimates a cartesian grid spacing of the spherical volume element (using dr and arc lengths along the other angular directions) -- since it's just used to define the stepping along the ray it doesn't have to be super exact. And things just seem to work now without the need for cranking up the sampling_factor!

Here's a projection of the ones field through a unit sphere

and one of the AMR images

In both cases, these look good down to low sample factors! and the artifacts from #159 (comment) seem to be fixed.

[chrishavlin]

ok, just took care of some redundant re-scaling of the domain that I don't think is necessary. also had to replace a deprecated usage of numpy.math in the new cython code here.

I think... this is good to go?? maybe I'll wait a day before merging...

[chrishavlin]

reading through the new docs one last time, found some typos

[chrishavlin]

hmm. had a shader compilation failure just now with this code locally, not sure why the tests didnt catch it.

[chrishavlin]

(so dont merge yet)

[chrishavlin]

oh i see. the test images are incorrect, it's just that tests dont fail on compilation errors. that might be nice to change.