Enzotools Blog

There are a couple other, more important things to consider when making the big volume renderings.  For starters, the entire structure of volume rendering in yt was not really created to generate a series of images -- only a single image.  The idea was that you would prepare a specific image, make it, and move on.  However, for this project, I want to do a zoomin, or possibly a more complicated camera path.

Additionally, one of the first things that we did with the volume rendering was silly: we applied no normalization to the output images.  That was a mistake, I see now.  Part of the reason for this was uncertainty in the correct normalization -- the bias that the user wanted to apply may not be the natural bias from the image.  But more than that, because the rendering algorithm itself was some what holistically settled upon (the original implementation, which we used for shell-style renderings, was not a "correct" implementation of alpha blending) a natural mechanism for scaling did not immediately present itself.  One likely exists, possibly dependent on the field of view, I simply do not yet know it.  This will have to be rectified, because the mechanism used for scaling a set of images will have to be different than the mechanism for scaling an image in isolation, or else frames will jump in brightness during the movie's course.

The final thing that I wanted to change was to add support for stereo rendering.  Rather than repeat any of the amazing discussion from Paul Bourke's website, I'll simply direct you there.  Everything you ever wanted to know about stereo rendering.  (When I was a first year grad student, we actually bought a copy of his site to use locally -- it was our way of showing support for him putting it online, and it also came with a bunch of source code for example applications.)  I first attempted to apply the correct method for stereo, where the view direction is parallel and the total view frustum is shifted.

This did not work.  In fact, it made me realize that all this time, the yt method for volume rendering is in fact ... not really a volume rendering method inasmuch as it is a planar ray-casting method.  Typically when doing volume rendering, there's a perspective applied to the image: the rays all emanate from a single place, creating a frustum.  But for yt, we actually set up a single plane of vectors at the back of the volume and advance that forward across the image.  This is good and bad; it's good in the sense that it's more clear precisely what is going on.  But it's bad in the sense that correct stereo is more difficult.  (Of course, on Bourke's page he has a workaround that may work for this, but I have not yet attempted it.)  Here's a rough depiction of the different between the two methods.

The upshot is that stereo doesn't seem to work unless you go with the "toe-in" method that can cause eyestrain after a long time and shows visible parallax at the edges.  I'm not sure if this is going to be a problem, but because I am not right now eager to rewrite the rendering backend, this is the way it is for the moment.

To set up the stereo rendering, I separated out the rendering mechanism from the objects to be rendered.  Previously, there was a single VolumeRendering object that you could create, raycast through, then discard.  I created a new camera object that accepted a homogenized volume and would call "traverse" on that volume, feeding a back and a front point.  The Volume is then responsible for passing off fixed-resolution grids to the camera, which accumulates an image buffer by calling the ray traversal functions.  The front and back points are essentially the only thing needed to know this order, but the camera also stores its three orientation vectors and its position that describe it in 3D space.  By separating out these two conceptual objects, we undo some of the "single, carefully constructed image" bias that was in the original volume renderer.  (And, we open ourselves up to being able to use the Camera with a hardware volume renderer, should that day ever come.)

So now we have a camera, and it makes images like this:

It's a little dim, but that's a task for another day.  The next step is taking that perspective and turning it into a set of stereo images.  To do that, I added a new class called StereoPairCamera.  It accepts a Camera object and turns it into two camera objects, where the final interocular distance is calculated relative to the image plane width.  As I mentioned above, this only operates via toe-in stereo, so it does this in the simplest manner possible: it moves each of the Left and Right cameras by half of the interocular distance away from the original location, and then recalculates a normal vector to point back at the original center.  Now we can generate left and right images:

Unfortunately, on my laptop (which is my primary/exclusive work computer) I don't have the ability to view these pairs.  To get around that, I wrote a simple stereo pair image viewer in OpenGL and imposed upon my friends that do have a stereo viz wall to test it out -- and after some fiddling with the interocular distance, we got what appeared to be workable stereo pairs.

The full code for generating camera paths as well as stereo pairs is already in yt, but the documentation is still being written; I might also clean up the interface a bit.  Additionally, at some point in the future, the issue of toe-in stereo versus correct parallel-frustum stereo will need to be dealt with; the last thing I really want to do is force people to only use a bad method for generating stereo pairs.  Hopefully that is something that can be dealt with at a later time.

Thanks to the wealth of resources out there for making this a relatively easy task: the aforementioned Paul Bourke website on stereo pairs, the PyOpenGL and PIL teams for making the image pair viewer easy, and everyone else whose work I've built on to make things like this.

Hi all, 

Matt just set up this blog to document some of the developments that are going on in yt, and he asked me to share some of the work we've recently done on a kD-Tree decomposition primarily used for the volume rendering, so here goes.

This semester I took a High Performance Scientific Computing course here at Boulder where a large fraction of our work revolved around a final project.  Luckily for me, I was able to work on a new parallel decomposition framework for the volume renderer inside yt.  What resulted was primarily a kD-tree decomposition of our AMR datasets.  While kD-tree are usually used to organize point-like objects, we modified it to separate regions of space.  

Our kD-tree has 2 primary objects: dividing nodes and leaf nodes.  Each dividing node has a position, a splitting axis, and pointers to left and right children.  The children are then either additional dividing nodes or leaf nodes.  The dividing nodes therefore recursively split up the domain until the "bottom" is composed of all leaf nodes.  Each leaf node represents a brick of data that has the same resolution and is covered by a single AMR grid.  A leaf doesn't have to (and rarely does) encompass an entire grid.  A leaf contains various information such as what AMR grid it resides in, what the start and end indices are in that grid, the estimated cost to render the leaf, and the left and right corners of the space it covers.

So what does this do for us?  Well, the kD-Tree decomposition allows a unique back-to-front ordering of all the data.  This ordering is determined on the fly by comparing the viewpoint with the dividing node position and splitting axis.  By comparing these two, we are able to navigate all the way to the "back" of the data and recursively step through the data.  This decomposition also provides at least 2 different ways of load-balancing the rendering.  By estimating the cost of each leaf, we can subdivide the back-to-front ordering of grids into Nprocs chunks and have each processor only work on their part.  This subdivision is viewpoint-dependent such that the first processor always gets the furthest 1/Nth of the data.  This subdivision is done on each processor and no communication between them is needed.  The second method to implement a "breadth-first" traversal of the tree where we first create Nprocs sub-trees that are not as well load balanced as the first method but are viewpoint-independent.  Each processor is then handed one of the sub-trees to completely render.  This time, when the viewpoint changes, the regions that each processor "owns" stays the same.  

Now one important part in both of these methods is how we composite the individual images back into a single image.  This has been studied a bunch in image processing and is what I think is known as an "over-operation".  For us, all this means is that we combine each pair of images using their alpha channels correctly to blend the images together.  This is essentially the same operation that is done during the actual ray-casting.  This composition is simple for the first decomposition method because we know exactly which order each processor's image goes in, and construct a binary tree reduction onto the root processor.  In the second case where there is a more arbitrary ordering of each of the images, we use the same knowledge from the kD-Tree to combine each pair of images in the correct order. 

For now, this implementation is living in the kd-render branch of the yt mercurial repository, but I imagine that after some more testing and trimming it will make its way back into the main trunk.  If you are interested in using it, you mainly have to change two calls in the rendering:

When instantiating the volume rendering, you need to add the kd=True and pf=pf keywords.  

Then, instead of 

vp.ray_cast()

we use either the back-to-front partitioning

vp.kd_ray_cast(memory_per_process=2**30) # 1GB per processor

or the breadth first 

There are two additional keywords shown above that somewhat useful for large datasets.  Because there can be several leafs that use the same vertex-centered-data of a single grid, it is important to save as much of that information in memory so that disk I/O doesn't hinder the performance.  To do this, we've implemented a rudimentary cacheing system that stores as much of the vertex centered data as possible.  The memory_factor keyword lets us know what fraction of the memory per processor to use, and memory_per_process tells us how much memory each processor has access to. This is fairly experimental and if you run out of memory, I'd suggest cutting one of these two options by a factor of 5-10.  

Anyways, that's about it for the discussion of the new kD-Tree decomposition for volume rendering.  As we make improvements some of this may change, so keep an eye out for those.  Let me just mention that this isn't the only thing this kD-Tree can be used for, so if you have any thoughts or questions, feel free to email me (sam.skillman@gmail.com) or perhaps comment on this post.  

Cheers,

Sam